简介潮流计算是电力系统规划、调度与安全校核的基础这份资源以一个9节点小型电网为对象提供了Python仿真练习包适合电力专业学生、能源行业工程师以及对电网分析感兴趣的入门者使用。压缩包内共三个文件核心是“9节点计算成品.py”脚本另附data.txt与data2.txt两个数据文件分别存储系统拓扑、线路阻抗、发电机出力和负荷需求等参数整体仅4KB轻巧易读方便直接运行和修改。脚本在建模时会包含发电机、负荷、变压器与线路等基本元件通过迭代求解非线性方程组输出各节点电压、相角和支路功率同时校验功率平衡与电压上限等约束。配合数据文件读者能够看清从数据解析到模型求解、再到结果分析的全过程并以此复现潮流计算实验或作为算法调试的基准。目前已有149人学习下载对于希望将Python应用到能源电力分析的学习者这套精简示例具有清晰的参考价值。1. 潮流计算到底在算什么做过配网规划或电网分析的工程师都清楚潮流计算不只是课本里的节点电压方程而是任何一张电网图纸从画出来到能运行之间必须跨过的那道门槛。给定发电机出力和负荷需求求解各母线电压幅值与相角、线路有功与无功分布听起来像线性代数作业但真实电网的约束条件是非线性的靠手算根本推不动几节点以上的网络。Python 在这个领域的价值在于NumPy 的矩阵运算让雅可比矩阵组装和稀疏求解变得极其顺手Pandas 能优雅地管理母线、支路和发电机三类表结构而 Matplotlib 又能把收敛过程直接可视化。这篇文章不打算介绍某个现成库的 API 用法而是从第一性原理出发用一套最小可运行的代码实现牛顿-拉夫逊法潮流计算并覆盖参数设置、收敛性判断、误差排查这些实际工程中必然遇到的问题。适合阅读本文的读者包括刚入门电力系统分析、想验证课本公式的在校学生需要快速搭建电网计算原型的算法工程师以及想知道自己拿到的潮流数据包里每一列到底在干什么的能源行业从业者。读完你应该能对一份陌生的潮流计算代码做出评估和修改而不仅仅是能跑。2. 潮流计算的数学模型与求解算法选型2.1 节点功率平衡方程的构建逻辑潮流计算的起点是每个节点的功率平衡约束。对于节点 i注入的有功功率和无功功率必须满足Pi Vi * Σ(Vj * (Gij * cosθij Bij * sinθij)) Qi Vi * Σ(Vj * (Gij * sinθij - Bij * cosθij))其中 Gij 和 Bij 是节点导纳矩阵 Y 的实部与虚部θij 是节点 i 和 j 的电压相角差。这个方程组里每个节点的未知量是电压幅值 Vi 和相角 θi而 Pi 和 Qi 在 PQ 节点上是已知量。整套方程组呈现高度非线性原因在于 Vi 与 Vj 的乘积项以及三角函数耦合。提示解析求解这个方程组是不现实的工程上全部采用迭代法逼近数值解。为什么不能像解线性方程组那样直接求逆因为雅可比矩阵本身依赖于当前迭代点的电压值每轮迭代都需要重新计算和分解。这意味着我们面对的不是一次性求解而是一个反复线性化-求解-更新的过程。理解了这一点就能明白为什么潮流计算的性能瓶颈往往不在迭代次数上而在每次迭代的矩阵组装和稀疏线性求解上。2.2 三种主流算法的适用边界电力系统潮流计算的算法谱系大致分为三类高斯-赛德尔法、牛顿-拉夫逊法和 PQ 分解法。高斯-赛德尔法实现最简单用节点电压方程直接迭代更新内存占用极低但收敛速度是线性的遇到重负荷系统时常需要几百次迭代甚至发散当前主要用于教学演示。牛顿-拉夫逊法是目前工业界的绝对主流以二次收敛速度著称。它把非线性方程组 F(x) 0 在每次迭代点做泰勒展开忽略二阶以上项得到一个线性修正方程组 J·Δx -F(x)。对于 1000 节点以下的中小型系统5 到 8 次迭代就能达到 1e-6 的精度。代价是需要每轮组装雅可比矩阵并做一次 LU 分解。PQ 分解法Fast Decoupled Load Flow利用了输电系统中有功功率-相角强耦合、无功功率-电压强耦合的物理特性将雅可比矩阵简化为两个常数矩阵 B 和 B只需在迭代开始前列式分解一次。3200 节点以上的大规模输电网中它的单次迭代耗时只有牛顿法的五分之一但迭代次数会增加到 15 到 25 次。对于配电网这类 R/X 比值较高的网络PQ 分解法的简化假设不再成立收敛性反而变差。2.3 为什么 Python 是搭建原型的合适选择有人会质疑 Python 的性能能否承担潮流计算任务。实际工程中我一般这样分层快速验证算法、教学演示、处理中小规模配电网数据用 Python 完全够用一旦到了在线调度、安全分析这类对单次计算时间苛刻的场景C 或 Fortran 仍是主力。不过Python 生态的真正优势体现在数据清洗和结果可视化上——用 Pandas 读变电站台账、合并负荷曲线、输出 Excel 报表这些是 C 需要写几百行代码才能达到的效果。追求性能时也无需立刻放弃 Python。可以先用纯 Python 跑通逻辑再用 Numba 的 jit 装饰器做即时编译或者用 scipy.sparse 库的稀疏矩阵存储导纳矩阵。对一个 118 节点的 IEEE 标准系统稀疏化之后的内存占用能降到全稠密的百分之一以下。这套先正确、再加速的路线恰恰是 Python 在电力计算领域被广泛接受的核心原因。3. 用 Python 从零实现牛顿-拉夫逊潮流的完整代码3.1 数据模型用 Pandas 定义母线、支路和发电机表动手写求解器之前先把数据模型定下来。电网数据本质上是三类表对象母线表描述节点参数支路表描述拓扑连接关系发电机表描述注入功率。我推荐直接仿照 MATPOWER 的案例格式来组织数据这能让你在后续对比验证时无缝衔接标准算例。import pandas as pd import numpy as np # 母线表bus_id, 类型(1PQ, 2PV, 3平衡), 有功负荷, 无功负荷 bus_data { bus_id: [1, 2, 3], type: [3, 2, 1], # 节点1为平衡节点节点2为PV节点节点3为PQ节点 Pd: [0.0, 0.0, 1.25], # 有功负荷单位pu标幺值 Qd: [0.0, 0.0, 0.55], # 无功负荷 V0: [1.0, 1.0, 1.0], # 电压幅值初值 theta0: [0.0, 0.0, 0.0] # 相角初值单位弧度 } bus pd.DataFrame(bus_data) # 支路表from, to, 电阻, 电抗, 电纳 branch_data { fbus: [1, 1, 2], tbus: [2, 3, 3], r: [0.02, 0.05, 0.035], x: [0.06, 0.15, 0.085], # 电抗远大于电阻符合输电网特征 b: [0.03, 0.025, 0.02] # 对地电纳单位S } branch pd.DataFrame(branch_data) # 发电机表母线号, 有功出力, 电压设定值 gen_data { bus_id: [1, 2], Pg: [0.65, 0.55], # 发电机有功出力单位pu Vg: [1.02, 1.01] # PV节点和平衡节点的电压幅值设定 } gen pd.DataFrame(gen_data)数据格式里最关键的是节点类型划分平衡节点type3承担功率缺额电压幅值和相角固定PV 节点type2固定有功出力和电压幅值无功出力待求PQ 节点type1的有功和无功负荷均为已知量。初值的选取对收敛速度有很大影响工程上一般取平启动即所有 PQ 节点电压幅值为 1.0相角为 0。3.2 导纳矩阵组装稀疏存储与向量化计算的配合导纳矩阵是后续所有计算的地基。其对角线元素 Yii 等于连接在节点 i 上的所有支路导纳之和非对角线元素 Yij 等于支路 i-j 的导纳取负。这里的导纳是复数实部为电导 G虚部为电纳 B。from scipy.sparse import lil_matrix def build_ybus(bus, branch): n len(bus) # 用LIL格式逐项填充最后转CSR提升计算效率 Y lil_matrix((n, n), dtypenp.complex128) for _, line in branch.iterrows(): f int(line[fbus]) - 1 t int(line[tbus]) - 1 # 支路导纳 y 1 / (r jx) y 1.0 / complex(line[r], line[x]) # 对地电纳作为附加导纳并联在两端节点上 ysh complex(0, line[b] / 2) # 更新自导纳和互导纳 Y[f, f] y ysh Y[t, t] y ysh Y[f, t] - y Y[t, f] - y return Y.tocsr() Ybus build_ybus(bus, branch) print(导纳矩阵实部 G:\n, Ybus.real.toarray()) print(导纳矩阵虚部 B:\n, Ybus.imag.toarray())注意这里用到了 scipy.sparse 的稀疏矩阵原因是真实电网中每个节点通常只与 2 到 5 条支路相连导纳矩阵的稀疏度超过 95%。代码中对地电纳 b 除以 2 是 π 型等值电路的常规处理方式支路两端各并联一半的电纳。如果你拿到的手工数据包里 b 已经是总电纳值务必做同样处理否则会导致无功计算明显偏离预期。3.3 牛顿-拉夫逊迭代主体从功率不平衡量到修正方程迭代的核心思路是反复计算当前电压状态下的功率注入与给定值做差得到不平衡量 ΔP 和 ΔQ再解修正方程得到电压幅值和相角的增量。实现时我习惯把功率计算和雅可比矩阵组装放在同一个函数里因为两者共享中间变量拆分反而会增加重复计算。def newton_raphson(bus, branch, gen, max_iter20, tol1e-8): Ybus build_ybus(bus, branch) n len(bus) # 初始化电压向量 V bus[V0].to_numpy(dtypenp.float64) * np.exp(1j * bus[theta0].to_numpy()) # 构建已知量数组平衡节点不参与迭代PV节点不参与无功方程 pq_idx bus[bus[type] 1].index.tolist() pv_idx bus[bus[type] 2].index.tolist() slack_idx bus[bus[type] 3].index.tolist() # 发电机有功注入按母线号映射 Pg np.zeros(n) gen_bus_map dict(zip(gen[bus_id] - 1, gen[Pg])) for idx, pg in gen_bus_map.items(): Pg[idx] pg # P和Q的给定值发电机注入减去负荷 P_spec Pg - bus[Pd].to_numpy() Q_spec -bus[Qd].to_numpy() # 无功给定初始视为零PV节点后续修正 for it in range(max_iter): # 计算当前注入功率 S V * np.conj(Ybus V) P_calc S.real Q_calc S.imag # 计算不平衡量PV节点和平衡节点跳过无功方程 dP P_spec - P_calc dQ Q_spec - Q_calc # 构建误差向量PQ节点含dP和dQPV节点只含dP mismatch np.concatenate([ dP[pq_idx pv_idx], dQ[pq_idx] ]) if np.max(np.abs(mismatch)) tol: print(f第 {it1} 次迭代收敛) break # 组装雅可比矩阵 J build_jacobian(Ybus, V, pq_idx, pv_idx) # 求解修正方程 dx np.linalg.solve(J, mismatch) # 分解修正量前三部分对应相角后面对应幅值 n_theta len(pq_idx) len(pv_idx) dtheta np.zeros(n) dtheta[pq_idx pv_idx] dx[:n_theta] dV_abs np.zeros(n) dV_abs[pq_idx] dx[n_theta:] # 更新电压PV节点幅值保持设定值 V V * np.exp(1j * dtheta) dV_abs * np.exp(1j * np.angle(V)) for g_idx in gen_bus_map: V[g_idx] abs(V[g_idx]) * np.exp(1j * np.angle(V[g_idx])) return V这段代码的数值处理有两点值得注意。第一相角修正量 dtheta 用在指数项里本质上是在做极坐标形式的更新比直接加法更稳定能避免大相角变化时因三角函数展开引入的误差第二PV 节点电压幅值在每轮迭代后被强制拉回设定值这个投影操作保证了最终结果满足恒定电压幅值的物理约束。如果不做这一步虚假的无功不平衡会导致 PV 节点电压逐步漂移。3.4 雅可比矩阵的解析表达式与组装细节雅可比矩阵的组装比功率计算稍微费力一些。根据极坐标下的功率方程对相角求偏导时非对角元素为 Vi·Vj·(Gij·sinθij - Bij·cosθij)对电压幅值求偏导时非对角元素为 Vi·(Gij·cosθij Bij·sinθij)。对角元素则涉及本节点自导纳和当前注入功率写法略有不同。def build_jacobian(Ybus, V, pq_idx, pv_idx): n len(V) n_pq len(pq_idx) n_pv len(pv_idx) n_theta n_pq n_pv J np.zeros((n_theta n_pq, n_theta n_pq)) # 提取导纳矩阵实虚部 G Ybus.real.toarray() B Ybus.imag.toarray() # 步骤1计算 H 子块 (dP/dθ) # 步骤2计算 N 子块 (dP/dV) # 步骤3计算 M 子块 (dQ/dθ) # 步骤4计算 L 子块 (dQ/dV) for i in range(n): for j in range(n): if i j: continue # 非对角元素的通用公式 dP_dtheta V[i] * V[j] * (G[i][j] * np.sin(np.angle(V[i]) - np.angle(V[j])) - B[i][j] * np.cos(np.angle(V[i]) - np.angle(V[j]))) # 按行映射到压缩索引 if i in pq_idx pv_idx: row_theta (pq_idx pv_idx).index(i) if j in pq_idx pv_idx: col_theta (pq_idx pv_idx).index(j) J[row_theta, col_theta] dP_dtheta if j in pq_idx: col_v n_theta pq_idx.index(j) J[row_theta, col_v] V[i] * (G[i][j] * np.cos(np.angle(V[i]) - np.angle(V[j])) B[i][j] * np.sin(np.angle(V[i]) - np.angle(V[j]))) # 对角元素需要加上本节点注入功率项 for i in pq_idx pv_idx: row (pq_idx pv_idx).index(i) Si V[i] * np.conj(np.sum(Ybus[i, :] V)) J[row, row] -V[i]**2 * B[i][i] - Si.imag # H对角 return J上述组装逻辑是完整的但循环方式在节点规模较大时效率偏低。实际工程中更推荐逐支路扫描的组装方式先初始化全零稠密矩阵遍历每条支路把四个端点的雅可比贡献一次性填入。这样做的好处是每轮迭代只扫描 O(branches) 次而非 O(n²)且与导纳矩阵的构建方式天然契合。代码中的双层循环仅为展示数学结构便于理解偏导关系的来源。4. 数据准备、精度控制与收敛性调优4.1 标幺值体系从有名值到无量纲的转换电力系统计算中最容易出错的一步不是算法选型而是单位换算。工程实际中的数据来自多种渠道变电站台账给出有名值的阻抗欧姆负荷预测系统给出兆瓦和兆乏而发电机参数表可能是以设备自身容量为基准的标幺值。在进入潮流计算前必须统一折算到同一个基准容量和基准电压之下。# 有名值转标幺值的示例 S_base 100e6 # 基准容量100MVA V_base 110e3 # 基准电压110kV Z_base V_base**2 / S_base # 基准阻抗 # 线路阻抗有名值欧姆 r_ohm 2.5 x_ohm 8.2 # 转标幺值 r_pu r_ohm / Z_base x_pu x_ohm / Z_base # 负荷从MW/MVar转pu P_pu 25e6 / S_base # 25MW负荷 Q_pu 9e6 / S_base # 9MVar无功负荷 print(f基准阻抗: {Z_base:.2f} 欧姆) print(f线路电阻标幺值: {r_pu:.4f}) print(f无功负荷标幺值: {Q_pu:.4f})使用标幺值的好处是数值范围适中雅可比矩阵的条件数较小收敛性也就更好。如果全部用有名值参与计算不同电压等级下的数值可能横跨 6 个数量级浮点舍入误差会在迭代中逐步放大。基准容量的选择以 100MVA 最常见但在新能源场站并网分析中以单台风机逆变器容量为基准的情况也不少见关键是整个计算域内始终保持统一。4.2 收敛判据的选取与误收敛的识别潮流计算的收敛判断并非只有一个标准答案。绝对误差判据检查不平衡量最大值对大规模系统过于严格相对误差判据又可能在轻载时产生过宽放松。我一般同时监控两个指标最大不平衡功率不超过 1e-6 pu 的绝对阈值同时检查电压幅值更新量的二范数。判断发散的经验法则是迭代过程出现振荡即误差不再单调递减或者直接跳到 NaN。此时应从四个方向排查。第一检查初值将平启动改为热启动利用前一次相近工况的结果通常能解决问题第二检查导纳矩阵尤其要注意变压器变比是否折算进 Y 矩阵第三降低负荷水平从 50% 负荷起步逐步上调定位引起不收敛的临界状态第四检查支路参数中是否存在 R 远大于 X 的异常数据这在配电网中常见需要考虑改用牛顿法内核。4.3 配电网高 R/X 比场景下的参数调整策略经典牛顿-拉夫逊法在输电网中性能优异但面对配电网时线路电阻与电抗之比常常大于 1此时 P-Q 解耦假设失效PQ 分解法基本不可用而牛顿法本身仍然是收敛的只是收敛域变小对初值更敏感。实用的参数调整方法包括对初值电压幅值采用从电源侧向末端逐级外推的方式替代统一平启动将雅可比矩阵的条件数打印出来超过 1e12 时考虑用 LU 分解加主元缩放。针对配电网调压器、分布式光伏接入等场景更可靠的做法是采用前推回代法。该算法利用辐射状拓扑的自底向上递推结构无需组装雅可比矩阵单次迭代耗时极短对 R/X 比值不敏感。如果你的对象是拓扑接近树状的低压配电网不妨将前推回代法作为备选方案只需把导纳矩阵求解替换为逐层功率流递推即可收敛条件从二次收敛降为线性但单步的吞吐量完全弥补了这一点。5. 从收敛曲线到工程应用的验证技巧5.1 绘制误差下降曲线识别收敛瓶颈收敛过程的误差变化轨迹包含大量诊断信息。理想的牛顿法应当呈现二次收敛特征前两轮迭代误差从 1e-2 压到 1e-4再掉到 1e-8 以下。如果误差曲线呈现线性下降说明问题出在初值质量太差如果误差在某个量级徘徊超过 5 轮大概率是 PV 节点无功越限未处理。import matplotlib.pyplot as plt # 修改迭代函数记录每轮的最大误差 errors [] # 在每次迭代计算mismatch后追加errors.append(np.max(np.abs(mismatch))) # 循环结束后绘制 plt.semilogy(errors, o-) plt.xlabel(iterations) plt.ylabel(max mismatch (pu)) plt.title(Newton-Raphson Convergence Curve) plt.grid(True, whichboth, linestyle--, alpha0.7) plt.show()如果收敛曲线在前两三轮正常后面突然变缓甚至上翘排查方向应转向无功越限。实际电力系统的发电机无功出力是有上下限的当某台 PV 节点的无功力需求超出其上限时该节点应转换为 PQ 节点带上 Qmax 作为固定无功注入重新迭代。不处理这个逻辑雅可比矩阵会在错误的节点类型假设下组装导致末端电压失真。5.2 前后校验功率平衡与线路载流量双检查一份可信的潮流结果至少要经过两道校验。第一道校验是全网功率平衡发电机总出力减去总负荷与总网损的差值应为零误差通常在 1e-6 pu 量级。第二道校验是逐条线路的视在功率不超过其容量限值。这两个校验都做过了结果才能放进电网规划报告。功率平衡校验的额外价值在于帮助发现数据录入错误比如某个变压器的分接头挡位填错或者某条双回线路被漏掉一条。很多情况下潮流软件不报错但网损数据离奇偏大回头检查往往就是支路参数录入问题。5.3 迷你算例的数字推演与复盘拿本文用的 3 节点系统来说IEEE 标准的 3 节点算例平衡节点出力约为 0.65 pu我们这里的发电机表给出的 0.55 pu 实际留有调节余量最终计算结果中平衡节点出力会因网损而略高于发电机表中的初设值。看结果时应当关注的不只是电压幅值还有相角差——正常输电网中相邻节点相角差很少超过 15 度如果你算出的相角差达到 30 度以上先检查支路电阻电抗是否发生量级混淆。电流热效应I²R导致的网损通常占总负荷的 2% 到 5%偏离这个区间说明模型的等值阻抗或负荷分布与实际情况存在偏差。当调试方向不再局限于让程序跑通时你会发现潮流计算的真正工程价值在于场景分析反复修改支路投切状态、负荷水平、发电机出力计划观察哪些节点电压逼近安全边界哪些线路面临过载风险。这套 Python 实现的优势是数据驱动——把 SCADA 系统导出的 CSV 替换掉代码里的静态表格你的计算就能跟随电网运行方式的变化自动更新这远比手工维护一份 Excel 计算表更符合业务需求。本文还有配套的精品资源点击获取