配电网潮流计算：基于IEEE33节点与牛顿-拉夫逊法的Python实现

简介面向电力系统专业学生与从业人员该MATLAB脚本实现了基于牛顿-拉弗森法NR法的IEEE33节点潮流计算可求解配电网各节点电压幅值与相角、支路功率及网损适用于教学演示、课程设计和科研算法验证。资源包内仅有1个m文件压缩包约2KB代码结构简洁、无额外依赖便于直接阅读和二次开发。脚本完整涵盖数据输入、初始状态设置、雅可比矩阵形成、迭代修正、收敛判断与结果输出等核心步骤并生成各节点电压和线路潮流数据能够直观呈现NR法求解非线性方程组的全过程。目前已有547人学习下载适合需要掌握电网潮流计算原理、以IEEE33节点为基准算例开展仿真与算法对比的读者。通过该脚本读者还可自行修改节点参数或网络拓扑进一步探索不同负荷、线路阻抗条件下的潮流特性。1. 为什么潮流计算的默认考题是IEEE33节点又为什么默认解法是NR法做配电网分析从IEEE33节点系统开始最省心。这个33节点配电网算例规模小、标准结果齐全无论是算网损、看电压分布还是做光伏接入和故障重构它都是默认的公共试验场。潮流计算的核心任务是在给定负荷下求出网络电压和功率分布NR法Newton-Raphson牛顿-拉夫逊法用雅可比矩阵逐轮逼近非线性功率方程在IEEE33这类小型系统上通常4~8次迭代就能收敛。但配电网线路R/X比大NR法的收敛半径问题不能忽视工程里经常要配合潮流计算最优因子法才能保证不翻车。下面直接围绕IEEE33 NR法这条线把建模、代码、收敛控制与结果校验一次讲透适合正在做配电网规划、微电网和分布式电源接入的从业者。2. 把IEEE33节点变成NR法能吃的输入拓扑、参数与极坐标方程2.1 节点-支路表是第一步IEEE33的拓扑与两套编号IEEE33节点系统由1个平衡根节点、32个PQ节点和32条运行支路组成基准电压12.66kV基准功率10MVA。公开资料里通常还列出5条联络开关支路默认断开所以潮流计算时网络是辐射状。全部负荷集中在节点上总负荷3.715MWj2.3Mvar这个总量可以用来核对数据抄写是否正确任何一张标准IEEE33数据集负荷相加都应该落在这一组数字附近。做输入前先确定编号规则。有一版常用数据把根节点编号为0从0到32共33个节点另一版从1号开始编号。两种编号在网络拓扑表里的支路起止节点完全不一样拿到数据后要先把根节点、负荷节点标注出来。我一般先画一张单线图再在程序里把索引写死避免后面对应错位。如果只拿一串CSV就开始跑最容易犯的错就是“看结果有点像但细看全错”。IEEE33的32条运行支路可以按四条路径来记主线0-1-2-...-17共17条从节点1伸出去的左侧分支1-18-19-20-21共4条从节点2伸出去的右侧分支2-22-23-24共3条从节点4伸出去的末端分支4-25-26-27-28-29-30-31-32共8条。合起来正好32条这也是为什么很多支路表看起来像四块拼起来而不是一条直线。分支路径起止节点支路数主线0-1717左侧分支1-214右侧分支2-243末端分支4-328这个分支结构有两个作用一是手填支路时可以按路径逐段检查二是后面做分布式电源接入时可以快速判断某个节点属于哪条馈线分支。注意联络开关支路不在上表里通常文献会额外列出5条编号在不同资料里有差异默认都断开。2.2 极坐标功率方程和雅可比矩阵NR法计算潮流的三个关键量NR法在电网潮流里常用极坐标形式。对于节点i注入功率的实部P_i和虚部Q_i可以写成P_i V_i * sum_j V_j (G_ij cosθ_ij B_ij sinθ_ij)Q_i V_i * sum_j V_j (G_ij sinθ_ij - B_ij cosθ_ij)其中θ_ij θ_i - θ_jG和B是节点导纳矩阵的实部和虚部。IEEE33系统里平衡节点0的电压幅值固定为1.0p.u.、相角固定为0其余32个PQ节点的状态量是电压幅值V和相角θ总共64个未知量。NR迭代要解的是[ΔP; ΔQ] J [Δθ; ΔV]J是2×2分块雅可比矩阵四块分别是有功对相角、有功对幅值、无功对相角、无功对幅值的偏导。写代码时按节点填充对角和非对角的公式要严格区分H_ii -Q_i - B_ii V_i^2H_ij V_i V_j (G_ij sinθ_ij - B_ij cosθ_ij)N_ii P_i/V_i V_i G_iiN_ij V_i (G_ij cosθ_ij B_ij sinθ_ij)K_ii P_i - G_ii V_i^2K_ij -V_i V_j (G_ij cosθ_ij B_ij sinθ_ij)L_ii Q_i/V_i - V_i B_iiL_ij V_i (G_ij sinθ_ij - B_ij cosθ_ij)这些公式就是后面代码里雅可比矩阵的填充依据不要凭印象写符号错一个电压结果就会往错误方向跑。另外要注意IEEE33原始系统全是PQ节点没有PV节点如果后续把某个节点改成分布式光伏并网点就需要在雅可比矩阵里删去该节点的无功失配行补上电压幅值约束行这一步NR法比前推回代法灵活得多。2.3 为什么NR法仍然值得在IEEE33上用配电网潮流还有一种更常见的前推回代法利用辐射状拓扑从末端回推电流、从根节点前推电压一次往返就能更新一轮编程简单、收敛也快。但它有个硬前提网络必须是辐射状很难处理PV节点和多源网络。NR法没有这个限制任意拓扑都能处理雅可比矩阵还天然提供灵敏度信息后面做DG接入、N-1校核或者电压无功优化都可以直接复用这个矩阵。代价是NR法的雅可比矩阵每轮都要重算33节点规模不明显到几百节点时计算成本会上升。但对IEEE33这个量级一次矩阵求解在毫秒级运算时间完全不是瓶颈这也是大家愿意拿它当教学算例的原因。真正要花心思的不是“能不能算”而是“怎么让它在重负荷下稳定收敛”这一点放到第4章专门说。3. 用Python跑通33节点潮流计算一个只用numpy的NR实现3.1 先准备IEEE33数据支路阻抗和节点负荷下面是一份可复制的数据准备代码把支路表整理成[i, j, R, X]四列负荷表按节点索引存放P(kW)和Q(kvar)。根节点编号用0这是最常见的一版编号规则。import numpy as np # IEEE33节点系统根节点编号0 # 支路: [送端, 受端, R(ohm), X(ohm)] branch 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], [4, 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] ]) # 负荷: [P(kW), Q(kvar)]节点0为根节点 load np.array([ [0, 0], [100, 60], [90, 40], [120, 80], [60, 30], [60, 20], [200, 100], [200, 100], [60, 20], [60, 20], [45, 30], [60, 35], [60, 35], [120, 80], [60, 10], [60, 20], [60, 20], [90, 40], [90, 40], [90, 40], [90, 40], [90, 40], [90, 50], [420, 200], [420, 200], [60, 25], [60, 25], [60, 20], [120, 70], [200, 600], [150, 70], [210, 100], [60, 40] ])这段代码里的支路表覆盖了主线、左侧分支、右侧分支和末端分支一共32条运行支路5条联络开关支路没有放进来默认开环。负荷数组第0行是根节点负荷为0后面每一行对应一个PQ节点的有功和无功单位是kW和kvar这些值会在下一小节换算成标幺值。如果你拿到的数据是1号起编记得先把根节点移到下标0再填负荷数组否则后面的雅可比矩阵全是错位。3.2 构建节点导纳矩阵和功率函数继续写代码把网络参数变成标幺值并构建节点导纳矩阵Y然后定义功率计算函数。Vbase 12.66 # kV Sbase 10.0 # MVA Zbase Vbase**2 / Sbase # 欧姆 n 33 # 节点导纳矩阵单位是标幺值 Y np.zeros((n, n), dtypecomplex) for f, t, R, X in branch: y 1.0 / ((R 1j*X) / Zbase) Y[f, f] y Y[t, t] y Y[f, t] - y Y[t, f] - y G Y.real B Y.imag def calc_pq(V, theta): dtheta theta[:, None] - theta[None, :] P V * (G * np.cos(dtheta) B * np.sin(dtheta)) V Q V * (G * np.sin(dtheta) - B * np.cos(dtheta)) V return P, Q逻辑说明支路阻抗先除以Zbase得到标幺值再取倒数得到支路导纳Y矩阵对角线累加自导纳非对角线累加互导纳这是标准节点导纳矩阵规则。calc_pq函数用dtheta构造出所有节点相角差的矩阵然后用矩阵乘法和电压数组完成全网络有功、无功计算比用双重循环逐节点累加更短也更容易发现维度错误。参数说明Zbase约等于16.03Ω也就是12.66kV的平方除以10MVA。Sbase用10.0而不是100因为IEEE33标准算例指定的就是10MVA如果按一般输电网习惯用100MVA所有标幺值都会缩小10倍结果看起来会像“没收敛”。写程序时把Vbase、Sbase放在数据准备区后续换算都是同一个常量避免在calc_pq里临时除。3.3 NR法主循环失配量、雅可比矩阵和电压更新下面是最核心的NR迭代部分这里先给出不带阻尼的原始版本方便对比后面加最优因子法的效果。# 给定功率负荷取负值 P_spec np.zeros(n) Q_spec np.zeros(n) P_spec[1:] -load[1:, 0] / (Sbase * 1000.0) Q_spec[1:] -load[1:, 1] / (Sbase * 1000.0) V np.ones(n) # 平启动单位p.u. theta np.zeros(n) # 相角单位rad tol 1e-8 # 标幺值失配功率阈值 max_iter 20 m n - 1 for k in range(max_iter): P, Q calc_pq(V, theta) dP P_spec[1:] - P[1:] dQ Q_spec[1:] - Q[1:] mis np.concatenate([dP, dQ]) if np.max(np.abs(mis)) tol: print(converged, iter, k 1) break J np.zeros((2 * m, 2 * m)) dtheta theta[:, None] - theta[None, :] for i in range(m): ni i 1 for j in range(m): nj j 1 if ni nj: J[i, j] -Q[ni] - B[ni, ni] * V[ni]**2 J[i, m j] P[ni] / V[ni] V[ni] * G[ni, ni] J[m i, j] P[ni] - G[ni, ni] * V[ni]**2 J[m i, m j] Q[ni] / V[ni] - V[ni] * B[ni, ni] else: s np.sin(theta[ni] - theta[nj]) c np.cos(theta[ni] - theta[nj]) J[i, j] V[ni] * V[nj] * (G[ni, nj] * s - B[ni, nj] * c) J[i, m j] V[ni] * (G[ni, nj] * c B[ni, nj] * s) J[m i, j] -V[ni] * V[nj] * (G[ni, nj] * c B[ni, nj] * s) J[m i, m j] V[ni] * (G[ni, nj] * s - B[ni, nj] * c) dx np.linalg.solve(J, mis) theta[1:] dx[:m] V[1:] dx[m:] print(Vmin, V.min(), Vmax, V.max()) print(P_loss(MW), (P[0] P_spec.sum()) * Sbase)这段代码的收敛判据是max(|ΔP|, |ΔQ|)小于1e-8这是比较严格的要求实际工程用1e-6已经足够因为标幺值1e-6对应的有名功率只有约10W。雅可比矩阵直接用解析偏导填充64×64规模很小np.linalg.solve的耗时可以忽略。如果代码正确IEEE33原始负荷下落点电压最低值通常在0.90p.u.附近不会低于0.85总网损在0.2MW上下。如果你跑出来的最低电压只有0.6p.u.大概率是负荷没除以Sbase而不是NR法本身的问题。原始NR法在平启动下通常能收敛但一旦修改支路参数或加重负荷就可能出现失配量反复震荡这正好引出下一章的最优因子法。4. 收敛不行就上潮流计算最优因子法原理与代码改造4.1 为什么配电网会让NR法翻车R/X比与初值敏感NR法的收敛性对初值非常敏感。输电网电抗远大于电阻雅可比矩阵条件数好平启动就能稳稳收敛配电网线路R/X比高IEEE33里很多支路电阻大于电抗矩阵对角占优性变差。这时候如果末端负荷很重电压真实解已经偏离1.0p.u.很远全步长修正很容易越过真实解表现为失配量先变小再变大电压出现负值或超过1.1p.u.。我实际跑过不少33节点改造模型把某条支路R放大或把末端负荷加到2倍原始NR法就开始发散。这种发散不是代码bug而是牛顿法收敛半径问题。遇到这种情况第一反应不要是去改收敛阈值那只会得到假收敛应该给修正量加一个自适应步长这就是潮流计算最优因子法要解决的场景。4.2 潮流计算最优因子法是什么给修正量乘一个标量μ最优因子法的思想很直观。NR法的修正量Δx来自线性化方程全步长xΔx只对线性模型最优非线性模型下这个步长可能太大。最优因子法把更新改成x^{k1} x^k μ Δx其中μ取某个标量让失配量函数F(μ) || f(x^k μΔx) ||^2尽量小。当μ1时就是原始NR法μ1时相当于阻尼步长μ1时可以加速收敛。精确求μ需要解一个高阶多项式工程上很少这么做更常见的做法是用回溯线搜索先试μ1如果更新后的失配量范数比更新前更大就把μ折半直到失配量下降或达到折半上限。4.3 代码改造加入最优因子后的NR主循环把上一章的更新部分替换成带线搜索的版本其余代码不变。# 在每次求得dx后执行替换 theta[1:] dx[:m] 那两行 max_mu_try 12 mu 1.0 old_norm np.linalg.norm(mis) for t in range(max_mu_try): V_can V.copy() theta_can theta.copy() V_can[1:] mu * dx[m:] theta_can[1:] mu * dx[:m] P_can, Q_can calc_pq(V_can, theta_can) new_mis np.concatenate([P_spec[1:] - P_can[1:], Q_spec[1:] - Q_can[1:]]) if np.linalg.norm(new_mis) old_norm: V, theta V_can, theta_can break mu * 0.5 else: print(line search failed at iter, k 1) break说明这段代码先尝试全步长只有失配量范数确实下降才接受否则把步长μ不断减半。max_mu_try12意味着最小步长约为2^-12足以应对绝大多数重负荷场景。如果连续多次线搜索失败或μ长期小于0.2说明问题大概率在数据本身而不是算法继续压步长只会让收敛变慢。和原始NR法相比这个改造增加的计算量在最坏情况下是每次迭代多算12次功率方程但对IEEE33来说每次calc_pq都是64维矩阵运算开销完全可以接受。收敛轮数可能会从4轮增加到6轮但换来的是不再随机发散这比那一点额外计算重要得多。4.4 最优因子法的三个使用注意第一不要无脑把所有场景都开成强阻尼。正常负荷下原始NR法4轮就收敛加了线搜索后如果μ一路保持1.0不影响结果但如果初始电压给得特别差μ连续多轮折半收敛速度会明显下降。这时可以改用上一轮潮流解作为初值而不是每次都从1.0p.u.平启动。第二最优因子法不能修复雅可比矩阵奇异。如果节点编号断链、支路阻抗填成0或者某个节点孤立Ybus对应行全为0np.linalg.solve依然会报singular matrix。线搜索只处理“方向对但步长太大”处理不了“方向本身错误”。第三μ折半上限要跟收敛阈值匹配。max_mu_try12对应最小步长约0.000244如果这个步长下失配量仍然不降那就不是步长问题应该停止迭代并检查数据。盲目把上限提高到50只会多算几十次calc_pq对结果没有帮助。5. IEEE33潮流计算避坑清单编号、基准值、收敛条件与拓扑5.1 节点编号从0还是1起编最容易引起结果对不上现象跑出来的电压分布趋势对但具体节点电压和公开参考值差一格越到末端越乱。原因同一份IEEE33单线图有人把根节点标0有人标1从1起编的数据里“支路1-2”对应0起编的“0-1”错位一个节点。解决写代码前先把源数据的路径画出来明确根节点索引我一般直接在支路表注释里写“0是变压器低压侧出口”再让程序打印每个节点电压肉眼核对首末端。注意负荷数组也要同步平移曾经有人只改了支路编号负荷还留在1号起编的位置结果总负荷分布完全错了。5.2 基准值选错12.66kV/10MVA背后的单位换算现象所有电压算出来只有1e-5量级或者网损几乎为0。原因把负荷kW直接当成标幺值或者支路阻抗没有除以Zbase。IEEE33的基准电压是12.66kV基准功率常用10MVAZbase约等于16.03Ω。负荷换算要先除以10MVA支路阻抗要先除以Zbase。解决在代码里保留Vbase、Sbase常量换算写在数组前面不要在calc_pq函数里临时处理单位。另一个常见误用是把Sbase设成100MVA结果电压剖面看起来“偏低但不离谱”这种错误最难排查建议每次打印总负荷标幺值确认它等于0.3715j0.2300。5.3 收敛条件只看ΔP不看ΔQ会得到假收敛现象迭代3轮就达到收敛但Q失配量还停在1e-2量级无功和网损明显不对。原因收敛判据只取了有功失配而PQ节点是P和Q两个方程同时要满足。解决用max(|ΔP|, |ΔQ|)合并判断且对标幺值取1e-6以下如果还关心电压稳定性再要求最大电压修正量小于1e-5。这一条在重负荷场景下特别重要因为无功失配往往比有功失配收敛得慢只看有功会把未收敛的解当成结果。5.4 5条联络开关不加区分33节点就变成闭环网络现象直接用文献里的完整33节点数据把5条联络开关也放进Ybus潮流结果与标准辐射状算例相差很大。原因IEEE33节点系统原始数据包含5条常开联络支路默认开环运行只有做网络重构时才逐个合上。解决支路表先只放32条运行支路把5条联络支路单独放一个数组需要时再追加。如果代码里Ybus是循环遍历branch构建的那就很容易把联络开关误算进去建议构建前打印branch.shape确认是32行而不是37行。5.5 雅可比矩阵奇异或迭代爆炸先查数据再怀疑算法现象np.linalg.solve报singular matrix或电压更新后出现负值。原因常见是节点编号断链、支路R或X填成0、某个节点负荷大到让该节点电压无解而不是NR法本身的锅。解决先打印Ybus的行非零个数确认每个节点都挂在网络上再用一个把各节点负荷乘0.1的轻载算例跑通确认算法稳定后再恢复全负荷。这样可以把数据问题和算法问题快速切开不用靠猜。轻载算例如果也发散那就是编号或雅可比矩阵代码问题轻载收敛但满载发散再去考虑最优因子法和初值。6. 结果可信度验证功率不平衡检查与前推回代交叉验证6.1 三个自检指标收敛后不要急着把V数组写进报告先用三个量做自检。第一个是根节点功率与全系统负荷的差值这个值就是总网损IEEE33在原始负荷下网损大约在0.2MW量级如果算出来是几MW要么是收敛判据太松要么是负荷单位错。第二个是节点电压范围正常应在0.90~1.0p.u.之间最低压节点通常在长分支末端任意节点超过1.05就说明潮流解有问题。第三个是最后一次迭代的失配量最大值打印max(|dP|, |dQ|)应该低于1e-6。这三项可以直接写成一个check函数每次跑完自动打印避免人工翻日志。6.2 用前推回代法做交叉验证NR法实现了不代表雅可比矩阵没有填错最稳妥的验证是用另一个算法交叉对比。IEEE33是辐射状网络前推回代法实现起来很简单从末端节点开始用节点功率和当前电压回代支路电流再从根节点向前推算各节点电压反复迭代到电压差小于阈值。把两种算法算出的末端电压放在一张表里对比偏差小于1e-4 p.u.就可以基本确认实现正确。我第一次调通NR法时就是靠前推回代交叉验证发现了支路表里一条R/X数据填反比对着参考值排查快得多。6.3 把最优因子法推广到DG接入场景如果后面要做分布式电源接入把某个负荷节点改成PV节点NR法的雅可比矩阵会多一行电压幅值约束这时最优因子法的价值更明显DG出力从0逐步增加到额定值每一步用上一轮的解做初值配合μ阻尼一般不会发散。我现在的习惯是跑任何配电网潮流都默认带最优因子开关先试一次μ1再允许线搜索介入这比反复调初值省事也少了很多玄学调参。希望帮到你。本文还有配套的精品资源点击获取