Python潮流计算实战：牛顿-拉夫逊法从源码到IEEE算例验证

简介：这份资源面向电力系统方向的工程师、研究人员及高校学生，聚焦如何用Python完成潮流计算这一电力网络分析的核心任务。内容围绕节点、支路、发电机与负荷建模，结合基尔霍夫定律构建非线性方程组，并借助NumPy、SciPy、Pandas等库实现求解与结果可视化，帮助读者理解电压分布与功率流动的评估方法。压缩包共27个文件，以png结果图、txt说明、py脚本为主，另含pyc、sln、pyproj等工程与项目配置文件，整体约687KB，结构紧凑便于按模块查阅。目前已有1332人学习下载。通过其中的示例脚本与案例数据，读者可对照代码走通建模、求解、迭代收敛到结果分析的完整流程，并进一步接触牛顿-拉弗森法、高斯-塞德尔迭代及无功补偿等进阶内容，为电力网络运行性能的优化提供可复用的实践参考。

1. 从一份 Python 潮流计算源码说起：它到底能算什么

如果你手头有一份电力系统的接线图，节点数从 3 个到几百个不等，想快速知道每个节点的电压幅值、相角，以及每条支路上流过的功率，那这份 Python 潮流计算资源就是干这个的。它不依赖商业仿真软件，用纯 Python 加 NumPy 就能把牛顿-拉夫逊法、PQ 分解法跑通，输入是节点参数和支路参数，输出是收敛后的电压向量和功率分布。适合电气工程专业的学生做课程设计，也适合刚转行进电力算法岗的工程师拿来当参考实现。我见过太多人卡在“公式看得懂、代码写不出”这一步，这份东西的价值就在于把雅可比矩阵的组装、稀疏线性求解、收敛判据这些环节全摊开了。下面我按“怎么装、怎么跑、怎么改、怎么排错”的顺序拆一遍，中间会穿插我实际调试时踩过的坑。

2. 环境准备与数据格式：把节点和支路参数喂进去

2.1 安装依赖与验证 NumPy 版本

这份源码的核心计算依赖 NumPy 的矩阵运算，尤其是numpy.linalg.solve和稀疏矩阵模块。如果你用的是 Python 3.8 及以上版本，直接 pip 安装即可。我一般会先建一个干净的虚拟环境，避免和系统里的其他包冲突。

python -m venv powerflow_env source powerflow_env/bin/activate # Windows 下用 powerflow_env\Scripts\activate pip install numpy scipy matplotlib

装完之后别急着跑主程序，先验证一下 NumPy 能不能正常做复数矩阵运算，因为潮流计算里电压和导纳都是复数。

import numpy as np # 构造一个简单的复数矩阵，验证 linalg.solve 是否支持 complex128 Y = np.array([[1+2j, 3-1j], [2+0j, 4+3j]], dtype=complex) V = np.array([1+0j, 1+0j], dtype=complex) I = Y @ V print(I) print(np.linalg.solve(Y, I))

这段代码先算Y @ V得到电流向量，再用solve反解电压，如果输出和初始的V一致，说明复数线性求解没问题。参数上注意dtype=complex必须显式指定，否则 NumPy 默认用 float64，遇到复数会直接报TypeError。我见过有人用python安装numpy库的方法搜到的老教程，装的是 1.19 之前的版本，linalg.solve对复数矩阵的支持有 bug，建议至少用 1.21 以上。

2.2 节点数据与支路数据的表格结构

潮流计算的数据分两张表：节点表（bus data）和支路表（branch data）。节点表每行描述一个节点，列包括节点编号、节点类型（平衡节点、PV 节点、PQ 节点）、电压幅值初值、相角初值、有功注入、无功注入。支路表每行描述一条支路，列包括首端节点、末端节点、电阻、电抗、对地电纳、变压器变比。

支路表里r和x是标幺值，b是对地总充电电纳，tap是变压器非标准变比，如果不是变压器支路就填 1.0。常见做法是把这两张表存成 CSV 或者直接写在 Python 脚本里作为列表。我一般会写一个load_case函数，把 CSV 读成 NumPy 数组，这样换算例的时候不用改代码。

import numpy as np def load_case(bus_file, branch_file): # 跳过表头，按逗号分隔读取 bus = np.loadtxt(bus_file, delimiter=',', skiprows=1) branch = np.loadtxt(branch_file, delimiter=',', skiprows=1) # bus 列顺序: id, type, Vm, Va, Pg, Qg, Pl, Ql # branch 列顺序: from, to, r, x, b, tap return bus, branch bus_data, branch_data = load_case('bus.csv', 'branch.csv') print('节点数:', bus_data.shape[0], '支路数:', branch_data.shape[0])

