国内做配电网研究的同学几乎都绕不开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 MVA32条支路参数全部为标幺值总负荷为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级别的精度。对比一下主流的几种潮流方法方法收敛速度代码复杂度处理PV节点/DG扩展配电网适用性高斯-赛德尔法慢线性收敛低较差教学为主牛顿-拉夫逊法快二次收敛中好天然支持PQ/PV/平衡节点中小规模系统很稳前推回代法快低较差处理PV和环网需改造大规模辐射网效率高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法的入口节点类型已知量待求量典型位置平衡节点Vθ节点电压幅值、相角有功、无功注入变电站母线、无穷大电源PQ节点有功、无功注入电压幅值、相角负荷节点、恒功率DGPV节点有功、电压幅值无功注入、相角恒电压控制DG、发电机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^2N_ii P_i G_ii * V_i^2K_ii P_i - G_ii * V_i^2L_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节点系统支路参数标幺值基准值 SB100 MVA, UB12.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 MVA100000 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/(RX)这种实数除法那整个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的电压下降得最快k1.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 四类常见问题速查表现象可能原因排查与处理第一轮迭代就发散负荷功率或阻抗量纲错误检查P/Q是否除以标幺基准阻抗是否从Ω换算到p.u.在真实解附近来回振荡雅可比矩阵符号错误或Q方程符号问题用数值差分验证J矩阵检查Hii/Lii正负号收敛但电压结果偏低支路数据版本不同或负荷数据有差异核对总负荷是否3.715 MW j2.3 Mvar节点30无功是否600 kvar迭代次数偏多20次负荷水平过高、初值不良或接近电压极限检查是否在负荷扫描模式下尝试放宽damping或更换初值这张表基本覆盖了我自己调试NR程序时遇到过的所有问题类型。如果你卡住了先对照这张表过一遍比盲目改代码效率高得多。5.4 收敛容差与迭代步长的经验值最后聊一下收敛容差和步长控制。收敛容差tol取1e-6到1e-8都是合理范围我一般取1e-8因为多迭代两次对33节点来说成本几乎可以忽略。但如果你要做成千上万次潮流的优化迭代建议放宽到1e-6能省不少时间。如果遇到NR法在重负荷下不收敛可以引入一个简单的阻尼策略在更新修正量时乘一个阻尼系数α比如α0.8或者0.5alpha 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节点、加光伏出力、做负荷水平扫描都只需要改输入数据算法代码完全不用动。这个封装习惯让我后续开发省了无数时间强烈建议你也这么做。