电力系统潮流计算通用程序：四表建模、三种解法与结果校验

简介：面向电力系统分析与潮流计算学习者的C++工具包，聚焦稳态网络电压、电流与功率分布的求解，适合课程设计、科研复现及工程入门者参考。压缩包共46个文件、约2.23MB，以18个txt算例与结果数据、3个xlsx表格、2个cpp源文件和2个h头文件为主，另含Visual Studio工程及编译中间文件，可直接编译运行并对照输入输出。程序采用牛顿法迭代求解非线性潮流方程，并借助tinny2稀疏矩阵技术降低大规模节点网络的存储与计算开销，提供IEEE14至300等标准系统算例，便于验证收敛性与效率。随附报告涵盖设计思路、算法描述与结果分析，已有840人学习。读者可据此理解潮流计算从方程线性化、雅可比矩阵组装到稀疏求解的完整实现路径，也能借助多规模算例开展性能比较与二次开发。

1. 一个“收敛了但结果不对”的 30 节点算例

去年帮一位做配网规划的朋友看一个 30 节点算例，程序收敛得干干净净，可现场拿到的母线电压与程序结果差了 0.08 p.u.，查了一整天才发现：变压器变比填在了 to 端，代码却按 from 端读。潮流计算（powerflow）最麻烦的情况往往是“算出来了但你不知道能不能信”，并非“算不出来”。这套电力系统潮流计算通用程序针对的就是这一环节：节点、支路、变压器、发电机拆成四张标准表，内置牛顿-拉夫逊、快速解耦、高斯-赛德尔三种解法，支持 PV/PQ/平衡节点混排、无功越限自动切换、批量算例扫描，输出母线电压、支路潮流、网损三张报表。想快速验证网架方案的可以直接用，想把求解器拆开换算法、接自己数据格式的也能改。

2. 架构与数据模型：四张表撑起整个潮流计算

这套程序最核心的设计取舍，是把“数据”和“算法”彻底切开。求解器不认得 bus 表长什么样，它只吃三样东西：节点导纳矩阵Ybus、注入功率向量Sbus、以及平衡节点/PV 节点/PQ 节点的索引数组。这条边界划清楚了，换算法就是换一个函数，换数据格式就是换一个 reader，两边互不干扰；边界没划清楚，改一个 CSV 列名就能让整个牛顿迭代崩掉，这类返工我见过不止一次。

2.1 目录结构：io、model、solver、report 四层

拿到源码包先别急着python run_pf.py，把目录看一遍比跑十次都值。这套程序的目录基本长这样：

powerflow/ ├── run_pf.py # 命令行入口：参数解析 + 调度 ├── config/ │ └── default.yaml # 默认容差、迭代上限、输出开关 ├── pf/ │ ├── model.py # Bus/Branch/Gen/Load 四张表 + 字段校验 │ ├── ybus.py # 节点导纳矩阵组装（含变压器、对地并联） │ ├── caseio.py # JSON/CSV 读写，标幺化统一在这一层做 │ ├── qlimit.py # 无功越限检测与 PV→PQ 切换 │ ├── solver/ │ │ ├── registry.py # 求解器注册表，二次开发的挂载点 │ │ ├── newton.py # 牛顿-拉夫逊 │ │ ├── fdpf.py # 快速解耦 │ │ └── gauss.py # 高斯-赛德尔 │ └── report.py # 母线/支路/网损三张报表 ├── cases/ │ ├── case14.json │ └── case30.json ├── ref/ │ └── case14_voltage_ref.csv └── tests/ └── test_newton_case14.py

四层的职责很明确。caseio负责把外部文件（JSON、CSV、甚至别的工具导出的文本）读成内存里的四张表，并且一次性完成标幺化——电压除以 baseKV、功率除以 baseMVA、阻抗除以 Zbase。很多新手版本把标幺化散落在各个模块里，最后一定会出现“线路阻抗是标幺值、负荷功率是有名值”的混合状态，算出来的电压能到 0.3 p.u.。

model只做校验不做计算：检查有没有孤立的节点编号、支路两端是否都在 bus 表里、PV 节点的 Vm 是否越界、平衡节点是不是只有一个。solver是纯算法层。report把求解器吐出来的复数电压向量还原成人能看的 kV、MW、MVar。

