简介针对电力系统网损偏高、电压质量不佳等运行优化需求该资源提供基于粒子群算法的无功优化实现代码采用IEEE30节点标准测试系统作为算例适合电力专业研究生、科研人员与算法爱好者参考学习。代码完整覆盖粒子群寻优流程每个粒子代表一组无功控制变量通过适应度评估、个体最优与全局最优迭代更新并在速度位置调整中考虑惯性权重和学习因子同时结合切线逼近法处理非线性约束部分脚本内涉及内点法的雅可比、海森矩阵生成与主计算模块划分便于逐段对照理解。源码包共14个文件全部为m脚本压缩后仅11KB规模小巧但逻辑完整适合快速阅读与二次开发。目前已有3357人学习运行该代码可观察粒子群算法在三十节点系统上的收敛过程理解参数调整对电压水平及网损的影响并为移植到其他算例提供可直接修改的框架。1. 粒子群无功优化是什么先用IEEE30跑通一套MATLAB实现夜间值班或做规划方案时最让人头疼的场景之一系统能运行但部分节点电压偏低、网损偏高翻来覆去试了几组电容器补偿方案也不见起色。这种问题一般交给无功优化来解——调整发电机电压、变压器分接头、电容器组三样手段在满足潮流约束的前提下把有功网损压下去。这套基于粒子群算法的MATLAB源代码正是把整套流程落在IEEE30标准测试节点系统上的可复现实现从节点数据和支路数据读入、潮流计算到粒子群迭代搜索和结果输出代码包里都串好了。适合电力系统方向的研究生、电网运行分析刚入门的工程师以及想快速拿一个标准算例验证PSO改进策略的算法研究者。下面从模型、编码、离散处理和调试排错几个层面把这份源代码的实际用法讲透。2. 无功优化模型与IEEE30节点数据结构目标函数、约束与潮流接口2.1 无功优化的数学模型网损最小与罚函数电力系统的无功优化本质上是一个带大量等式约束和不等式约束的非线性规划。目标函数最常见的是有功网损最小也可以写成网损加电压偏差的加权形式。数学模型上若以 Ploss 表示全网有功损耗则min Ploss sum_{(i,j)∈Lines} G_ij (V_i^2 V_j^2 - 2 V_i V_j cos(θ_i - θ_j))其中 G_ij 是支路电导V_i、V_j 是支路两端节点电压幅值θ 是相角。等式约束是各节点的有功、无功注入功率必须满足潮流方程不等式约束则集中在四类发电机节点电压幅值、可调变压器变比、并联电容器组投入容量、负荷节点电压幅值。把这套约束写进粒子群通常不使用专门的约束函数库而是采用罚函数。也就是先算潮流看哪些状态变量越限把越限量平方后乘以一个惩罚系数加到网损上去。代码里适应度函数的一种典型骨架如下function f objFun(x, data) % x: 控制变量向量顺序为发电机电压、变压器变比、电容器容量 Ug x(1:data.nG); tap round(x(data.nG1:data.nGdata.nT)); Qc round(x(data.nGdata.nT1:end)); [V, converged] powerFlow(data.bus, data.branch, Ug, tap, Qc); if ~converged f 1e8; % 潮流发散直接给一个极大罚值 return; end Ploss calcLoss(V, data.bus, data.branch); Vmax 1.05; Vmin 0.95; ov sum(max(0, V - Vmax).^2) sum(max(0, Vmin - V).^2); f Ploss 5000 * ov; % 权重按网损量级调整 end这个函数里有一个容易被新手忽略的地方惩罚系数不是越大越好。IEEE30系统优化前的网损一般在5~8MW标幺值下是0.05~0.08p.u.如果惩罚系数设为50000电压越限0.01p.u.就会产生5单位的罚值目标函数完全被电压项主导粒子会只顾抬电压而不管网损如果设为50电压越限一点也可以被网损降低抵消优化结果可能长期停留在越限状态。我一般会在正式优化前把初始潮流跑一遍统计电压越限水平再按“越限罚值约为正常网损的3~5倍”来定系数。适应度函数中的 round 是对离散变量的量化。发电机电压是连续变量不取整变压器变比和电容器容量必须取整。这一行看起来简单但它关系到粒子群能不能收敛到实际可执行的档位后面第4章会专门展开。另外要提一句潮流不收敛时返回1e8的处理方式对粒子群来说并不完全安全。粒子群在更新个体最优时会拿当前适应度和历史最优比较如果某个粒子恰好飞到了潮流发散区域适应度变成1e8这个值可能比某些差解的适应度大但不会被误认为最优。可一旦全局最优还没被初始化好比如第一代恰好全是发散解1e8也可能成为gbestFit后面所有正常粒子都会被这个错误最优吸引一段时间。更稳的做法是让objFun在发散时返回NaN并在主循环更新最优时做isnan判断。2.2 IEEE30节点的数据组织方式与控制变量索引IEEE30节点系统是电力系统优化算法里最常用的标准算例之一规模为30个节点、41条支路、6台发电机带4个左右的可调变压器支路和若干并联电容器组。发电机的挂接节点通常是1、2、5、8、11、13这个信息会直接影响控制变量的维度。打开代码包的数据文件时第一件事不是运行主程序而是先核对数据矩阵的列顺序。比如bus矩阵如果按Matpower格式组织各列含义如下列位置字段说明1节点编号1~302节点类型1代表PQ2代表PV3代表平衡节点3有功负荷单位MW4无功负荷单位Mvar5~8并联导纳电容器补偿通常与此相关11电压幅值下限标幺值12电压幅值上限标幺值branch矩阵里第9列是变压器变比非1的那几条支路就是可调变压器。很多翻车都是因为索引错位把节点类型那一列当成电压上下限或者把线路首末端编号当成控制变量结果优化出来的结果完全没法解释。我的习惯是先写一小段脚本打印bus(1:5, :)、branch(:, 9)的非1位置确认数据语义后再往下调主函数。控制变量在粒子群中统一编码成一个行向量x。前nG维是发电机电压中间nT维是变压器变比档位编号后面nC维是电容器投切组数。这样编码的好处是pbest、gbest、速度矩阵都是同构的更新公式不需要区分变量类型离散物理量在评估前再处理搜索空间仍是连续空间。2.3 潮流计算与优化循环的配合方式无功优化里的PSO每评估一个粒子都要解一次潮流粒子群60个、迭代100代就是6000次潮流计算。IEEE30规模小牛顿法单次只需几毫秒但若每次都从平直初值开始累计耗时就会变得明显如果你后续要把这套代码迁移到更大系统这就是第一个需要改的地方。资源代码里常用的提速方式是把上一次收敛的潮流结果作为下一次评估的初值相邻两次评估之间控制变量变化很小初值越接近真解牛顿法迭代次数就越少。如果代码包自带的潮流函数没有做这个优化我建议你自己改。在objFun内部维护一个持久的V0、theta0初值变量每次收敛后更新它若潮流发散再把初值重置回平启动值。这属于很小但收益很直接的改动通常能让总耗时下降20%~40%。数据接入部分常见写法是这样mpc loadcase(case30_ieee); bus mpc.bus; branch mpc.branch; gen mpc.gen; genBus gen(:, 1); % 发电机挂接节点 tapIdx find(branch(:, 9) ~ 0); % 可调变压器支路 capNodes find(bus(:, 5) 0); % 有并联补偿的节点这里如果数据文件不是Matpower格式就按资源里的read_data函数实际返回的字段来。无论什么格式都需要确保三个索引的维度加起来等于控制变量向量的总维度。有人把并联电容补偿节点加到了变压器变比索引里导致维度多了一位程序直接下标越界这个错很隐蔽因为数据文件里数字都对就是索引错位。3. 粒子群算法在MATLAB中的实现框架初始化、速度更新与边界处理3.1 粒子结构与初始化一组控制变量就是一个粒子粒子群算法的思路很直观一群粒子在控制变量空间里飞行每个粒子代表一组调节方案例如6个发电机电压、4个变压器档位、3个电容器组容量合并成一个13维向量。每次迭代时粒子根据自己历史最好的位置和整个种群历史最好的位置来调整速度再移动到自己新的位置。下面这段代码是这套资源里PSO初始化的典型结构nPop 60; dim data.nG data.nT data.nC; lb [data.ug_min, data.tap_min*ones(1,data.nT), data.qc_min*ones(1,data.nC)]; ub [data.ug_max, data.tap_max*ones(1,data.nT), data.qc_max*ones(1,data.nC)]; pop lb (ub - lb) .* rand(nPop, dim); vel zeros(nPop, dim); fitP inf(nPop, 1); pbest pop; gbest pop(1, :); gbestFit inf;这里的lb和ub按控制变量顺序拼好顺序必须和objFun里解包的顺序完全一致。粒子数60是经验值维度13的情况下够用如果你想缩小种群做实验最低建议保持在维度数的3倍以上否则搜索空间太稀疏很容易掉进局部最优。3.2 速度与位置更新惯性权重降速模型标准PSO的更新公式是v wv c1r1*(pbest - x) c2r2(gbest - x) x x vw是惯性权重c1控制粒子向自身历史最优飞的水平c2控制向全局最优飞的水平。代码里实现为线性递减是很常见的思路早期w0.9粒子保持比较大的探索步长后期w降到0.4让粒子在最优解附近做精细搜索。for it 1:maxIt w 0.9 - (0.9 - 0.4) * it / maxIt; r1 rand(nPop, dim); r2 rand(nPop, dim); vel w .* vel c1 .* r1 .* (pbest - pop) c2 .* r2 .* (gbest - pop); vel max(min(vel, Vmax), -Vmax); % 速度限制 pop pop vel; pop max(min(pop, ub), lb); % 边界截断 for i 1:nPop f objFun(pop(i, :), data); if ~isnan(f) f fitP(i) pbest(i, :) pop(i, :); fitP(i) f; end if ~isnan(f) f gbestFit gbestFit f; gbest pop(i, :); end end end这段代码里有两个细节值得注意。第一速度限制解耦了更新公式和物理边界如果没有Vmax粒子可能在两个边界之间反复弹跳第二isnan判断非常关键。前面2.1节说了发散粒子返回NaN更稳妥这里直接用isnan把发散结果排除在最优更新之外比返回1e8再比较更稳。c1、c2在代码里一般取2.0这个值是标准PSO论文里常见的默认值如果你发现收敛太慢可以把c2提到2.2如果发现振荡就把两者都降到1.5左右。3.3 边界约束的三种处理方式与离散变量的坑边界截断pop max(min(pop, ub), lb)是最简单的处理方式对发电机电压这样的连续变量很有效越界就拉回边界速度方向保留下次迭代可能继续飞出界但因为适应度惩罚会把它推回来。但对变压器变比和电容器组这类离散变量简单截断会带来一个问题粒子位置越界后总是停留在边界档位整群粒子在边界处堆积离散档位的多样性下降。更讲究一点的实现是边界反射法如果粒子往上越界就把位置反射回边界内侧同时把速度方向反转。对连续变量这个做法能减少边界堆积对离散变量效果也不差。但对大多数IEEE30无功优化场景我建议先用截断法因为它简单可靠真正需要操心的是适应度函数中round带来的不连续性。3.4 Vmax取值与初值设计Vmax的取值这里再展开说一下。控制变量的物理范围不一样发电机电压是0.95~1.10跨度0.15变压器变比0.90~1.10跨度0.20电容器投切组数可能是0~10。如果Vmax在所有维度上都设成同一个值那么电压维度可能飞得太快档位维度又飞得太慢。比较合理的做法是按各维度范围的比例设置比如通用做法是取范围的10%~15%。用代码实现就是range ub - lb; Vmax 0.15 * range; % 每个维度独立的速度上限这样每个维度的最大步长约为该变量允许范围的15%粒子一步最多跨过约7档变压器不会出现一步从0挡飞到10挡的夸张轨迹。初值方面uniform随机采样在边界内生成即可不要对初始种群做过度的可行化筛选因为无功优化里初始解越限是常态罚函数会把粒子引导回可行域。过度筛选反而会降低初期的探索性。4. 变压器分接头与电容器组的离散处理混合变量的量化与评估4.1 变压器变比档位索引与连续量的映射变压器分接头调节是无功优化最重要的离散变量实际设备并不是连续可调而是在一个范围内按固定档距跳变。IEEE30标准数据里可调变压器变比一般取0.90~1.10常见分成17档或33档步长0.0125或0.00625。粒子群产生的连续位置比如1.0563物理上并不存在必须映射回最近的可用档位。常见做法是先把连续位置映射到档位序号再反算变比nTaps 17; % 分接头总档数 tapIdx round((tapPos - data.tap_min) / (data.tap_max - data.tap_min) * (nTaps - 1)) 1; tapIdx min(max(tapIdx, 1), nTaps); % 防止四舍五入后越界 tapVal data.tap_min (tapIdx - 1) * (data.tap_max - data.tap_min) / (nTaps - 1);这段代码逻辑是先把连续的tapPos线性映射到[1, nTaps]区间round成整数档位再反算回实际变比。最后的min/max很关键否则当tapPos恰好略低于下限时round后的索引可能是0。如果数据文件里分接头不是均匀档距就不能用线性公式要改成查表tapVal tapTable(tapIdx)表里每个档位对应实际变比。4.2 电容器组按组投切的整数化并联电容器补偿容量通常是按组投切每一组容量固定比如5Mvar一组最多投入10组那么总补偿容量就是0、5、10……50Mvar的整数倍。在粒子位置里这个维度保存的可以是“组数”的连续量评估前同样取整qcPos x(data.nG data.nT j); % 第j个补偿节点对应的连续位置 qcGroup round(qcPos); qcGroup max(min(qcGroup, data.qc_max_group), 0); qcMvar qcGroup * data.qc_step; % 实际投入到潮流的补偿容量注意这里qcPos的物理单位是“组数”而不是Mvar。有些实现会让粒子直接搜索Mvar值评估时除以组步长再取整效果类似但容易因为单位不一致导致速度限制难调。我建议统一把位置变量的单位定义成“组数”这样上下限、Vmax都是无量纲整数附近的连续值参数好理解。4.3 为什么不在粒子更新阶段直接取整很多初用混合变量PSO的人会问反正最终要取整为什么不能每次更新位置后马上round实际上这会破坏粒子群算法的基本机制。粒子群的速度更新依赖pbest和当前位置之间的差向量如果当前位置已经被round成档位那么这个差向量丢失了亚档位的梯度信息粒子只知道“自己在第5档”但不知道“自己其实偏向第6档”后续速度更新会在同一个档位内反复振荡收敛速度明显变差。正确做法是把取整延迟到评估函数里。粒子的位置、速度、pbest、gbest全部保持浮点只有调用objFun时才对变压器档位和电容器组做round。这样一来粒子在连续空间正常飞行评估结果却能反映离散模型。这也是这套源码在IEEE30系统上稳定工作的核心设计之一。4.4 个体最优的比较基准浮点位置还是离散方案还有一个容易被忽略的细节pbest和gbest保存的应该是浮点控制向量还是取整后的离散方案我的习惯是保存浮点向量但每一次适应度比较时都基于round之后的结果。换句话说pbest更新不能直接拿objFun返回的f去比较否则会出现一个尴尬情况连续位置是5.9的时候适应度好存进pbest后下一次粒子飞向5.9但实际可执行的档位是65.9和6的差异在连续空间里被当成“还有0.1的提升空间”实际评估却总是被round成6。解决的办法是修改个体最优的比较逻辑评估粒子时先把当前粒子的离散版本算出来然后以离散版本去更新pbest。这样pbest中保存的连续值虽然仍是浮点但它对应的离散方案一定是历史上真实评估过的方案。代码上可以写成xDiscrete data.roundMap(x); % 对离散维度统一取整 f objFun(xDiscrete, data); % 用离散方案评估 if f fitP(i) pbest(i, :) xDiscrete; % 个体最优也存离散方案 fitP(i) f; endroundMap这个函数就是把上一节里的变压器档位和电容器组取整逻辑封装起来统一调用避免在主循环里重复写round。这样改完之后收敛曲线后段的抖动通常会明显减少因为向量与档位之间不再存在“永远差半步”的假象。5. 避坑收敛抖动、发散与结果不稳的排查与修复5.1 迭代中段网损曲线突然冲高又回落现象迭代到第70代全局最优网损从5.3MW瞬间跳到60MW再过几代又回到5.3MW收敛曲线出现一根尖刺。很多人第一反应是数据错了或者算法崩了。原因某个粒子在位置更新后飞到了潮流不收敛区域objFun返回了一个非常大但有限的值比如1e8。如果主循环里没有对NaN做处理只是简单比较数值大小那么这个粒子可能在某次比较中被误选为新的gbest污染了全局最优其他粒子被短暂吸引后又被正常适应度拉回来曲线上就出现尖峰。解决把潮流发散粒子的适应度定义为NaN并且在更新pbest和gbest时都加上isnan过滤。代码写法在第3.2节里已经给出两处更新前都判断~isnan(f)。虽然只改两行却是粒子群实现里最常见的修正方式。5.2 变压器分接头在两个档位之间反复横跳现象多次运行后某个变压器分接头的最优结果一直在第9档和第10档之间变化网损相差不到0.01MW收敛曲线后期几乎是一条平线但两端结果不一致。原因这两个档位在目标函数上几乎是等效的粒子群在平坦区域里无法形成有效的梯度方向。这不是算法实现错误而是目标函数本身存在多个近似等优的解。解决可以先判断档位步长是否过细。如果0.01MW的网损差异对应一个档位实际工程中完全可以通过减少可调档数来处理如果不改模型则在输出阶段做等效解聚类把最优网损落在同一小范围内的方案归为一类再按“调节动作最少”或“变压器动作次数最少”的原则选取最终方案。这个思路比强行让算法收敛到唯一解更符合实际调度。5.3 发电机电压顶到上限负荷节点电压仍然偏低现象优化结果里发电机电压几乎全部贴在上限1.10p.u.但部分远端负荷节点电压依然低于0.95p.u.。网损虽然降下去了电压质量却不达标。原因控制变量的灵敏度分布不合理。发电机电压对远端节点的抬升作用有限而电容器补偿节点可能没有覆盖到那些薄弱节点或者补偿容量上限给得不够。这属于建模层面的问题不是粒子群本身的问题。解决先用一次潮流做灵敏度分析找出对薄弱节点电压灵敏度最高的补偿节点把电容器组挂到这些节点上重跑。如果补偿节点已经覆盖但容量上限不够就把上限调大如果模型没得改就在目标函数里加入电压偏差权重让算法在网损和电压质量之间做权衡。代码里常见做法是把objFun的惩罚部分扩展成带权重的组合voltageLoss sum(max(0, V - Vmax).^2) sum(max(0, Vmin - V).^2); f Ploss lambda1 * voltageLoss lambda2 * flowPenalty;lambda1 按电压偏差量级和网损量级的比例来取我一般从0.1开始试观察两者的变化趋势再微调。5.4 固定数据多次运行结果差异大到不可接受现象同一份数据同一台机器连续跑三遍网损结果分别为5.18MW、5.26MW、5.11MW最优方案中变压器档位甚至差两档。如果拿单次结果去写报告很难自圆其说。原因粒子群初始化、速度更新、位置更新都依赖rand产生的随机序列没有固定随机种子的话每次运行都是一次新的随机搜索。解决在main函数开头加rng(2024)固定种子。然后做一轮多运行统计每次用不同种子跑共20次记录每次最优网损取最小值作为最终方案同时看平均值和最大值。平均值和最优值差距在0.1MW以内说明算法稳定性可接受如果差距太大优先检查Vmax和粒子数——这两项是影响稳定性的主要参数。固定种子还带来一个附加好处你可以把每次的实验参数完整存档后续所有调试都有后悔药。5.5 Matpower版本差异导致数据读取报错现象代码包在某台机器上运行正常换到另一台机器直接报错提示Undefined function runpf或者提示结构体字段读取失败。原因这套源码依赖Matpower的潮流函数而Matpower不同版本对runpf的输入输出定义有差别。老版本返回多个输出参数新版本往往返回一个mpc结构体如果代码按旧接口写的新版本里第一步就挂了。解决先执行addpath把Matpower路径加进来再用which runpf确认调用到的确实是目标版本接着打开代码里的潮流封装函数看它期望的调用方式。如果资源自带了潮流函数并兼容多个版本那就优先用自带封装不要直接改调runpf如果你用的是新版本Matpower而代码里还是[baseMVA, bus, gen, branch] runpf(mpc)这种老式写法改成mpc2 runpf(mpc)后从mpc2.bus、mpc2.branch取数即可。6. 验证与进阶固定种子复现、性能对比与迁移到实际系统拿到这套代码跑通只是第一步真正能写进报告的是验证后的结果。我一般按三件事收尾。第一件事固定随机种子做多轮复现在main开头加rng(42)跑20次统计最优、平均、最差网损和每次找到的变压器档位方案确认结果在可接受区间内。第二件事和基准对比把发电机电压设为1.0、变压器变比设为1.0、电容器全部不投算一次纯初始潮流得到基准网损再和优化后的网损对比算出下降百分比。如果优化后网损只比基准低不到1%大概率是参数没调开或约束设置太紧而不是粒子群本身不行。迁移到自己的系统时最容易翻车的地方是控制变量索引。IEEE30的发电机节点是1、2、5、8、11、13换成实际系统的数据后发电机数量、节点编号排序、支路中变压器位置都会变。如果沿用原代码里写死的索引轻则优化结果不符合预期重则直接数组越界。我的习惯是写一个独立的extractControls()函数只负责从bus、branch、gen三个矩阵里提取控制变量索引并在主函数外单独打印一遍核对数量与objFun解包顺序一致后再交给PSO。进阶方向有两个比较实用。一是把PSO的最终结果作为局部搜索的初值再调用一次梯度类优化器或枚举搜索对离散变量做微调通常能把网损再压0.1%~0.3%二是针对分接头和电容器组做定向扰动优化后在最优解附近逐档枚举观察是否存在等效解这个方法对运行方式分析特别有用——用它来判断变压器是否真的需要动作。以前我不固定随机种子跑到一个可用的结果就当作结论回头看才发现其中不少是运气从那以后我每一次跑无功优化都会强制走一遍固定种子、多轮统计、保留收敛曲线这个流程这个习惯让我避开了很多“假最优”带来的回头返工。希望帮到你。本文还有配套的精品资源点击获取