潮流计算在电力系统里属于那种“看起来简单、写起来全是细节”的东西。很多教材把公式推导梳理得很漂亮但一到 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 - θ_jH_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) 有功注入 无功注入] % 类型1PQ, 2PV, 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(np1: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 - PspdQ Qc - Qsp。很多教材和网上的代码用的是另一种 “给定值减计算值”那么 B1、B2 前面的符号就要相应翻转。我自己调试时最大的体会是先固定一套残差符号别在代码里混用不然很容易得到一套“看起来在迭代、实际上完全走反方向”的结果。4. 三节点算例实测两种算法的收敛过程与结果对比4.1 同一套数据下的最终计算结果我用上面这段三节点数据实际跑下来牛顿-拉夫逊法在第 4 到第 5 次迭代时残差降到 1e-8 以下P-Q 分解法则需要大约 13 到 17 次迭代。两者的最终解应该完全一致我拿到的节点电压和相角大致是节点电压幅值 p.u.相角 deg节点类型11.0000.0平衡节点2约 0.987约 -4.6PQ节点31.000约 -2.1PV节点节点 2 带的是纯负荷电压略低于 1.0 符合直觉节点 3 是 PV 节点电压被锁定在 1.0。有功从节点 1 和节点 3 分别流向节点 2所以节点 2 的相角最低这也符合电力系统里“相角从发电机端往负荷端递减”的经验。4.2 收敛速度与单次迭代成本的双重对比很多人看到 P-Q 分解法迭代次数比牛顿法多好几倍就误以为它没有优势这是一个非常常见的误解。判断一个算法实际快不快要看总耗时而不是迭代次数。牛顿法每轮迭代都要重新计算雅可比矩阵的所有元素再对变化的高阶矩阵做一次三角分解而 P-Q 分解法的 B、B 矩阵是常数只需要在最开始做一次三角分解后续每轮只是用这个分解好的因子去回代求解。指标牛顿-拉夫逊法P-Q分解法迭代次数4~513~17单次迭代计算量大需组装雅可比并三角分解小常数矩阵回代收敛特性二次收敛近似线性收敛内存占用每轮都要存新矩阵只需存两个常数矩阵编程难度雅可比矩阵符号容易出错矩阵含义更直观从表格能看出牛顿法适合对精度要求高、系统规模不算特别大的场景P-Q 分解法则在在线调度、状态估计这类需要反复批量求解的大规模场景中更有优势。真实系统里P-Q 分解法一次有功修正和一次无功修正交替进行配合稀疏技术速度优势会进一步拉大。4.3 收敛过程的中间数据如何阅读我在调试时会打印每一轮的残差范数观察下降趋势。牛顿法最典型的表现是前三轮下降平缓第四轮开始残差突然从 1e-3 掉到 1e-10P-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/VP-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 标准节点系统会顺手很多。
