1. 为什么无功优化是电力系统里一块值得啃的硬骨头如果你刚接了一个课题导师只丢给你一句话用粒子群算法做电力系统无功优化测试系统用IEEE14节点平台用Matlab。剩下的全靠自己。这时候最茫然的不是怎么调PSO参数而是先搞清楚我们到底在优化什么为什么偏偏选这个算法、这个系统先说无功优化本身。电力系统里节点电压要想维持在允许范围内光有有功功率是不够的系统里必须按时供给足够多的无功功率。无功不足电压往下掉无功过多电压往上冲而无功在输电线路上来回流动又会产生额外损耗和电压降落。所谓无功优化就是在满足各种运行约束的前提下合理调节发电机机端电压、有载调压变压器的分接头位置、并联电容器等无功补偿装置的投切容量让某些运行指标达到最优——最常见的目标就是全网有功损耗最小同时兼顾电压质量。这个问题难在哪里难在它是一个大规模、非线性的混合整数规划问题。潮流方程本身是非线性的变压器分接头又是离散变量无功补偿容量既有连续特性又有离散投切逻辑再加上一大堆安全约束。传统的线性规划、非线性规划方法要么需要求导、对初值敏感要么处理离散变量很费劲。这时候启发式算法就派上用场了而粒子群算法因为结构简单、参数少、收敛快成了最热门的选型之一。它不需要目标函数的梯度信息只要你能给每个候选解算出一个适应度值它就能迭代搜索这和无功优化的适配度很高。至于IEEE14节点那是国际通用的标准测试系统14个节点、5台发电机、20条支路规模不大不小既能体现算法效果又不至于在调试时被数据量淹没。对初学者来说在这个系统上把PSO和无功优化的耦合逻辑跑通后面移植到IEEE30、IEEE118甚至实际电网数据只是改数据文件的事。这篇文章我不打算讲空泛的理论而是把我实际跑通这个课题的完整思路、代码结构和踩坑记录都摊开来讲。不论你是本科生做课程设计、研究生入门科研还是刚接触无功优化的工程师按这个路线去实现应该能省掉不少瞎摸索的时间。2. IEEE14节点系统的数学建模与Matlab数据准备2.1 系统的物理构成与控制变量盘点在做任何优化之前先要把IEEE14节点的家底弄清楚。这个系统虽然叫14节点但麻雀虽小五脏俱全节点1是平衡节点同时也是发电节点节点2、3、6、8是PV节点各有一台发电机节点9附近挂有并联电容器作为无功补偿装置。支路里面4-7、4-9、5-6这三条是带可调分接头的变压器支路其余是普通输电线路。控制变量一共9个我习惯用九维向量来表示一个粒子变量类别数量取值范围标幺值作用发电机机端电压5个0.95 ~ 1.10励磁系统直接可调变压器分接头变比3个0.90 ~ 1.10改变无功分布走向补偿电容容量1个0 ~ 0.30就地供给无功这里有个容易踩的坑很多同学会把发电机有功出力和负荷大小也当作控制变量这是不对的。无功优化的前提是典型运行方式下有功潮流已经确定优化只作用于无功相关的控制手段发电机有功出力保持初始设定值不动。状态变量则是那些不直接可控、但必须监视的量包括所有PQ节点的电压幅值、发电机的无功出力、部分支路的无功潮流。PSO只会搜索控制变量空间状态变量通过潮流计算得到如果越限就靠惩罚函数拉回可行域。2.2 Matlab中组织系统参数的核心数据结构写代码的第一步不是写PSO而是先把系统数据在Matlab里组织好。Matpower用的是矩阵结构bus矩阵、branch矩阵、gen矩阵但我自己写代码时更习惯用结构体因为可读性更高% 系统基准容量 baseMVA 100; % 节点数据节点编号, 节点类型(1平衡,2PV,3PQ), 有功负荷, 无功负荷, 电压幅值初值, 电压上限, 电压下限 bus [ 1 1 0.000 0.000 1.060 1.10 0.95; 2 2 0.217 0.127 1.045 1.10 0.95; 3 2 0.942 0.190 1.010 1.10 0.95; % ... 其余节点略 ]; % 发电机数据所在节点, 有功出力, 无功出力下限, 无功出力上限, 机端电压设定值 gen [ 1 0.000 0 inf 1.060; 2 0.400 -0.4 0.5 1.045; % ... 其余发电机略 ]; % 支路数据起点, 终点, 电阻, 电抗, 对地导纳, 变比0表示非变压器支路 branch [ 1 2 0.01938 0.05917 0.0528 0; 4 7 0.00000 0.20912 0.0000 0.978; % ... 其余支路略 ];数值有两点要注意。第一变压器支路的变比在潮流计算里通常放到起点侧或者终点侧方向不同初始值写法也不同建议以自己用的潮流函数说明为准或者直接用Matpower里IEEE14的标准数据文件。第二补偿电容器不是简单往bus矩阵的负荷里加一个负的无功它是并联导纳无功出力会随节点电压的平方变化不是恒定值。这两点如果搞错后面潮流计算怎么调都不对劲。2.3 潮流计算引擎的选择粒子群每迭代一次每个粒子都要跑一遍潮流计算来评价适应度。以种群40、迭代80次为例一共要跑3200次潮流。这时候潮流求解器的稳定性和速度就很关键了。我的建议是整个项目早期的验证阶段直接用Matpower的runpf函数做潮流引擎把精力集中在PSO和无功优化的耦合逻辑上等算法逻辑全部验证通过再考虑自己写牛顿-拉夫逊法代码替换。原因很简单手写潮流涉及节点导纳矩阵、雅可比矩阵、稀疏分解等一大堆细节这些环节本身就会引入bug如果优化结果不对你很难分清是PSO的问题还是潮流计算的问题。如果你决定自己写潮流流程是根据bus和branch数据构造节点导纳矩阵Y把发电机有功、无功注入看作是节点注入功率先把PV节点和无功初值折算进去迭代计算不平衡量求解修正方程更新电压幅值和相角收敛后从电压分布反推所有支路功率和全网网损。这个流程说起来就四步但每一步都有无数细节所以我还是推荐先用Matpower跑通闭环再来做替换。3. PSO原理与面向无功优化的算法设计3.1 粒子群算法的核心迭代逻辑粒子群算法的灵感来自鸟群觅食。每只鸟粒子代表解空间里一个候选解整个鸟群通过个体经验和群体经验不断调整飞行方向。在无功优化问题里一个粒子就是一组9维的控制变量也就是一套发电机电压变压器变比补偿容量方案。速度更新和位置更新是算法的两条腿v(i,:) w * v(i,:) c1 * rand * (pbest(i,:) - x(i,:)) c2 * rand * (gbest - x(i,:)); x(i,:) x(i,:) v(i,:);这组公式看着简单里面的门道不少。速度项w*v(i,:)是惯性部分让粒子保留之前的运动趋势不至于每代都彻底转向pbest部分是个体记忆把粒子拉向自己历史上最好的位置gbest部分是群体共享把粒子拉向全群最优位置。三个力合成新的速度粒子在解空间里飞。初学的时候我犯过一个典型的错误把位置更新和速度更新混在一起先更新位置然后拿新位置去算速度导致收敛曲线一路乱跳。正确做法是先用旧位置算出新速度再用新速度更新位置两者必须严格先后顺序。3.2 关键参数怎么定种群、迭代、学习因子与惯性权重PSO参数选择是有经验可循的我实际跑下来推荐这样一套初始配置参数推荐值调整思路种群规模30~509个变量40个粒子足够最大迭代次数50~100先跑80次看收敛趋势学习因子c1, c22.0, 2.0经典的对称配置兼顾个体与群体惯性权重w0.9递减至0.4前期全局探索后期局部精细搜索惯性权重的线性递减是最实用的一种方案。迭代早期w大粒子速度快有利于跳出局部最优迭代后期w小粒子速度慢有利于在最优解附近精细搜索。递减公式写成w w_max - (w_max - w_min) * iter / max_iter;有个说法是学习因子c1和c2也可以随迭代阶段变化前期c1大、后期c2大但在9维这种小规模问题上收益不明显我试过差异不大用对称2.0就够了。真正影响结果的反而是下面要说的边界处理和归一化。3.3 位置速度边界与归一化技巧9个控制变量的取值范围差异很大发电机电压在0.95到1.10之间晃变压器变比在0.90到1.10之间微调补偿容量上限0.30。如果直接在这三个不同尺度上搜索速度更新的量级会失衡——电压方向稍微一动补偿容量方向上可能已经飞出去十万八千里了。我推荐的做法是归一化。所有控制变量先映射到[0,1]区间PSO在统一的单位超立方体里搜索评价函数内部再还原成真实物理量% 解码从归一化值 [0,1] 还原为物理量 orig lb x_norm .* (ub - lb);这样做有两个好处一是不同变量对位置扰动的敏感度一致算法不用纠结每个维度的速度上限二是边界处理统一——位置越界直接镜像反射或钳制到边界不需要每个变量单独判断。速度上限也要设通常取变量区间宽度的10%~20%防止粒子一步飞出太远。比如归一化后变量范围是0到1速度上限就设0.2左右。4. 目标函数、约束条件与越限处理设计4.1 网损最小化的目标函数建模无功优化的目标函数可以写成网损最小也可以用最小网损加电压偏移最小化做多目标。入门阶段先用单目标——全网有功损耗最小这个目标物理意义清晰也最容易验证算法正确性。全网有功损耗可以用两种方式获得一种是把所有支路的有功损耗求和一种是用所有发电机有功出力之和减去所有负荷有功之和。在Matlab里我用支路功率的代数关系计算function ploss calc_loss(voltage) % 根据节点电压结果计算全网有功损耗 ploss sum(sum(branch_loss_matrix)); end具体实现里如果调用了Matpower的runpf返回的结果结构体里会有branch矩阵第14列或对应的列是支路有功损耗求个和就行。如果自己写潮流那就从节点电压和导纳矩阵出发逐条支路算。4.2 等式约束与不等式约束哪些约束必须处理无功优化的约束分两类等式约束是无功潮流方程本身。每个节点必须满足注入无功等于流出无功这个约束在潮流计算内部已经自动满足了不需要在PSO里额外处理。这也是为什么粒子群算法和无功优化这么好结合——潮流计算天然保证了等式约束。不等式约束才是PSO要操心的包括控制变量的上下限如发电机电压不越限、变压器变比范围状态变量的上下限如PQ节点电压在[0.95,1.05]之间发电机无功出力不超过上下限部分情况下还要盯支路无功潮流上限。需要注意的是控制变量的上下限靠边界处理解决状态变量越限靠惩罚函数解决。这两个思路不能搞混。控制变量的边界你一钳制它就乖了状态变量你不能直接钳制因为你改不了它只能靠目标函数里的惩罚项把粒子往安全区域引导。4.3 惩罚函数的设计经验惩罚函数是最容易翻车的地方。我记得第一次做的时候惩罚系数取得极大觉得越限就该狠狠惩罚结果算法宁可在墙角里蹲着也不敢往边界附近探索收敛结果惨不忍睹。后来才琢磨明白惩罚系数太大惩罚项直接吞掉了网损的梯度信息适应度函数变成一片大平原PSO根本找不到方向惩罚系数太小电压越限的粒子又大摇大摆被认为是最优解。我实用的方案是动态惩罚越限程度越严重惩罚越重越限不大则温和引导penalty lambda * sum(max(0, state_var - upper) .^ 2 ... % 电压越上限部分 max(0, lower - state_var) .^ 2 ... % 电压越下限部分 max(0, qgen - qmax) .^ 2 ... % 发电机无功越上限 max(0, qmin - qgen) .^ 2); % 发电机无功越下限 fitness ploss penalty;lambda建议初期取100量级后期可以视收敛情况加大到1000。注意电压越限平方项里用的是状态变量解码后的实际值而不是归一化值否则越限0.1个标幺值电压和越限0.1个标幺值无功的惩罚量级会失真。我还发现一个小技巧可以把状态变量越限分为越上限和越下限分别统计输出到日志里。调试时一眼就能看出某一代是电压普遍偏高还是无功过载针对性调整约束范围或初始粒子分布。5. 主程序实现PSO与潮流计算的耦合细节5.1 完整的算法执行流程整个主程序我习惯按下面这个顺序组织逻辑清晰也好调试加载IEEE14节点系统数据定义基准容量初始化PSO参数种群规模、迭代次数、学习因子、惯性权重上下限、变量边界初始化粒子群每个粒子的9维位置随机生成在[0,1]区间速度随机初始化为小区间进入迭代循环对每个粒子解码控制变量还原为发电机电压、变压器变比、补偿容量更新系统参数调用潮流计算提取潮流结果计算网损检查状态变量越限情况计算带惩罚的适应度值更新个体最优与全局最优所有粒子评价完毕后更新速度和位置记录本代最优适应度、平均适应度、最优解对应的控制变量迭代结束后输出结果绘制收敛曲线。这里有个细节全局最优gbest的更新是整个算法的心脏它决定收敛方向。我每次更新gbest时都会顺手把对应的控制变量存一份解码后的真实值避免最后要从归一化坐标再解一次码省得出错。5.2 核心matlab代码结构主循环代码骨架如下风格是Matlab原生矩阵化尽量少用巨型for循环套娃% 初始化 nVar 9; % 控制变量维度 popSize 40; % 种群规模 maxIter 80; % 最大迭代次数 x rand(popSize, nVar); % 归一化位置 v 0.1 * randn(popSize, nVar); % 归一化速度 pbest x; pbestFit inf(popSize, 1); gbest x(1, :); gbestFit inf; for iter 1:maxIter w 0.9 - (0.9 - 0.4) * iter / maxIter; % 惯性权重线性递减 for i 1:popSize % 解码并评价 ctrl decode(x(i,:)); % 还原成真实控制变量 [ploss, violation] evaluate(ctrl); % 潮流计算约束检查 fit ploss penalty(violation); % 适应度 % 更新个体最优 if fit pbestFit(i) pbestFit(i) fit; pbest(i, :) x(i, :); end % 更新全局最优 if fit gbestFit gbestFit fit; gbest x(i, :); end end % 更新速度与位置 for i 1:popSize v(i, :) w * v(i, :) ... 2.0 * rand * (pbest(i, :) - x(i, :)) ... 2.0 * rand * (gbest - x(i, :)); v(i, :) bound_velocity(v(i, :), 0.2); x(i, :) x(i, :) v(i, :); x(i, :) bound_position(x(i, :)); % 越界镜像或钳制 end records(iter) gbestFit; end这段代码的evaluate函数是整个程序的绝对核心也是出错率最高的地方。我建议把它独立成一个函数文件传入控制变量解码后的物理值返回网损和越限信息。5.3 评价函数的编写细节评价函数内部按下面几步组织第一步把控制变量写入系统数据。发电机节点电压直接覆盖gen矩阵的电压设定值列变压器变比写入branch矩阵对应变压器支路的变比列补偿容量写入bus数据里的并联导纳字段注意要用实际电纳值不是直接填容量。第二步调用潮流计算。用Matpower的话mpc.gen(:, VG_col) ctrl(1:5); % 发电机端电压 mpc.branch(:, TAP_col) ctrl(6:8); % 变压器变比 mpc.bus(:, BS_col) ctrl(9); % 补偿电纳 result runpf(mpc);第三步从result里提取网损。Matpower返回的result.branch矩阵里有一个损耗列直接求和result.bus里的状态变量用于越限判断。如果潮流不收敛通常返回的result里会有收敛失败的标志此时直接把适应度设为一个很大的值。第四步输出越限信息。电压越限的量、发电机无功越限的量都存下来供惩罚函数调用。我特别提醒一点控制变量解码后要检查潮流是否收敛。PSO初期粒子是随机生成的组合出来的电压和变比很可能离谱导致潮流压根不收敛。如果不加判断程序会在计算网损时报错或返回一个垃圾值被当成优秀解。我的处理是evaluate函数里对result.success 1做判断不收敛直接返回ploss 1e10。这个数值要足够大大到任何正常解都不会超过它保证收敛失败的粒子被自然淘汰。6. 仿真结果与调参排坑实录6.1 收敛曲线与优化前后对比有一次我设置种群40、迭代80次粒子群适应度收敛曲线大致是这样一个规律前15到20代曲线下降特别陡网损从初始的0.138标幺值基准容量100MVA对应的13.8MW快速往下降20到50代进入缓慢下降期算法在小范围内精调50代以后基本趋于平稳再迭代只是微小波动。最终网损优化到了0.126左右下降了约8.7%。优化前后对比如下指标优化前优化后全网有功损耗标幺值约0.138约0.126节点最低电压标幺值0.940.99电压越限节点数3个0个这个结果的意义不只是网损降了百分之几更重要的是把电压拉回了安全范围。多数工程场景里电压质量约束往往比网损降低更刚需算法能同时兼顾这两点就已经有实用价值了。收敛曲线的绘制很简单迭代次数做横轴、最优适应度做纵轴画出来观察形态。如果曲线是一条笔直下降拉到最后说明算法还在探索期迭代次数可以加大如果曲线在低处水平抖动说明收敛到平台期了继续迭代收益有限。6.2 我在调试中踩过的几个坑第一个坑是初始粒子生成太随意导致收敛速度慢。虽然理论说PSO对初值不敏感但实际中如果粒子初始位置全在可行域边缘一堆潮流不收敛的粒子会把gbest引导到奇怪的地方。我的解决办法是先跑一次初始潮流以初始运行点附近的控制变量为基础加上小幅随机扰动生成初始种群。这样做之后第一代粒子几乎全部可行收敛曲线明显变稳。第二个坑是变压器变比的编码方向。IEEE14节点系统的变压器支路在Matpower标准数据里变比设置在特定一侧如果解码后没有按数据约定写入优化结果会显示网损下降但电压分布乱七八糟因为变压器的档位效果完全反了。碰到这种问题建议先固定PSO不迭代只验证一个特定粒子比如全取1.0变比的潮流结果与初始潮流是否一致通了再跑优化。第三个坑是惩罚系数的一刀切问题。不同状态变量越限的严重程度不一样PQ节点电压越限0.05和发电机无功越限0.1对电网安全的危害不是一个量级。可给电压越限和高一级的惩罚权重。我实际用的权重是电压越限乘2、无功越限乘1效果比权重全等要好。第四个坑是结果随机性带来的复现问题。粒子群算法是随机算法每次跑的gbest都可能不同。我在实验记录时会把随机种子固定下来保证结果可复现rng(2024);调参阶段用固定的随机种子对比不同参数配置才公平否则你很难判断网损的改善是参数变好了还是运气变好了。最后一个心得PSO不是万能的不要迷信一次就能找到全局最优。我跑了很多组实验不同参数下最终网损有细微差异这是正常的。真正重要的是算法能不能稳定收敛、约束能不能满足、结果是否可解释。这三点做到了这个课题就算真正吃透了。后续如果想深入还可以把单目标扩展到网损电压偏移多目标优化用MOPSO实现帕累托前沿那将是另一个有意思的方向。
