简介面向电力系统稳态分析的C潮流计算通用程序围绕牛顿法求解大规模节点网络这一核心引入tinny2稀疏矩阵技术降低存储与计算开销支持txt格式数据导入导出适合电力专业学生、科研人员与工程技术人员学习算法实现、验证标准算例。压缩包共46个文件、2.23MB包含cpp/h源码、可执行程序、详细设计报告doc以及14、30、57、118、300节点等标准系统的txt输入数据与xlsx结果表格代码、数据、文档分块清晰便于对照调用。已有840人学习。通过阅读源码可理解牛顿迭代与稀疏矩阵存储的具体配合方式直接运行算例可检验不同规模网络的收敛效果配套报告还整理了设计思路、算法细节与结果分析能帮助读者从数学原理走到工程实现是一份适合教学自学与项目起步的完整参考资料。1. 电力系统潮流计算通用程序为什么要自己造这个轮子拿到一张电网节点数据表很多人的第一反应是丢进现成软件点一下等结果出来直接抄数。但真到自己写规划算法、批量改数据、对接非标准格式时现成工具就变成了黑匣子——你不知道它内部怎么编号、怎么处理变压器变比、收敛判据是多少换一套数据就提心吊胆。电力系统潮流计算通用程序要解决的就是这个问题把数据解析、网络建模、非线性求解、结果输出这条链路全部握在自己手里无论换节点规模、换拓扑还是换数据来源都能稳定给出电压幅值、相角和支路功率分布。潮流计算本质上是在给定网络参数和各类节点注入条件下求解一组非线性功率平衡方程。所谓“通用”意味着程序不绑定某一张特定网架能自适应 PQ、PV、平衡节点的混合、变压器支路、环网和放射状结构。“通用”二字也正是一道分水岭脚本化的专用程序容易写但换一张图就崩真正通用的程序难在数据模型和收敛控制而不在牛顿法本身。这篇文章适合三类人电力系统方向的学生想把手动推导变成可复现代码做调度、规划、可靠性分析的工程师需要一个能嵌进自己程序的求解内核以及想摆脱商业软件约束、自己掌握计算细节的团队。我会沿着原理选型、数据格式、代码实现、避坑、工程化的顺序把一套能跑的通用程序拆开讲清楚。2. 牛顿-拉夫逊法通用程序默认引擎的选型理由与极坐标修正方程2.1 三类节点已知量、求解量与统一编号约定潮流计算的第一步不是写迭代而是把每个节点按“已知什么、待求什么”分成三类。这个分类直接决定方程组的维度和雅可比矩阵的结构。PQ 节点已知有功和无功待求电压幅值和相角典型角色是负荷节点也可能含固定无功出力的发电节点PV 节点已知有功和电压幅值待求无功和相角典型角色是有调压能力的发电节点需要注意的是它带有一组无功上、下限平衡节点习惯上也叫 Vθ 节点已知电压幅值和相角待求有功和无功全网必须且只能有一个它负责吸收整个系统的不平衡功率包括网损的差额。节点类型已知量待求量典型设备PQP、QV、δ负荷、无功固定的发电机PVP、VQ、δ有调压能力的发电机平衡/VθV、δP、Q平衡机组全网唯一个这套分类是几乎所有通用潮流程序的内核约定。工程上还有一个容易忽略的约束PV 节点的无功出力不能超过发电机上下限一旦越限节点会在迭代中从 PV 转为 PQ这一点我会在第 5 章和第 6 章分别讲坑和自动处理。程序内部统一用编号 0/1/2 表示 PV/PQ/Vθ 也可以但输入文件里建议直接用可读字符串解析时再映射后面的数据结构会体现这个思路。2.2 极坐标修正方程四个分块雅可比矩阵怎么填通用程序为什么默认选牛顿-拉夫逊法而不是高斯-赛德尔两条硬理由一是收敛速度快平启动下通常 5 到 8 次迭代就能达到 1e-6 精度而高斯法在节点多时收敛慢且容易发散二是扩展性好PQ 分解法、最优潮流、灵敏度计算几乎都是复用牛顿法的雅可比矩阵先写牛顿法等于给后续工程化铺路。极坐标下每个节点的电压表达为 V∠δ功率偏差方程写成F(x) 0, x [δ; V]在每次迭代中把 F 在当前点做泰勒展开保留一阶项就得到修正方程 J Δx −F。其中雅可比矩阵 J 是四个分块矩阵的组合本文采用 ΔV/V 作为电压未知量这是工程程序里的常见选择公式最干净分块含义非对角元i≠j对角元H∂P/∂δViVj(Gij sinθij − Bij cosθij)−Qi − BiiVi²N(∂P/∂V)·VViVj(Gij cosθij Bij sinθij)Pi GiiVi²J∂Q/∂δ−ViVj(Gij cosθij Bij sinθij)Pi − GiiVi²L(∂Q/∂V)·VViVj(Gij sinθij − Bij cosθij)Qi − BiiVi²这里 G、B 分别是导纳矩阵的实部和虚部θij θi − θj。我把未知量取成 ΔV/V更新时 V_new V_old ΔV其中 ΔV V_old × ΔV/V代码里直接用增量乘以 V 再叠加避免把两种未知量约定混在一起。这是写程序时最容易翻车的细节同一个公式在不同教材里有 ΔV 和 ΔV/V 两套版本混着抄必然收敛不了。2.3 收敛判据与迭代参数的务实选择收敛判据决定了程序“什么时候停下来”。我一般用功率偏差的无穷范数所有非平衡节点的有功偏差、所有 PQ 节点的无功偏差最大绝对值小于 tol 就判定收敛默认 tol 1e-6标幺值。这个精度对绝大多数工程分析足够需要和外部程序对账时改成 1e-8但要注意浮点精度和迭代次数的配合。最大迭代次数默认给 30 次牛顿法正常迭代不会超过 10 次。如果超过 20 次还不收敛通常是数据或代码有问题而不是“再多给几次就能行”盲目提上限只会掩盖错误。初值方面行业里的常见做法叫平启动所有 PQ 节点电压幅值取 1.0、相角取 0PV 节点用给定电压和相角 0平衡节点不动。这套初值对大多数输电网都稳定比“按经验给任意值”可靠得多而冷启动给很低电压初值只在极少数高 R/X 比配电系统里有优势通用程序默认平启动即可。3. 数据模型先行节点支路文件格式与解析代码通用性的第一道坎3.1 节点数据与支路数据一个可直接运行的 6 节点算例通用程序能不能换张图就运行关键在输入数据模型。我见过不少代码把节点和支路硬编码在脚本里换系统等于改代码这根本不叫通用。常见做法是仿照工程潮流卡的风格用文本文件分块描述节点和支路节点块里每行一个节点支路块里每行一条支路以中括号行作为分节标志。下面是一个 6 节点算例的完整输入文件节点 1 是平衡节点节点 2、3 是 PV其余是 PQ支路里特意混入一条带变比的变压器支路和一个环网结构# 节点数据: 编号 类型(平衡/VTH, 发电机/PV, 负荷/PQ) 有功注入 无功注入 电压幅值 相角(度) [NODE] 1 VTH 0.0 0.0 1.05 0.0 2 PV 0.6 0.0 1.00 0.0 3 PV 0.3 0.0 1.00 0.0 4 PQ -0.2 -0.10 1.0 0.0 5 PQ -0.1 -0.05 1.0 0.0 6 PQ -0.15 -0.08 1.0 0.0 # 支路数据: 首端节点 末端节点 电阻R 电抗X 对地电纳B/2 变比tap # 变比tap仅变压器支路非1tap定义在i侧k高压侧电压/低压侧电压标幺值 [BRANCH] 1 2 0.019 0.072 0.026 1.0 1 4 0.013 0.048 0.019 1.0 2 3 0.003 0.024 0.0 0.96 2 5 0.015 0.062 0.020 1.0 3 6 0.009 0.042 0.015 1.0 4 5 0.011 0.043 0.017 1.0这份文件是全程序的第一份“契约”。注意三个约定第一功率注入以标幺值表示发电为正、负荷为负和工程习惯一致第二节点编号可以不是从 1 开始的连续整数但解析时要做内部重编号原始编号要保留下来方便和外部结果对照第三变比 tap 是支路首端i 侧归算的标幺值物理含义是 i 侧电压除以 j 侧电压记住这个方向定义第 5 章会专门讲它的坑。3.2 文本解析要点注释、空行与类型码映射解析文本的目标不是“把文件读进来”而是“把外部可读格式转成程序内部的数据结构”。我习惯用按行解析加字符串分割不用复杂正则理由很简单潮流数据文件是给人和程序共同看的格式保持简单反而最不容易出错。下面这段解析代码覆盖了注释、空行、分节和类型映射import numpy as np def parse_case_file(path): buses [] # 每个元素是 dict描述一个节点 branches [] # 每个元素是元组描述一条支路 node_map {} # 原始编号 - 内部连续编号 section None with open(path, r, encodingutf-8) as f: for raw_line in f: line raw_line.split(#, 1)[0].strip() # 去掉注释忽略空行 if not line: continue if line.startswith([): section line[1:-1].strip().upper() continue parts line.split() if section NODE: orig_idx int(parts[0]) internal_idx len(buses) node_map[orig_idx] internal_idx buses.append({ type: parts[1], # VTH / PV / PQ Pspec: float(parts[2]), Qspec: float(parts[3]), V: float(parts[4]), theta: float(parts[5]) * np.pi / 180.0 }) elif section BRANCH: i_orig int(parts[0]) j_orig int(parts[1]) branches.append(( node_map[i_orig], node_map[j_orig], # 映射为内部连续编号 float(parts[2]), # R float(parts[3]), # X float(parts[4]), # B/2 float(parts[5]) # tap )) return buses, branches这段代码有两个值得留意的参数行为。一是split(#, 1)[0]会把每行井号之后的内容全部丢弃所以注释可以写在行尾不需要单独处理二是节点编号映射发生在 NODE 部分解析时而 BRANCH 部分引用的编号一定是前面已经出现过的节点否则node_map会抛 KeyError这种“报错即提示数据乱序”的行为其实是好事能尽早暴露数据问题。内部编号从 0 开始连续递增正好匹配 NumPy 数组下标后面构建导纳矩阵时不再碰原始编号。3.3 为什么不用数据库可读、可 diff、可被外部程序生成有人会问输入数据为什么不用数据库或者 JSON我的习惯是数据库适合管理海量历史数据但不适合作为潮流程序的直接输入JSON 嵌套结构反而不如清晰的表格文本直观。文本格式最大的三个好处第一是人和程序都能读拿着文件就能手算核对某条支路第二是能直接丢进版本管理工具做 diff数据改动在哪里一目了然第三是外部程序生成非常容易不管上游是 Excel 导出、数据库查询还是优化算法输出拼成这个约定格式只需循环打印。等到程序真正工程化需要处理成百上千个算例时再在这个文本格式外面包一层数据管理也不迟。4. 用 Python 手写极坐标牛顿-拉夫逊法导纳矩阵、迭代主循环与参数调整4.1 构建导纳矩阵 Ybus普通支路与变压器的统一处理导纳矩阵是潮流计算的骨架所有功率计算和雅可比矩阵都建立在它之上。这里采用统一的两端口模型普通线路的 tap 为 1变压器支路通过变比 tap 修正导纳值。变比定义在 i 侧的物理模型是i 侧电压经过理想变压器折算后与 j 侧相连串联导纳 y 1/(RjX)对地电纳 B/2 分摊到两端。构建函数如下def build_ybus(buses, branches): 根据节点和支路数据构建复数导纳矩阵 Ybus。 n len(buses) Y np.zeros((n, n), dtypecomplex) for i, j, r, x, b_half, tap in branches: y 1.0 / (r 1j * x) # 串联导纳 b 1j * b_half # 对地电纳的 jB/2 形式 # 理想变压器在 i 侧tap 为标幺变比 Y[i, i] y / (tap * tap) b Y[j, j] y b Y[i, j] -y / tap Y[j, i] -y / tap return Y这段代码的关键在变压器处理Y[i,i]除以 tap²、Y[i,j]和Y[j,i]除以 tapY[j,j]不变。如果把 tap 定义到 j 侧所有公式要反过来后文第 5 章会单独讲这个坑。另一个参数细节是b的写法B/2在数据文件里是实数电纳程序里要乘虚数单位变成jB/2如果解析时已经把它当复数这里就不该再乘数据格式与代码实现必须严格对应这是我自己踩过的血泪经验。4.2 迭代主循环功率偏差、雅可比组装与状态更新有了 Ybus迭代主循环就清晰了。程序先按平启动给全部节点赋初值然后在每轮迭代里完成四件事计算当前电压下的节点注入功率、比较与给定值的偏差、组装雅可比矩阵、解线性方程并更新状态。下面这段是完整可运行的核心循环为了让新手能对照 2.2 节的公式阅读我把雅可比矩阵拆成 H、N、J、L 四个分块分别填充def calc_s(Y, V, theta): 由当前电压向量计算各节点注入有功、无功。 Vc V * np.exp(1j * theta) Ic Y Vc Sc Vc * np.conj(Ic) return Sc.real, Sc.imag def assemble_jacobian(Y, V, theta, P, Q, non_slack, pq_idx): 组装极坐标牛顿法雅可比矩阵未知量取 Δδ 与 ΔV/V。 n len(V) G, B Y.real, Y.imag m len(non_slack) # 相角待求量个数 k len(pq_idx) # 电压幅值待求量个数 H np.zeros((m, m)); N np.zeros((m, k)) J np.zeros((k, m)); L np.zeros((k, k)) # H 分块∂P/∂δ for a, i in enumerate(non_slack): for b, j in enumerate(non_slack): if i j: H[a, b] -Q[i] - B[i, i] * V[i] ** 2 else: dij theta[i] - theta[j] H[a, b] V[i] * V[j] * (G[i, j] * np.sin(dij) - B[i, j] * np.cos(dij)) # N 分块(∂P/∂V)·V for a, i in enumerate(non_slack): for b, j in enumerate(pq_idx): if i j: N[a, b] P[i] G[i, i] * V[i] ** 2 else: dij theta[i] - theta[j] N[a, b] V[i] * V[j] * (G[i, j] * np.cos(dij) B[i, j] * np.sin(dij)) # J 分块∂Q/∂δ for a, i in enumerate(pq_idx): for b, j in enumerate(non_slack): if i j: J[a, b] P[i] - G[i, i] * V[i] ** 2 else: dij theta[i] - theta[j] J[a, b] -V[i] * V[j] * (G[i, j] * np.cos(dij) B[i, j] * np.sin(dij)) # L 分块(∂Q/∂V)·V for a, i in enumerate(pq_idx): for b, j in enumerate(pq_idx): if i j: L[a, b] Q[i] - B[i, i] * V[i] ** 2 else: dij theta[i] - theta[j] L[a, b] V[i] * V[j] * (G[i, j] * np.sin(dij) - B[i, j] * np.cos(dij)) return np.block([[H, N], [J, L]]) def run_nr(buses, branches, tol1e-6, max_iter30): 牛顿-拉夫逊法潮流主循环返回最终电压、相角与迭代信息。 n len(buses) Y build_ybus(buses, branches) # 分类索引 slack_idx [i for i, b in enumerate(buses) if b[type] VTH] pv_idx [i for i, b in enumerate(buses) if b[type] PV] pq_idx [i for i, b in enumerate(buses) if b[type] PQ] non_slack pv_idx pq_idx # 平启动初值 V np.array([b[V] if b[type] ! PQ else 1.0 for b in buses], dtypefloat) theta np.array([b[theta] if b[type] VTH else 0.0 for b in buses], dtypefloat) Pspec np.array([b[Pspec] for b in buses], dtypefloat) Qspec np.array([b[Qspec] for b in buses], dtypefloat) for it in range(max_iter): P, Q calc_s(Y, V, theta) # 功率偏差用 calc - spec求解时对右侧取负 res_P P[non_slack] - Pspec[non_slack] res_Q Q[pq_idx] - Qspec[pq_idx] rhs -np.concatenate([res_P, res_Q]) if np.max(np.abs(rhs)) tol: print(fconverged in {it1} iterations) return V, theta, it 1 J assemble_jacobian(Y, V, theta, P, Q, non_slack, pq_idx) dx np.linalg.solve(J, rhs) m len(non_slack) theta[non_slack] dx[:m] V[pq_idx] * (1.0 dx[m:]) raise RuntimeError(not converged within max_iter iterations)代码里的功率偏差有意写成“计算值减给定值”而修正方程右边取负这是为了保证雅可比矩阵各分块的符号与 2.2 节公式完全一致。很多初学实现翻车就是因为这里混用了“给定值减计算值”的写法和教材里的雅可比符号最后整个矩阵差一个负号表现为迭代发散或振荡。np.linalg.solve用的是稠密矩阵直接分解对 500 节点以内的系统完全够用再往上走要换稀疏求解器那是第 6 章的内容。4.3 三个必调参数初值、收敛精度与最大迭代次数运行主循环之前有三个参数值得花时间理解因为它们决定了程序在不同系统上的表现。第一个是初值代码里采用的平启动对输电网非常稳健但如果你拿到的是一个重负荷配电系统平启动可能让迭代卡在局部振荡区这时可以尝试冷启动把 PQ 节点电压初值改成 0.9 再跑第二个是收敛精度 tol1e-6 是默认值适合大多数工程分析如果结果要和外部程序严格对齐改到 1e-8 或 1e-10 都行但要同时放宽最大迭代次数否则可能在精度和迭代上限之间卡死第三个是最大迭代次数30 次足够如果 10 次以内没收敛先别急着加次数优先检查导纳矩阵的正负号、节点类型设置和变比方向。调试时还有一个很有用的习惯每轮迭代打印rhs的最大绝对值。正常情况它应该每轮下降几倍到几十倍呈单调收敛如果它先降后升说明雅可比矩阵有符号错误如果它下降缓慢说明数据里可能有病态的高 R/X 支路。这套“看中间量”的方法比只看最终报错信息高效得多。5. 潮流计算避坑排查最容易让通用程序翻车的 5 个细节5.1 现象节点编号不连续导致奇异矩阵现象是程序在第一次解方程时抛出奇异矩阵错误或者干脆算到一半结果离谱。原因很典型某开发者在解析数据时直接用原始设备编号作为数组下标比如节点编号是 1、3、7、9但程序只创建了 4 个节点的数组下标 7 越界或者留下大量全零行雅可比矩阵自然奇异。解决方法是程序内部一律使用连续编号原始编号只保留在映射字典里也就是第 3 章代码中的node_map方案。需要强调的是这个编号映射不是可选项而是通用程序的必需品——你永远不知道下一个算例的节点编号长什么样。5.2 现象变压器变比搞反收敛了结果也不对现象比奇异矩阵更隐蔽程序正常收敛但算出来的节点电压和支路功率与外部校验数据系统性偏大或偏小尤其是变压器两侧的电压比例不对。原因几乎都出在变比定义方向上。数据文件里 tap 0.96 时如果程序把它当成 j 侧归算而实际约定是 i 侧归算算出的高压侧电压会差一个 0.96 的倍数。解决方法分两步第一步程序里只保留“tap 定义在 i 侧”这一种约定文件头部用注释写明第二步解析时把 tap 的数值原样传入build_ybus不要在解析函数里做任何符号调整这样一旦数据文件发生方向变化只需要改文件约定程序代码不用跟着改。这是我能给出的最稳妥的后悔药方案。5.3 现象PV 节点无功越限迭代在几个值之间来回跳现象是程序既不收敛也不发散每轮迭代的功率偏差在某个水平来回振荡电压幅值在一个节点附近小幅抖动。原因基本可以锁定为 PV 节点的无功出力超过了发电机实际能力但程序仍然强行维持电压幅值不变导致数学上有解而物理上无解。解决方法是用“越限转 PQ”策略每轮迭代后检查所有 PV 节点的无功注入如果超过上限或低于下限就把该节点类型改成 PQ并把无功给定值锁定在越限边界值然后继续迭代。要注意的是一旦转成 PQ整个迭代过程中不要再转回 PV否则容易在两种状态之间反复切换这就是振荡的来源。第 6 章会给出自动切换的具体实现思路。5.4 现象判据设成 1e-12程序“卡死”在迭代上限现象是迭代次数到了 30 次还没收敛但打印出来的偏差已经小到 1e-10 级别肉眼看上去结果完全可用了。原因很简单浮点数的有限精度让偏差很难稳定压到 1e-12尤其在节点数较多时求和运算带来的舍入误差就足以让残差在 1e-11 左右波动。解决方法是把收敛判据和迭代次数看作一对组合tol 设为 1e-6 或 1e-8max_iter 保持 30 次如果确实需要更高精度用双精度浮点并把 max_iter 加到 100 次但这种情况在工程分析中几乎不存在。把判据设得过严不是严谨而是给自己找麻烦。5.5 现象换成 PQ 分解法后重负荷系统不收敛现象是同一个算例牛顿法 5 次迭代就收敛换成快速分解法PQ 分解后要么不收敛要么收敛精度差一截。原因在于 PQ 分解法把雅可比矩阵里的 N、J 分块置零并使用固定的 B′、B″ 矩阵近似这些近似只有在输电网络高电压比、低 R/X 比、轻负荷条件下才成立一旦系统重负荷或存在大量电缆线路近似误差就会让迭代发散。解决方法是通用程序默认使用牛顿-拉夫逊法PQ 分解只作为“可选加速开关”保留并且至少要加一个 R/X 比检查当某条支路 R/X 大于 0.3 时自动退回牛顿法。这是工程经验和数学边界共同给出的判断不要迷信“更快的方法一定更好”。6. 从通用走向实用稀疏求解、PV 自动转 PQ 与灵敏度验证6.1 用稀疏矩阵与节点重编号处理大规模网络当节点数超过一两千稠密矩阵的内存和求解时间都会迅速失控这时要把np.linalg.solve换成稀疏求解器。常见做法是先用scipy.sparse.lil_matrix组装 Ybus 和雅可比矩阵再转成csc_matrix配合spsolve求解如果能用节点重编号把矩阵带宽压下来比如按电气距离做参考排序则迭代稳定性和求解速度都会明显改善。雅可比矩阵的稀疏模式和 Ybus 几乎一致沿用同一套索引即可。这一步能让程序从“能跑”变成“能用”。6.2 PV 转 PQ 的自动切换逻辑工程级的通用程序不会因为 PV 无功越限就罢工。做法是在第 4 章迭代循环里加一段“越限检查”每轮更新完状态后对当前仍标记为 PV 的节点计算 Q如果 Q 大于上限或小于下限就把该节点改为 PQ 类型把 Qspec 固定为越限值然后立即重建雅可比矩阵并继续迭代。这里的关键参数是切换的触发时机每轮都检查但一旦切换就保持 PQ 到底不做回切直到收敛后输出时再提醒“哪些节点发生了类型切换”。这样既保证收敛也保留了排查轨迹。6.3 用雅可比矩阵校验结果与获得灵敏度一次收敛的潮流结果不仅要“电压合理”还要满足全网功率平衡所有节点的注入功率之和近似等于网损平衡节点的出力落在合理范围内。雅可比矩阵在这里还有第二重用途——它本身就是电压对注入功率的灵敏度信息收敛后求一次稀疏逆就能得到 dV/dP、dV/dQ 的近似值这是静态安全分析和电压稳定性评估的常用输入。我自己最早写的版本就是栽在符号约定上后来把雅可比矩阵拆成 H、N、J、L 四个小函数对照公式翻车率大幅下降。通用潮流程序没有太多玄学把数据契约定好、符号约定写清楚、收敛过程看仔细它就能成为你后续所有电力系统分析工具的可靠底座。希望帮到你。本文还有配套的精品资源点击获取