提示：先看config/default.yaml，里面决定了容差和迭代上限；默认容差如果被改成1e-4，牛顿法通常 3 次就“收敛”，但结果只能信到小数点后两位。

2.2 四张核心表：bus、branch、gen、load 的字段设计

数据模型是一套潮流程序的骨架。这套程序用的是电力系统里最通用的“四表结构”，字段命名也和主流工具对齐，方便你把别处的算例直接导进来。

gen表和load表是 bus 表的补充：一个节点可以挂多台机组、多个负荷，caseio在加载阶段把它们按bus_id聚合回Pg/Qg/Pd/Qd。这种“明细表 + 聚合”的写法比直接在 bus 表里堆字段灵活，尤其适合后续要接时序数据的场景。

常见做法是在model.py里定义一个轻量的Case数据类，把四张表装进去，顺便缓存一次索引：

from dataclasses import dataclass import numpy as np @dataclass class Case: bus: np.ndarray # (nb, ncol) 浮点/整型混合，用列号访问 branch: np.ndarray # (nl, ncol) gen: np.ndarray load: np.ndarray baseMVA: float = 100.0 def __post_init__(self): # 节点类型索引：求解器只认这三个数组，不认 bus 表 types = self.bus[:, 1].astype(int) self.ref = np.where(types == 3)[0] # 平衡节点，必须恰好 1 个 self.pv = np.where(types == 2)[0] self.pq = np.where(types == 1)[0] if len(self.ref) != 1: raise ValueError(f"平衡节点数量应为 1，实际 {len(self.ref)}")

__post_init__里做索引缓存是这套代码里最划算的一处优化：牛顿迭代每轮都要用pvpq去切雅可比矩阵，如果每轮都重新where一遍，30 节点看不出来，118 节点、上千次批量扫描时的开销就很明显了。ref必须恰好一个这条校验也放在这里，因为平衡节点缺失是“雅可比矩阵奇异”最常见的诱因，早报错比让np.linalg.solve抛LinAlgError友好得多。

2.3 Ybus 组装：稀疏矩阵、变压器变比与对地并联

节点导纳矩阵是整个程序的地基，写错了后面全是白算。标准支路（含变压器）的四项导纳贡献是这样的：

import numpy as np from scipy.sparse import lil_matrix def build_ybus(case, include_shunt=True): nb = case.bus.shape[0] Y = lil_matrix((nb, nb), dtype=complex) for k in range(case.branch.shape[0]): br = case.branch[k] if int(br[10]) == 0: # status=0，停运支路直接跳过 continue f, t = int(br[0]) - 1, int(br[1]) - 1 r, x, b = br[2], br[3], br[4] tap, shift = br[8], br[9] if tap == 0: tap = 1.0 ys = 1.0 / complex(r, x) # 串联导纳 tap_c = tap * np.exp(1j * np.deg2rad(shift)) # 复数变比 # 变压器按 f 端建模：f 侧除以 |tap|^2，t 侧不变 Y[f, f] += (ys + 1j * b / 2) / (tap_c * np.conj(tap_c)) Y[t, t] += (ys + 1j * b / 2) Y[f, t] += -ys / np.conj(tap_c) Y[t, f] += -ys / tap_c if include_shunt: for i in range(nb): Y[i, i] += complex(case.bus[i, 9], case.bus[i, 10]) # Gs + jBs return Y.tocsr()

这段代码有三个地方值得反复看。第一是对地电纳要除以 2：b在 branch 表里存的是全长充电电纳，等分到两端，所以(ys + 1j*b/2)；如果算例给的是半长值，这里就会多出一倍，表现为轻载时电压偏高、网损偏小，属于典型的“结果看起来挺合理”的错误。

第二是变比的共轭位置。Y[f,t] = -ys / conj(tap_c)、Y[t,f] = -ys / tap_c，这两个共轭一旦写反，变压器两侧的有功方向就整体反了。我见过有人在调试时把这两行对调“试试看”，结果潮流方向对了、网损不对，因为对调只是把变比挪到了 t 侧，而表里存的是 f 侧。

第三是稀疏结构的选择。lil_matrix适合逐元素累加，组装完立刻tocsr()转成行压缩格式，后面Ybus @ V的矩阵向量乘法就快了。用稠密numpy数组在 30 节点没问题，上万节点的配网会直接把内存吃满。

