配电网潮流计算入门：基于Matlab的IEEE 33节点牛顿-拉夫逊法实现与验证

刚入行做配电网分析那会儿，我的导师给我丢了第一个任务：用Matlab把IEEE 33节点系统的潮流算出来。当时我手里只有一段不知道从哪儿抄来的牛顿-拉夫逊法脚本，跑出来的结果和文献对不上，电压曲线七拐八拐，网损差了一倍多，整个人处于一种"代码在跑、但我完全不知道它在算什么"的状态。后来老老实实把NR法的数学原理啃了一遍，雅可比矩阵一行一行推导，再用IEEE 33节点这个经典算例反复对照，才算真正把潮流计算这件事搞透。

时至今日，IEEE 33节点系统依然是配电网研究里出场率最高的基准算例，NR法也是各种潮流算法里最通用、最值得手写一遍的算法。这篇文章就把当年踩过的坑、验证过有效的实现思路、可以直接抄作业的Matlab代码结构全部梳理出来，顺便讨论几个实际应用时会遇到的扩展问题。适合刚接触电力系统仿真的学生，也适合正在做配电网规划、分布式电源接入、主动配电网分析的工程师参考。

1. IEEE 33节点系统解析：为什么它成了配电网的"标准考场"

1.1 拓扑结构与基准参数

IEEE 33节点系统最早是针对配电网馈线优化问题提出的标准算例，结构是一条从变电站母线引出的12.66kV辐射状馈线，包含33个节点、32条支路，另外还有5条联络开关支路，合计37条支路。节点1是平衡节点，代表变电站二次侧母线，电压恒定为额定值；其余节点大多数是PQ节点，负荷按三相平衡的恒功率负荷处理。

系统的标准参数如下：

这里要特别注意阻抗基值的计算：Zbase = U²/S = 12.66² / 10 = 16.03Ω。算例中给出的支路参数，通常直接是每段线路的电阻电抗有名值，单位是欧姆，使用时需要除以Zbase归算成标幺值。我在初学阶段就是在这里吃过亏——用有名值直接堆进标幺公式里，结果导纳矩阵差了16倍，算出来的电压、网损全都对不上号。

1.2 参数口径与版本陷阱

网上能找到的IEEE 33节点数据其实有多个版本，这个坑很少被教材提到，但实操中非常要命。不同版本之间的差异主要体现在三个地方：

第一个是单位。大多数论文和开源代码用的是有名值电阻电抗，单位是Ω，但有些教材或代码为了方便，直接给出标幺值支路参数。如果你拿到的数据没有注明单位，一定要先通过Zbase换算验证一下。

第二个是联络开关支路。IEEE 33节点系统在原始定义里包含5条联络开关，默认是断开状态，也就是说正常运行方式是辐射状开环。有的复现代码为了研究网络重构，会把联络开关全部闭合，变成弱环网运行。这两种工况下潮流结果完全不同，网损和最低电压都有明显差异。

第三个是负荷单位。标准的33节点系统总负荷是3715kW+j2300kvar，但有些文献把某个节点的负荷做了微调，或者增加了一个"DG接入点"版本。做结果对比时，别只看电压曲线形状，要把总负荷、基准容量、联络开关状态全部核对一致，再下结论。

1.3 为什么配电网研究偏爱这个算例

IEEE 33节点系统之所以在配电网论文里高频出现，原因很实在：规模适中、数据完整、结论公认。相比3节点9节点系统，它更接近真实馈线的拓扑特征——线路长度几十米到几百米不等，负荷沿馈线分散分布，末端电压跌落明显；相比IEEE 123节点这种大规模系统，它又足够轻量，Matlab脚本几秒钟就能跑完，非常适合作为算法验证平台。

更关键的是，这个算例的基准结果在文献里高度一致。系统总有功网损大约是202~211kW，最低电压出现在节点18附近，电压幅值大概在0.903~0.913p.u.之间。不同文献因为收敛精度和参数口径略有差异，但都在这个范围内。这个"共识结果"就是天然的对拍基准，你写完了NR法程序，算出来落在这个区间，基本可以确定核心逻辑是对的；如果偏差很大，那就可以直接进入排查模式了。

2. NR法原理解读：不是套公式，而是理解雅可比矩阵怎么来的

2.1 潮流方程的极坐标写法与失配量