这里skiprows=1是跳过表头，如果你的 CSV 没有表头就去掉这个参数。delimiter=','对应逗号分隔，用 Excel 导出的 CSV 默认就是逗号。注意节点编号必须从 1 开始连续，因为后面组装雅可比矩阵时用编号直接做索引，中间缺号会导致矩阵维度对不上。这个坑我在第一次跑 IEEE 14 节点算例时就遇到了，当时手改数据删了一行，结果程序报IndexError，查了半天才发现是编号不连续。

2.3 导纳矩阵的组装逻辑

导纳矩阵 Y 是潮流计算的核心数据结构，它是一个 n×n 的复数方阵，n 是节点数。对角元素 Y[i][i] 是节点 i 的自导纳，等于所有连在 i 上的支路导纳之和加上对地电纳；非对角元素 Y[i][j] 是节点 i 和 j 之间的互导纳，等于它们之间支路导纳的负值。如果支路有变压器变比，还要在导纳上除以变比或者乘以变比的共轭。

def build_ybus(bus_data, branch_data): n = bus_data.shape[0] Y = np.zeros((n, n), dtype=complex) for br in branch_data: f = int(br[0]) - 1 # 转成 0 基索引 t = int(br[1]) - 1 r, x, b, tap = br[2], br[3], br[4], br[5] z = r + 1j * x y = 1.0 / z # 变压器支路处理：非标准变比在首端 if tap != 1.0: y = y / tap Y[f][f] += y / tap # 首端自导纳要除以 tap 的平方 else: Y[f][f] += y Y[t][t] += y + 1j * b / 2 Y[f][t] -= y Y[t][f] -= y return Y

这段代码里tap的处理是最容易翻车的地方。标准变比在首端时，首端自导纳要除以tap的平方，互导纳除以tap，末端自导纳不变。我见过有人把tap直接乘上去，结果算出来的电压偏差特别大，还以为是算法不收敛。另外b是对地总电纳，要分一半到首端、一半到末端，如果支路表里给的是单侧电纳，这里就要改成1j * b而不是1j * b / 2。参数含义一定要和填表的人对齐，否则后面所有结果都是错的。

3. 牛顿-拉夫逊法核心：雅可比矩阵与迭代收敛

3.1 功率不平衡量的计算

牛顿-拉夫逊法的思路是：先给所有节点一个电压初值（通常是平启动，即 PQ 节点电压设为 1.0∠0°，PV 节点电压幅值固定、相角初值 0），然后计算每个节点的注入功率，和给定值比较得到不平衡量 ΔP 和 ΔQ。对于 PQ 节点，ΔP 和 ΔQ 都要算；对于 PV 节点，只算 ΔP，因为电压幅值已经给定；平衡节点不参与方程。

def calc_power(Y, V, bus_data): n = len(V) S = V * np.conj(Y @ V) # 注入功率 S = V * (YV)* P = S.real Q = S.imag dP = np.zeros(n) dQ = np.zeros(n) for i in range(n): typ = int(bus_data[i][1]) if typ == 1: # 平衡节点 continue Pg, Qg = bus_data[i][4], bus_data[i][5] Pl, Ql = bus_data[i][6], bus_data[i][7] dP[i] = Pg - Pl - P[i] if typ == 3: # PQ 节点才算无功不平衡 dQ[i] = Qg - Ql - Q[i] return dP, dQ

这里S = V * np.conj(Y @ V)是标准的注入功率公式，np.conj取共轭。注意dP和dQ对平衡节点和 PV 节点的处理不同：平衡节点直接跳过，PV 节点的dQ保持 0。如果这里写错，把 PV 节点的无功不平衡也算进去，雅可比矩阵的维度就会多出来，导致求解失败。我一般会在循环里加一个断言，检查dP和dQ中非零元素的个数是否等于方程数。

3.2 雅可比矩阵的四个子块

雅可比矩阵 J 是一个 (n_pq + n_pv) × (n_pq + n_pv) 的实数矩阵，分四个子块：H = ∂P/∂θ，N = ∂P/∂V，M = ∂Q/∂θ，L = ∂Q/∂V。对于 i ≠ j 的情况，H[i][j] = -V[i]V[j](G[i][j]*sin(θ[i]-θ[j]) - B[i][j]cos(θ[i]-θ[j]))，N[i][j] = -V[i](G[i][j]*cos(θ[i]-θ[j]) + B[i][j]*sin(θ[i]-θ[j]))，M 和 L 类似。对角元素要额外加上注入功率项。

