简介：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节说明）：
- 忽略支路电阻对电纳的影响，令r=0，支路导纳的实部为零，只有纯电纳。
- 认为电压幅值变化主要由无功功率决定，相角变化主要由有功功率决定，于是雅可比矩阵中∂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、Bpp（Bp需要将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对象代替lu，MATLAB会按矩阵稀疏结构自动选择排序策略。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-6，14节点系统迭代次数通常从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/R>3才实用。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分解法最大的回报。
本文还有配套的精品资源，点击获取
