简介这份资源面向电力系统工程师、研究人员及电气类学生聚焦用Python实现潮流计算这一电力网络分析的核心技术帮助读者掌握节点电压、支路功率与变压器功率分配的求解方法。压缩包共27个文件约687KB以py脚本、txt说明、png结果图为主另含pyc、sln、pyproj等工程与项目配置文件便于直接运行与二次开发。内容围绕网络建模、求解器选择、KCL与KVL方程构建、迭代收敛及结果可视化展开并涉及牛顿-拉弗森法、高斯-塞德尔迭代等算法思路以及无功补偿、动态稳定性等进阶方向。已有1332人学习读者可借助示例脚本与案例数据理解NumPy、SciPy、Pandas在电力计算中的配合方式快速搭建可复用的潮流计算流程为后续优化与科研实践提供参考。1. 用 Python 把电力系统潮流计算跑起来从“能算”到“算得对”的分水岭很多做电气仿真的朋友第一次用 Python 写潮流计算都会经历同一个心理落差代码明明收敛了输出却和商业软件对不上。问题往往不在编程而在对模型和算法的理解。潮流计算本质上是求解一组非线性代数方程核心是节点功率平衡。Python 的价值在于把这份“黑匣子”打开让你能改导纳矩阵、换迭代方法、看每一步残差。它适合三类人想搞懂算法细节的学生、需要批量扫描工况的工程师、以及要把潮流嵌入优化或机器学习流程的开发者。这一篇不讲空泛概念只讲怎么用 NumPy 和 SciPy 把牛顿-拉夫逊法写对以及那些让结果“差一点”的坑。2. 从节点导纳矩阵到功率方程潮流计算的数学骨架2.1 节点导纳矩阵到底怎么组装电力网络的拓扑结构最终都浓缩到一个复数矩阵 ( Y_{bus} ) 里。对于 ( n ) 个节点的系统它是一个 ( n \times n ) 的对称复数矩阵。对角元素 ( Y_{ii} ) 是连接到节点 ( i ) 的所有支路导纳之和非对角元素 ( Y_{ij} ) 是节点 ( i ) 和 ( j ) 之间支路导纳的负值。如果节点之间有变压器还要考虑变比带来的非对称修正。我一般用三个列表来描述网络支路首端、支路末端、支路阻抗。这样组装起来最不容易出错。下面这段代码演示了从支路表生成 ( Y_{bus} ) 的最小实现。import numpy as np def build_ybus(branch_from, branch_to, branch_z, shunt_yNone): branch_from, branch_to: 支路首末端节点编号从0开始 branch_z: 支路复阻抗 shunt_y: 节点对地并联导纳可选 n max(max(branch_from), max(branch_to)) 1 ybus np.zeros((n, n), dtypecomplex) y_series 1.0 / np.array(branch_z) # 支路串联导纳 for k in range(len(branch_from)): i branch_from[k] j branch_to[k] ybus[i, i] y_series[k] ybus[j, j] y_series[k] ybus[i, j] - y_series[k] ybus[j, i] - y_series[k] if shunt_y is not None: for i, y in enumerate(shunt_y): ybus[i, i] y return ybus逻辑说明先根据最大节点编号确定矩阵维度然后逐条支路累加。串联导纳加到对角负值加到非对角。并联导纳直接加在对角。参数上要注意branch_z必须是复数数组如果只给电阻会丢掉电抗导致结果完全不对。节点编号从 0 开始如果原始数据从 1 开始记得先减 1。2.2 功率方程与雅可比矩阵的构造节点功率方程是 ( S_i V_i \sum_{j} Y_{ij}^* V_j^* )。牛顿-拉夫逊法需要雅可比矩阵它由四个子块组成( H \partial P / \partial \theta )( N \partial P / \partial V )( M \partial Q / \partial \theta )( L \partial Q / \partial V )。很多手写代码在这里翻车因为极坐标下的偏导公式容易抄错。我习惯用极坐标形式把电压写成 ( V_i V_i \angle \theta_i )。雅可比矩阵的维度是 ( 2n_{PQ} n_{PV} )。下面给出一个经过验证的雅可比构造片段。def build_jacobian(Ybus, V, theta, P_spec, Q_spec, pq, pv): Ybus: 节点导纳矩阵 V, theta: 电压幅值和相角弧度 P_spec, Q_spec: 指定有功和无功 pq: PQ节点索引列表 pv: PV节点索引列表 n len(V) G Ybus.real B Ybus.imag dP_dtheta np.zeros((n, n)) dP_dV np.zeros((n, n)) dQ_dtheta np.zeros((n, n)) dQ_dV np.zeros((n, n)) for i in range(n): for j in range(n): if i j: # 对角元素 dP_dtheta[i, i] -Q_spec[i] - B[i, i] * V[i]**2 dP_dV[i, i] P_spec[i] / V[i] G[i, i] * V[i] dQ_dtheta[i, i] P_spec[i] - G[i, i] * V[i]**2 dQ_dV[i, i] Q_spec[i] / V[i] - B[i, i] * V[i] else: dP_dtheta[i, j] V[i] * V[j] * (G[i, j] * np.sin(theta[i]-theta[j]) - B[i, j] * np.cos(theta[i]-theta[j])) dP_dV[i, j] V[i] * (G[i, j] * np.cos(theta[i]-theta[j]) B[i, j] * np.sin(theta[i]-theta[j])) dQ_dtheta[i, j] -V[i] * V[j] * (G[i, j] * np.cos(theta[i]-theta[j]) B[i, j] * np.sin(theta[i]-theta[j])) dQ_dV[i, j] V[i] * (G[i, j] * np.sin(theta[i]-theta[j]) - B[i, j] * np.cos(theta[i]-theta[j])) # 按 PQ 和 PV 节点组装 pq_pv sorted(pq pv) J11 dP_dtheta[np.ix_(pq_pv, pq_pv)] J12 dP_dV[np.ix_(pq_pv, pq)] J21 dQ_dtheta[np.ix_(pq, pq_pv)] J22 dQ_dV[np.ix_(pq, pq)] return np.block([[J11, J12], [J21, J22]])参数说明pq和pv是节点类型列表平衡节点不参与。P_spec和Q_spec是净注入功率发电机出力减负荷。注意对角元素的公式里用到了当前功率所以每次迭代都要重新计算。如果雅可比矩阵条件数很大说明系统接近电压崩溃需要检查数据。2.3 牛顿-拉夫逊迭代的收敛判据迭代过程就是不断求解修正量 ( \Delta x -J^{-1} \Delta y )然后更新电压幅值和相角。收敛判据一般看最大功率不平衡量阈值取 ( 10^{-8} ) 到 ( 10^{-10} ) 之间。我见过有人用 ( 10^{-4} ) 就停结果和商业软件差在小数点后三位做灵敏度分析时完全不能用。def newton_raphson(Ybus, Sbus, V0, theta0, pq, pv, tol1e-10, max_iter20): V V0.copy() theta theta0.copy() for it in range(max_iter): # 计算注入功率 V_complex V * np.exp(1j * theta) S_calc V_complex * np.conj(Ybus V_complex) dP Sbus.real - S_calc.real dQ Sbus.imag - S_calc.imag # 不平衡量向量 mismatch np.concatenate([dP[pq pv], dQ[pq]]) if np.max(np.abs(mismatch)) tol: print(f收敛于第 {it} 次迭代) return V, theta, S_calc J build_jacobian(Ybus, V, theta, S_calc.real, S_calc.imag, pq, pv) dx np.linalg.solve(J, -mismatch) # 更新相角和幅值 n_pq_pv len(pq pv) theta[pq pv] dx[:n_pq_pv] V[pq] dx[n_pq_pv:] raise RuntimeError(潮流计算未收敛检查数据或增加迭代次数)逻辑说明每次迭代先算当前功率再算不平衡量。如果最大不平衡量小于容差就返回。否则构造雅可比矩阵解线性方程组得到修正量。注意更新时只改 PQ 和 PV 节点的相角只改 PQ 节点的幅值。PV 节点幅值保持不变平衡节点完全不更新。如果 20 次还不收敛基本可以断定是数据问题比如阻抗为 0 或者负荷过大。3. 用 SciPy 稀疏矩阵加速从 5 节点到 300 节点的跨越3.1 稀疏存储为什么能快一个数量级节点数超过 100 以后稠密矩阵的存储和求逆会变得非常吃力。一个 300 节点的系统雅可比矩阵是 600×600稠密存储约 2.8 MB求逆的浮点运算量在 ( 10^8 ) 量级。而实际电力网络是稀疏的每个节点平均只连 2 到 3 条支路。用scipy.sparse存储非零元素只有几千个求解速度能提升 10 倍以上。我一般用lil_matrix组装然后转成csr_matrix求解。下面是一个稀疏雅可比的构造示例。from scipy.sparse import lil_matrix, csr_matrix from scipy.sparse.linalg import spsolve def build_jacobian_sparse(Ybus, V, theta, P_spec, Q_spec, pq, pv): n len(V) pq_pv sorted(pq pv) n_pq_pv len(pq_pv) n_pq len(pq) dim n_pq_pv n_pq J lil_matrix((dim, dim), dtypefloat) G Ybus.real B Ybus.imag # 填充 H 块 for idx_i, i in enumerate(pq_pv): for idx_j, j in enumerate(pq_pv): if i j: J[idx_i, idx_j] -Q_spec[i] - B[i, i] * V[i]**2 else: J[idx_i, idx_j] V[i] * V[j] * (G[i, j] * np.sin(theta[i]-theta[j]) - B[i, j] * np.cos(theta[i]-theta[j])) # 填充 N 块 for idx_i, i in enumerate(pq_pv): for idx_j, j in enumerate(pq): if i j: J[idx_i, n_pq_pv idx_j] P_spec[i] / V[i] G[i, i] * V[i] else: J[idx_i, n_pq_pv idx_j] V[i] * (G[i, j] * np.cos(theta[i]-theta[j]) B[i, j] * np.sin(theta[i]-theta[j])) # 填充 M 块 for idx_i, i in enumerate(pq): for idx_j, j in enumerate(pq_pv): if i j: J[n_pq_pv idx_i, idx_j] P_spec[i] - G[i, i] * V[i]**2 else: J[n_pq_pv idx_i, idx_j] -V[i] * V[j] * (G[i, j] * np.cos(theta[i]-theta[j]) B[i, j] * np.sin(theta[i]-theta[j])) # 填充 L 块 for idx_i, i in enumerate(pq): for idx_j, j in enumerate(pq): if i j: J[n_pq_pv idx_i, n_pq_pv idx_j] Q_spec[i] / V[i] - B[i, i] * V[i] else: J[n_pq_pv idx_i, n_pq_pv idx_j] V[i] * (G[i, j] * np.sin(theta[i]-theta[j]) - B[i, j] * np.cos(theta[i]-theta[j])) return J.tocsr()参数说明lil_matrix适合逐元素填充转成csr_matrix后求解效率最高。注意索引映射前n_pq_pv行对应相角修正后n_pq行对应幅值修正。列索引同理。如果节点数少于 50稀疏矩阵的优势不明显直接用稠密矩阵更简单。3.2 用 spsolve 替代求逆np.linalg.solve内部会做 LU 分解但对于稀疏矩阵scipy.sparse.linalg.spsolve会利用稀疏结构速度快很多。下面是对比测试的代码片段。import time # 稠密求解 J_dense J.toarray() t0 time.time() dx_dense np.linalg.solve(J_dense, -mismatch) t1 time.time() # 稀疏求解 t2 time.time() dx_sparse spsolve(J, -mismatch) t3 time.time() print(f稠密求解耗时: {t1-t0:.4f} 秒) print(f稀疏求解耗时: {t3-t2:.4f} 秒) print(f结果差异: {np.max(np.abs(dx_dense - dx_sparse)):.2e})逻辑说明两种方法的结果应该几乎一致差异在 ( 10^{-12} ) 量级。如果差异很大说明稀疏矩阵组装有误。实际测试中300 节点系统稀疏求解能快 5 到 8 倍。注意spsolve返回的是稀疏解如果需要稠密数组用.toarray()转换。3.3 节点数超过 500 时的内存与速度权衡节点数到 500 以上雅可比矩阵维度超过 1000。这时候即使稀疏存储LU 分解的填充元也会增加。我一般会做三件事一是用scipy.sparse.linalg.splu做预分解二是检查节点编号顺序用reverse_cuthill_mckee减少带宽三是考虑用pypardiso这类并行求解器。下面是一个带宽优化的例子。from scipy.sparse.csgraph import reverse_cuthill_mckee def reorder_buses(Ybus): 按 RCM 算法重排节点减少稀疏矩阵带宽 graph (Ybus ! 0).astype(int) perm reverse_cuthill_mckee(graph, symmetric_modeTrue) return perm # 使用示例 perm reorder_buses(Ybus) Ybus_reordered Ybus[perm][:, perm] # 后续计算用重排后的矩阵最后再映射回原顺序参数说明reverse_cuthill_mckee返回的是排列索引。重排后导纳矩阵的稀疏结构更紧凑LU 分解的填充元更少。注意重排后节点类型列表也要相应调整。这个技巧在 1000 节点以上效果明显小系统可以忽略。4. 避坑与排查那些让潮流结果“差一点”的常见问题4.1 收敛了但结果不对检查平衡节点功率现象迭代 5 次就收敛但平衡节点的有功出力比预期大很多。原因平衡节点承担了所有网损和功率不平衡如果线路电阻设得太大网损会异常。解决检查支路电阻是否用了标幺值基准功率是否统一。我见过有人把欧姆值直接填进去结果网损差了 100 倍。4.2 雅可比矩阵奇异PV 节点转 PQ 节点现象迭代过程中出现LinAlgError: Singular matrix。原因某个 PV 节点的无功出力越限但代码没有把它转成 PQ 节点。解决每次迭代后检查 PV 节点的无功是否在限值内如果越限就转为 PQ并重新构造雅可比矩阵。这个逻辑在商业软件里是自动的手写代码容易漏。4.3 迭代震荡不收敛负荷模型与初值现象不平衡量在 ( 10^{-2} ) 附近来回跳20 次不收敛。原因初值太平所有电压设为 1.0 且相角为 0对于重载系统太乐观。解决用直流潮流给一个初始相角或者把电压初值设为 0.95 到 1.05 之间的随机值。另外检查负荷是否用了恒功率模型重载时改成恒阻抗模型能改善收敛性。4.4 标幺值系统不统一变压器变比翻车现象变压器两侧电压等级不同但潮流结果功率不守恒。原因变比没有折算到统一基准。解决所有阻抗和导纳都转成标幺值变压器的非标准变比要单独处理。我一般会在支路数据里加一列tap组装导纳矩阵时对变压器支路做修正。4.5 稀疏矩阵组装慢循环填充的代价现象300 节点系统组装雅可比矩阵花了 2 秒比求解还慢。原因用lil_matrix逐元素循环填充Python 循环开销大。解决用coo_matrix预先构造行、列、数据三个列表一次性生成。或者用 NumPy 向量化计算所有非对角元素再筛选非零项。5. 进阶技巧用自动微分和 GPU 把潮流计算推到更大规模5.1 用 JAX 自动微分替代手写雅可比手写雅可比矩阵容易出错尤其是极坐标下的偏导公式。JAX 可以自动求导把功率方程写成一个函数然后用jax.jacfwd直接得到雅可比矩阵。这样代码量减少一半而且不会抄错公式。import jax.numpy as jnp from jax import jacfwd def power_equations(V, theta, G, B): 计算所有节点的有功和无功注入 V_complex V * jnp.exp(1j * theta) S V_complex * jnp.conj((G 1j*B) V_complex) return jnp.concatenate([S.real, S.imag]) # 自动微分得到雅可比 jac_fn jacfwd(power_equations, argnums(0, 1))参数说明jacfwd返回的是前向模式雅可比对于潮流计算这种输入维度远小于输出维度的场景前向模式效率很高。注意 JAX 默认用 32 位浮点做潮流计算要开 64 位jax.config.update(jax_enable_x64, True)。这个方案在 1000 节点以内比手写雅可比快而且代码更短。5.2 批量工况扫描用 vmap 一次算 100 个场景做新能源接入分析时经常要扫描几百个出力场景。用 JAX 的vmap可以把牛顿迭代向量化一次算完所有场景。下面是一个简化示例。from jax import vmap def solve_power_flow(V0, theta0, P_spec, Q_spec): # 简化的牛顿迭代实际需要循环 # 这里只展示 vmap 的用法 return V0, theta0 # 批量扫描 batch_solve vmap(solve_power_flow, in_axes(None, None, 0, 0)) V_batch, theta_batch batch_solve(V0, theta0, P_batch, Q_batch)逻辑说明in_axes指定哪些参数需要批处理。P_batch和Q_batch的形状是(n_scenarios, n_buses)。vmap会自动把函数映射到每个场景上。注意牛顿迭代内部有while循环JAX 的vmap对动态控制流支持有限需要用jax.lax.while_loop重写。这个方案适合场景数多但节点数少的批量分析。5.3 结果验证和商业软件对标的三个检查点写完潮流程序怎么知道对不对我一般做三个检查。第一平衡节点功率是否等于所有发电机出力减去总负荷和网损。第二PV 节点的电压幅值是否严格等于设定值。第三把结果和商业软件对比电压幅值差异应在 ( 10^{-6} ) 以内相角差异在 ( 10^{-5} ) 弧度以内。如果差得多先检查标幺值基准再检查变压器变比最后检查雅可比矩阵的对角元素。我自己的习惯是每写一个新版本先用 5 节点系统跑通再上 14 节点最后才碰 300 节点。中间任何一步对不上绝不往下走。这个习惯帮我省了无数个调试的夜晚。希望帮到你。本文还有配套的精品资源点击获取