注意：status字段的位置（第 10 列）在不同算例文件里可能不一致，导入前一定对着表头核一遍。索引错一列，程序不会报错，只会安静地算错。

3. 核心算法拆解：一次牛顿迭代到底发生了什么

牛顿-拉夫逊法之所以是潮流计算的主力，是因为它在正常运行方式下二次收敛——误差从 1e-2 掉到 1e-4，再到 1e-8，通常 4 到 5 次迭代就够。但“二次收敛”是有前提的：初值不能太离谱，雅可比矩阵不能奇异，不平衡量的符号不能搞反。这一章就沿着newton_pf这条链路，把一次迭代拆成三段看。

3.1 不平衡量：从 Sbus 到 F 向量

迭代开始前要先拼出指定注入功率Sbus。规则很朴素：PQ 节点取-(Pd + jQd)，PV 节点取(Pg - Pd) + j(-Qd)，平衡节点的有功无功都当成未知量、不参与方程。

def make_sbus(case): nb = case.bus.shape[0] Sbus = np.zeros(nb, dtype=complex) for i in range(nb): Pd, Qd = case.bus[i, 4], case.bus[i, 5] Pg, Qg = case.bus[i, 6], case.bus[i, 7] if i in case.pq: Sbus[i] = complex(-Pd, -Qd) # PQ：注入 = -负荷 elif i in case.pv: Sbus[i] = complex(Pg - Pd, -Qd) # PV：Q 是未知量，置 0 else: Sbus[i] = complex(Pg - Pd, Qg - Qd) # 平衡节点，实际不参与方程 return Sbus

PV 节点的无功为什么要置 0？因为它的无功出力是待求量，方程里只保留有功平衡，无功方程直接不列。等迭代收敛后，再用Q = imag(V * conj(Ybus @ V))反算出这台机组的实际无功，去跟Qmax/Qmin比。这一步在qlimit.py里做，越限就把该节点从 PV 改成 PQ，无功固定在上限或下限，重新迭代。

有了Sbus，不平衡量就是一行：

S = V * np.conj(Ybus @ V) mis = np.r_[(Sbus.real - S.real)[pvpq], (Sbus.imag - S.imag)[pq]]

注意这里是指定值减计算值。符号反过来，迭代的方向就反了，表现为电压越算越离谱、两三次之后直接发散。判断依据很简单：如果第一次迭代的不平衡量是正数、第二次变成更大的正数，八成是符号搞反了。

pvpq是np.r_[case.pv, case.pq]，顺序无所谓，但整个程序里必须保持一致——雅可比矩阵的行列、状态量的回写位置都依赖这个顺序。中途有人把pvpq改成np.sort(...)，结果收敛了但电压对不上，就是顺序不一致导致的错位。

3.2 雅可比矩阵的四个分块

雅可比矩阵是[dP/dVa, dP/dVm; dQ/dVa, dQ/dVm]的四个分块，规模是(npvpq+npq) × (npvpq+npq)。手推公式容易出错，实践中更稳的写法是先算复数域的dS/dVa和dS/dVm，再取实部虚部。

def dSbus_dV(Ybus, V): """返回 dS/dVa 与 dS/dVm，均为 nb×nb 复矩阵""" Y = Ybus.toarray() if hasattr(Ybus, "toarray") else np.asarray(Ybus) Ibus = Y @ V diagV, diagI = np.diag(V), np.diag(Ibus) diagVn = np.diag(V / np.abs(V)) # 单位相量 dS_dVm = diagV @ np.conj(Y @ diagVn) + np.conj(diagI) @ diagVn dS_dVa = 1j * diagV @ np.conj(diagI - Y @ diagV) return dS_dVa, dS_dVm

这两行公式看着唬人，来源其实很清楚：S = diag(V) · conj(Y·V)，对Va求导时V的幅值不变、只转角度，所以多了个1j；对Vm求导时，V = Vm · e^{jθ}，幅值变化等价于乘单位相量。diagI - Y @ diagV这一项是链式法则带出来的，漏掉它，雅可比矩阵就退化成一个近似版本，迭代次数从 4 次变成十几次，节点一多直接不收敛。

分块组装：

