尧图网站建设 
简介IEEE标准14节点PQ分解法MATLAB程序.m是一份面向电力系统专业学生、研究人员与工程初学者的仿真算法源码用于在14节点标准算例上进行快速潮流计算与稳态分析。程序基于PQ分解思想将潮流方程组拆分为P-θ与Q-V两类子问题通过交替迭代求解节点的电压幅值与相角从而获得支路潮流分布。压缩包共1个文件为可直接运行的MATLAB脚本约5KB代码从节点参数、线路阻抗输入开始依次完成导纳矩阵生成、迭代初值设置、收敛判据判断与结果输出逻辑完整、注释清晰适合对照教材逐段学习或嵌入更复杂的电力系统研究项目。已有879人学习下载。使用该程序可改变节点负荷或发电机出力设置观察电网电压与潮流的相应变化帮助理解PQ分解法的收敛特性和IEEE标准节点系统的结构特点也为后续开展最优潮流、故障分析等拓展研究提供了可复用的计算入口。1. PQ分解法在IEEE 14节点上算什么为什么40年后还要动手写拿到一个IEEE 14节点标准算例第一件事不是打开MATLAB敲代码而是想清楚要算什么。这个系统包含14条母线其中1个平衡节点、4个PV节点发电机节点、9个PQ节点支路里既有变压器又有线路对地电容。工程上在这个系统做潮流计算一般优先选PQ分解法也叫快速解耦法它在牛顿-拉夫逊法的基础上利用输电网有功—相角、无功—电压的弱耦合特性把一个2N阶的修正方程组拆成两个N阶方程迭代时只需对两个常数矩阵各做一次LU分解。速度比牛拉法快2到4倍内存占用也更低。1974年Stott和Alsac提出这个方法后至今仍是能量管理系统EMS里在线潮流的主力算法之一。这篇内容就顺着“14节点数据准备 → B和B矩阵形成 → MATLAB主程序 → 调参与排错 → N-1分析扩展”这条路把程序从零写通。2. 从牛拉法到PQ分解法B、B矩阵怎么形成14节点数据去哪拿2.1 为什么PQ分解法能省这么多计算量极坐标形式的牛顿-拉夫逊法每轮迭代要解一个维度为(2N-1-Npv)的雅可比矩阵方程雅可比矩阵元素依赖当前电压幅值和相角所以每轮都要重新计算、重新分解。网络规模到几千个节点时分解高维稀疏矩阵的耗时占总计算时间的绝大部分。PQ分解法做两个近似这里以BX方案的取法为准它与XB方案的差别在第2.2节说明忽略支路电阻对电纳的影响令r0支路导纳的实部为零只有纯电纳。认为电压幅值变化主要由无功功率决定相角变化主要由有功功率决定于是雅可比矩阵中∂P/∂V和∂Q/∂θ两个子块直接置零。基于这两条近似原本耦合的修正方程变成两个独立方程组[ΔP/V] B * [Δθ] [ΔQ/V] B * [ΔV]B和B都是常数矩阵迭代开始前各做一次LU分解之后每轮只做前代回代。这是PQ分解法速度快的根本原因。14节点系统看不出太大优势到上百节点规模时总耗时可能只相当于牛拉法的三分之一而内存省得更多。2.2 B、B矩阵分别取哪些支路参数这是程序里最容易被写错的环节。B矩阵用于有功—相角修正只由支路电抗构成电阻置零。对角元是所有连接到该节点的1/x_ij之和非对角元是-1/x_ij。不包含对地电容和变压器非标准变比。B矩阵用于无功—电压修正由支路导纳的虚部保持原有b构成考虑变压器和线路充电电容。对角元是Σ(1/x_ij b_ij/2)非对角元是导纳虚部的负值。由于14节点输电线路的X/R比普遍大于5忽略电阻引入的误差非常小整个迭代过程能保持稳定收敛。如果网络含有大量低X/R比支路配电网常见PQ分解法可能发散这一点在第4.3节专门展开。另一个实现细节是在标准算法里B维度是去掉平衡节点后的(N-1)×(N-1)B维度进一步去掉PV节点是NPQ×NPQ。手写代码时频繁做矩阵索引容易出错。我通常的做法是保留全尺寸矩阵然后把这些不参与修正的节点对应的行、列对角元置一个大数比如1e10非对角元置零。这样解出来的对应修正量近似为零效果与缩维一致代码却简单得多。以IEEE 14节点为例从支路数据形成Y矩阵以及B、B的核心代码如下Nbus size(bus, 1); Y zeros(Nbus, Nbus); for k 1:size(branch, 1) fb branch(k, 1); tb branch(k, 2); r branch(k, 3); x branch(k, 4); b branch(k, 5); z r 1i*x; y 1/z; ysp b/2; if branch(k, 6) 0 % 非变压器支路两端各加半条对地电纳 Y(fb, fb) Y(fb, fb) y 1i*ysp; Y(tb, tb) Y(tb, tb) y 1i*ysp; Y(fb, tb) Y(fb, tb) - y; Y(tb, fb) Y(tb, fb) - y; else % 变压器支路变比折算到送端 kt branch(k, 6); Y(fb, fb) Y(fb, fb) y/(kt^2); Y(tb, tb) Y(tb, tb) y; Y(fb, tb) Y(fb, tb) - y/kt; Y(tb, fb) Y(tb, fb) - y/kt; end end这段代码先建立节点导纳矩阵Y。变压器变比的处理是关键标准IEEE数据里变比一般填在送端侧导纳折算到送端时要除以k²互导纳除以k。漏掉这个折算会让B矩阵出现明显偏差潮流算出来的电压分布不对。形成B矩阵时需要把上面循环里的电阻r全部置零、对地电纳b置零然后重新走一遍循环取虚部。这样得到的矩阵纯由1/x构成满足B的定义。2.3 14节点的数据从哪来怎么装成表格IEEE 14节点系统数据有几个常见来源IEEE经典测试系统文档上世纪60年代公开的节点电气参数、MATPOWER工具包自带的case14文件、各类电力系统教材附录。数据量小到可以手工录入但不建议手工输入——14条支路参数弄错一个程序排错要花半天。我一般用MATPOWER做数据源它还包含发电机出力上限、电压上下限等附加信息后续做N-1分析会用到。MATPOWER的case14结构体里有两个关键数组bus和branch各列含义按官方文档统一约定。bus矩阵的type列标识节点类型1为PQ节点2为PV节点3为平衡节点。branch矩阵前六列是送端母线号、受端母线号、电阻r、电抗x、对地电纳b、变压器变比非变压器支路填0。读取方式很简单mpc loadcase(case14); bus mpc.bus; branch mpc.branch; gen mpc.gen;不想依赖MATPOWER时可以把bus和branch直接硬编码成m文件。两种方式我都试过最后选择了独立数据脚本的设计把ieee14_data.m单独放一个文件程序主体不掺数据换IEEE 30节点时只改数据源主函数一行不用动。3. 用MATLAB写PQ分解法主程序从节点数据读入到潮流收敛3.1 主函数框架和初值选择主函数设计成接收bus、branch、gen三个矩阵返回电压幅值、相角、迭代次数和收敛标志。这样后续做批量工况、N-1扫描时换数据文件即可函数体不改。初值用平启动所有PQ节点电压幅值取1.0 p.u.相角取0。以下是一个可以直接运行的完整主程序骨架function [V, theta, iter, flag] pq14_flow(bus, branch, gen, baseMVA, eps, maxIter) % PQ分解法潮流计算 % 输入: % bus : 每行 [母线号, 类型, Pd, Qd, Vm, Va, ...] % branch : 每行 [送端, 受端, r, x, b, 变比, ...] % gen : 每行 [母线号, Pg, Qg, ...] % baseMVA: 基准容量 % eps : 收敛阈值 % maxIter: 最大迭代次数 % 输出: % V : 电压幅值 % theta : 相角(弧度) % iter : 实际迭代次数 % flag : 1收敛 0不收敛 Nbus size(bus, 1); % 节点类型 typeList bus(:, 2); pqIdx find(typeList 1); pvIdx find(typeList 2); slackIdx find(typeList 3); % 净注入功率标幺化 netP zeros(Nbus, 1); netQ zeros(Nbus, 1); for k 1:size(gen, 1) gi find(bus(:,1) gen(k, 1)); netP(gi) gen(k, 2) / baseMVA; netQ(gi) gen(k, 3) / baseMVA; end netP netP - bus(:, 3) / baseMVA; netQ netQ - bus(:, 4) / baseMVA; % 形成Y、Bp、BppBp需要将r和b置零后重新走2.2节循环 Y formY(bus, branch); Bp formBp(bus, branch); Bpp imag(Y); % PV和平衡节点不参与电压修正对角元置大数 for k [pvIdx; slackIdx] Bp(k, :) 0; Bp(:, k) 0; Bp(k, k) 1e10; Bpp(k, :) 0; Bpp(:, k) 0; Bpp(k, k) 1e10; end % Bp、Bpp在迭代前各做一次LU分解 [Lp, Up] lu(Bp); [Lpp, Upp] lu(Bpp); % 平启动初值 Vm ones(Nbus, 1); Va zeros(Nbus, 1); flag 0; for iter 1:maxIter % 计算节点注入功率 Vc Vm .* exp(1i * Va); S Vc .* conj(Y * Vc); Pcal real(S); Qcal imag(S); % 有功不平衡平衡节点不参与 dP (netP - Pcal) ./ Vm; dP(slackIdx) 0; % 解 B * dTheta dP/V dVa Up \ (Lp \ dP); Va Va dVa; % 重新计算用最大不平衡量判断是否进入无功环 Vc Vm .* exp(1i * Va); S Vc .* conj(Y * Vc); Pcal real(S); dP (netP - Pcal) ./ Vm; if max(abs(dP)) eps % 无功不平衡PV和平衡节点不参与 dQ (netQ - Qcal) ./ Vm; dQ(pvIdx) 0; dQ(slackIdx) 0; % 解 B * dV dQ/V dVm Upp \ (Lpp \ dQ); Vm Vm dVm; Vc Vm .* exp(1i * Va); S Vc .* conj(Y * Vc); Qcal imag(S); dQ (netQ - Qcal) ./ Vm; dQ(pvIdx) 0; dQ(slackIdx) 0; if max(abs(dQ)) eps flag 1; break; end end end V Vm; theta Va; end提一个参数取舍eps取1e-4标幺值即最大不平衡量小于0.0001 p.u.14节点系统通常6到9次迭代收敛。maxIter取30足够。代码里dP、dQ都除以了Vm这与方程形式ΔP/V BΔθ完全对应。如果漏掉这个除法收敛判据会随着电压水平漂移电压偏低时可能提前判定收敛电压偏高时又过于严苛。3.2 为什么先分解Bp、Bpp再进入循环这段代码与牛拉法最本质的区别就在这两行LU分解出现在迭代之前而且只执行一次。每轮循环里使用的都是同一个Lp、Up、Lpp、Upp迭代只做前代回代。牛拉法则每轮要重新计算雅可比矩阵矩阵元素随电压变化必须重新分解。所以PQ分解法虽然迭代次数通常比牛拉法多50%到100%但单轮开销只有几分之一综合耗时反而更少。扩展到几百节点时这个优势就从“可以感觉到”变成“一眼可见”。如果想进一步压耗时可以对Bp和Bpp用decomposition对象代替luMATLAB会按矩阵稀疏结构自动选择排序策略。14节点规模差异不明显到1000节点以上时能差出几倍。另一个容易忽略的点是Bp和Bpp中的1e10大数会略微破坏矩阵的条件数求解时可能会多一些舍入误差但相对1e-4的收敛阈值来说影响可以忽略。3.3 计算净注入功率时的符号约定IEEE 14节点系统里负荷消耗功率发电机注入功率。潮流方程S V·conj(Y·V)得到的实部、虚部都是“注入”方向。如果直接把负荷功率填成正数而不叠加发电机出力所有节点净注入都是负值潮流不可能收敛。因此netP和netQ的计算必须是发电机出力标幺化减去负荷功率标幺化。在MATPOWER的数据里mpc.gen的PG、QG列是发电机注入mpc.bus的Pd、Qd列是负荷两者正好一正一负。上面的代码已经按这个约定装配。自己在写数据文件时最容易犯的错是把某台发电机的出力漏掉导致该节点注入功率偏低电压偏低PV节点电压又拉不回来。排查方法很简单看平衡节点出力——如果算出来的平衡节点输出远超合理范围14节点系统正常情况下约在20~80MW区间先查净注入装配。4. PQ分解法的收敛判据与三个必调参数松弛因子、ε、X/R比4.1 阈值ε怎么取取错了会怎样收敛阈值控制的是潮流计算停在哪个精度上。工程实践建议用途eps建议值说明初值计算、系统粗扫1e-2 ~ 1e-3几十毫秒给出可用初值常规潮流分析1e-4校验过载和电压越限足够状态估计/高精度研究1e-6迭代次数增加10%左右不建议更低阈值与迭代次数近似对数关系从1e-4收严到1e-614节点系统迭代次数通常从6次增加到10次左右。反过来阈值放太松会让电压误差超过0.1%判断电压越限时可能误判。我一般固定用1e-4不会为了省几次迭代去动它。还有一个维度容易被忽略dP和dQ是标幺功率和基准容量有关。基准容量取100 MVA时1e-4对应0.01 MW改成1000 MVA时同样1e-4对应0.1 MW等效精度变差。更换baseMVA时务必重新审视这个阈值。4.2 松弛因子α哪些场景值得改PQ分解法迭代公式可以加松弛因子Va_new Va_old α * ΔVa Vm_new Vm_old α * ΔVmα1是标准PQ分解法。重负荷、电压偏低导致收敛慢时把α调到1.2到1.5能加速收敛超过1.6很容易振荡甚至发散。更稳的做法是动态松弛连续两轮最大不平衡量递减时把α增大反弹时减半。不过14节点这类输电网算例固定α1通常就够改反而不安全。对配电网或辐射状网X/R比低收敛性先天不好α1.3常能把振荡迭代拉回正轨。但要注意α只改变迭代路径不改变收敛结果。负荷超过发电机极限、物理上无解时改α不会让它收敛只会让残差停在某个下界附近。4.3 X/R比是PQ分解法收敛的命门PQ分解法成立的前提是支路X/R比足够大工程上一般认为X/R3才实用。IEEE 14节点输电系统的支路X/R比多在5以上所以收敛表现很好。如果把同一套程序直接拿去算某条10kV馈线X/R比经常小于1最常见的故障现象是相角迭代震荡残差降不下去需要远超正常水平的迭代次数才到1e-4极端情况下完全发散但牛拉法算同一网络却收敛正常遇到这类现象第一件事不是加maxIter而是查支路参数里是否混入了低X/R比支路。确认是配电网类型后两个可行方向一是改用牛拉法或保留雅可比矩阵的分块分解二是如果必须用PQ分解法可以考虑对低X/R比支路做串联补偿近似但这会引入额外误差工程上很少用。4.4 结果校核怎么判断程序算对没有收敛后的结果不能直接信先做三个快速校验。第一所有节点注入功率之和加上网络损耗应等于零平衡节点出力应落在合理区间第二电压幅值应全部在0.95到1.06 p.u.之间低于0.9或高于1.1说明数据装配错了第三用MATPOWER的runpf对同一数据跑一遍牛拉法对比两种方法的电压结果误差一般应小于1e-4。如果误差偏大多半是B矩阵漏了某条支路的对地电纳或者变压器变比折算方向写反。5. 扩展用同一套PQ分解法程序做N-1静态安全分析和负荷增长扫描5.1 支路开断模拟与结果记录潮流计算函数稳定后N-1静态安全分析只剩下循环。把branch矩阵中某条支路的状态改为停运重新跑一次潮流检查是否有支路越限、电压越限。停运模拟在MATPOWER数据里可以把该支路的status列置0或者像我这样把变比列强行置0效果等价但要注意别残留导通状态的对地电容。function N1_report run_n1(bus, branch, gen, baseMVA) % 对每条支路做N-1开断扫描 nb size(branch, 1); N1_report zeros(nb, 3); for k 1:nb br_test branch; br_test(k, 6) 0; % 变比置0表示停运注意不要残留b值 [V, ~, ~, flag] pq14_flow(bus, br_test, gen, baseMVA, 1e-4, 30); if flag ~ 1 N1_report(k, :) [-1, -1, -1]; % 不收敛单独标记 continue; end N1_report(k, 1) min(V); % 最低电压 N1_report(k, 2) max(V); % 最高电压 N1_report(k, 3) check_overload(V, br_test, baseMVA); % 最严重负载率 end endN-1分析里有一个容易误读的现象原工况接近极限时开断一条关键支路后潮流不收敛这在物理上代表电压崩溃或失稳不是程序bug。把这批开断单独列出来往往就是系统最薄弱的环节。报告里三项指标齐备后按最低电压排序很快就能定位到最关键的1到2条支路。5.2 负荷增长扫描从14节点换到30节点配电网和输电网规划里最常问的问题是“负荷涨20%后哪个节点先越限”。同一套程序外层循环改一下bus矩阵中的Pd、Qd列即可。比如用0.9、1.0、1.1、1.2四档负荷倍率扫描每档跑一次潮流记录关键母线电压变化画出来就是一条负荷-电压关系曲线。代码上没有任何新难点核心收获是把潮流计算封装成输入输出清晰的小函数后所有分析任务都变成调用它效率提升非常明显。从14节点扩展到IEEE 30节点只需要把数据文件替换成case30。MATPOWER的case30和case14数据结构完全一致读入后直接传给pq14_flow就行。需要注意的一点是case30里变压器支路的变比列格式与14节点一致但发电机节点分布更密集PV节点更多B矩阵中置大数的行数也更多不影响正确性。第一次跑30节点时建议把maxIter临时调到50观察收敛过程确认迭代次数在合理区间后再调回30。这套程序一旦在14节点上调通换数据就能直接用在更大的标准算例上这也是动手写一遍PQ分解法最大的回报。本文还有配套的精品资源点击获取 用MATLAB实现PQ分解法潮流计算：从IEEE 14节点到N-1分析
简介IEEE标准14节点PQ分解法MATLAB程序.m是一份面向电力系统专业学生、研究人员与工程初学者的仿真算法源码用于在14节点标准算例上进行快速潮流计算与稳态分析。程序基于PQ分解思想将潮流方程组拆分为P-θ与Q-V两类子问题通过交替迭代求解节点的电压幅值与相角从而获得支路潮流分布。压缩包共1个文件为可直接运行的MATLAB脚本约5KB代码从节点参数、线路阻抗输入开始依次完成导纳矩阵生成、迭代初值设置、收敛判据判断与结果输出逻辑完整、注释清晰适合对照教材逐段学习或嵌入更复杂的电力系统研究项目。已有879人学习下载。使用该程序可改变节点负荷或发电机出力设置观察电网电压与潮流的相应变化帮助理解PQ分解法的收敛特性和IEEE标准节点系统的结构特点也为后续开展最优潮流、故障分析等拓展研究提供了可复用的计算入口。1. PQ分解法在IEEE 14节点上算什么为什么40年后还要动手写拿到一个IEEE 14节点标准算例第一件事不是打开MATLAB敲代码而是想清楚要算什么。这个系统包含14条母线其中1个平衡节点、4个PV节点发电机节点、9个PQ节点支路里既有变压器又有线路对地电容。工程上在这个系统做潮流计算一般优先选PQ分解法也叫快速解耦法它在牛顿-拉夫逊法的基础上利用输电网有功—相角、无功—电压的弱耦合特性把一个2N阶的修正方程组拆成两个N阶方程迭代时只需对两个常数矩阵各做一次LU分解。速度比牛拉法快2到4倍内存占用也更低。1974年Stott和Alsac提出这个方法后至今仍是能量管理系统EMS里在线潮流的主力算法之一。这篇内容就顺着“14节点数据准备 → B和B矩阵形成 → MATLAB主程序 → 调参与排错 → N-1分析扩展”这条路把程序从零写通。2. 从牛拉法到PQ分解法B、B矩阵怎么形成14节点数据去哪拿2.1 为什么PQ分解法能省这么多计算量极坐标形式的牛顿-拉夫逊法每轮迭代要解一个维度为(2N-1-Npv)的雅可比矩阵方程雅可比矩阵元素依赖当前电压幅值和相角所以每轮都要重新计算、重新分解。网络规模到几千个节点时分解高维稀疏矩阵的耗时占总计算时间的绝大部分。PQ分解法做两个近似这里以BX方案的取法为准它与XB方案的差别在第2.2节说明忽略支路电阻对电纳的影响令r0支路导纳的实部为零只有纯电纳。认为电压幅值变化主要由无功功率决定相角变化主要由有功功率决定于是雅可比矩阵中∂P/∂V和∂Q/∂θ两个子块直接置零。基于这两条近似原本耦合的修正方程变成两个独立方程组[ΔP/V] B * [Δθ] [ΔQ/V] B * [ΔV]B和B都是常数矩阵迭代开始前各做一次LU分解之后每轮只做前代回代。这是PQ分解法速度快的根本原因。14节点系统看不出太大优势到上百节点规模时总耗时可能只相当于牛拉法的三分之一而内存省得更多。2.2 B、B矩阵分别取哪些支路参数这是程序里最容易被写错的环节。B矩阵用于有功—相角修正只由支路电抗构成电阻置零。对角元是所有连接到该节点的1/x_ij之和非对角元是-1/x_ij。不包含对地电容和变压器非标准变比。B矩阵用于无功—电压修正由支路导纳的虚部保持原有b构成考虑变压器和线路充电电容。对角元是Σ(1/x_ij b_ij/2)非对角元是导纳虚部的负值。由于14节点输电线路的X/R比普遍大于5忽略电阻引入的误差非常小整个迭代过程能保持稳定收敛。如果网络含有大量低X/R比支路配电网常见PQ分解法可能发散这一点在第4.3节专门展开。另一个实现细节是在标准算法里B维度是去掉平衡节点后的(N-1)×(N-1)B维度进一步去掉PV节点是NPQ×NPQ。手写代码时频繁做矩阵索引容易出错。我通常的做法是保留全尺寸矩阵然后把这些不参与修正的节点对应的行、列对角元置一个大数比如1e10非对角元置零。这样解出来的对应修正量近似为零效果与缩维一致代码却简单得多。以IEEE 14节点为例从支路数据形成Y矩阵以及B、B的核心代码如下Nbus size(bus, 1); Y zeros(Nbus, Nbus); for k 1:size(branch, 1) fb branch(k, 1); tb branch(k, 2); r branch(k, 3); x branch(k, 4); b branch(k, 5); z r 1i*x; y 1/z; ysp b/2; if branch(k, 6) 0 % 非变压器支路两端各加半条对地电纳 Y(fb, fb) Y(fb, fb) y 1i*ysp; Y(tb, tb) Y(tb, tb) y 1i*ysp; Y(fb, tb) Y(fb, tb) - y; Y(tb, fb) Y(tb, fb) - y; else % 变压器支路变比折算到送端 kt branch(k, 6); Y(fb, fb) Y(fb, fb) y/(kt^2); Y(tb, tb) Y(tb, tb) y; Y(fb, tb) Y(fb, tb) - y/kt; Y(tb, fb) Y(tb, fb) - y/kt; end end这段代码先建立节点导纳矩阵Y。变压器变比的处理是关键标准IEEE数据里变比一般填在送端侧导纳折算到送端时要除以k²互导纳除以k。漏掉这个折算会让B矩阵出现明显偏差潮流算出来的电压分布不对。形成B矩阵时需要把上面循环里的电阻r全部置零、对地电纳b置零然后重新走一遍循环取虚部。这样得到的矩阵纯由1/x构成满足B的定义。2.3 14节点的数据从哪来怎么装成表格IEEE 14节点系统数据有几个常见来源IEEE经典测试系统文档上世纪60年代公开的节点电气参数、MATPOWER工具包自带的case14文件、各类电力系统教材附录。数据量小到可以手工录入但不建议手工输入——14条支路参数弄错一个程序排错要花半天。我一般用MATPOWER做数据源它还包含发电机出力上限、电压上下限等附加信息后续做N-1分析会用到。MATPOWER的case14结构体里有两个关键数组bus和branch各列含义按官方文档统一约定。bus矩阵的type列标识节点类型1为PQ节点2为PV节点3为平衡节点。branch矩阵前六列是送端母线号、受端母线号、电阻r、电抗x、对地电纳b、变压器变比非变压器支路填0。读取方式很简单mpc loadcase(case14); bus mpc.bus; branch mpc.branch; gen mpc.gen;不想依赖MATPOWER时可以把bus和branch直接硬编码成m文件。两种方式我都试过最后选择了独立数据脚本的设计把ieee14_data.m单独放一个文件程序主体不掺数据换IEEE 30节点时只改数据源主函数一行不用动。3. 用MATLAB写PQ分解法主程序从节点数据读入到潮流收敛3.1 主函数框架和初值选择主函数设计成接收bus、branch、gen三个矩阵返回电压幅值、相角、迭代次数和收敛标志。这样后续做批量工况、N-1扫描时换数据文件即可函数体不改。初值用平启动所有PQ节点电压幅值取1.0 p.u.相角取0。以下是一个可以直接运行的完整主程序骨架function [V, theta, iter, flag] pq14_flow(bus, branch, gen, baseMVA, eps, maxIter) % PQ分解法潮流计算 % 输入: % bus : 每行 [母线号, 类型, Pd, Qd, Vm, Va, ...] % branch : 每行 [送端, 受端, r, x, b, 变比, ...] % gen : 每行 [母线号, Pg, Qg, ...] % baseMVA: 基准容量 % eps : 收敛阈值 % maxIter: 最大迭代次数 % 输出: % V : 电压幅值 % theta : 相角(弧度) % iter : 实际迭代次数 % flag : 1收敛 0不收敛 Nbus size(bus, 1); % 节点类型 typeList bus(:, 2); pqIdx find(typeList 1); pvIdx find(typeList 2); slackIdx find(typeList 3); % 净注入功率标幺化 netP zeros(Nbus, 1); netQ zeros(Nbus, 1); for k 1:size(gen, 1) gi find(bus(:,1) gen(k, 1)); netP(gi) gen(k, 2) / baseMVA; netQ(gi) gen(k, 3) / baseMVA; end netP netP - bus(:, 3) / baseMVA; netQ netQ - bus(:, 4) / baseMVA; % 形成Y、Bp、BppBp需要将r和b置零后重新走2.2节循环 Y formY(bus, branch); Bp formBp(bus, branch); Bpp imag(Y); % PV和平衡节点不参与电压修正对角元置大数 for k [pvIdx; slackIdx] Bp(k, :) 0; Bp(:, k) 0; Bp(k, k) 1e10; Bpp(k, :) 0; Bpp(:, k) 0; Bpp(k, k) 1e10; end % Bp、Bpp在迭代前各做一次LU分解 [Lp, Up] lu(Bp); [Lpp, Upp] lu(Bpp); % 平启动初值 Vm ones(Nbus, 1); Va zeros(Nbus, 1); flag 0; for iter 1:maxIter % 计算节点注入功率 Vc Vm .* exp(1i * Va); S Vc .* conj(Y * Vc); Pcal real(S); Qcal imag(S); % 有功不平衡平衡节点不参与 dP (netP - Pcal) ./ Vm; dP(slackIdx) 0; % 解 B * dTheta dP/V dVa Up \ (Lp \ dP); Va Va dVa; % 重新计算用最大不平衡量判断是否进入无功环 Vc Vm .* exp(1i * Va); S Vc .* conj(Y * Vc); Pcal real(S); dP (netP - Pcal) ./ Vm; if max(abs(dP)) eps % 无功不平衡PV和平衡节点不参与 dQ (netQ - Qcal) ./ Vm; dQ(pvIdx) 0; dQ(slackIdx) 0; % 解 B * dV dQ/V dVm Upp \ (Lpp \ dQ); Vm Vm dVm; Vc Vm .* exp(1i * Va); S Vc .* conj(Y * Vc); Qcal imag(S); dQ (netQ - Qcal) ./ Vm; dQ(pvIdx) 0; dQ(slackIdx) 0; if max(abs(dQ)) eps flag 1; break; end end end V Vm; theta Va; end提一个参数取舍eps取1e-4标幺值即最大不平衡量小于0.0001 p.u.14节点系统通常6到9次迭代收敛。maxIter取30足够。代码里dP、dQ都除以了Vm这与方程形式ΔP/V BΔθ完全对应。如果漏掉这个除法收敛判据会随着电压水平漂移电压偏低时可能提前判定收敛电压偏高时又过于严苛。3.2 为什么先分解Bp、Bpp再进入循环这段代码与牛拉法最本质的区别就在这两行LU分解出现在迭代之前而且只执行一次。每轮循环里使用的都是同一个Lp、Up、Lpp、Upp迭代只做前代回代。牛拉法则每轮要重新计算雅可比矩阵矩阵元素随电压变化必须重新分解。所以PQ分解法虽然迭代次数通常比牛拉法多50%到100%但单轮开销只有几分之一综合耗时反而更少。扩展到几百节点时这个优势就从“可以感觉到”变成“一眼可见”。如果想进一步压耗时可以对Bp和Bpp用decomposition对象代替luMATLAB会按矩阵稀疏结构自动选择排序策略。14节点规模差异不明显到1000节点以上时能差出几倍。另一个容易忽略的点是Bp和Bpp中的1e10大数会略微破坏矩阵的条件数求解时可能会多一些舍入误差但相对1e-4的收敛阈值来说影响可以忽略。3.3 计算净注入功率时的符号约定IEEE 14节点系统里负荷消耗功率发电机注入功率。潮流方程S V·conj(Y·V)得到的实部、虚部都是“注入”方向。如果直接把负荷功率填成正数而不叠加发电机出力所有节点净注入都是负值潮流不可能收敛。因此netP和netQ的计算必须是发电机出力标幺化减去负荷功率标幺化。在MATPOWER的数据里mpc.gen的PG、QG列是发电机注入mpc.bus的Pd、Qd列是负荷两者正好一正一负。上面的代码已经按这个约定装配。自己在写数据文件时最容易犯的错是把某台发电机的出力漏掉导致该节点注入功率偏低电压偏低PV节点电压又拉不回来。排查方法很简单看平衡节点出力——如果算出来的平衡节点输出远超合理范围14节点系统正常情况下约在20~80MW区间先查净注入装配。4. PQ分解法的收敛判据与三个必调参数松弛因子、ε、X/R比4.1 阈值ε怎么取取错了会怎样收敛阈值控制的是潮流计算停在哪个精度上。工程实践建议用途eps建议值说明初值计算、系统粗扫1e-2 ~ 1e-3几十毫秒给出可用初值常规潮流分析1e-4校验过载和电压越限足够状态估计/高精度研究1e-6迭代次数增加10%左右不建议更低阈值与迭代次数近似对数关系从1e-4收严到1e-614节点系统迭代次数通常从6次增加到10次左右。反过来阈值放太松会让电压误差超过0.1%判断电压越限时可能误判。我一般固定用1e-4不会为了省几次迭代去动它。还有一个维度容易被忽略dP和dQ是标幺功率和基准容量有关。基准容量取100 MVA时1e-4对应0.01 MW改成1000 MVA时同样1e-4对应0.1 MW等效精度变差。更换baseMVA时务必重新审视这个阈值。4.2 松弛因子α哪些场景值得改PQ分解法迭代公式可以加松弛因子Va_new Va_old α * ΔVa Vm_new Vm_old α * ΔVmα1是标准PQ分解法。重负荷、电压偏低导致收敛慢时把α调到1.2到1.5能加速收敛超过1.6很容易振荡甚至发散。更稳的做法是动态松弛连续两轮最大不平衡量递减时把α增大反弹时减半。不过14节点这类输电网算例固定α1通常就够改反而不安全。对配电网或辐射状网X/R比低收敛性先天不好α1.3常能把振荡迭代拉回正轨。但要注意α只改变迭代路径不改变收敛结果。负荷超过发电机极限、物理上无解时改α不会让它收敛只会让残差停在某个下界附近。4.3 X/R比是PQ分解法收敛的命门PQ分解法成立的前提是支路X/R比足够大工程上一般认为X/R3才实用。IEEE 14节点输电系统的支路X/R比多在5以上所以收敛表现很好。如果把同一套程序直接拿去算某条10kV馈线X/R比经常小于1最常见的故障现象是相角迭代震荡残差降不下去需要远超正常水平的迭代次数才到1e-4极端情况下完全发散但牛拉法算同一网络却收敛正常遇到这类现象第一件事不是加maxIter而是查支路参数里是否混入了低X/R比支路。确认是配电网类型后两个可行方向一是改用牛拉法或保留雅可比矩阵的分块分解二是如果必须用PQ分解法可以考虑对低X/R比支路做串联补偿近似但这会引入额外误差工程上很少用。4.4 结果校核怎么判断程序算对没有收敛后的结果不能直接信先做三个快速校验。第一所有节点注入功率之和加上网络损耗应等于零平衡节点出力应落在合理区间第二电压幅值应全部在0.95到1.06 p.u.之间低于0.9或高于1.1说明数据装配错了第三用MATPOWER的runpf对同一数据跑一遍牛拉法对比两种方法的电压结果误差一般应小于1e-4。如果误差偏大多半是B矩阵漏了某条支路的对地电纳或者变压器变比折算方向写反。5. 扩展用同一套PQ分解法程序做N-1静态安全分析和负荷增长扫描5.1 支路开断模拟与结果记录潮流计算函数稳定后N-1静态安全分析只剩下循环。把branch矩阵中某条支路的状态改为停运重新跑一次潮流检查是否有支路越限、电压越限。停运模拟在MATPOWER数据里可以把该支路的status列置0或者像我这样把变比列强行置0效果等价但要注意别残留导通状态的对地电容。function N1_report run_n1(bus, branch, gen, baseMVA) % 对每条支路做N-1开断扫描 nb size(branch, 1); N1_report zeros(nb, 3); for k 1:nb br_test branch; br_test(k, 6) 0; % 变比置0表示停运注意不要残留b值 [V, ~, ~, flag] pq14_flow(bus, br_test, gen, baseMVA, 1e-4, 30); if flag ~ 1 N1_report(k, :) [-1, -1, -1]; % 不收敛单独标记 continue; end N1_report(k, 1) min(V); % 最低电压 N1_report(k, 2) max(V); % 最高电压 N1_report(k, 3) check_overload(V, br_test, baseMVA); % 最严重负载率 end endN-1分析里有一个容易误读的现象原工况接近极限时开断一条关键支路后潮流不收敛这在物理上代表电压崩溃或失稳不是程序bug。把这批开断单独列出来往往就是系统最薄弱的环节。报告里三项指标齐备后按最低电压排序很快就能定位到最关键的1到2条支路。5.2 负荷增长扫描从14节点换到30节点配电网和输电网规划里最常问的问题是“负荷涨20%后哪个节点先越限”。同一套程序外层循环改一下bus矩阵中的Pd、Qd列即可。比如用0.9、1.0、1.1、1.2四档负荷倍率扫描每档跑一次潮流记录关键母线电压变化画出来就是一条负荷-电压关系曲线。代码上没有任何新难点核心收获是把潮流计算封装成输入输出清晰的小函数后所有分析任务都变成调用它效率提升非常明显。从14节点扩展到IEEE 30节点只需要把数据文件替换成case30。MATPOWER的case30和case14数据结构完全一致读入后直接传给pq14_flow就行。需要注意的一点是case30里变压器支路的变比列格式与14节点一致但发电机节点分布更密集PV节点更多B矩阵中置大数的行数也更多不影响正确性。第一次跑30节点时建议把maxIter临时调到50观察收敛过程确认迭代次数在合理区间后再调回30。这套程序一旦在14节点上调通换数据就能直接用在更大的标准算例上这也是动手写一遍PQ分解法最大的回报。本文还有配套的精品资源点击获取