潮流计算本质上是在求解一个非线性方程组：在给定的网络拓扑和负荷条件下，找到一组节点电压幅值和相角，使得每个节点注入的复功率满足基尔霍夫定律。用公式表达，节点注入电流等于导纳矩阵乘以节点电压向量，而注入功率等于电压乘以注入电流的共轭。

把节点电压写成极坐标形式U∠θ，导纳矩阵拆成电导G和电纳B，那么节点i的注入有功和无功可以写为：

P_i = U_i Σ U_j (G_ij cosθ_ij + B_ij sinθ_ij)

Q_i = U_i Σ U_j (G_ij sinθ_ij - B_ij cosθ_ij)

其中θ_ij = θ_i - θ_j。

所谓的失配量（mismatch），就是给定状态量计算出来的注入功率与实际注入功率之差：ΔP = P_sp - P_calc，ΔQ = Q_sp - Q_calc。当电压解收敛时，失配量趋向于零。NR法的迭代过程，本质上就是反复计算失配量，并通过雅可比矩阵求出状态量的修正方向，逐步把失配量压到阈值以下。

2.2 雅可比矩阵的四个分块与向量化构造

NR法比高斯-赛德尔法收敛快得多，核心就在于用雅可比矩阵给出了准确的修正方向。在极坐标形式下，雅可比矩阵划分为四个分块：H、N、M、L。其中H对应P对相角θ的偏导，N对应P对电压幅值U的偏导（配合ΔU/U的修正变量形式），M对应Q对相角的偏导，L对应Q对电压幅值的偏导。

这里有一个特别实用的编程技巧：不需要一上来就用二阶循环一个个填雅可比矩阵元素，而是可以利用Matlab的向量化运算，一次把四个分块全部算出来。

% 节点相角差矩阵和电压外积矩阵 SinT = sin(theta - theta.'); CosT = cos(theta - theta.'); U_out = U * U.'; % 电压外积，U是列向量 % 先填所有节点对的非对角项（包括i==j时的自导纳项） H = U_out .* (G .* SinT - B .* CosT); N = U_out .* (G .* CosT + B .* SinT); M = -N; L = H; % 再补对角修正项 H = H - diag(Q); N = N + diag(P); M = M + diag(P); L = L + diag(Q); % 组装完整雅可比，变量顺序为 [theta_all; U_all] J = [H N; M L];

这个写法的逻辑是：先按统一的公式填满所有元素，再利用对角修正项把偏导中的额外项补齐。这样写出来的代码简洁、没有下标循环，而且不容易在i=j和i≠j的边界条件上出错。

2.3 修正方程的组装与状态量更新

完整的雅可比矩阵是2n×2n，n为节点数，但并不是所有变量都参与迭代。平衡节点电压幅值和相角都已知，不参与修正；PV节点电压幅值已知，无功方程不参与迭代。所以最后真正求解的修正方程是：

J_red · Δx = ΔS_red

其中J_red是删除冗余行列后的雅可比，ΔS_red是保留节点的有功和无功失配量。修正变量Δx包含两部分：所有非平衡节点的相角修正量Δθ，以及所有PQ节点的电压幅值相对修正量ΔU/U。

Matlab里做这个删除操作时，我建议先构造完整的J和失配量向量，再用索引一次性删掉不需要的行列，逻辑最清晰：

ref = 1; % 平衡节点编号 pq = find(type == 0); % PQ节点 pv = find(type == 2); % PV节点 % 待求解的变量索引：非平衡节点的θ，PQ节点的U/U tVar = (1:nb)'; tVar(ref) = []; uVar = pq; dS_red = [dP(tVar); dQ(uVar)]; J_red = J([tVar; nb + uVar], [tVar; nb + uVar]); dx = J_red \ dS_red; % 回代更新状态量 nt = length(tVar); theta(tVar) = theta(tVar) + dx(1:nt); U(uVar) = U(uVar) + U(uVar) .* dx(nt+1:end);

注意这里修正变量用的是ΔU/U而不是ΔU，所以更新时要用U点乘修正量。这种写法数值上更平稳，也是Matpower等成熟工具采用的变量形式。

3. Matlab实现全流程：从零把NR法跑通

3.1 数据准备：33节点参数怎么录入

写NR法程序，第一步是把网络参数组织成程序能读的数据结构。我的习惯是仿照Matpower的bus和branch矩阵格式，用两个矩阵分别保存节点数据和支路数据。节点矩阵每一行是：节点编号、有功负荷、无功负荷、有功发电、无功发电、节点类型。支路矩阵每一行是：首端节点、末端节点、电阻、电抗、变比标志。

