简介：这份资源面向电力工程专业学生与研究人员，提供基于MATLAB的电力系统常规潮流计算完整实现，算例采用美国西部9节点系统（WSCC9）。内容同时覆盖牛顿迭代法与PQ分解法两条技术路线：牛顿法通过构建雅可比矩阵迭代求解非线性方程组，收敛快但计算量大；PQ分解法将节点划分为PQ与PV节点并分解雅可比矩阵，降低解空间复杂度，更适合WSCC9这类小型系统。压缩包共19个文件，约1.98MB，以13个m脚本和3个mlx实时脚本为核心，辅以pdf说明文档、license与md说明，脚本涵盖主流程、雅可比矩阵求解、功率不平衡量计算、节点电压与相角修正、导纳矩阵生成及批量测试等模块，便于读者对照代码理解算法细节并复现计算结果。目前已有109人学习下载，适合作为电力系统分析课程设计、毕业设计及MATLAB编程练习的参考工具，帮助读者掌握潮流计算建模、收敛判断与结果解读的完整思路。
1. 从 WSCC9 算例说起：为什么潮流计算值得用两套算法对照着做
如果你手头有一个 WSCC9 节点系统的算例压缩包，里面放着线路参数、发电机出力和负荷数据，想用 MATLAB 把潮流算出来，那么真正的问题从来不是「能不能算」，而是「用哪种方法算、算出来的电压幅值和相角能不能信」。WSCC9 是美国西部系统协调委员会给出的一个 9 节点、3 台发电机、3 个负荷的标准测试系统，规模小到可以在笔记本上秒级求解，却又完整保留了平衡节点、PV 节点、PQ 节点的全部特征，是验证潮流程序的经典试金石。牛顿迭代法收敛快、迭代次数少，但初值敏感、雅可比矩阵推导繁琐；PQ 分解法利用高压电网有功与电压相角、无功与电压幅值之间的弱耦合关系，把 2n 维方程拆成两个低维方程交替求解，单次迭代便宜但迭代次数多。把两套方法放在同一个算例上跑，既能互相校验结果，又能直观看到收敛特性的差异。这篇笔记面向的是正在做电力系统课程设计、仿真复现或算法对照的工程师和学生，从数据组织、Y 阵生成、雅可比矩阵构造一路写到收敛判据和结果核对，中间会点出几个我实际调试时翻过车的地方。
2. 牛顿迭代法潮流计算：从节点导纳矩阵到雅可比矩阵的完整链路
2.1 节点分类与功率方程组的建立
潮流计算的核心是一组非线性代数方程。对每个节点，注入功率必须等于网络吸收的功率，写成复数形式就是节点功率方程。WSCC9 系统里，节点 1 通常设为平衡节点（Slack），电压幅值和相角给定，功率待求；节点 2、3 是 PV 节点，电压幅值和有功给定，无功和相角待求；其余节点是 PQ 节点，有功和无功给定，电压幅值和相角待求。这个分类直接决定了未知量的个数和雅可比矩阵的维数。
把节点功率方程按实部虚部展开，得到有功和无功的注入方程。对 PQ 节点，有功和无功都是已知量，方程数为 2；对 PV 节点，有功已知、电压幅值已知，无功和相角待求，方程数为 1。设系统有 n 个节点，其中 PQ 节点 m 个，PV 节点 r 个，那么待求的相角有 n-1 个（平衡节点相角为参考），待求的电压幅值有 m 个。牛顿法把这组方程在当前点做泰勒展开，忽略高阶项，得到修正量方程，雅可比矩阵就是功率方程对相角和电压幅值的偏导数矩阵。
这里有一个容易被忽略的点：雅可比矩阵的元素不是常数，每次迭代都要重新计算。很多初学者为了省事，在第一次迭代后固定雅可比矩阵，结果收敛速度从平方降为线性，迭代次数翻倍。WSCC9 规模小，重算雅可比的开销可以忽略，但在大系统上是否重算就需要权衡。
2.2 用 MATLAB 生成 Y 阵并组装雅可比矩阵
节点导纳矩阵是潮流计算的入口数据。WSCC9 的线路参数通常以 R、X、B 的形式给出，变压器支路还有变比。下面这段代码从支路表生成 Y 阵，支路表按「首节点 末节点 R X B 变比」组织，变比为 0 表示普通线路。
% branch: [from to R X B tap] % tap = 0 表示无变压器，按 1 处理 function Y = buildYbus(branch, nbus) Y = zeros(nbus, nbus); for k = 1:size(branch, 1) f = branch(k, 1); t = branch(k, 2); R = branch(k, 3); X = branch(k, 4); B = branch(k, 5); tap = branch(k, 6); if tap == 0 tap = 1; end Z = R + 1j * X; y = 1 / Z; % 对角元素累加自导纳和充电电容 Y(f, f) = Y(f, f) + y / tap^2 + 1j * B / 2; Y(t, t) = Y(t, t) + y + 1j * B / 2; % 非对角元素为互导纳 Y(f, t) = Y(f, t) - y / tap; Y(t, f) = Y(t, f) - y / tap; end end这段代码的关键在于变比的处理。变压器支路的变比影响自导纳和互导纳，如果变比放在首节点侧，首节点自导纳要除以变比平方，互导纳除以变比。WSCC9 里变压器支路的变比通常是非标准值，抄数据时容易把变比位置搞反，导致 Y 阵不对称。检查方法是看 Y 阵是否满足 Y(f,t) 与 Y(t,f) 在无变比时相等，有变比时按变比关系缩放。
雅可比矩阵的组装是牛顿法最繁琐的部分。下面给出一个按极坐标形式组装的框架，H、N、M、L 四个子块分别对应有功对相角、有功对幅值、无功对相角、无功对幅值的偏导。
% 计算雅可比矩阵，V 为电压幅值向量，theta 为相角向量 function J = buildJacobian(Y, V, theta, P, Q, pq, pv) nbus = length(V); npv = length(pv); npq = length(pq); % 未知量排列：先非平衡节点相角，再 PQ 节点幅值 ntheta = nbus - 1; J = zeros(ntheta + npq, ntheta + npq); for i = 1:nbus for j = 1:nbus if i ~= j % 互导纳相关偏导 Gij = real(Y(i,j)); Bij = imag(Y(i,j)); dtheta = theta(i) - theta(j); % H = dP/dtheta, N = dP/dV H = V(i)*V(j)*(Gij*sin(dtheta) - Bij*cos(dtheta)); N = V(i)*(Gij*cos(dtheta) + Bij*sin(dtheta)); M = -V(i)*V(j)*(Gij*cos(dtheta) + Bij*sin(dtheta)); L = V(i)*(Gij*sin(dtheta) - Bij*cos(dtheta)); else % 对角元素用求和形式，避免重复推导 % 此处省略具体累加，按标准公式填充 end % 按节点类型把 H/N/M/L 填入 J 的对应位置 end end end实际写的时候，对角元素用求和公式更稳妥：H(i,i) = -Q(i) - B(i,i)*V(i)^2，N(i,i) = P(i)/V(i) + G(i,i)*V(i)，M(i,i) = P(i) - G(i,i)*V(i)^2，L(i,i) = Q(i)/V(i) - B(i,i)*V(i)。这些公式在多数教材里都有，但符号容易记混。我的习惯是先用一个 2 节点系统验证雅可比矩阵，确认修正量方向正确后再上 WSCC9。
2.3 迭代流程与收敛判据的参数设置
牛顿法的迭代流程可以概括为：给定初值，计算功率不平衡量，判断是否收敛，不收敛则组装雅可比矩阵、求解修正量、更新电压，重复直到收敛。初值一般取平启动，即所有 PQ 节点电压幅值为 1.0，相角为 0，PV 节点幅值取给定值。WSCC9 用平启动通常 4 到 5 次迭代收敛。
收敛判据用最大功率不平衡量的无穷范数，阈值一般取 1e-6 到 1e-8。阈值太松，结果精度不够；太紧，可能因为浮点误差在最后几次迭代反复震荡。我一般取 1e-8，配合最大迭代次数 20 作为保护。如果 20 次还没收敛，多半是数据有问题或者雅可比矩阵推导错了，继续迭代没有意义。
tol = 1e-8; maxiter = 20; for iter = 1:maxiter % 计算注入功率和不平衡量 [Pcalc, Qcalc] = calcPower(Y, V, theta); dP = Pspec - Pcalc; dQ = Qspec - Qcalc; mismatch = max(abs([dP(2:end); dQ(pq)])); if mismatch < tol fprintf('收敛于第 %d 次迭代，最大不平衡量 %.2e\n', iter, mismatch); break; end % 组装雅可比并求解修正量 J = buildJacobian(Y, V, theta, Pcalc, Qcalc, pq, pv); dx = J \ [dP(2:end); dQ(pq)]; % 更新相角和幅值 theta(2:end) = theta(2:end) + dx(1:nbus-1); V(pq) = V(pq) + dx(nbus:end); end这段代码里，dP 去掉了平衡节点的有功不平衡量，因为平衡节点的有功不受约束；dQ 只取 PQ 节点的无功不平衡量，PV 节点的无功是待求量。修正量的排列顺序必须和雅可比矩阵的行列顺序一致，否则更新时会把相角修正量加到电压上，这种错误在调试时表现为电压越限或者迭代发散。
3. PQ 分解法：把 2n 维方程拆成两个低维问题的工程取舍
3.1 有功-相角与无功-幅值的解耦依据
PQ 分解法的出发点是高压输电网的一个物理事实：线路的电抗远大于电阻，有功功率的传输主要受电压相角影响，无功功率的传输主要受电压幅值影响。把牛顿法雅可比矩阵中的耦合项 N 和 M 忽略，H 子块只与相角有关，L 子块只与幅值有关，于是 2n 维的修正方程拆成两个 n 维方程，交替求解。
这个近似在高压电网（X/R 比大于 3）成立得很好，在低压配电网或者电缆线路居多的场景下会变差。WSCC9 是高压输电网，X/R 比足够大，PQ 分解法适用。如果拿它去算一个低压配电系统，迭代次数可能从十几次涨到几十次甚至不收敛，这时候应该回到牛顿法或者改用其他方法。
解耦之后，H 子块可以进一步写成 B' 矩阵的形式，L 子块写成 B'' 矩阵。B' 和 B'' 在忽略电阻后就是电纳矩阵的虚部，且 B' 不含平衡节点对应的行列，B'' 只含 PQ 节点对应的行列。这两个矩阵在迭代过程中保持不变，可以预先做三角分解，每次迭代只做前代和回代，单次迭代的计算量大幅下降。
3.2 B' 和 B'' 矩阵的构造与因子分解
B' 矩阵的构造规则是：对非平衡节点，取 Y 阵的虚部，去掉平衡节点对应的行和列。B'' 矩阵只取 PQ 节点对应的行和列。下面代码展示从 Y 阵提取 B' 和 B'' 的过程。
% 从 Y 阵提取 B' 和 B'' function [Bp, Bpp] = buildBmatrices(Y, slack, pq) nbus = size(Y, 1); B = imag(Y); % B' 去掉平衡节点 keep = setdiff(1:nbus, slack); Bp = B(keep, keep); % B'' 只保留 PQ 节点 Bpp = B(pq, pq); end提取之后对 B' 和 B'' 做 LU 分解，MATLAB 里直接用 lu 函数。注意 B' 和 B'' 是对称矩阵，但 lu 分解不利用对称性，如果想省内存可以用 LDL 分解。WSCC9 规模小，直接 lu 就够了。
[Lp, Up] = lu(Bp); [Lpp, Upp] = lu(Bpp);迭代时，有功修正量除以电压幅值得到相角修正方程，无功修正量除以电压幅值得到幅值修正方程。具体来说，dP/V 作为右端项，解 B' * dtheta = dP/V；dQ/V 作为右端项，解 B'' * dV = dQ/V。这里的除法是逐元素除，不是矩阵除。
3.3 迭代顺序与收敛速度的实测对比
PQ 分解法的迭代顺序有两种：先解相角再解幅值，或者先解幅值再解相角。多数实现采用先相角后幅值，因为相角修正对有功不平衡更敏感。每次迭代做一次相角修正和一次幅值修正，算作一轮。WSCC9 用 PQ 分解法通常需要 8 到 12 轮收敛，比牛顿法的 4 到 5 次多，但每轮的计算量小。
下面是一个简化的 PQ 分解迭代框架。
for iter = 1:maxiter [Pcalc, Qcalc] = calcPower(Y, V, theta); dP = (Pspec - Pcalc) ./ V; dQ = (Qspec - Qcalc) ./ V; mismatch = max(abs([dP(2:end); dQ(pq)])); if mismatch < tol break; end % 解相角修正 dtheta = Up \ (Lp \ dP(2:end)); theta(2:end) = theta(2:end) + dtheta; % 解幅值修正 dV = Upp \ (Lpp \ dQ(pq)); V(pq) = V(pq) + dV; end实测对比时，我会记录两种方法的迭代次数、单次迭代耗时和总耗时。WSCC9 上牛顿法总耗时通常更短，因为迭代次数少；PQ 分解法在节点数增加到几百上千时优势才明显。这个对比结果本身就是课程设计里值得写进报告的内容。
4. 避坑与排查：WSCC9 潮流计算里最容易翻车的五个地方
4.1 现象：迭代发散，电压幅值越算越大
原因通常是雅可比矩阵符号错误或者修正量更新方向反了。牛顿法里，修正量方程是 J * dx = mismatch，更新时 theta = theta + dx。如果误写成 theta = theta - dx，相当于往错误方向修正，电压会迅速发散。排查方法是打印第一次迭代的 dx，看相角修正量是否在合理范围（通常零点几弧度以内），如果出现几弧度甚至几十弧度的修正量，基本可以确定符号或排列有问题。
4.2 现象：收敛后平衡节点功率与预期不符
平衡节点的有功和无功是潮流结果的一部分，不是给定值。如果算出来的平衡节点功率和参考值差很多，先检查负荷和发电机的有功无功是否填对。WSCC9 的数据里，负荷和发电机的单位可能是 MW 和 MVar，而程序内部用标幺值，忘记除以基准功率会导致功率差 100 倍。另一个常见原因是 PV 节点的电压给定值填错，导致无功出力异常。
4.3 现象：PQ 分解法迭代次数远多于预期
如果 WSCC9 上 PQ 分解法迭代超过 30 次还不收敛，检查 B' 和 B'' 是否包含了不该包含的节点。B' 要去掉平衡节点，B'' 只保留 PQ 节点，这两个规则搞混会导致方程维数和未知量不匹配。另外，如果线路电阻填得过大，解耦假设不成立，迭代次数也会增加。WSCC9 的线路电阻是标准值，一般不会出这个问题，但如果抄数据时把 R 和 X 的位置对调，就会人为制造一个低 X/R 比的系统。
4.4 现象：Y 阵不对称或者对角元素异常
Y 阵应该是对称矩阵（无变比时），如果不对称，检查支路表的首末节点是否填反、变比是否只加在一侧。对角元素异常小，可能是充电电容 B 的符号搞反，容性充电电纳在 Y 阵里是正虚部，如果写成负的，会导致无功分布异常。WSCC9 的线路充电电纳不大，符号错误不会立刻导致发散，但会影响无功和电压结果。
4.5 现象：两种方法结果不一致
牛顿法和 PQ 分解法收敛后，电压幅值和相角的差异应该在收敛阈值量级，比如 1e-6 以内。如果差异达到 1e-3 甚至更大，说明其中一种方法没有真正收敛，或者收敛判据用的不是同一个量。检查方法是把两种方法的最终功率不平衡量都打印出来，确认都小于阈值。如果都小于阈值但结果仍有差异，可能是 PQ 分解法的解耦近似导致的固有偏差，WSCC9 上这个偏差通常很小，但如果线路参数特殊，偏差会放大。
5. 用结果核对与灵敏度分析验证你的潮流程序
写完程序、跑出结果之后，怎么确认它是对的？最直接的办法是用已知的 WSCC9 标准结果对照。WSCC9 的电压幅值和相角在多数电力系统教材和论文里都有给出，把算出来的结果和标准值逐节点对比，误差在 1e-4 以内基本可以认为程序正确。如果没有标准结果，可以用功率平衡来验证：所有发电机有功之和减去所有负荷有功之和，应该等于网络损耗；无功同理。网络损耗可以通过 Y 阵和电压向量算出来，和功率差额对比。
进阶一点的做法是做灵敏度分析。比如把某个负荷的有功增加 10%，看平衡节点的有功变化和关键节点电压的变化，验证是否符合物理直觉。或者把 PV 节点的电压给定值调高，看无功出力的变化方向。这些测试不需要额外数据，改几行输入就能跑，但能暴露很多隐藏的错误。
我自己的习惯是，每改一次雅可比矩阵或者 B' 矩阵的构造逻辑，都先用 2 节点系统跑一遍，确认修正量方向正确，再上 WSCC9。2 节点系统可以手算验证，是最便宜的后悔药。WSCC9 虽然小，但节点类型齐全，一旦出错，排查起来比 2 节点麻烦得多。另一个习惯是把每次迭代的电压幅值和相角打印出来，观察收敛过程是否单调。牛顿法在接近收敛时应该是平方收敛，不平衡量每次迭代下降两个数量级；如果下降很慢或者震荡，说明雅可比矩阵有问题。这些打印语句在调试完成后可以注释掉，但调试阶段不要省。
最后说一个具体技巧：MATLAB 的稀疏矩阵。WSCC9 只有 9 个节点，稠密矩阵足够，但如果想把程序扩展到几十上百个节点，Y 阵和雅可比矩阵都应该用 sparse 存储。sparse 的构造方式是在组装前用 spalloc 预分配，或者用三元组 i、j、v 最后一次性 sparse(i,j,v,n,n)。稠密转稀疏不会自动利用稀疏性，必须从构造阶段就用稀疏格式。这个习惯在 WSCC9 上看不出差别，但到了 IEEE 118 节点系统上，稀疏矩阵能把求解时间从秒级降到毫秒级。希望帮到你。
本文还有配套的精品资源，点击获取
