简介：面向电力系统专业学生与研究人员的牛顿拉夫逊法潮流计算MATLAB程序包，针对潮流计算中非线性方程组的迭代求解问题，给出了从数学模型到代码实现的完整方案。程序涵盖节点导纳矩阵构建、雅可比矩阵计算、功率不平衡量求解及迭代收敛判断等关键模块，并附带使用说明文档与课本例题解答，便于对照理论逐步验证结果。压缩包共包含18个文件，包括14个m脚本、2个docx文档、1个pdf参考资料及1个txt备注，m脚本中既有主运行程序，也有独立的子函数与数据输入模块，结构清晰、便于二次开发。压缩包整体大小仅为565KB，资源紧凑但功能完整，便于下载和部署。目前已有941人学习下载，适合刚接触电力系统计算的新手，也适合需要快速搭建潮流计算工具的开发者参考借鉴。 研究生上电力系统分析那会儿，第一次自己写潮流计算程序，对着教材上的公式硬啃了一个多星期，反复看书上推导过程，再看代码，好不容易才能把三节点的算例跑通。这几年做电网仿真和教学辅导，用MATLAB写过好几个版本的潮流程序，从三节点教学算例到上百节点的配电网模型都折腾过。回头看，牛顿拉夫逊法这个算法本身并不神秘，关键是把极坐标下的功率方程、雅可比矩阵的物理含义和迭代更新的逻辑捋顺。这篇就专门讲清楚这些问题，给一套能直接跑通的MATLAB实现，再把我踩过的坑和调试经验一并交代。
这篇文章适合三类读者：正在学电力系统分析、需要交潮流计算大作业的学生；刚接触电网仿真、想用MATLAB快速算一个算例的工程师；以及想弄懂牛顿拉夫逊法内部推导、不满足于直接调用工具箱的爱好者。看完你不仅能把程序跑起来，还能知道每一步背后的原因，遇到不收敛、矩阵奇异这类问题也知道从哪里排查。
1. 牛顿拉夫逊法在潮流计算中的位置与选择理由
1.1 潮流计算到底要解决什么问题
潮流计算，简单说就是给定电网的拓扑结构、线路参数和部分节点的功率注入，求每个节点的电压幅值和相角，进而算出每条支路的功率分布。电力系统是一个高度非线性的网络，节点注入功率和节点电压之间不是线性关系，所以必须求解非线性方程组。
工程上通常把节点分成三类：PQ节点，只给定有功和无功注入，待求电压幅值和相角；PV节点，给定有功注入和电压幅值，待求无功和相角；平衡节点，一般只有一个，给定电压幅值和相角，用来平衡全网功率，待求有功和无功注入。节点类型的划分不是一个数学技巧问题，它对应的是电力系统实际的物理运行状态。比如普通负荷母线就是PQ节点，发电机的机端母线往往是PV节点，而系统里需要有一台机组承担频率和电压基准，它就是平衡节点。
把这三类节点弄清楚，后面程序里哪些方程该写、哪些变量该更新，就一目了然了。
1.2 为什么选牛顿拉夫逊而不是高斯赛德尔
早期潮流计算常用高斯赛德尔法，它实现简单，迭代格式也很直观，但收敛速度慢，病态系统下还容易发散。牛顿拉夫逊法最大的优势是二阶收敛，迭代次数和系统规模关系不大，通常三五次就能把不平衡量压到1e-6以下，而且收敛特性好，对初值不太敏感，只要初值不偏离得太离谱，基本都能收敛。
代价是每次迭代都要重新组装雅可比矩阵，并且求解一个线性方程组。对于中型系统，这个计算量完全可接受；对于超大规模电网，可以配合稀疏矩阵技术和因子表分解来提速，这就留到后面扩展部分再说。
2. 推导极坐标功率方程与雅可比矩阵
2.1 极坐标下的节点功率方程
潮流计算的数学模型是节点功率平衡方程。在极坐标下，节点i的注入功率可以写成：
P_i = V_i Σ_{j=1}^{n} V_j (G_ij cos(θ_i - θ_j) + B_ij sin(θ_i - θ_j))
Q_i = V_i Σ_{j=1}^{n} V_j (G_ij sin(θ_i - θ_j) - B_ij cos(θ_i - θ_j))
其中G_ij和B_ij分别是节点导纳矩阵的实部和虚部，θ_i是节点i的电压相角。这两个方程看着长，其实每一项的物理含义很明确：节点电压通过导纳矩阵和相邻节点耦合，功率是电压乘积与导纳的三角函数组合。
程序里计算Pcal和Qcal时，我习惯直接把公式写成两层循环，避免调用向量化函数导致逻辑不透明。虽然效率不如矩阵运算，但对理解算法和调试非常有帮助，算三节点、五节点这种规模根本不在乎这点循环开销。
2.2 雅可比矩阵四块结构与物理含义
对功率方程求偏导，可以得到雅可比矩阵，它分成四块：
- H = ∂P/∂θ，对应相角对有功的影响；
- N = V * ∂P/∂V，对应电压幅值对有功的影响；
- K = ∂Q/∂θ，对应相角对无功的影响；
- L = V * ∂Q/∂V，对应电压幅值对无功的影响。
这四个子矩阵组装起来就是完整的雅可比矩阵J。PV节点不参与无功方程，所以对应K和L的行列要剔除；平衡节点完全不需要参与迭代，它不产生不平衡量。
非对角元素和对角元素的公式需要区分：
非对角（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)
K_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)
对角（i = j）：
H_ii = -Q_i - B_ii V_i²
N_ii = P_i + G_ii V_i²
K_ii = P_i - G_ii V_i²
L_ii = Q_i - B_ii V_i²
这里特别说明一下N和L矩阵，公式里乘了电压幅值V_j，目的就是让修正方程里对应的未知量是ΔV/V而不是ΔV。这样处理的好处是雅可比矩阵各元素的量纲更一致，数值上也更稳定。这个细节在编程时非常关键，如果忘记乘V_j，修正量的物理意义就变了，更新电压时直接加ΔV反而会出错。
2.3 迭代修正方程与收敛判据
每次迭代都要解这样一个修正方程：
J * [Δθ; ΔV/V] = [ΔP; ΔQ]
其中ΔP和ΔQ是给定注入功率和当前电压计算得到的注入功率之差。求出来的Δθ和ΔV/V，分别加到相角和电压幅值上，就得到新一轮的电压估计值。
收敛判据我一般用inf范数，就是取所有不平衡量的最大绝对值，小于容差就认为收敛。容差设多少合适？教学场景1e-6足够了，工程计算可以设1e-8甚至更小，但对应的时间成本也会上升。程序里我会设置最大迭代次数，比如20次，防止某些病态情况下无限循环。
3. MATLAB完整实现与核心代码解读
3.1 数据结构和节点导纳矩阵构建
先给出一套完整的MATLAB实现，用的算例是三节点系统：节点1是PQ节点（负荷），节点2是PV节点（发电机），节点3是平衡节点。程序从零构建节点导纳矩阵，不依赖任何内置工具箱函数，任何版本的MATLAB都能运行。
clear; clc; %% 输入数据 % 节点数组: [类型, P, Q, V, theta] % 类型: 1=PQ, 2=PV, 3=平衡 % P和Q是节点注入功率，负荷取负值 node = [ 1, -0.5, -0.2, 1.0, 0; 2, 0.6, 0, 1.05, 0; 3, 0, 0, 1.0, 0 ]; % 支路数组: [节点i, 节点j, R, X, B/2] % 这里B/2是线路对地电纳的一半，单位标幺值 branch = [ 1, 2, 0.010, 0.050, 0.020; 1, 3, 0.015, 0.060, 0.025; 2, 3, 0.012, 0.055, 0.020 ]; n = size(node, 1); %% 构建节点导纳矩阵 Y = zeros(n, n); for k = 1:size(branch, 1) i = branch(k, 1); j = branch(k, 2); z = branch(k, 3) + 1j * branch(k, 4); % 支路阻抗 y = 1 / z; % 支路导纳 b = branch(k, 5); % 对地电纳 % 自导纳叠加 Y(i, i) = Y(i, i) + y + 1j * b; Y(j, j) = Y(j, j) + y + 1j * b; % 互导纳 Y(i, j) = Y(i, j) - y; Y(j, i) = Y(j, i) - y; end G = real(Y); B = imag(Y);节点导纳矩阵是潮流计算的基石，它的构造逻辑其实很朴素：每条支路都对两端节点的自导纳有贡献，同时对两端节点的互导纳产生负的贡献。变压器支路、线路对地电容都在这个基础上扩充，理解了这个框架，后续加任何元件都不会慌。
3.2 主迭代循环与雅可比矩阵组装
接下来是核心迭代部分。程序里先根据节点类型把参与P方程和Q方程的节点列表整理出来，然后进入迭代循环，每次迭代都重算功率不平衡量和雅可比矩阵。
%% 初始化 V = node(:, 4); theta = node(:, 5); type = node(:, 1); Psp = node(:, 2); % 给定有功 Qsp = node(:, 3); % 给定无功 pq = find(type == 1); % PQ节点列表 pv = find(type == 2); % PV节点列表 slack = find(type == 3); % 平衡节点列表 np = [pq; pv]; % 参与P方程的所有节点 nq = pq; % 参与Q方程的节点，只有PQ节点 tol = 1e-6; maxiter = 20; for iter = 1:maxiter % 计算当前电压下的节点功率 Pcal = zeros(n, 1); Qcal = zeros(n, 1); for i = 1:n for j = 1:n theta_ij = theta(i) - theta(j); Pcal(i) = Pcal(i) + V(i) * V(j) * ... (G(i,j) * cos(theta_ij) + B(i,j) * sin(theta_ij)); Qcal(i) = Qcal(i) + V(i) * V(j) * ... (G(i,j) * sin(theta_ij) - B(i,j) * cos(theta_ij)); end end % 不平衡量 dP = Psp - Pcal; dQ = Qsp - Qcal; dF = [dP(np); dQ(nq)]; if max(abs(dF)) < tol fprintf('收敛，迭代次数: %d\n', iter); break; end if iter == maxiter error('达到最大迭代次数，潮流不收敛'); end %% 组装雅可比矩阵 H = zeros(n, n); N = zeros(n, n); K = zeros(n, n); L = zeros(n, n); for i = 1:n for j = 1:n if i == j H(i, i) = -Qcal(i) - B(i, i) * V(i)^2; N(i, i) = Pcal(i) + G(i, i) * V(i)^2; K(i, i) = Pcal(i) - G(i, i) * V(i)^2; L(i, i) = Qcal(i) - B(i, i) * V(i)^2; else theta_ij = theta(i) - theta(j); H(i, j) = V(i) * V(j) * ... (G(i,j) * sin(theta_ij) - B(i,j) * cos(theta_ij)); N(i, j) = V(i) * V(j) * ... (G(i,j) * cos(theta_ij) + B(i,j) * sin(theta_ij)); K(i, j) = -V(i) * V(j) * ... (G(i,j) * cos(theta_ij) + B(i,j) * sin(theta_ij)); L(i, j) = V(i) * V(j) * ... (G(i,j) * sin(theta_ij) - B(i,j) * cos(theta_ij)); end end end % 按参与节点做行列缩减 J = [H(np, np), N(np, nq); K(nq, np), L(nq, nq)]; % 解修正方程 dX = J \ dF; % 更新电压 nnp = length(np); theta(np) = theta(np) + dX(1:nnp); V(nq) = V(nq) + V(nq) .* dX(nnp+1:end); % 输出每次迭代的最大不平衡量，方便观察收敛过程 fprintf('迭代第 %d 次，最大不平衡量: %.2e\n', iter, max(abs(dF))); end %% 结果输出 fprintf('\n===== 潮流结果 =====\n'); for i = 1:n fprintf('节点 %d: V = %.6f, theta = %.6f rad, P = %.6f, Q = %.6f\n', ... i, V(i), theta(i), Pcal(i), Qcal(i)); end3.3 修正量更新时容易踩的坑
上面这段代码里最容易被忽略的就是更新电压那一行。因为雅可比矩阵用的是N和L，对应的是∂P/∂V乘上V，解出来的修正量后半部分是ΔV/V，不是ΔV。所以更新时不能直接写V(nq) = V(nq) + dX(nnp+1:end)，必须乘上原来的电压V(nq)。我第一次写的时候就栽在这里，结果迭代过程剧烈振荡，始终不收敛。花了大半天排查，最后一行一行对照教材公式才发现问题。
另一个常见的坑是参与方程的节点列表。PV节点虽然待求无功，但在极坐标牛顿法中它的Q方程不参与迭代，因为Q不是给定的，没有不平衡量。如果错把PV节点也放进nq列表，雅可比矩阵会多出几行几列，方程数大于未知量数，线性方程组直接变成超定问题，\运算虽然不会报错，但算出的结果毫无意义，潮流也收敛不到正确值。
4. 调试实录：收敛失败与数值问题的排查
4.1 潮流不收敛的常见根因
程序写完后第一件事不是看结果，而是看收敛曲线。如果最大不平衡量一路下降，很快进到1e-6以内，说明算法没问题。如果出现振荡、卡住不动甚至发散，就要按顺序排查这几个地方。
先检查节点导纳矩阵。拿一个只有两个节点的简单系统手算Y矩阵，跟程序输出的对比，这是最有效的验证手段。很多错误出在支路对地电纳的叠加方式上，我见过有人把B/2只加在一个端点上，差点把对称性丢了。Y矩阵不对称，雅可比矩阵算出来必然有问题。
再检查功率计算的符号。负荷功率在程序里是负值，发电机是正值，这是电力系统的约定俗成。如果符号弄反，物理上相当于负荷在发电，发电机在用电，潮流不可能收敛到有意义的解。我给学生改作业时发现这个错误出现频率相当高。
如果Y矩阵和功率计算都对，还是不收敛，就要看初值。极坐标牛顿法对初值不算敏感，但也不是完全无所谓。我习惯用平启动：所有PQ节点电压设1.0∠0°，PV节点设给定的电压幅值、相角0°。绝大多数算例都能从平启动收敛。如果遇到重负荷系统或者高R/X比的配电网，平启动可能失效，这时可以对电压幅值赋一个略高于1.0的初值，比如1.05，往往能改善收敛。
4.2 雅可比矩阵程序bug的典型特征
雅可比矩阵如果组装错了，最典型的现象是收敛速度非常慢，或者迭代几步后不平衡量跳到另一条曲线上。正常牛顿法在不平衡量比较大的时候，每步迭代数值下降得非常猛，降不了说明方向错了。
一个有效的自查办法是数值差分验证。对某个节点电压加一个微小扰动，比如1e-6，重新算功率不平衡量，看它和解析雅可比矩阵对应元素是否一致。这个方法虽然笨，但验证角角落落非常可靠。我每次把雅可比公式记不太清的时候，就写个临时脚本做差分校验，几分钟就能确认公式写没写对。
另一个隐蔽问题是MATLAB矩阵索引搞混。J矩阵的组装顺序是：[H(np,np), N(np,nq); K(nq,np), L(nq,nq)]，对应的未知量顺序是[Δθ(np); ΔV/V(nq)]。如果行列顺序不一致，矩阵维度虽然对得上，但方程对应的物理量张冠李戴，结果完全不可信。这个bug不报错，最难发现，建议在更新前输出size(dX)和对应列表长度，确认一致再继续。
5. 从三节点到实际电网：扩展与性能建议
5.1 大规模系统的改进方向
上面这套程序跑通三节点、五节点毫无压力，但拿到上百节点的真实电网模型就会暴露出性能问题。根源主要在两个方面：每次迭代都要重新组装完整雅可比矩阵并求解稠密线性方程组，规模一大，计算量以二次方甚至三次方的速度增长。
工程上的标准做法是把雅可比矩阵当稀疏矩阵处理，用sparse来存，配合MATLAB的稀疏线性方程组求解器，速度和内存占用都会明显改善。更进一步可以用因子表法，在迭代前先对雅可比矩阵做LU分解，因为牛顿法接近收敛时雅可比矩阵变化不大，可以用同一个因子表迭代多次再重新分解，这就是所谓的"保留非线性"或"快速解耦"思路。
快速解耦法（P-Q分解法）是另一个实用改进方向。利用高压输电网络有功主要取决于相角、无功主要取决于电压幅值的物理特性，把雅可比矩阵简化为常数矩阵B'和B''，每次迭代只需要一次前代回代，速度能提升一个量级。很多商业潮流软件在牛顿法收敛困难时会自动切换到P-Q分解法，两者结合使用非常稳健。
5.2 配电网潮流计算中的适用性
牛顿拉夫逊法是输电网潮流计算的事实标准，但在配电网里要小心。配电网线路电阻和电抗之比很高，R/X比经常大于1，传统牛顿法的收敛特性会明显恶化，甚至直接发散。配电网更主流的方法是前推回代法，它利用辐射状结构的层次关系，每一步从末端往电源推功率，再从电源往末端推电压，计算量小、收敛性对R/X不敏感。
不过现在很多配电网网架越来越复杂，出现了联络开关形成的弱环结构，纯前推回代法处理起来很麻烦，这时牛顿法反而能处理。我的经验是：先判断网络结构，辐射状用前推回代，弱环网用牛顿法加环网处理，这样既稳又快。
最后再说一个折腾过很多次的经验：MATLAB代码里尽量别把所有东西都写到主脚本里，把Y矩阵构建、雅可比组装、核心迭代拆成独立函数，输入输出接口留清楚。这样调试一个模块的时候不用每次都从头跑，改成新算例也只要换数据文件。这个习惯让我后来从三节点换到IEEE标准算例时，改动量小到只需要替换输入数据表而已。
本文还有配套的精品资源，点击获取