% bus: [节点编号 Pd(kW) Qd(kvar) Pg(kW) Qg(kvar) 类型] % 类型：1=平衡节点，0=PQ节点，2=PV节点 bus = [ 1 0 0 0 0 1; 2 100 60 0 0 0; 3 90 40 0 0 0; 4 120 80 0 0 0; 5 60 30 0 0 0; 6 60 20 0 0 0; % ... 完整数据见IEEE 33节点标准算例 ]; % branch: [首端 末端 R(Ω) X(Ω) 变比] branch = [ 1 2 0.0922 0.0470 1; 2 3 0.4930 0.2511 1; 3 4 0.3660 0.1864 1; % ... 完整数据见IEEE 33节点标准算例 ];

注意负荷单位是kW和kvar，计算时统一除以基准容量10MVA转换成标幺值。完整的33节点参数表在IEEE官方文档和大量论文附录里都能找到，GitHub上也有许多开源代码仓库，建议下载后先核对总负荷是否3715kW+j2300kvar，再继续往下走。

3.2 导纳矩阵与标幺值基准

导纳矩阵是潮流的骨架，所有计算都建立在Ybus之上。构造方法很直接：遍历每条支路，把支路导纳加到对应的自导纳和互导纳位置。配电网没有变压器（或者用变比=1简化处理），所以不需要处理非标准变比支路。

Vbase = 12.66e3; % 基准电压 12.66kV Sbase = 10e6; % 基准容量 10MVA Zbase = Vbase^2 / Sbase; % 阻抗基值 = 16.03Ω nb = size(bus, 1); nl = size(branch, 1); Y = zeros(nb, nb); for k = 1:nl i = branch(k, 1); j = branch(k, 2); z = (branch(k, 3) + 1j * branch(k, 4)) / Zbase; % 支路阻抗归算 y = 1 / z; 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);

这里有个细节值得多说一句：33节点系统支路数只有32条，用这种双层循环完全没性能压力。但如果以后需要算几百上千节点的配电网，建议用稀疏矩阵构造Ybus，MATLAB里把Y初始化成sparse(nb,nb)即可，速度能快一个数量级。

3.3 完整迭代主循环与关键代码

完成数据准备和Ybus构造后，接下来就是NR法的核心——迭代主循环。我把一个能跑的框架放在这里，关键步骤都有注释，你可以直接在自己的环境里运行并观察中间结果。

