MATLAB实现PQ分解法潮流计算：从原理到源码解析

简介这是一份基于PQ分解法的电力系统潮流计算MATLAB源码面向电力工程专业学生、科研人员以及需要掌握数值计算技巧的MATLAB学习者。程序通过区分PQ节点与PV节点结合基尔霍夫定律构建代数方程组并采用牛顿-拉弗森迭代进行求解最终输出各节点电压、功率等关键参数。核心文件PQ.m完整覆盖了网络拓扑定义、节点分类、方程组构建、迭代求解与结果输出等环节共1个m文件压缩包仅2KB结构精炼适合逐行研读与二次开发。目前已有160人学习下载。读者借助该源码可深入理解潮流计算原理掌握MATLAB矩阵运算、迭代收敛判断及程序调试方法也可在此基础上扩展实现牛顿法、高斯-塞德尔法等其他算法对提升电力系统分析与MATLAB编程综合能力有明显帮助。1. 为什么快速潮流计算都绕不开PQ分解法电力系统调度员在做N-1校核或在线安全分析时经常要在一两秒内跑完上千节点的潮流。完整牛顿法每次迭代都要更新和分解雅可比矩阵计算量随节点数近似三阶增长实际在线场景往往扛不住。PQ分解法正是针对这个问题出现的它利用有功功率主要取决于相角、无功功率主要取决于电压幅值这一物理特性把二维修正方程拆成两个一维方程而且系数矩阵只要在迭代前算一次。这份matlab源码下载里只包含一个PQ.m文件却完整覆盖了节点分类、导纳矩阵、BB矩阵、迭代收敛和结果输出全链路。对电力专业学生、MATLAB编程初学者以及需要快速理解潮流算法实现的工程师都是一个很适合拆解的样本。2. PQ分解法如何从牛顿法里“拆”出效率2.1 从极坐标牛顿法到常数系数矩阵潮流计算的本质是求解一组非线性功率方程。对任意节点i注入功率的极坐标形式为Pi Ui * sum_{j} Uj * (Gij cosθij Bij sinθij) Qi Ui * sum_{j} Uj * (Gij sinθij - Bij cosθij)这里θij θi - θj。写成牛顿法的分块修正方程[ΔP; ΔQ] [H N; J L] [Δθ; ΔU/U]H、N、J、L是雅可比矩阵的四个分块。PQ分解法的第一步是忽略N和J。为什么能这样忽略因为高压输电网的支路电抗远大于电阻正常运行节点间的相角差一般在10度以内因此Gij sinθij很小而Bij cosθij占主导于是ΔP几乎只对Δθ敏感ΔQ几乎只对ΔU敏感。第二步假设各节点电压都接近1.0 p.u.那么H和L矩阵可以近似为常数只用节点导纳矩阵的虚部就能构成。经过标幺值换算后得到两个线性修正方程B Δθ ΔP / U B ΔU ΔQ / U这才是PQ分解法真正高效的原因雅可比矩阵被两个常系数矩阵B和B取代迭代过程中只需要做LU分解一次之后每次迭代都是两次前代回代而不是每轮都重做矩阵分解。对于大型稀疏网络节省的开销往往能达到一个数量级。这里需要说明上面公式里的U是当前迭代电压幅值向量ΔP/U和ΔQ/U是代数量源码中通常用“dP./V”这样逐元素除来实现。这里需要提醒的是上述推导建立在rx的假设上。对电缆线路或配电网电阻较大PQ分解法的收敛性会变差这时更稳妥的做法是使用保留N和J部分的快速解耦法或者直接回到牛顿法。MATLAB源码里如果只按经典公式实现拿到配电网算例就很容易不收敛这不是程序bug而是算法适用边界的问题。使用者要在写代码前判断自己的网络是否满足高压输电这一前提。2.2 节点分类PQ、PV与平衡节点在实现PQ分解法之前必须把网络节点分类清楚因为B和B的维度由节点类型决定。表格如下。节点类型已知量未知量典型设备PQ节点有功P、无功Q电压幅值V、相角θ负荷、电容器、电抗器PV节点有功P、电压幅值V无功Q、相角θ发电机、调相机平衡节点电压幅值V、相角θ有功P、无功Q等值大电网在MATLAB程序里节点类型通常用数字1、2、3表示。第一次阅读PQ.m时先看清楚类型定义然后找到下面这类索引语句pq find(bus(:,2) 1); pv find(bus(:,2) 2); slack find(bus(:,2) 3);这三行代码决定了后续B、B怎么切分也决定了修正方程组的节点顺序。如果节点编号不连续或bus矩阵没有按编号排序用find取出索引后后续所有运算都要保持引用一致否则会出现“索引错位”问题这是排错过程中最隐蔽的一个点。另外还要注意平衡节点一般不参与B和B的任何一维求解但并不代表它不需要功率方程。平衡节点的电压和相角固定它的作用是在每次迭代后吸收全网的功率不平衡量所以程序里必须有单独一段代码计算它的P和Q。2.3 为什么用MATLAB实现最有性价比PQ分解法大量涉及矩阵行/列筛选、稀疏矩阵分解和迭代向量更新MATLAB对这类操作提供了非常自然的语法。比如Bp \ dP就能完成求解不用手写高斯消去稀疏矩阵的sparse()函数可以大幅降低内存占用让千节点计算瞬发完成。更重要的是MATLAB工作区可以直接观察迭代中每一个矩阵方便在课堂和实验室复原算法细节。相对而言用C/C写同样功能需要额外维护稀疏矩阵库用Python则要谨慎处理numpy的细粒度内存复制容易因为深拷贝导致性能劣化。如果只是研究PQ分解法本身MATLAB这份源码的难度曲线最平缓。之后想迁移到生产环境再按相同逻辑改写成C或Python都不迟。3. PQ.m源码结构与核心参数设置3.1 数据输入节点矩阵与支路矩阵运行一个潮流程序最先面对的问题是“怎么把网络描述清楚”。PQ.m常见做法是定义两个矩阵bus和branch。bus的每行描述一个节点列含义为编号、类型、电压幅值初值、相角初值、负荷有功、负荷无功、发电机有功、发电机无功branch的每行描述一条支路列含义为首端节点、末端节点、电阻、电抗、充电电纳。为了便于理解我通常直接在脚本开头给出如下示例% bus: [节点 类型 幅值 相角 负荷P 负荷Q 发电P 发电Q] bus [ 1 3 1.05 0 0 0 0 0; 2 2 1.00 0 0 0 1.2 0; 3 1 1.00 0 0.8 0.4 0 0; ]; % branch: [首端 末端 电阻R 电抗X 充电电纳B] branch [ 1 2 0.02 0.10 0.02; 1 3 0.04 0.20 0.02; 2 3 0.03 0.15 0.02; ];这里所有电气量都采用标幺值p.u.。如果原始数据是实际值需要先除以基准容量。例如算例基准容量S_base100 MVA那么80 MW负荷写为0.8。初值一般让PQ节点电压幅值为1.0PV节点用发电机给定的电压幅值相角统一设为0这样不至于让迭代一开始就偏离解太远。branch里的充电电纳B一般指总充电电纳程序中要除以2再分别加到两端节点上你阅读源码时注意看它有没有做这一步。3.2 从导纳矩阵生成B和B构建节点导纳矩阵Y是所有潮流计算的第一步。普通线路采用π型等值串联电纳和充电电纳都要累加。下面这段代码是PQ.m中常见做法Y complex(zeros(nb, nb)); for k 1:size(branch,1) i branch(k,1); j branch(k,2); r branch(k,3); x branch(k,4); b 2*branch(k,5); ys 1/(r 1i*x); Y(i,i) Y(i,i) ys 1i*b; Y(i,j) Y(i,j) - ys; Y(j,i) Y(j,i) - ys; Y(j,j) Y(j,j) ys 1i*b; end B imag(Y);注意这里的b 2*branch(k,5)因为输入的是总充电电纳需要平均分到线路两端。随后从B生成B和Bpq find(bus(:,2) 1); pv find(bus(:,2) 2); slack find(bus(:,2) 3); Bp B([pq; pv], [pq; pv]); % 相角修正矩阵 Bpp B(pq, pq); % 电压修正矩阵B包含PQ和PV两种节点因为这两类节点的相角都是未知的B只包含PQ节点因为PV节点的电压幅值固定不需要对它求解。很多初学者在这里踩坑把PV节点也放进了B导致矩阵维度变大、迭代结果震荡甚至直接奇异。判断维度是否正确的办法很简单B的行数必须等于PQ节点个数B的行数必须等于PQ节点数加PV节点数。3.3 迭代求解修正量如何更新状态核心迭代部分可以写成下面的结构。先把当前电压V和相角theta代入功率方程计算节点注入功率的实部和虚部然后计算有功不平衡量ΔP和无功不平衡量ΔQ。for iter 1:iter_max [Pcal, Qcal] calc_power(bus, Y, theta, V); dP (bus(:,7) - bus(:,5)) - Pcal; dQ (bus(:,8) - bus(:,6)) - Qcal; dTheta(valid) Bp \ (dP(valid) ./ V(valid)); theta(valid) theta(valid) dTheta(valid); dV(pq) Bpp \ (dQ(pq) ./ V(pq)); V(pq) V(pq) dV(pq); if max(abs(dTheta)) tol max(abs(dV)) tol break; end end上面的valid表示PQPV节点的索引pq表示PQ节点索引。MATLAB的反斜杠运算符会为稠密矩阵做LU分解如果是稀疏矩阵则自动切换到稀疏分解路径这里不需要手动指定。每次迭代只解两个低阶方程组所以即使迭代次数比牛顿法多总时间仍然占优势。值得关注的是PV节点不参与dQ的修正所以Bpp只针对PQ节点计算dP时也不能让PV节点的电压不定而要在每次迭代后把PV节点的V恢复为初值。3.4 关键参数表和调整建议参数默认值调整建议收敛精度tol1e-6教学场景可放宽到1e-4在线分析建议1e-6最大迭代次数iter_max30超过20次不收敛先检查BB组装基准容量S_base100 MVA与全网数据保持一致否则功率错一个数量级电压初值1.0 p.u.PQ节点固定1.0PV节点用发电机给定电压相角初值0重载网络可改为平启动后再用牛顿法预热这些参数直接影响收敛速度。收敛精度从1e-6改成1e-4迭代次数通常会减少两三次但结果精度也下降如果只是为了课程实验够用即可。最大迭代次数设成30已经覆盖绝大多数正常情况如果30步还不收敛说明问题不在参数而在建模。另外如果发现迭代前期正常、后期反复震荡可以适当减小电压阻尼因子比如把dV乘以0.8这属于很实用的小技巧。4. 运行PQ.m的完整流程与收敛判据4.1 准备一个可复现的三节点算例为了验证源码是否正常我建议先用一个最小系统跑通再换自己的大网络。三节点系统包含一个平衡节点、一个发电机节点和一个负荷节点节点1作为平衡节点2作为PV节点3作为PQ。输入矩阵如下S_base 100; % 基准容量单位MVA % 编号 类型 幅值 相角 负荷P 负荷Q 发电P 发电Q bus [ 1 3 1.05 0 0.00 0.00 0.00 0.00; 2 2 1.00 0 0.00 0.00 0.50 0.00; 3 1 1.00 0 0.70 0.30 0.00 0.00; ]; branch [ 1 2 0.02 0.10 0.02; 1 3 0.04 0.20 0.02; 2 3 0.03 0.15 0.02; ];这里负荷为70 MW 30 Mvar发电机出力0.5 p.u.不足部分由平衡节点补足。把这段代码复制到PQ.m的开头运行后观察循环迭代次数。如果一切正常应该看到迭代次数在3到6次之间每个节点电压幅值都在0.95到1.1之间。如果超出这个范围优先检查支路电抗和充电电纳是否写反。4.2 看懂输出结果源码运行结束后一般会在命令行打印节点电压、相角、发电机出力和负荷功率。常见的输出格式是迭代次数4 节点 类型 电压幅值(p.u.) 相角(deg) 发电有功 发电无功 1 3 1.0500 0.0000 0.5934 0.1202 2 2 1.0000 -0.4821 0.5000 -0.0284 3 1 0.9824 -1.2510 0.0000 0.0000节点1的发电有功0.5934意味着平衡节点掏出了约59.34 MW的有功节点2的发电无功是负值说明它实际在吸收无功。这些都是正常现象。重点看节点3的电压幅值是否太低如果低于0.95说明负荷过重或无功补偿不足。相角差如果超过15度也要回头检查支路参数。4.3 功率平衡校验结果对了才算跑通潮流计算最后必须满足节点功率平衡。校验公式是总发电总负荷总网损。可以用下面这段代码在脚本末尾检查total_P_load sum(bus(:,5)); total_Q_load sum(bus(:,6)); total_P_gen sum(bus(:,7)) P_slack; total_Q_gen sum(bus(:,8)) Q_slack; P_loss total_P_gen - total_P_load; Q_loss total_Q_gen - total_Q_load; fprintf(网损: %.6f p.u., 无功损耗: %.6f p.u.\n, P_loss, Q_loss);P_loss应该略大于0通常在0.001到0.02之间。如果P_loss是负的说明回路中某条支路方向定义反了如果P_loss为0系统里可能有零阻抗支路物理上不合理。这里P_slack和Q_slack可以从源码内部变量里直接取不同版本命名不一样可以搜索slack找到。4.4 常见运行错误与排查表错误现象可能原因排查方法索引超出矩阵维度branch里的节点编号大于bus行数用max(max(branch(:,1:2)))检查B奇异或行列式为0B包含了PV节点确认Bpp的行列数等于PQ节点个数迭代发散出现NaN支路电抗x为0或初值离解太远检查branch第4列是否为0降低负荷相角结果反复震荡bus编号不连续或重复对bus按编号排序并检查重复项结果与MATPOWER差距大充电电纳是否除以2检查branch第5列处理方式遇到NaN时不要先怀疑算法先用dbstop if nan设置断点停在第一个NaN位置然后查看迭代次数和当前变量。绝大多数NaN来自dP或dQ除以零也就是某个节点电压幅值在迭代中变成0。这种情况往往是因为B或B组装错误导致更新步长过大。5. PQ.m的进阶动作节点越限转换与结果导出5.1 让PV节点越限后自动转成PQ节点在实际电网里发电机无功出力有上下限。如果PQ.m计算出某台发电机的无功超过限额说明当前运行点不可行需要把这个节点从PV改成PQ并令其无功固定为限额值。实现方式是在主循环外面再套一层外层循环Qgen calc_Q_from_V(bus, Y, theta, V); % 由结果反推无功 pv find(bus(:,2) 2); for k pv if Qgen(k) bus(k,9) || Qgen(k) bus(k,10) bus(k,2) 1; % PV转PQ bus(k,8) clamp(Qgen(k), bus(k,10), bus(k,9)); fprintf(节点%d转成PQ节点\n, bus(k,1)); end end注意这里bus矩阵需要增加两列存放Q上限和Q下限。转成PQ节点后B的维度会变化因为该节点的电压幅值从已知变成了未知。因此每次更新节点类型后要重新调用Bpp B(pq, pq)不能沿用旧矩阵。如果忽略这一步源码会报维度错误或者干脆用错误维度算出错误结果。5.2 把节点结果导出成带表头的CSV通过源码计算完之后如果想交给其他工具做报告用writetable是最省事的T table(bus(:,1), bus(:,2), V, theta_deg, Pgen, Qgen, ... VariableNames, {NodeID,Type,Voltage,Angle_deg,P_gen,Q_gen}); writetable(T, pq_result.csv);这个CSV文件可以用Excel打开也可以用Python的pandas继续处理。注意输出给外部系统时最好把角度从弧度换算成度并保留6位小数避免下游去重烦恼。5.3 用MATPOWER结果交叉验证我拿到任何修改过的PQ.m都会用MATPOWER的runpf跑同一个系统做背靠背对比。具体步骤是用loadcase(case3.m)载入三节点算例调用runpf得到节点电压然后和PQ.m结果逐项相减观察最大绝对误差。误差小于1e-5说明B、B组装正确迭代顺序也没有错误。这样做的好处是能把矩阵索引、节点类型转换这些细节问题快速暴露出来而不需要自己手动算一遍潮流方程。本文还有配套的精品资源点击获取