dS_dVa, dS_dVm = dSbus_dV(Ybus, V) J11 = dS_dVa[np.ix_(pvpq, pvpq)].real J12 = dS_dVm[np.ix_(pvpq, pq)].real J21 = dS_dVa[np.ix_(pq, pvpq)].imag J22 = dS_dVm[np.ix_(pq, pq)].imag J = np.block([[J11, J12], [J21, J22]])

np.ix_是这里的常用技巧：行列同时切片，拿到的是子矩阵；如果写成dS_dVa[pvpq][:, pvpq]，两次索引会先复制再切，节点上千时性能差一个量级，而且中间结果占内存。J12的列是pq（只有 PQ 节点有电压幅值方程），J21的行是pq（只有 PQ 节点有无功方程），这两个索引最容易被统一写成pvpq，一写就错。

提示：调试时可以检查J的对角元素。正常情况下J11的对角元是负数（有功对相角），J22的对角元也是负数（无功对幅值），如果出现大面积正值，基本可以断定Ybus的符号写反了。

3.3 修正方程求解与状态量回写

把上面几段拼起来就是完整的迭代循环：

def newton_pf(Ybus, Sbus, V0, case, tol=1e-8, max_iter=20): V = V0.astype(complex).copy() pvpq = np.r_[case.pv, case.pq] npvpq, npq = len(pvpq), len(case.pq) history = [] for it in range(max_iter): S = V * np.conj(Ybus @ V) mis = np.r_[(Sbus.real - S.real)[pvpq], (Sbus.imag - S.imag)[case.pq]] err = np.max(np.abs(mis)) history.append(err) if err < tol: return V, it, history # 收敛，返回迭代次数 dS_dVa, dS_dVm = dSbus_dV(Ybus, V) J = np.block([ [dS_dVa[np.ix_(pvpq, pvpq)].real, dS_dVm[np.ix_(pvpq, case.pq)].real], [dS_dVa[np.ix_(case.pq, pvpq)].imag, dS_dVm[np.ix_(case.pq, case.pq)].imag], ]) dx = np.linalg.solve(J, mis) # F = Sspec - Scalc，F' = -J，故 dx = J⁻¹·mis Va, Vm = np.angle(V), np.abs(V) Va[pvpq] += dx[:npvpq] Vm[case.pq] += dx[npvpq:] V = Vm * np.exp(1j * Va) raise RuntimeError(f"牛顿法 {max_iter} 次未收敛，末次不平衡量 {history[-1]:.3e}")

dx = solve(J, mis)这个符号关系值得单独说。设F(x) = Sspec - Scalc(x)，我们要F=0，牛顿步是x_{k+1} = x_k - F/F'；而F' = -J（因为J是Scalc的导数），所以dx = J⁻¹ · F，直接加而不是减。很多实现里写成dx = -solve(J, mis)，然后不平衡量又定义成Scalc - Sspec，两次负号抵消，看起来也能跑，但只要有人改动其中一个，整套就崩。

回写只更新Va[pvpq]和Vm[pq]：平衡节点的幅值相角都固定，PV 节点的幅值固定、只更新相角。写成Vm[:] += ...是最常见的越界错误，PV 节点电压被推着走，表现为这台机组的无功越算越大、最后报越限。

max_iter默认 20 已经足够宽松，正常算例 5 次内收敛。如果 20 次还不收敛，别急着调大——先看history的衰减速度：误差不降反升是发散，误差在 1e-3 附近反复横跳是震荡，两种情况的原因完全不同，前者多半是初值或数据问题，后者常见于重载或接近电压崩溃点的算例。

3.4 快速解耦与高斯-赛德尔：什么场景该降级

牛顿法精度高，但每轮都要重新组装并分解雅可比矩阵，118 节点以下的算例感觉不到，上千节点、批量扫描几百个方式时就慢了。这套程序提供了两种降级方案。

快速解耦（PQ 解耦法）利用高压网里R << X的特点，把雅可比矩阵近似成两个常数矩阵，迭代中不再更新：