% 初始化 U = ones(nb, 1); % 电压幅值初值，平启动 theta = zeros(nb, 1); % 相角初值 Psp = (bus(:,4) - bus(:,2)) / Sbase; % 节点有功给定值 Qsp = (bus(:,5) - bus(:,3)) / Sbase; % 节点无功给定值 pq = find(bus(:,6) == 0); pv = find(bus(:,6) == 2); ref = find(bus(:,6) == 1); tVar = (1:nb)'; tVar(ref) = []; uVar = pq; tol = 1e-8; % 收敛阈值 maxIter = 50; for iter = 1:maxIter % 计算当前状态下的注入功率 V = U .* exp(1j * theta); I = Y * V; S = V .* conj(I); P = real(S); Q = imag(S); dP = Psp - P; dQ = Qsp - Q; % 收敛判据：最大失配量小于阈值 if max(abs([dP(tVar); dQ(uVar)])) < tol fprintf('收敛于第%d次迭代\n', iter); break; end % 构造雅可比矩阵 SinT = sin(theta - theta.'); CosT = cos(theta - theta.'); U_out = U * U.'; H = U_out .* (G .* SinT - B .* CosT); N = U_out .* (G .* CosT + B .* SinT); M = -N; L = H; H = H - diag(Q); N = N + diag(P); M = M + diag(P); L = L + diag(Q); J = [H N; M L]; J_red = J([tVar; nb + uVar], [tVar; nb + uVar]); dS_red = [dP(tVar); dQ(uVar)]; dx = J_red \ dS_red; % 更新状态量 nt = length(tVar); theta(tVar) = theta(tVar) + dx(1:nt); U(uVar) = U(uVar) + U(uVar) .* dx(nt+1:end); if iter == maxIter error('NR法未在最大迭代次数内收敛'); end end

这段代码里有几个容易踩的地方：一是初始值，配电网推荐用平启动，即所有节点电压幅值初值设为1.0p.u.、相角设为0，在这个算例上NR法通常5~8次迭代就能收敛。二是收敛判据要盯失配量的最大值，而不是盯修正量dx的大小——dx很小但失配量还很大，是可能出现的情况，说明修正被雅可比病态放大或缩小了。三是求解修正方程建议用左除号反斜杠，Matlab对稠密矩阵会自动选择合适的数值分解，比手算逆矩阵又快又稳。

3.4 结果验证：网损、电压分布与画图

跑完迭代不等于完成任务，还得把结果从标幺值还原成有名值，并和文献基准结果对拍。电压幅值直接乘以12.66kV就是有名值电压，节点功率乘以10MVA就是有名值功率。网损的计算方法是把所有支路损耗求和，或者用根节点注入功率减去总负荷。

% 支路潮流与网损 lossTotal = 0; for k = 1:nl i = branch(k, 1); j = branch(k, 2); yij = 1 / ((branch(k,3) + 1j*branch(k,4)) / Zbase); Vi = U(i) * exp(1j * theta(i)); Vj = U(j) * exp(1j * theta(j)); Iij = yij * (Vi - Vj); Sij = Vi * conj(Iij); Sji = Vj * conj(-Iij); lossTotal = lossTotal + real(Sij + Sji); end fprintf('系统总有功网损: %.4f MW (%.2f kW)\n', lossTotal * Sbase / 1e6, lossTotal * Sbase / 1e3);

可视化方面，我一般画出两个图：一个是节点电压幅值分布曲线，横轴是节点编号，纵轴是电压标幺值，可以直接看到最低电压点；另一个是支路有功功率分布，能看到潮流沿馈线递减的趋势。画图代码很简单：

figure; plot(1:nb, U, 'o-', 'LineWidth', 1.5); xlabel('节点编号'); ylabel('电压幅值 (p.u.)'); title('IEEE 33节点系统电压分布'); grid on;

画完图先别急着收工，看一眼最低电压点和网损值是否落在合理范围。如果你的最低电压在0.90p.u.附近、网损在200~211kW之间，基本说明程序跑通了。

4. 收敛性暗礁：配电网NR法的经典陷阱与排查方法

4.1 大R/X比值带来的麻烦

很多人把高压输电网的程序直接拿来算配电网，结果发现迭代发散或者收敛极慢，这是配电网潮流计算最经典的坑。原因在于配电线路的R/X比值远高于输电网，输电网动辄0.1以下的比值，配电网经常在0.5~2甚至更高。

这个特性对PQ分解法是致命的，因为PQ分解法在推导过程中为了解耦，忽略了电阻的影响，R/X比值越大，近似误差越大，甚至可能完全不收敛。NR法因为用的是完整雅可比矩阵，理论上不受R/X比值影响，但实际表现为收敛速度放慢、对初值更敏感。所以我一直建议，配电网仿真优先用NR法，别图省事用PQ分解法，后者在33节点这类系统上翻车概率非常高。

4.2 初值、重负荷与发散现象归类

迭代发散并不总是程序bug，也有可能是系统本身已经处于重负荷甚至无解状态。我总结了几种常见的收敛异常现象，可以对照排查：

有一个很实用的经验：当程序不收敛时，先把所有负荷乘以0.5运行，如果轻载下能收敛，说明核心逻辑大概率没问题，问题出在运行点本身。这个方法帮我在调试中快速区分"程序错了"和"负荷太重了"。

4.3 定位雅可比符号错误的差分校验法

NR法程序里最隐蔽的错误就是雅可比矩阵符号写反，或者i和j的偏导方向搞混，导致迭代第一步失配量不减反增。要快速定位这类错误，可以用有限差分法逐列校验雅可比矩阵。

做法是：在某个电压状态点（比如平启动点），分别在第j个变量上加一个微小扰动ε=1e-6，重新计算注入功率，用 (P(x+εe_j)-P(x))/ε 作为偏导的数值近似，和你解析算出的雅可比矩阵第j列对比。用Matlab写一个校验函数，几行代码的事：

% 差分校验雅可比，以P对theta的偏导为例 eps = 1e-6; H_num = zeros(nb, nb); for j = 1:nb theta_pert = theta; theta_pert(j) = theta_pert(j) + eps; P_pert = powerCalc(Y, U, theta_pert); H_num(:, j) = (P_pert - P) / eps; end disp(max(abs(H_num - H(:,:)))); % 接近0说明解析雅可比正确

注意无功和有功的注入功率计算公式都要校验。函数的实现需要与主程序共用同一个Ybus和功率计算函数。这个技巧我在接手别人的老代码时用过无数次，效率极高，比一行行对公式盯半天靠谱得多。

4.4 收敛性改进的实用手段

如果NR法在重负荷工况下始终不收敛，可以尝试几个基本手段。第一是加阻尼因子，迭代更新时把修正量乘一个0.5~0.9的系数，牺牲部分收敛速度来换取稳定性，尤其适合接近电压崩溃点的运行状态。第二是改用平启动加电压幅值下限约束，提前把启动电压设置为0.95p.u.左右，避免迭代过程中电压跌到负值导致计算崩溃。第三是使用连续潮流法，从轻载开始逐步增加负荷，每一级用上一级的解作为初值，这种方法能算出完整的PV曲线，也是分析配电网极限供电能力的有力工具。

还有一个小细节经常被忽略：Matlab默认的double类型精度足够高，但收敛阈值不要设置得太苛刻。对于33节点这类系统，收敛阈值1e-8已经非常安全，继续往下降不会明显改善结果精度，反而可能因为数值噪声导致无法达到收敛判据。

5. 从标准算例到实用扩展：DG接入与程序模块化

5.1 标准算例结果对照

把上面的程序跑通后，用IEEE 33节点标准参数算一遍，正常情况应该得到以下结果（以10MVA为基准，收敛阈值1e-8）：

最低电压出现的节点是18，这在配电网拓扑上很有代表性：节点18位于最长的一条分支线路末端，距离电源点最远，线路阻抗累积最大，压降也最明显。这也是为什么33节点系统常被用来做电压无功控制研究——它的"弱点"非常清晰，改进措施的效果也容易量化。

如果你的计算结果落在这个区间，恭喜你，NR法主程序已经验证通过。这之后才能真正进入"应用"阶段，而不是停留在复现代码。

5.2 分布式电源接入：从PQ到PV节点的改造

配电网研究中经常要模拟分布式光伏或风电接入。分布式电源在潮流计算里有两种建模方式：一种是恒功率因数控制，等效为PQ节点，把负的有功负荷叠加到对应节点即可，这种最简单；另一种是电压控制模式，等效为PV节点，需要在NR法里增加PV节点处理逻辑。

前面第二章已经提到了PV节点的行列删除规则，这里补充一下PV节点的完整处理流程。给定PV节点有功和电压幅值后，NR迭代中PV节点不参与无功失配的计算，但每次迭代完成后需要计算该节点的无功注入Q，并检查是否超过限值。如果Q越限，说明该分布式电源已经无法维持电压，需要从PV节点退化为PQ节点，以限值作为无功给定值，然后重新迭代。

MATLAB中的实现要点就是：在组装修正方程之前，把PV节点的Q方程行和U列删掉；在每次循环计算P、Q时，保留PV节点的Q值用于越限判断。这个改造在我的实际项目里几乎必做，因为分布式电源接入研究避不开电压越限和功率越限问题。

5.3 把潮流程序封装成可复用的分析工具

等NR法调通以后，我建议立即做一件事：把整个程序从脚本改造成函数。输入是bus数据和branch数据，输出是节点电压、支路潮流、网损和收敛信息。这样后续做配电网重构、网损优化、分布式电源选址定容时，只需要循环调用这个潮流函数，不需要每次复制粘贴代码。

我自己的工程模板大致是这样的：

function [U, theta, P, Q, loss, iter] = runNR(bus, branch, options) % 输入：bus节点数据，branch支路数据，options结构体（收敛阈值等） % 输出：电压幅值、相角、注入功率、网损、迭代次数 % 内部：Ybus构造、NR迭代主循环、结果后处理

有了这个函数后，做DG接入分析时，只需要修改bus矩阵中对应节点的发电功率，再调用runNR即可，几分钟就能得到一组新场景的潮流结果。

我在实际做配电网规划项目时，经常需要在十几个方案之间反复比较电压分布和网损指标，这时候模块化的潮流函数就是整个分析平台的基石。再往后做优化算法时也只是在runNR外面套一层寻优循环，每次评估适应度时调用一次潮流，一套代码通吃多种研究场景。

如果你也是刚把NR法跑通，我强烈建议下一步手边放一份Matpower的源码，不一定要调用它，而是对照着它的makeJac函数读一遍，你会发现自己刚写出来的向量化雅可比矩阵和成熟工具的思路高度一致。那才是真正摸到潮流计算门道的感觉。