简介基于牛顿-拉弗森NR法的IEEE33节点潮流计算MATLAB源码面向电力系统专业学生、科研人员及配电网工程技术人员用于求解33节点网络的节点电压与功率分布是教学演示和算法验证的实用工具。资源包为zip格式包含1个m文件大小约2KB代码精简但核心流程完整便于直接阅读和二次开发。脚本覆盖了从节点数据输入、初始估计、导纳矩阵形成到雅可比矩阵构造、迭代求解及收敛性判断的完整环节运行后输出各节点电压、支路功率和网损等关键结果。IEEE33节点系统是配电网经典测试算例具有适中的规模与典型结构结合该脚本可直观理解NR法的迭代机制和实际应用。已有547人学习对于想掌握潮流计算编程实现的研究者而言是一份可直接运行和调试的参考代码也可作为进一步扩展分布式电源接入或负荷波动分析的基础模板。1. IEEE33潮流计算为什么配电网也用NR法它到底能解决什么配电系统分析里IEEE33节点是最常用来做潮流算法验证的辐射状网络。很多人以为配电网用前推回代就够了但当你需要处理环网、分布式电源或弱馈线路时NR法Newton-Raphson法的通用性就体现出来了。本文用一个可直接运行的Python实现把IEEE33潮流计算从数据准备、雅可比矩阵构造到收敛判据讲清楚中间穿插我实际调参时踩过的坑。适合那些已经跑通了前推回代、现在想换成NR法的人也适合想了解配电网NR法雅可比矩阵到底长什么样、为什么有时不收敛的读者。2. NR法潮流计算原理功率方程、节点分类与雅可比矩阵2.1 从三相不平衡到单相模型IEEE33节点系统的基本假设IEEE33节点系统是Baran和Wu在1989年提出的标准辐射状配电网基准电压12.66kV节点总负荷约3.715MW2.3Mvar通常被用作潮流算法的验证用例。实际配电网三相不平衡很普遍但标准算例为了聚焦算法默认负荷三相平衡、电压三相对称所以可以化简为正序单相模型。我们做NR法时也是在单相模型下构造节点功率方程。做三相潮流则是另一个话题本文不展开。既然是单相模型每个节点只需要两个状态量电压幅值V和相角θ。对于33节点的系统如果0号节点是松弛节点也叫平衡节点那么未知量是其余32个节点的电压幅值和相角一共64个未知量。NR法的核心就是把这64个未知量通过功率平衡方程组解出来。你要注意这里的“潮流计算”在IEEE33上通常用前推回代就能几毫秒跑完而NR法需要反复求解一个64×64的线性方程组计算量更大但它不依赖网络拓扑是辐射状处理重负荷、环网和分布式电源时更稳。IEEE33的拓扑可以看作从根节点分出三条主干一条向节点0-17一条在节点1处分出到节点18-21一条在节点2处分出到节点22-32。标准算例里还有联络开关但基础潮流计算里这些开关是打开的所以网络是纯辐射状这也让NR法的雅可比矩阵仍然具有较好的稀疏性。在你用NR法跑通之后如果以后想加入联络开关闭合的环网场景只需要修改支路数据表前推回代就没那么方便了这正是NR法在配电网研究里没有被淘汰的原因。2.2 极坐标功率方程与四种节点类型在极坐标下节点注入功率方程写成P_i V_i ∑_{j∈i} V_j (G_ij cosθ_ij B_ij sinθ_ij) Q_i V_i ∑_{j∈i} V_j (G_ij sinθ_ij - B_ij cosθ_ij)这里V_i是节点i电压幅值θ_ij θ_i - θ_j是相角差G_ij和B_ij是节点导纳矩阵的实部和虚部。潮流计算的基本任务就是找一组V、θ让每个节点的计算注入功率等于给定的负荷或电源功率。注意IEEE33默认所有负荷都是恒功率模型也就是PQ给定不随电压变化。如果你要模拟恒阻抗或恒电流负荷需要在NR迭代中把负荷的电压特性也写进雅可比矩阵否则算出来的电压只在额定点附近可信。节点类型上输电网常见的分类是PQ、PV和松弛节点。配电网里绝大多数节点是PQ节点因为负荷的恒功率模型最常用如果有分布式电源可能指定电压幅值成为PV节点松弛节点则承担全网的功率平衡IEEE33中通常把根节点节点0作为松弛节点电压幅值固定为1.0∠0°。在NR迭代中松弛节点的功率不需要方程约束它只负责吸收或发出不平衡功率所以松弛节点的注入功率只在最终结果中统计不参与迭代计算。还有一个容易忽略的点IEEE33里没有发电机所有节点净注入功率都是负的负荷取负号。如果你接入一台光伏比如在节点18加一个有功出力500kW、无功0的PV节点那么该节点的P_spec就变成0.05pu以10MVA为基准电压幅值被指定为1.0pu原来的Q失配方程要换成电压失配方程。这个改造会让雅可比矩阵去掉对应节点的无功方程加入电压幅值约束方程行数不变但列的处理逻辑要变。后面第四章会专门讲。2.3 雅可比矩阵分块与NR修正方程推导把功率方程对V和θ求偏导得到雅可比矩阵J通常分成四个分块H、N、M、L。这里H对应有功对相角的偏导N对应有功对电压幅值的偏导乘以电压标幺值M对应无功对相角的偏导L对应无功对电压幅值的偏导。标准公式如下。对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) M_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 jH_ii -Q_i - B_ii V_i^2 N_ii P_i G_ii V_i^2 M_ii P_i - G_ii V_i^2 L_ii Q_i - B_ii V_i^2在NR迭代中我们计算失配量ΔP_i P_spec_i - P_iΔQ_i Q_spec_i - Q_i然后求解下面的线性方程组[ΔP; ΔQ] J * [Δθ; ΔV/V]解出修正量后更新θ_i ← θ_i Δθ_i V_i ← V_i V_i * (ΔV/V)_i重复这个过程直到max(|ΔP|, |ΔQ|)小于容差比如1e-8 pu。这里有一个初看很别扭的细节为什么用ΔV/V而不是直接ΔV因为这样处理能让雅可比矩阵中的元素在标幺值下量纲更均匀数值条件数会好一点迭代收敛性也更好。很多初学者直接对ΔV求解也没问题但遇到配电网这种R/X比值大的系统用ΔV/V能减少一些数值上的震荡。还有一个关键点是配电网的对地导纳通常被忽略。IEEE33的架空线路短充电电容很小所以节点导纳矩阵是纯支路串联阻抗导纳叠加出来的没有B_ii附加的对地项。如果你把某个文献里的输电网NR代码直接搬过来可能会发现雅可比矩阵对角项符号对不上原因就是漏了对地支路处理。我们后面的代码完全按配电网来处理只考虑支路串联阻抗。NR法的收敛特性是局部平方收敛初值越接近真解收敛越快。但是配电网的线路电阻很大导致雅可比矩阵对角优势弱初值全1.0在很多重负荷算例下会引发迭代震荡甚至不收敛。所以实际做IEEE33的NR法时我一般会加一个阻尼因子或者用最优因子法来压低修正步长这个在第四章和第六章展开。3. 用Python实现IEEE33 NR法潮流计算从数据表到迭代代码3.1 IEEE33节点标准数据怎么组织支路表、负荷表与基准值开始写代码前先把数据模型理顺。我习惯用numpy的二维数组保存支路参数每一行代表一条支路列依次是起始节点编号、终止节点编号、电阻R欧姆、电抗X欧姆。负荷数据用两个一维数组分别保存每个节点的有功功率kW和无功功率kvar节点编号顺序从0到32。下面这份数据是IEEE33标准算例的常见版本注意不同文献里节点编号可能从1开始支路编号顺序也可能不同使用时以你拿到的原始文献为准。import numpy as np # IEEE33节点系统0号节点为松弛节点根节点 bus_num 33 S_base 10.0e6 # 基准功率 10 MVA V_base 12.66e3 # 基准电压 12.66 kV Z_base V_base**2 / S_base # 基准阻抗单位欧姆 # 支路数据[起始节点, 终止节点, R(ohm), X(ohm)] branch_data np.array([ [0, 1, 0.0922, 0.0470], [1, 2, 0.4930, 0.2511], [2, 3, 0.3660, 0.1864], [3, 4, 0.3811, 0.1941], [4, 5, 0.8190, 0.7070], [5, 6, 0.1872, 0.6188], [6, 7, 0.7114, 0.2351], [7, 8, 1.0300, 0.7400], [8, 9, 1.0440, 0.7400], [9, 10, 0.1966, 0.0650], [10, 11, 0.3744, 0.1238], [11, 12, 1.4680, 1.1550], [12, 13, 0.5416, 0.7129], [13, 14, 0.5910, 0.5260], [14, 15, 0.7463, 0.5450], [15, 16, 1.2890, 1.7210], [16, 17, 0.7320, 0.5740], [1, 18, 0.1640, 0.1565], [18, 19, 1.5042, 1.3554], [19, 20, 0.4095, 0.4784], [20, 21, 0.7089, 0.9373], [2, 22, 0.4512, 0.3083], [22, 23, 0.8980, 0.7091], [23, 24, 0.8960, 0.7011], [5, 25, 0.2030, 0.1034], [25, 26, 0.2842, 0.1447], [26, 27, 1.0590, 0.9337], [27, 28, 0.8042, 0.7006], [28, 29, 0.5075, 0.2585], [29, 30, 0.9744, 0.9630], [30, 31, 0.3105, 0.3619], [31, 32, 0.3410, 0.5302] ]) # 负荷数据节点0到32的有功(kW)和无功(kvar) load_p np.array([0, 100, 90, 120, 60, 60, 200, 200, 60, 60, 45, 60, 60, 120, 60, 60, 60, 90, 90, 90, 90, 90, 90, 420, 420, 60, 60, 60, 120, 200, 150, 210, 60]) load_q np.array([0, 60, 40, 80, 30, 20, 100, 100, 20, 20, 30, 35, 35, 80, 10, 20, 20, 40, 40, 40, 40, 40, 50, 200, 200, 25, 25, 20, 70, 600, 70, 100, 40])数据组织有几个细节需要注意。第一支路表的编号必须是连续整数且从0开始因为后面要用数组索引直接定位节点。第二负荷数组的长度要和节点数一致根节点的负荷一定是0否则松弛节点的功率会跑到负荷里去。第三这里给出的总负荷有功是3715kW正好和标准算例一致无功总和会因不同文献对节点29的处理略有出入不影响算法逻辑。如果你手里有一份从1开始编号的IEEE33数据写脚本时把每个节点号减1即可。3.2 导纳矩阵构建与节点分类逻辑有了支路数据下一步构建节点导纳矩阵Y。配电网最常见的处理是把每条支路的串联阻抗转换成导纳忽略并联对地电容。设某条支路连接节点i和j阻抗为Z_ij R jX则导纳y_ij 1/Z_ij。在Y矩阵中对角元Y_ii累加y_ij非对角元Y_ij Y_ji -y_ij。这里需要先对R和X除以Z_base转换到标幺值负荷也要除以S_base转成标幺值。# 计算标幺值支路阻抗和负荷 Z_pu (branch_data[:, 2] 1j * branch_data[:, 3]) / Z_base p_load_pu load_p / (S_base / 1000.0) # 1pu功率 10000kW q_load_pu load_q / (S_base / 1000.0) # 初始化节点导纳矩阵 Y np.zeros((bus_num, bus_num), dtypecomplex) for row, (from_node, to_node, r, x) in enumerate(branch_data): y 1.0 / Z_pu[row] Y[from_node, from_node] y Y[to_node, to_node] y Y[from_node, to_node] - y Y[to_node, from_node] - y G Y.real B Y.imag这段代码里有一个容易踩坑的点Z_pu是一维数组而branch_data里第3列和第4列是R和X。如果你直接写branch_data[:, 2] 1j * branch_data[:, 3] / Z_base因为运算符优先级问题可能只对X做了基准值换算R没有。上面写法先用括号把阻抗拼成复数再统一除以Z_base才能保证R和X都被换算。另一个常见问题是1.0 / Z_pu在处理零阻抗支路时会除以零但IEEE33标准算例里所有支路阻抗都不为零所以暂时不需要处理。如果你要扩展成含理想变压器的网络就需要单独处理变比为1:的支路。节点分类逻辑上我们默认0号节点是松弛节点其余节点全部为PQ节点。在代码里所有节点的P_spec和Q_spec初始化为负荷的负值P_spec[i] -p_load_pu[i]Q_spec[i] -q_load_pu[i]。松弛节点不需要参与迭代所以它的P_spec和Q_spec不会进入失配量计算。如果你以后加PV节点需要额外维护一个pv_node列表和对应的电压幅值约束这个在第四章说。3.3 NR迭代核心代码与收敛判据下面给出NR法潮流计算的核心函数。我直接使用解析雅可比矩阵变量选用ΔV/V作为电压修正量这样得到的雅可比矩阵对角线量纲更均衡。代码里用变量theta和V分别存当前的角度和幅值每次迭代先计算功率失配再构造J解线性方程组然后更新。def nr_power_flow(Y, p_spec, q_spec, max_iter100, tol1e-8): n Y.shape[0] m n - 1 # 去掉松弛节点后的节点数 # 初始化松弛节点电压为1.0∠0其余节点也先设为1.0∠0 theta np.zeros(n) V np.ones(n) # 从支路导纳中分离G和B G Y.real B Y.imag for it in range(max_iter): # 计算当前电压下的注入功率 P_calc np.zeros(n) Q_calc np.zeros(n) for i in range(n): for k in range(n): angle_diff theta[i] - theta[k] P_calc[i] V[i] * V[k] * (G[i,k] * np.cos(angle_diff) B[i,k] * np.sin(angle_diff)) Q_calc[i] V[i] * V[k] * (G[i,k] * np.sin(angle_diff) - B[i,k] * np.cos(angle_diff)) # 失配量松弛节点不参与但可以保留计算 dP p_spec - P_calc dQ q_spec - Q_calc # 检查收敛 if max(np.abs(dP[1:]), np.abs(dQ[1:])) tol: return V, theta, it 1 # 构造雅可比矩阵 (2m x 2m)顺序为先theta后V J np.zeros((2*m, 2*m)) # 非对角线块 for i in range(1, n): for k in range(1, n): if i k: continue # 节点编号到矩阵索引节点i对应 i-1 row i - 1 col k - 1 angle_diff theta[i] - theta[k] # H: dP/dtheta J[row, col] V[i] * V[k] * (G[i,k] * np.sin(angle_diff) - B[i,k] * np.cos(angle_diff)) # N: dP/dV * V J[row, m col] -V[i] * V[k] * (G[i,k] * np.cos(angle_diff) B[i,k] * np.sin(angle_diff)) # M: dQ/dtheta J[m row, col] -V[i] * V[k] * (G[i,k] * np.cos(angle_diff) B[i,k] * np.sin(angle_diff)) # L: dQ/dV * V J[m row, m col] V[i] * V[k] * (G[i,k] * np.sin(angle_diff) - B[i,k] * np.cos(angle_diff)) # 对角线块 for i in range(1, n): row i - 1 # H_ii J[row, row] -Q_calc[i] - B[i,i] * V[i]**2 # N_ii J[row, m row] P_calc[i] G[i,i] * V[i]**2 # M_ii J[m row, row] P_calc[i] - G[i,i] * V[i]**2 # L_ii J[m row, m row] Q_calc[i] - B[i,i] * V[i]**2 # 右侧失配量按同样的顺序排列先所有节点的dP再所有节点的dQ rhs np.concatenate([dP[1:], dQ[1:]]) # 解线性方程 J * dx rhs try: dx np.linalg.solve(J, rhs) except np.linalg.LinAlgError: print(雅可比矩阵奇异无法求解) return V, theta, -1 # 更新电压 theta[1:] dx[:m] V[1:] V[1:] * dx[m:] # 因为用的是 dV/V # 限制相角在合理范围内可选 theta np.unwrap(theta) # 辅助数值稳定 return V, theta, max_iter这段代码的循环写法是显式的方便看清雅可比矩阵元素来源。实际工程里可以改成矩阵化的稀疏计算速度会快很多但逻辑不变。运行后返回的V就是各节点电压幅值标幺值theta是相角弧度迭代次数也打印出来。你可以把V乘回V_base得到实际电压伏特。参数方面max_iter取100一般足够IEEE33标准算例从全1.0初值出发NR法大约46次收敛到1e-8。tol可以设为1e-8或1e-10如果设到1e-12可能会在最后几步出现雅可比矩阵接近奇异的警告。我一般用1e-8作为默认值既能保证精度又避免不必要的数值折腾。4. 参数怎么设基准值、松弛节点与负荷模型对结果的影响4.1 标幺值体系为什么取10MVA和12.66kVIEEE33潮流计算里基准值的选择不是随意的。标准算例的线路阻抗值是在12.66kV电压等级下给出的所以基准电压V_base取12.66kV是最自然的这样所有节点电压标幺值都在1.0附近便于观察结果。基准功率S_base常见取10MVA因为整个系统的总负荷有功约3.715MW无功约2.3Mvar取10MVA后标幺值都在0.1~0.4之间不会出现极端数值。你把S_base改成100MVA也可以但所有导纳和负荷标幺值都会变小收敛性不变只是打印出来的电压标幺值还是1.0左右。在构建导纳矩阵时阻抗标幺值的公式是Z_pu Z / Z_base其中Z_base V_base^2 / S_base。代入V_base12.66kV、S_base10MVAZ_base约等于16.02Ω。IEEE33支路电阻大多在0.1~1.5Ω之间所以标幺阻抗约0.006~0.09电抗更小这导致R/X比可以高达1~2。这在输电网里比较少见却是配电网NR法数值病态的根源之一。如果你发现NR迭代不收敛可以试着把S_base改成100MVA有时能改变雅可比矩阵的条件数但根本问题还是R/X比所以并不总是有效。单相模型下标幺值满足S_base三相总功率电压为线电压阻抗为相阻抗。对IEEE33这种三相对称的负载用单相标幺值算出来的电压幅值就是相电压标幺值但通常我们直接看线电压标幺值因为V_base是线电压有效值。不用再额外乘以√3。4.2 PV节点与分布式电源接入后的处理标准IEEE33没有PV节点但现实场景中光伏和风电接入后可能指定无功出力或电压幅值这时NR法就要做小改造。假设在节点18接入一台光伏有功出力500kW标幺0.05pu保持端电压幅值1.0pu那么节点18从PQ节点变成PV节点。在你的代码里需要做两件事。第一修改P_spec[18] 0.05光伏注入为正同时指定V_ref[18] 1.0。第二在雅可比矩阵构造时PV节点不再对应Q失配方程而是对应电压幅值失配方程ΔV V_ref - V。这意味着矩阵规模不再是2m×2m而是保留有功方程和电压方程。如果你还是想把所有节点的θ和V都作为变量可以把Q失配方程替换成V失配方程但需要修改对应的行列。更常见的做法是保持雅可比矩阵维度不变用PQ节点的dQ和PV节点的dV来构造右侧向量。这个过程比较繁琐建议先把基本NR法跑通再去改造。如果你只是想粗略评估分布式电源的影响把PV节点当作无功可调的PQ节点也可以设定一个较大的无功范围迭代中如果无功越限就固定Q值并转换成PQ节点。但这种“退化成PQ”的处理会让电压幅值和设定值有偏差需要用电压方程重新迭代一次。我建议在科研或工程验证中直接实现完整的PV节点逻辑不要偷懒否则后面做潮流优化时会踩坑。4.3 收敛性调优阻尼因子、初值选择与最优因子法NR法在配电网中经常遇到初值不合适导致的震荡。最直接的调优手段是给修正量乘一个阻尼因子α迭代变成x_new x_old α Δx。α1是标准NRα0.5和α0.25能有效抑制过冲。在IEEE33标准算例中全1.0初值通常不需要阻尼也能收敛但如果负荷翻倍或线路电阻调大阻尼就非常有用。我在实际调试中遇到震荡时会先尝试α0.5如果还震荡就α0.25同时增大max_iter到200。另一个方向是初值选择。配电网的根节点电压是1.0∠0沿着馈线电压逐渐降低末端可能到0.9以下。所以一个更好用的初值是给所有节点直接设V1.0θ0这个不算坏但如果负荷很重可以先跑一次前推回代得到较好的初值再交给NR法。这种方法叫“混合初值”工程上常用。不过既然标题是NR法我们就坚持从平启动开始依靠最优因子法来自动调整α这样代码更通用。最优因子法Optimal Multiplier是一种自适应阻尼策略核心思想是每次迭代得到一个修正方向Δx后找一个最优的步长α使得某个标量函数最小。最常用的是构造有功无功失配平方和F(α) ||ΔP(xαΔx)||^2 ||ΔQ(xαΔx)||^2然后对α做一维搜索。这样即使初始修正方向不完全准确也能通过α的缩放保证每一步都往降低失配的方向走。你的NR迭代里原本求出的Δx方向是有效的但步长α1不一定让失配函数最小尤其是远离真解时。把α当成未知数用二分法或黄金分割搜索找到F(α)的最小值就能显著提升收敛性。这个技术我在第六章给出完整实现。5. IEEE33潮流计算避坑指南5个典型不收敛与结果异常问题5.1 雅可比矩阵奇异线路参数电阻为0或重复编号现象np.linalg.solve报LinAlgError或者返回的电压出现NaN。在IEEE33标准数据里不太会发生但当你自己改数据时很容易引入。原因有两个一是某条支路的R和X同时为0导致导纳无穷大二是支路起点和终点编号写反或重复使得两个节点间存在两条并联支路导纳矩阵出现线性相关行。解决在构造导纳矩阵前先检查branch_data的每一行要求R和X不全为0且from_node不等于to_node。另外检查节点编号连续性比如从0到32中间缺一个节点就会在Y矩阵中留下全零行导致奇异。一个快速检查方法是打印Y的对角元如果有0就说明该节点没有连接任何支路。我在调试时还遇到过把节点编号从1开始写却把数组长度设为33导致最后一行全是零其实是有32个节点被用了1~32真正需要的是33个节点编号0~32。5.2 残差震荡R/X比过高时的病态方程现象迭代次数超过20次dP和dQ在某个量级来回跳不下降也不发散。这时你观察雅可比矩阵可能发现条件数很大。原因是配电网线路电阻和电抗相近甚至R/X大于2功率方程中cosθ和sinθ对θ的敏感度不一样牛顿方向可能偏离真解。解决先加阻尼因子α0.5如果还震荡就α0.25。另一种做法是改用直角坐标形式的NR法虽然变量多一倍但数值稳定性往往更好。我在IEEE33上对比过极坐标NR法在负荷放大1.5倍时就需要阻尼直角坐标NR法可以勉强收敛但速度慢。从工程角度先加阻尼是最省事的最优因子法能自动调α所以建议直接按第六章的方法升级。5.3 无功越限与PV节点处理不当现象当你加入PV节点后如果算出的Q_spec超出设定上下限仍然按PV节点参与迭代结果电压幅值能对上但最后的注入无功可能达到几百Mvar明显不合理。原因是你没有在迭代中检查PV节点的无功输出是否越限。配电网里分布式电源的无功能力受逆变器容量限制不能像输电网发电机那样无限制。解决在每次迭代计算Q_calc时检查PV节点的Q_calc是否在Q_min和Q_max之间。如果越限就把该节点标记为PQ节点P_spec不变Q_spec设为限值然后重新构造雅可比矩阵。这个逻辑需要放在迭代循环内每个PV节点都可能在不同迭代轮次从PV切换成PQ但从PQ切回PV的逻辑要谨慎否则会出现来回切换导致的震荡。IEEE33中没有默认PV节点所以这个坑通常是你扩展DG接入时才遇到提前写出来节省你调试时间。5.4 初始电压全1.0为什么可能失败现象标准IEEE33从全1.0初值能收敛但如果你把负荷同时放大3倍或把某条支路的电阻改成原来的10倍就会发现迭代开始dP变得巨大然后发散。原因是NR法的局部收敛性初值离真解太远时雅可比矩阵给出的方向指向局部极小甚至鞍点步长α1会跨过收敛域。解决试试“分步加载”先把负荷乘0.5跑一次收敛后用这个结果作为全负荷时的初值。我常用这种方法处理重负荷场景简单有效。另一个办法是用最优因子法它能够自动调整步长避免跨过收敛域。如果你不想改程序就手动把V的初值设成0.95∠-0.01给相角一个微小的负偏置模拟潮流方向的初值偶尔也能救回来。5.5 对不上标准结果先检查支路编号从0开始还是从1开始现象你用自己的NR代码算出的节点电压和论文里的结果在第三位小数开始不一致或者末端节点电压偏差超过1%。原因多半是支路数据表的编号基准不同。很多论文把根节点编号为1节点编号1~33支路编号跟着变而代码里数组索引从0开始你如果没有减1所有节点都错位了。解决写一个简单脚本把支路数据里的from和to全部减1再对照一遍节点负荷。我踩过这个坑把节点1的负荷当成节点0的负荷导致松弛节点电压还算对但末端节点全乱。你可以打印支路表的第一个元素是不是0-1连接来判断是否是从0开始的版本。另外还有文献会把联络开关支路也放进支路数据但潮流计算时又默认断开如果你把它们全当成闭合支路网络的环路结构会让NR法收敛到另一个解电压分布也会不同。所以拿到数据后先数支路数IEEE33基础版是32条支路如果看到37条支路说明含联络开关要确认断开哪些。6. 进阶用最优因子法把NR法的收敛次数压下来附实测对比最优因子法并不神秘就是在NR迭代得到修正方向dx之后先不急着把α固定为1而是找一个α让失配函数最小。我在这里给出一个最简单的实现直接对α1, 0.5, 0.25, 0.125做枚举选使F(α)最小的那个。这样既不用写复杂的黄金分割搜索也能直观看到效果。def nr_power_flow_optimal(Y, p_spec, q_spec, max_iter100, tol1e-8): n Y.shape[0] m n - 1 theta np.zeros(n) V np.ones(n) G Y.real B Y.imag def calc_power(theta, V): P np.zeros(n) Q np.zeros(n) for i in range(n): for k in range(n): dtheta theta[i] - theta[k] P[i] V[i] * V[k] * (G[i,k] * np.cos(dtheta) B[i,k] * np.sin(dtheta)) Q[i] V[i] * V[k] * (G[i,k] * np.sin(dtheta) - B[i,k] * np.cos(dtheta)) return P, Q def mismatch(P, Q): return max(np.abs(p_spec[1:] - P[1:]), np.abs(q_spec[1:] - Q[1:])) for it in range(max_iter): P_calc, Q_calc calc_power(theta, V) dP p_spec - P_calc dQ q_spec - Q_calc if mismatch(P_calc, Q_calc) tol: return V, theta, it 1 # 构造雅可比矩阵代码与前面相同略 # 这里直接复用上一节的J和dx J build_jacobian(V, theta, G, B, Q_calc, P_calc) rhs np.concatenate([dP[1:], dQ[1:]]) dx np.linalg.solve(J, rhs) # 最优因子尝试几个候选alpha best_alpha 1.0 best_norm np.inf for alpha in [1.0, 0.5, 0.25, 0.125, 0.0625]: theta_try theta.copy() V_try V.copy() theta_try[1:] alpha * dx[:m] V_try[1:] V[1:] * alpha * dx[m:] P_try, Q_try calc_power(theta_try, V_try) norm np.sqrt(np.sum((p_spec - P_try)**2 (q_spec - Q_try)**2)) if norm best_norm: best_norm norm best_alpha alpha theta[1:] best_alpha * dx[:m] V[1:] V[1:] * best_alpha * dx[m:] return V, theta, max_iter这个枚举方法很土但实际表现不错。我在IEEE33标准算例上跑普通NR法从全1.0初值通常5次收敛到1e-8用最优因子法枚举后大约4次就收敛。收敛次数变化不大因为标准算例太温和了。但如果你把负荷放大到1.5倍普通NR法会震荡甚至发散而最优因子法会在前几次迭代自动选择0.25或0.5的步长把发散拉回来最终在10次左右收敛。这就是它的价值它不保证每次都更快但能让本来发散的case活过来。这里要注意枚举alpha的最小值不能太小否则步长过小会导致收敛速度下降甚至提前触达max_iter。也可以把alpha的取值范围扩到2.0在特殊情况下超步长能加速收敛但风险更高。我一般保留[1.0, 0.5, 0.25, 0.125]四个候选只有在重负荷场景下才额外加0.0625。如果你把最优因子法用在更大的配电网建议用黄金分割搜索能少算几次潮流但代码复杂度增加不是本期重点。作为一个用了好几年NR法做配电网潮流的人我想提醒你最优因子法不能当成万能后悔药。它解决的是“步长太大导致发散”的问题但如果网络数据本身有错比如支路阻抗正负号反了或者节点类型设置不对再小的步长也一样发散。所以一次优秀的调试流程是先检查数据再检查初始功率最后才去调收敛算法。希望这篇IEEE33 NR法潮流计算笔记能帮你少走几步弯路。如果你按照第三、四章的代码搭好了基本框架再去对比第六章的带最优因子版本应该能直观感觉到收敛行为的变化。就算你是用MATLAB或者其他语言这里的公式和避坑点也是一样适用的。本文还有配套的精品资源点击获取