def build_jacobian(Y, V, bus_data): n = len(V) G = Y.real B = Y.imag theta = np.angle(V) Vm = np.abs(V) # 确定方程顺序：先所有 PQ 和 PV 节点的 dP，再所有 PQ 节点的 dQ pq = [i for i in range(n) if int(bus_data[i][1]) == 3] pv = [i for i in range(n) if int(bus_data[i][1]) == 2] p_buses = pq + pv # 有功方程对应的节点 q_buses = pq # 无功方程对应的节点 np_ = len(p_buses) nq_ = len(q_buses) J = np.zeros((np_ + nq_, np_ + nq_)) # 填充 H, N, M, L 的代码略，核心是按上述偏导公式逐元素计算 return J, p_buses, q_buses

雅可比矩阵的组装是整份代码里最长的部分，但逻辑很直接：遍历所有节点对，按公式算偏导。我建议先把p_buses和q_buses两个列表确定好，后面填矩阵的时候直接用它们的索引，不要用原始节点编号，否则维度会乱。常见做法是用np.ix_做花式索引，但新手容易搞混，老老实实写双重循环反而更稳。参数上注意G和B是从Y里直接取的实部和虚部，不要重新算。

3.3 迭代求解与收敛判据

每一轮迭代的流程是：算不平衡量 → 如果最大不平衡量小于阈值就退出 → 组装雅可比矩阵 → 解线性方程组 J * Δx = Δy → 更新电压幅值和相角。收敛判据一般取 1e-6 到 1e-8，单位是标幺值。我一般用 1e-8，因为 IEEE 标准算例的参考结果通常给到小数点后 6 位。

def newton_raphson(Y, bus_data, tol=1e-8, max_iter=20): n = len(bus_data) V = np.ones(n, dtype=complex) for i in range(n): V[i] = bus_data[i][2] * np.exp(1j * np.radians(bus_data[i][3])) for it in range(max_iter): dP, dQ = calc_power(Y, V, bus_data) mismatch = np.concatenate([dP[[i for i in range(n) if int(bus_data[i][1]) != 1]], dQ[[i for i in range(n) if int(bus_data[i][1]) == 3]]]) if np.max(np.abs(mismatch)) < tol: print(f'第 {it} 次迭代收敛') return V J, p_buses, q_buses = build_jacobian(Y, V, bus_data) dx = np.linalg.solve(J, mismatch) # 更新相角和幅值，代码略 raise RuntimeError('达到最大迭代次数仍未收敛')

这里max_iter=20是经验值，牛顿法一般 3 到 5 次就收敛，超过 10 次还没收敛基本就是数据有问题。mismatch的拼接顺序必须和雅可比矩阵的行顺序一致，否则solve出来的dx对应不上。我见过有人把dP和dQ的顺序搞反，结果电压越迭代越离谱，最后发散。如果遇到不收敛，先把tol放宽到 1e-4 看看能不能出结果，能出就说明是数据精度问题，不能出就是雅可比矩阵组装错了。

4. 避坑与排查：那些让程序跑不通的细节

4.1 现象：迭代次数正常但电压结果明显偏离参考值

原因通常是支路参数的单位不统一。节点注入功率用的是标幺值，但支路表里的r和x如果填的是欧姆值，导纳就会差好几个数量级。解决方法是检查基准功率和基准电压是否一致，一般基准功率取 100 MVA，基准电压取各电压等级的额定值，支路阻抗要按Z_base = V_base^2 / S_base换算成标幺值。我一般会在load_case之后加一行打印，看看Y矩阵对角元素的模值是不是在 1 到 100 之间，如果出现 0.001 或者 10000 这种极端值，基本就是单位错了。

4.2 现象：程序报LinAlgError: Singular matrix

雅可比矩阵奇异，说明方程数不够或者有冗余。最常见的原因是节点类型设置错了，比如把平衡节点设成了 PQ 节点，导致没有参考相角，整个系统缺少一个约束。另一个原因是孤岛节点，某个节点没有任何支路连接，它的导纳矩阵对角元素为 0，雅可比矩阵对应行全零。解决方法是检查节点表里type列是否恰好有一个 1（平衡节点），以及每个节点是否至少出现在一条支路里。我一般会写一个check_topology函数，用邻接矩阵的连通性判断有没有孤岛。

4.3 现象：PV 节点的无功出力超出限值但程序不处理

标准牛顿-拉夫逊法不处理 PV 节点的无功限值，如果算出来的无功超过发电机能力，电压幅值应该从 PV 转成 PQ。这份源码默认不做这个转换，所以如果你算的算例里 PV 节点无功越限，结果会偏乐观。解决方法是手动加一个外层循环：每次收敛后检查 PV 节点的无功，如果越限就把该节点改成 PQ，无功固定为限值，重新跑一遍。这个逻辑在python电力系统潮流计算的进阶用法里很常见，但基础版本一般不带。

4.4 现象：用python构建邻接矩阵时节点编号不连续导致索引错位

