简介面向电力系统分析与潮流计算学习者的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 的字段设计数据模型是一套潮流程序的骨架。这套程序用的是电力系统里最通用的“四表结构”字段命名也和主流工具对齐方便你把别处的算例直接导进来。bus 表字段类型含义说明bus_idint节点编号从 1 开始内部索引统一减 1typeint节点类型1PQ2PV3平衡节点Vmfloat电压幅值p.u.初值一般给 1.0Vafloat电压相角度求解器内部转弧度Pd/Qdfloat负荷p.u.Pg/Qgfloat发电机出力p.u.PQ 节点通常为 0Qmax/Qminfloat无功上下限越限切换用baseKVfloat基准电压报表还原用branch 表字段含义备注f, t首端、末端节点号变压器变比默认挂在 f 端r, x, b电阻、电抗、对地电纳p.u.b 是全长充电电纳rateA长期允许容量MVA报表里标越限tap变比0 视为 1.0shift移相角度status投运状态0 表示退出组装 Ybus 时跳过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_shuntTrue): nb case.bus.shape[0] Y lil_matrix((nb, nb), dtypecomplex) for k in range(case.branch.shape[0]): br case.branch[k] if int(br[10]) 0: # status0停运支路直接跳过 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|^2t 侧不变 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()这段代码有三个地方值得反复看。第一是对地电纳要除以 2b在 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, dtypecomplex) 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) # PVQ 是未知量置 0 else: Sbus[i] complex(Pg - Pd, Qg - Qd) # 平衡节点实际不参与方程 return SbusPV 节点的无功为什么要置 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]的四个分块规模是(npvpqnpq) × (npvpqnpq)。手推公式容易出错实践中更稳的写法是先算复数域的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, tol1e-8, max_iter20): 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 - ScalcF -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)我们要F0牛顿步是x_{k1} 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, tol1e-8, max_iter30): 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, tol1e-8, max_iter500): 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 依赖安装与最小可跑路径依赖很轻numpyscipy是核心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的参数不多但每个都影响结果值得对着表过一遍参数默认值说明--case无单个算例文件路径JSON/CSV--case-dir无批量模式目录下所有算例逐个求解--methodnrnr/fdpf/gs--tol1e-8功率不平衡量收敛阈值p.u.--max-iter20迭代上限--sweep-load无负荷倍率扫描格式起:止:步长--outresults输出目录--verbosefalse打印每轮最大不平衡量单算例是最常用的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里已经按baseMVA100把功率折算成 p.u.而run_pf.py在传给求解器本文还有配套的精品资源点击获取