def fast_decoupled_pf(Ybus, Sbus, V0, case, tol=1e-8, max_iter=30): Y = Ybus.toarray() if hasattr(Ybus, "toarray") else Ybus V = V0.astype(complex).copy() pvpq, pq = np.r_[case.pv, case.pq], case.pq Bp = -Y.imag[np.ix_(pvpq, pvpq)] # B'：有功-相角 Bpp = -Y.imag[np.ix_(pq, pq)] # B''：无功-幅值 for it in range(max_iter): S = V * np.conj(Y @ V) dP = (Sbus.real - S.real)[pvpq] dQ = (Sbus.imag - S.imag)[pq] if max(np.max(np.abs(dP)), np.max(np.abs(dQ))) < tol: return V, it Va, Vm = np.angle(V), np.abs(V) Va[pvpq] += np.linalg.solve(Bp, dP / Vm[pvpq]) Vm[pq] += np.linalg.solve(Bpp, dQ / Vm[pq]) V = Vm * np.exp(1j * Va) raise RuntimeError("快速解耦未收敛")

Bp、Bpp是常数矩阵，可以在循环外做一次 LU 分解，之后每轮只做前代回代，这就是它在批量扫描里快的原因。代价是收敛变成线性，迭代次数通常 8 到 15 次。使用时有一条边界条件：线路的 R/X 比值不能太大。城市配网电缆的 R/X 能到 1 以上，解耦前提不成立，这时快速解耦会震荡甚至发散，应换回牛顿法。这也是为什么程序默认--method nr而不默认 fdpf。

高斯-赛德尔法只作为兜底和教学用，收敛慢但每轮代价极低、内存占用小，适合几千节点的稀疏配网先跑个粗略初值，再把结果作为牛顿法的热启动。这套代码里它就是一个简单循环：

def gauss_seidel_pf(Ybus, Sbus, V0, case, tol=1e-8, max_iter=500): Y = Ybus.toarray() V = V0.astype(complex).copy() for it in range(max_iter): for i in range(len(V)): if i in case.ref or i in case.pv: continue # 把第 i 行的自导纳项单独解出来 Ii = Sbus[i] / np.conj(V[i]) V[i] = (Ii - Y[i, :] @ V + Y[i, i] * V[i]) / Y[i, i] if np.max(np.abs(V - V0)) < tol: return V, it V0 = V.copy() raise RuntimeError("高斯-赛德尔未收敛")

三种解法的选择经验是：小于 300 节点、要精确结果，用nr；高压网、上千节点、要批量扫方式，用fdpf；数据质量差、初值离谱、只是想先看看电压大概形态，用gs热启动再切nr。程序里--method就是这三个字符串。

4. 跑起来：环境、命令行参数与结果校验

算法讲清楚之后，得先让它跑出第一组数，否则后面全是纸上谈兵。这一章按“装环境 → 跑单个算例 → 批量扫 → 对参考解”的顺序来，每步都能独立验证。

4.1 依赖安装与最小可跑路径

依赖很轻，numpy+scipy是核心，pandas只用于报表导出，pyyaml读配置，matplotlib可选（画收敛曲线）。建议开虚拟环境，别装到系统 Python 里。

# 1. 建环境（Python 3.9 以上都行，3.10/3.11 最稳） python -m venv .venv source .venv/bin/activate # Windows: .venv\Scripts\activate # 2. 装依赖 pip install numpy scipy pandas pyyaml # 或者 pip install -r requirements.txt # 3. 自检：跑测试用例，先确认算法层是好的 python -m pytest tests/ -q # 4. 跑第一个算例 python run_pf.py --case cases/case14.json --method nr --verbose

第 3 步容易被跳过，但它是排错时的分水岭：测试通过说明算法和 Ybus 组装没问题，后面跑不通就是数据或参数的问题；测试不通过就先别碰算例文件。tests/test_newton_case14.py里通常会硬编码两三个节点的参考电压，跑一下就知道牛顿法有没有被改坏。

常见的新手卡点是ModuleNotFoundError: No module named 'pf'，原因是没在项目根目录执行，或者run_pf.py里的sys.path没把当前目录加进去。简单办法是在根目录执行，或者pip install -e .装成可编辑包。

4.2 三种入口：单算例、批量扫描、配置驱动

run_pf.py的参数不多，但每个都影响结果，值得对着表过一遍：

单算例是最常用的：

python run_pf.py --case cases/case30.json --method nr \ --tol 1e-8 --max-iter 20 \ --out results/case30 --verbose