很多人习惯用节点编号直接做矩阵下标，但实际电网的节点编号可能不连续，比如 1, 2, 5, 7。如果直接拿编号当索引，Y[5][5]就会越界。解决方法是先建一个映射字典，把原始编号映射到 0 到 n-1 的连续索引。我一般会在load_case里就做这个映射，后面所有代码都用映射后的索引，输出结果时再映射回原始编号。这个坑在python构建邻接矩阵的教程里经常被忽略，但实际工程数据里编号不连续是常态。

4.5 现象：收敛后功率不平衡量在 1e-6 附近震荡不下降

这是典型的数值精度问题，通常出现在节点数超过 100 的大算例里。原因是雅可比矩阵的条件数太大，np.linalg.solve的默认精度不够。解决方法是改用scipy.sparse.linalg.spsolve做稀疏求解，或者把tol放宽到 1e-6。如果还不行，检查一下有没有特别小的支路阻抗（比如小于 1e-5），这种支路会导致导纳矩阵元素极大，条件数爆炸。常见做法是给支路阻抗加一个下限，比如 1e-4，或者用标幺值重新换算。

5. 进阶用法：从固定算例到批量计算与结果验证

5.1 批量读取算例并自动生成结果表

实际工作中不会只算一个算例，通常要跑几十个不同负荷水平的场景。我一般会写一个批处理脚本，遍历一个文件夹下的所有 CSV，每个算例跑一遍潮流，把收敛后的电压幅值、相角、支路功率写到一个汇总表里。下面是一个简化版的批处理框架。

import os import pandas as pd def batch_run(case_dir, output_csv): results = [] for fname in os.listdir(case_dir): if not fname.endswith('_bus.csv'): continue prefix = fname.replace('_bus.csv', '') bus_file = os.path.join(case_dir, fname) branch_file = os.path.join(case_dir, prefix + '_branch.csv') bus_data, branch_data = load_case(bus_file, branch_file) Y = build_ybus(bus_data, branch_data) try: V = newton_raphson(Y, bus_data) for i, v in enumerate(V): results.append({ 'case': prefix, 'bus_id': int(bus_data[i][0]), 'Vm': abs(v), 'Va': np.degrees(np.angle(v)) }) except RuntimeError as e: print(f'{prefix} 不收敛: {e}') pd.DataFrame(results).to_csv(output_csv, index=False)

这里用pandas做结果汇总，case列区分不同算例，bus_id保留原始编号。注意newton_raphson里如果抛异常，这里要捕获并跳过，不能让一个坏算例中断整个批处理。参数上output_csv建议用绝对路径，避免相对路径在不同工作目录下找不到。我一般还会在结果表里加一列converged标记是否收敛，方便后续筛选。

5.2 用 IEEE 标准算例验证实现正确性

写完代码后必须验证，最直接的方法是用 IEEE 14 节点、30 节点、57 节点标准算例跑一遍，和公开的参考结果对比。参考结果里通常给出每个节点的电压幅值和相角，以及平衡节点的注入功率。我一般会算电压幅值的最大绝对误差，如果小于 1e-4 就认为实现正确。

如果误差在 1e-3 以上，优先检查雅可比矩阵的对角元素，尤其是N和L子块里有没有漏掉注入功率项。另一个常见错误是平衡节点的电压相角没有固定为 0，导致整个系统的相角参考漂移。我一般会在newton_raphson开头强制把平衡节点的Va设为 0，并在每次更新后重新置零。

5.3 一个具体技巧：用numpy的广播机制加速雅可比矩阵组装

基础版的雅可比矩阵组装用双重循环，节点数一多就慢。我后来改成用 NumPy 的广播机制一次性算完所有非对角元素。核心思路是把theta和Vm扩展成二维矩阵，用theta[:, None] - theta[None, :]得到所有节点对的相角差，然后直接套公式。

def build_jacobian_fast(Y, V, bus_data): n = len(V) G, B = Y.real, Y.imag theta = np.angle(V) Vm = np.abs(V) dtheta = theta[:, None] - theta[None, :] # H 非对角部分 H = -Vm[:, None] * Vm[None, :] * (G * np.sin(dtheta) - B * np.cos(dtheta)) # N 非对角部分 N = -Vm[:, None] * (G * np.cos(dtheta) + B * np.sin(dtheta)) # 对角元素修正和子块提取代码略 return H, N

广播写法比循环快 10 倍以上，但内存占用也大，节点数超过 1000 的时候要注意。我一般会在节点数小于 500 的时候用广播，超过 500 就换稀疏矩阵。这个技巧在python量化交易策略代码里也常见，本质都是用向量化替代循环。从那以后我每次写矩阵组装都先想能不能广播，实在不行才写循环。希望帮到你。

本文还有配套的精品资源，点击获取