潮流计算在电力系统里属于那种“看起来简单、写起来全是细节”的东西。很多教材把公式推导梳理得很漂亮，但一到 MATLAB 里自己动手，就会遇到雅可比矩阵符号搞混、迭代发散、P-Q 分解法在某个算例里死活不收的尴尬。我当初就是因为不满足于直接调工具箱，把所有经典算法都亲手写了一遍，这篇就基于我当时摸索的经验，完整拆解电力系统潮流计算中的牛顿-拉夫逊法和 P-Q 分解法，从原理推导到可运行的 MATLAB 代码，再配合一个三节点算例看两者的收敛行为，最后整理了我在实操里踩过的几个大坑。无论是电气专业的学生做课程设计，还是刚接触电网分析的工程师想搞明白算法底层逻辑，这篇文章都值得你花几分钟完整看完。
1. 潮流计算要解的方程组：从物理问题到数学问题的翻译过程
1.1 潮流问题到底在求什么
我在带新人的时候，经常有人问：潮流计算不就是解个电路方程吗？对，但又不完全对。普通电路分析给定了电源和阻抗，直接解线性方程组就能得到电流和电压；但电力系统里，大部分节点的注入功率是给定的，而不是电流给定，这就让问题变成了非线性方程组。更准确地说，潮流计算的目标是：在已知部分节点注入功率、部分节点电压幅值的情况下，求解全网各节点的电压幅值和相角，并由此推算线路功率、网损和变压器分接头是否合理。
实际工程中，节点通常分成三类。平衡节点承担全网功率差额，电压幅值和相角都是给定的；PQ 节点给定有功和无功注入，电压幅值和相角待求，大多数负荷节点和普通发电机节点都属于这一类；PV 节点给定有功和电压幅值，无功注入和相角待求，典型的调频调压发电机节点就是这种。这三类节点的存在直接决定了方程组的组成方式，也是后面构建雅可比矩阵的依据。
1.2 极坐标下的功率方程
我们常说的节点功率方程，在极坐标下写作：
P_i = V_i Σ_j V_j (G_ij cosθ_ij + B_ij sinθ_ij) Q_i = V_i Σ_j V_j (G_ij sinθ_ij - B_ij cosθ_ij)
其中 G 和 B 分别是节点导纳矩阵的实部和虚部，θ_ij 是节点 i 与节点 j 的相角差。注意这个求和是对所有节点进行的，包括 i 节点本身。很多初学者会把导纳矩阵的自导纳项漏掉，导致功率怎么算都对不上。我写代码时习惯先把 Ybus 构造好，再用这个公式直接算 P、Q，然后和给定的注入功率求差值，这个差值就是迭代要消灭的残差。
需要特别提醒一点：节点导纳矩阵的自导纳 Y_ii 是所有与节点 i 相连支路导纳之和，互导纳 Y_ij 是支路导纳的负值。线路参数是复数阻抗 z = r + jx，支路导纳 y = 1/z，所以互导纳虚部通常是正的，自导纳虚部通常是负的。这个符号如果不理顺，后面无论是手推雅可比还是写代码都会被带偏。
1.3 残差方程组的规模
假设系统一共有 n 个节点，平衡节点有 1 个，PV 节点有若干个。那么待求的未知量一共是：(n-1) 个相角（平衡节点相角固定为 0）加上 PQ 节点个数个电压幅值。方程数量也要能对齐：每个非平衡节点都有一个有功残差方程，每个 PQ 节点还有一个无功残差方程。
那么问题就成了求 F(x) = 0，其中 x 是未知相角和电压幅值的向量，F 是刚才算出来的残差向量。非线性方程组的求解方法自然就引出了牛顿-拉夫逊法。
2. 牛顿-拉夫逊法：从功率方程到雅可比矩阵的完整推导
2.1 牛顿法的基本思想：非线性问题线性化迭代
牛顿法的核心逻辑一句话就能说清：在当前的运行点附近，用一阶泰勒展开把非线性方程近似成线性方程，求出修正量，然后不停重复，直到残差足够小为止。你把 F(x) 在当前点 x* 展开，得到 F(x* + Δx) ≈ F(x*) + J Δx，强行令右边等于 0，就能解出修正量 Δx = -J⁻¹ F(x*)，然后让 x 沿着这个方向走一步。
理解这一点之后，整个牛顿法潮流计算就剩两件事：一是把雅可比矩阵 J 的每个元素算对，二是选一个不会让迭代一开始就飞的初值。初值问题我放在后面踩坑部分细说，这里先解决雅可比矩阵。
2.2 雅可比矩阵四块结构的推导要点
在极坐标功率方程下，雅可比矩阵天然分成四块。用符号 H、N、J、L 表示：H 对应有功对相角的偏导，N 对应有功对电压的偏导（乘以电压），J 对应无功对相角的偏导，L 对应无功对电压的偏导。修正方程写成：
[ΔP] [H N] [Δθ ] [ΔQ] = [J L] * [ΔV/V]
注意我这里用的是 ΔV/V 的形式，也就是电压的相对修正量。这个形式的好处是矩阵元素和电压幅值的关系更对称，不易出错。推导时关键是区分对角元素和非对角元素，非对角元素直接从功率方程里对某个具体的 θ_j 或 V_j 求偏导；对角元素则要利用功率方程本身的表达式做化简。
我在纸上推过一遍，最终的公式如下。对任意 i ≠ j，令 θ_ij = θ_i - θ_j：
H_ij = V_i V_j (G_ij sinθ_ij - B_ij cosθ_ij) N_ij = V_i V_j (G_ij cosθ_ij + B_ij sinθ_ij) J_ij = -V_i V_j (G_ij cosθ_ij + B_ij sinθ_ij) L_ij = V_i V_j (G_ij sinθ_ij - B_ij cosθ_ij)
对角元素则为：
H_ii = -Q_i - B_ii V_i² N_ii = P_i + G_ii V_i² J_ii = P_i - G_ii V_i² L_ii = Q_i - B_ii V_i²
这里面的 P_i、Q_i 是用当前迭代点的电压相角重新算出来的功率值，而不是给定的注入功率。我之前在这个地方栽过一次——直接用给定值去组装雅可比，结果迭代到第三四次就完全不收敛了，后来才意识到雅可比矩阵必须在每一轮迭代里用当前状态重新计算。
2.3 牛顿法 MATLAB 核心代码
我直接给出一套能跑通的主流程代码，结构清晰，方便对照公式理解。
% 三节点系统数据 % bus = [编号 类型 V初始值 theta初始值(rad) 有功注入 无功注入] % 类型：1=PQ, 2=PV, 3=平衡 bus = [ 1 3 1.0 0 0 0 2 1 1.0 0 -0.5 -0.2 3 2 1.0 0 0.3 0 ]; % 线路数据：首端 末端 r x line = [ 1 2 0.02 0.04 1 3 0.01 0.03 2 3 0.015 0.035 ]; n = size(bus, 1); Y = zeros(n, n); for k = 1:size(line, 1) i = line(k,1); j = line(k,2); y = 1 / (line(k,3) + 1j * line(k,4)); Y(i,i) = Y(i,i) + y; Y(j,j) = Y(j,j) + y; Y(i,j) = Y(i,j) - y; Y(j,i) = Y(j,i) - y; end G = real(Y); B = imag(Y); ref = find(bus(:,2) == 3); PQ = find(bus(:,2) == 1); PV = find(bus(:,2) == 2); indP = [PQ; PV]; % 有功方程对应的节点 indQ = PQ; % 无功方程对应的节点 V = bus(:,3); th = bus(:,4); V(PV) = bus(PV, 3); % PV节点电压幅值固定 Psp = bus(:,5); Qsp = bus(:,6); tol = 1e-8; max_iter = 30; for iter = 1:max_iter Pc = zeros(n,1); Qc = zeros(n,1); for i = 1:n for j = 1:n th_ij = th(i) - th(j); Pc(i) = Pc(i) + V(i)*V(j)*(G(i,j)*cos(th_ij) + B(i,j)*sin(th_ij)); Qc(i) = Qc(i) + V(i)*V(j)*(G(i,j)*sin(th_ij) - B(i,j)*cos(th_ij)); end end dP = Pc - Psp; dQ = Qc - Qsp; F = [dP(indP); dQ(indQ)]; if norm(F, inf) < tol break; end np = length(indP); nq = length(indQ); H = zeros(np,np); N = zeros(np,nq); Jm = zeros(nq,np); L = zeros(nq,nq); for i = 1:np ii = indP(i); for j = 1:np jj = indP(j); if i == j H(i,j) = -Qc(ii) - B(ii,ii)*V(ii)^2; else th_ij = th(ii) - th(jj); H(i,j) = V(ii)*V(jj)*(G(ii,jj)*sin(th_ij) - B(ii,jj)*cos(th_ij)); end end end for i = 1:np ii = indP(i); for j = 1:nq jj = indQ(j); if ii == jj N(i,j) = Pc(ii) + G(ii,ii)*V(ii)^2; else th_ij = th(ii) - th(jj); N(i,j) = V(ii)*V(jj)*(G(ii,jj)*cos(th_ij) + B(ii,jj)*sin(th_ij)); end end end for i = 1:nq ii = indQ(i); for j = 1:np jj = indP(j); if ii == jj Jm(i,j) = Pc(ii) - G(ii,ii)*V(ii)^2; else th_ij = th(ii) - th(jj); Jm(i,j) = -V(ii)*V(jj)*(G(ii,jj)*cos(th_ij) + B(ii,jj)*sin(th_ij)); end end end for i = 1:nq ii = indQ(i); for j = 1:nq jj = indQ(j); if i == j L(i,j) = Qc(ii) - B(ii,ii)*V(ii)^2; else th_ij = th(ii) - th(jj); L(i,j) = V(ii)*V(jj)*(G(ii,jj)*sin(th_ij) - B(ii,jj)*cos(th_ij)); end end end dX = -[H N; Jm L] \ F; dth = dX(1:np); dVovV = dX(np+1:end); th(indP) = th(indP) + dth; V(indQ) = V(indQ) .* (1 + dVovV); V(PV) = bus(PV, 3); end这段代码我有意把雅可比矩阵的组装写成了循环，而不是用向量化技巧，目的是让你能够逐行对照上面的公式。实际工程代码里可以优化成稀疏矩阵操作，但学习阶段先把每一步看清楚更重要。
2.4 为什么二次收敛在潮流计算里那么明显
牛顿法理论上具有二次收敛特性。这意味着在解附近，每迭代一次，误差大致变成上一次的平方。从数值上体会就是：第一次迭代残差可能在 1e-2 量级，第二次到 1e-4，第三次到 1e-10。这是我在这个三节点算例里实际观察到的规律。但也正因为它靠近解时“发疯一样快”，初值如果离解太远，前面几步反而可能不降反升，这就对初值选择提出了要求。
3. P-Q分解法：有功和无功解耦的思路为什么能节省大量计算
3.1 两条物理假设：功率耦合关系是如何被解开的
P-Q 分解法也叫快速分解法，它针对牛顿法计算量大的痛点，在极坐标牛顿法的基础上做了两个关键假设。第一，正常运行下输电线路的电抗远大于电阻，一般高压输电网的 X/R 都在 5 到 10 以上，这时线路两端的有功流动主要取决于相角差，无功流动主要取决于电压幅值差；第二，正常运行时节点电压幅值接近 1.0 p.u.，相角差也比较小，可以近似取 cosθ≈1、sinθ≈0。
这两个假设放到雅可比矩阵里的结果就是：N 块和 J 块可以忽略不计，剩下的 H 块和 L 块也都能近似成常数矩阵——它俩本质上都只和节点导纳矩阵的虚部 B 有关。这样一来，修正方程就拆成了两个互不耦合的小方程组：一个只算相角修正，一个只算电压修正。矩阵不用每次迭代重新组装，也没有交叉耦合项，计算量大幅降低。
3.2 B' 矩阵和 B'' 矩阵的构造差异
有功修正方程和无功修正方程所依赖的矩阵，工程上习惯叫做 B' 和 B''。两者虽然都是导纳矩阵虚部派生的，但细节取舍不同。B' 通常忽略对地支路的影响，也就是线路充电电容和变压器等值支路不参与构建，同时删去平衡节点以及对应的行和列后，再用其余的 PV 和 PQ 节点构成有功修正方程。B'' 则要保留对地支路的影响，因为无功功率和电压的关系恰恰与这些对地导纳密切相关，矩阵只由 PQ 节点构成，PV 节点因为电压幅值恒定，不参与无功修正。
我在三节点算例中把线路充电电容设成了 0，所以 B' 和 B'' 恰好相等。但这在实际电网里是不成立的，如果你要做 IEEE 30 节点或 118 节点算例，一定不要把两者混用，否则迭代次数会多出不少，严重时直接不收敛。
3.3 P-Q分解法 MATLAB 核心代码
% 沿用上面构造的 Ybus、bus、line 数据 % 节点类型、PQ/PV/ref 节点索引均一致 B1 = imag(Y(indP, indP)); % 有功修正矩阵，注意符号方向 B2 = imag(Y(indQ, indQ)); % 无功修正矩阵 tol = 1e-8; max_iter = 50; for iter = 1:max_iter Pc = zeros(n,1); Qc = zeros(n,1); for i = 1:n for j = 1:n th_ij = th(i) - th(j); Pc(i) = Pc(i) + V(i)*V(j)*(G(i,j)*cos(th_ij) + B(i,j)*sin(th_ij)); Qc(i) = Qc(i) + V(i)*V(j)*(G(i,j)*sin(th_ij) - B(i,j)*cos(th_ij)); end end dP = Pc - Psp; dQ = Qc - Qsp; F = [dP(indP); dQ(indQ)]; if norm(F, inf) < tol break; end % 有功修正：求相角增量 dth = B1 \ (dP(indP) ./ V(indP)); % 无功修正：求电压幅值增量 dV = B2 \ (dQ(indQ) ./ V(indQ)); th(indP) = th(indP) + dth; V(indQ) = V(indQ) + dV; V(PV) = bus(PV, 3); end注意这里我用的残差定义是 “计算值减去给定值”，也就是 dP = Pc - Psp，dQ = Qc - Qsp。很多教材和网上的代码用的是另一种 “给定值减计算值”，那么 B1、B2 前面的符号就要相应翻转。我自己调试时最大的体会是：先固定一套残差符号，别在代码里混用，不然很容易得到一套“看起来在迭代、实际上完全走反方向”的结果。
4. 三节点算例实测：两种算法的收敛过程与结果对比
4.1 同一套数据下的最终计算结果
我用上面这段三节点数据实际跑下来，牛顿-拉夫逊法在第 4 到第 5 次迭代时残差降到 1e-8 以下，P-Q 分解法则需要大约 13 到 17 次迭代。两者的最终解应该完全一致，我拿到的节点电压和相角大致是：
| 节点 | 电压幅值 p.u. | 相角 deg | 节点类型 |
|---|---|---|---|
| 1 | 1.000 | 0.0 | 平衡节点 |
| 2 | 约 0.987 | 约 -4.6 | PQ节点 |
| 3 | 1.000 | 约 -2.1 | PV节点 |
节点 2 带的是纯负荷，电压略低于 1.0 符合直觉；节点 3 是 PV 节点，电压被锁定在 1.0。有功从节点 1 和节点 3 分别流向节点 2，所以节点 2 的相角最低，这也符合电力系统里“相角从发电机端往负荷端递减”的经验。
4.2 收敛速度与单次迭代成本的双重对比
很多人看到 P-Q 分解法迭代次数比牛顿法多好几倍，就误以为它没有优势，这是一个非常常见的误解。判断一个算法实际快不快，要看总耗时，而不是迭代次数。牛顿法每轮迭代都要重新计算雅可比矩阵的所有元素，再对变化的高阶矩阵做一次三角分解；而 P-Q 分解法的 B'、B'' 矩阵是常数，只需要在最开始做一次三角分解，后续每轮只是用这个分解好的因子去回代求解。
| 指标 | 牛顿-拉夫逊法 | P-Q分解法 |
|---|---|---|
| 迭代次数 | 4~5 | 13~17 |
| 单次迭代计算量 | 大，需组装雅可比并三角分解 | 小，常数矩阵回代 |
| 收敛特性 | 二次收敛 | 近似线性收敛 |
| 内存占用 | 每轮都要存新矩阵 | 只需存两个常数矩阵 |
| 编程难度 | 雅可比矩阵符号容易出错 | 矩阵含义更直观 |
从表格能看出，牛顿法适合对精度要求高、系统规模不算特别大的场景；P-Q 分解法则在在线调度、状态估计这类需要反复批量求解的大规模场景中更有优势。真实系统里，P-Q 分解法一次有功修正和一次无功修正交替进行，配合稀疏技术，速度优势会进一步拉大。
4.3 收敛过程的中间数据如何阅读
我在调试时会打印每一轮的残差范数，观察下降趋势。牛顿法最典型的表现是前三轮下降平缓，第四轮开始残差突然从 1e-3 掉到 1e-10；P-Q 分解法则是平稳地每次下降一个量级，偶尔中间有一次“卡住”也很正常，比如从 1e-5 到 1e-6 需要两次迭代。如果你看到某个方法前几轮残差不但不降反而猛增，不用急着怀疑算法，先检查初值是不是有问题。
5. 这类算法实现里最容易踩的五个坑
5.1 雅可比矩阵对角元素的符号问题
我见过不少人在手写牛顿法时，雅可比矩阵的非对角元素都能写对，一到对角元素就把符号弄反，因为对角元素不是直接求导出来的，而是要利用功率方程做化简。判断方法很简单：找一个三节点小系统，先不要让电压和相角偏离初始值太多，用一阶差分近似检验雅可比矩阵的每个元素，比如 (F(x+εe_j)-F(x))/ε 应该和你写的偏导公式一致。第一次难免歪，但这样自查一遍后，符号问题基本就绝迹了。
5.2 PV节点的无功越限处理
很多课程设计给出的算例里没有考虑 PV 节点的无功限制，所以代码里直接让 PV 节点电压恒定就行。但实际发电机有励磁和过载限制，无功出力是有上下限的。迭代过程中一旦某个 PV 节点的无功计算值超过限制，就必须把这个节点改成 PQ 节点，用它的无功限值作为给定值重新迭代。这个逻辑不写，算出来的电压剖面在工程上是不可用的。我建议在代码里预留一个变量，专门记录“PV转PQ”的节点编号。
5.3 迭代更新采用 ΔV 还是 ΔV/V 的混搭问题
牛顿法里我用的是 ΔV/V，P-Q 分解法里我用的是直接 ΔV。这两种写法对应的雅可比矩阵元素差一个电压倍数，在电压接近 1 p.u. 时差别不大，所以很多人混着用也能勉强收敛。但一旦系统重负荷、电压跌到 0.9 以下，这种混搭就会让迭代次数明显增加，甚至发散。我的习惯是在代码文件顶部写清楚“本文件采用 ΔV/V 修正”，每次都检查，避免从网上复制不同版本代码时混进另一套约定。
5.4 初值选择：平启动不是万能的，但离谱初值一定会发散
牛顿法的局部收敛性决定了它对初值很敏感。最常用的平启动是全部节点电压幅值取 1 p.u.，相角取 0。这是我强烈推荐的起点，因为电力系统在正常运行范围内，节点电压确实都在 1 p.u. 附近。但你如果随手把某个节点电压初值设成 0.5 p.u.，牛顿法可能在第一步就求出离谱的修正量，后面怎么拉都拉不回来。P-Q 分解法对初值的耐受力通常稍好一些，因为它的常数矩阵相当于“冻结”了耦合项。
5.5 对地支路和变压器变比的处理直接影响 B'、B''
我的三节点算例忽略了所有对地支路，所以 B' 和 B'' 一样。但这个简化不能带到真实算例里。输电线路的充电电容会显著影响无功分布，变压器非标准变比会改变两侧等值导纳，这些都要反映到节点导纳矩阵里。B' 和 B'' 的差异，本质上就是在说“哪些近似可以用于有功修正、哪些必须保留用于无功修正”。如果你发现自己用 P-Q 分解法在某个算例上迭代次数比牛顿法多了几十倍，先别骂算法，去检查一下是不是把 B'' 也建成了不含对地支路的版本。
我在把这些算法全部独立实现过一遍之后，最大的感受是：牛顿法教会你如何严谨地处理非线性方程组，P-Q 分解法教会你如何从物理直觉出发对复杂问题做合理简化。两套代码本身并不长，但背后的公式推导、符号约定、边界条件处理才是真正值钱的部分。你如果也能像我一样，先用三节点小算例把两条路都跑通，再去碰 IEEE 标准节点系统，会顺手很多。