加--verbose后终端会打出类似iter 1: 2.31e-01 / iter 2: 8.42e-03 / iter 3: 1.9e-05 / iter 4: 6.1e-10的序列，这条曲线比任何日志都好用：正常应该是阶梯式下降，每次大约掉两个数量级；如果第 2 次就掉到 1e-9，八成是算例太简单或者容差被改过；如果一直卡在 1e-3，去看第 5 章的越限问题。

负荷扫描是规划场景里最实用的功能，一行命令能出一族结果：

python run_pf.py --case cases/case14.json --sweep-load 0.8:1.2:0.1 \ --method fdpf --out results/sweep14

它会把负荷倍率从 0.8 到 1.2 按 0.1 步长各算一次，每次结果单独存目录，最后在summary.json里汇总各方式的收敛情况和最低电压。做重载校核时，这条命令比手工改十遍文件靠谱得多。

批量模式适合回归验证：

python run_pf.py --case-dir cases/ --method nr --out results/batch # 任何一个算例不收敛，进程返回非零退出码，方便挂到 CI 上

4.3 输出报表与跟参考解对齐的方法

跑完之后，results/case30/下会有三张表和一份汇总：

results/case30/ ├── bus_result.csv # bus_id, Vm, Va_deg, Pg, Qg, Pd, Qd ├── branch_result.csv # f, t, Pf, Qf, Pt, Qt, Ploss, Qloss, loading ├── summary.json # 迭代次数、末次误差、总网损、最大负载率 └── converge.log # 每轮不平衡量，便于画图

branch_result.csv里的Ploss/Qloss是两端功率之差，loading是max(|Pf|,|Pt|)/rateA，这条列直接决定你要不要回头改网架。所有报表的功率单位是 MVA/MW，电压是 p.u. 加 kV 两列，省得自己乘baseKV。

结果到底可不可信，靠对参考解。ref/case14_voltage_ref.csv存了 14 节点测试系统的标准电压解，写一段校验脚本：

import json import numpy as np import pandas as pd ref = pd.read_csv("ref/case14_voltage_ref.csv") # 列：bus_id, vm_ref, va_ref cur = pd.read_csv("results/case14/bus_result.csv") err_v = np.max(np.abs(cur["Vm"].values - ref["vm_ref"].values)) err_a = np.max(np.abs(cur["Va_deg"].values - ref["va_ref"].values)) print(f"电压幅值最大偏差 {err_v:.3e} p.u.，相角最大偏差 {err_a:.3e} 度") assert err_v < 1e-6, "电压幅值对不上，优先查标幺化和 Ybus" assert err_a < 1e-4, "相角对不上，优先查变压器移相角和变比方向" # 功率平衡校验：总注入 = 总负荷 + 总网损（精度取决于 tol） loss = pd.read_csv("results/case14/branch_result.csv")[["Ploss", "Qloss"]].sum() summary = json.load(open("results/case14/summary.json")) gap = abs(summary["total_ploss"] - loss["Ploss"]) assert gap < 1e-4, f"网损汇总不一致，差 {gap:.3e} MW"

这两条断言是我每次改完代码都会跑的。第一条卡电压，第二条卡网损，两者同时通过，基本可以认定结果可信。相角容差放到1e-4度是因为不同实现的角度基准可能差一点点，幅值容差1e-6p.u. 则没什么商量余地。

注意：参考解文件里的节点顺序必须和算例一致。有些参考文件按bus_id排序，而算例里是乱序的，直接按行号比较会得到一堆莫名其妙的偏差。稳妥做法是先merge到bus_id上再比。

5. 避坑与排查：五个让我返工过的问题

潮流计算程序的坑有一个共同特征：程序不报错，结果看着也挺合理。下面五条都是我自己或者身边的人真实踩过的，按“现象 → 原因 → 解决”写，遇到了可以对号入座。

5.1 基准值不一致：电压算出来是 0.4 而不是 1.0

现象：第一个算例跑完，bus_result.csv里一半节点的Vm在 0.3 到 0.5 之间，但收敛正常、迭代 4 次、网损也在合理范围。

原因：标幺化做了两次。caseio里已经按baseMVA=100把功率折算成 p.u.，而run_pf.py在传给求解器

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