电力系统潮流与短路电流计算程序从零实现指南

发布时间:2026/9/18 17:08:50
电力系统潮流与短路电流计算程序从零实现指南 简介电力系统潮流与短路电流计算是电网规划、运行与控制中的关键环节。这份PDF以3机9节点系统为算例面向电力专业学生、考研备考者及一线工程师系统演示了从节点导纳矩阵形成、发电机与负荷等效建模、阻抗矩阵求逆到正常潮流计算与三相短路电流计算的完整流程。压缩包内为1个PDF文档整包仅132KB轻量便携便于下载后随时查阅。已有363人学习使用。文档在短路部分分别采用精确算法与近似算法并详细对比两种方法下故障点电流、节点电压及支路电流分布同时给出节点导纳矩阵Y、节点阻抗矩阵第4列、电压相角等关键中间结果内容层次清晰、步骤完整兼具理论推导与算例验证可帮助读者深入理解节点导纳法在潮流与短路计算中的实际应用提升电力系统分析与编程实践能力。1. 电力系统潮流及短路电流计算程序先看清那份 PDF 里有什么很多工程文件夹里都躺着一本《电力系统潮流及短路电流计算程序.pdf》它往往不是可以直接安装的软件而是一份把电力系统潮流计算和短路电流计算结果写进计算书的技术文档。PDF 里出现最多的是节点导纳矩阵、牛顿-拉夫逊迭代表、三相短路电流值以及一串手工验算步骤。对于第一次接触的人这份 PDF 最大的价值是描述了标准计算流程但它并不能帮你直接把下一个算例跑出来。这类 PDF 真正值得提取的部分有两块一块是把潮流方程和短路电流公式沉淀成可复用的程序函数另一块是给出了你后续调试程序时用来对结果的验算基准。实际开发时我会把电力系统潮流计算和短路电流计算拆成几个独立模块先建导纳矩阵再跑潮流然后从潮流结果计算故障前电压最后进入短路电流计算。下面按这条链路展开使用标幺值代码以 MATLAB 风格为主必要时给出 Python 版本。2. 潮流计算程序的基础节点导纳矩阵与牛顿-拉夫逊实现2.1 节点导纳矩阵是潮流和短路电流共用的输入潮流计算的起点不是迭代公式而是把网络拓扑变成矩阵。节点导纳矩阵 Y 同时服务于潮流计算和短路电流计算短路电流计算中的节点阻抗矩阵正是由 Y 求逆得到的。手算时自导纳 Yii 是连接在该节点上所有支路导纳之和互导纳 Yij 是两节点间支路导纳的负值。程序里不需要手写这些规则只需要遍历线路参数表function Y build_ybus(branch) % branch: 每行 [首端节点, 末端节点, R(pu), X(pu), 对地导纳B/2(pu)] nb max(max(branch(:, 1:2))); Y zeros(nb, nb); for k 1:size(branch, 1) f branch(k, 1); t branch(k, 2); z branch(k, 3) 1i * branch(k, 4); y 1 / z; Y(f, f) Y(f, f) y 1i * branch(k, 5); Y(t, t) Y(t, t) y 1i * branch(k, 5); Y(f, t) Y(f, t) - y; Y(t, f) Y(t, f) - y; end Y sparse(Y); end这段代码的关键有三点。第一参数必须是标幺值电阻电抗必须在同一基准容量下折算否则导纳矩阵整体偏移潮流和短路电流都会出错。第二对地导纳 B/2 通常由线路充电电容产生程序里要分别加到首端和末端节点。第三最后强制转成稀疏矩阵因为实际电网动辄几百上千节点全矩阵存储会让后面每一步都变慢。潮流计算里节点功率方程使用的是极坐标形式。对节点 i注入功率满足P_i V_i * sum_j V_j * (G_ij * cos(theta_i - theta_j) B_ij * sin(theta_i - theta_j)) Q_i V_i * sum_j V_j * (G_ij * sin(theta_i - theta_j) - B_ij * cos(theta_i - theta_j))牛顿-拉夫逊法的思路是在当前电压幅值 V 和相角 theta 附近把上述非线性方程做一阶泰勒展开迭代修正各节点的 V 和 theta。迭代修正量由雅可比矩阵 J 和功率不平衡量 dS 组成的线性方程组决定。2.2 用 MATLAB 或 Python 写牛顿-拉夫逊潮流计算程序的最小可运行版本下面给一个三节点系统的 Python 版本。节点 1 是平衡节点节点 2 和节点 3 都是 PQ 节点负荷作为负的注入功率。为了把注意力放在算法结构上这里用数值差分求雅可比矩阵避免手推十几行偏导公式时出现符号错位import numpy as np Ybus np.array([ [6-20j, -412j, -28j], [-412j, 8-22j, -410j], [-28j, -410j, 6-18j] ], dtypecomplex) # 节点注入功率发电为正负荷为负 S_spec np.array([00j, -1.5-0.3j, -0.6-0.1j], dtypecomplex) V np.array([1.0, 1.0, 1.0]) # 电压幅值 th np.array([0.0, 0.0, 0.0]) # 电压相角 def calc_S(V, th): ph V * np.exp(1j * th) return ph * np.conj(Ybus ph) def F(x): th[1] x[0] th[2] x[1] V[1] x[2] V[2] x[3] S calc_S(V, th) return np.array([ S[1].real - S_spec[1].real, S[2].real - S_spec[2].real, S[1].imag - S_spec[1].imag, S[2].imag - S_spec[2].imag ]) x np.array([0.0, 0.0, 1.0, 1.0]) tol 1e-8 for it in range(15): f F(x) if np.max(np.abs(f)) tol: print(converged at iter, it) break J np.zeros((4, 4)) h 1e-7 for j in range(4): xp x.copy() xm x.copy() xp[j] h xm[j] - h J[:, j] (F(xp) - F(xm)) / (2 * h) dx np.linalg.solve(J, -f) x x dx th[1] x[0] th[2] x[1] V[1] x[2] V[2] x[3] print(V , V) print(theta(rad) , th)运行后能看到节点 2 的电压幅值略低于 1.0节点 3 更低这是有功和无功负荷共同作用的结果。代码里的 F 函数返回四个残差节点 2 和节点 3 的有功不平衡量、无功不平衡量。数值差分雅可比矩阵虽然比解析雅可比慢但胜在维护简单适合做教学和测试。生产环境建议换成解析雅可比或调用稀疏矩阵求解器。提示这段代码里的 Ybus 是示意参数实际工程中应该由 build_ybus 根据线路参数自动生成不要手工填写。2.3 PQ 分解法什么时候用参数怎么设当系统规模超过 500 个节点并且输电线路的 R/X 比值较小时牛顿-拉夫逊法的雅可比矩阵每次迭代都要重新计算代价较大。PQ 分解法把有功和无功方程解耦用两个常数矩阵 B 和 B 代替雅可比矩阵迭代格式简化为delta_theta inv(B) * (delta_P / V) delta_V inv(B) * (delta_Q / V)写成程序时B 是节点电纳矩阵去掉平衡节点后的子矩阵B 则进一步忽略影响较小支路。收敛判据仍然使用功率不平衡量最大值通常取 1e-6 pu。PQ 分解法对重负荷系统和低压配电网可能不收敛当线路 R/X 大于 0.3 时我一般会切回牛顿-拉夫逊法。这个边界条件在 PDF 计算书里很少写但它决定了程序在哪些算例上可信。3. 短路电流计算程序的核心从等值阻抗到对称分量法3.1 先拿到故障点的戴维南等值阻抗短路电流程序不是从零开始的独立模块它需要潮流计算的结果作为输入。具体说短路发生前故障节点的电压 V_f 由潮流计算给出短路发生后故障点被看作一个电压源其源内阻抗就是系统戴维南等值阻抗 Z_th。求 Z_th 最直接的方法是构建节点阻抗矩阵Z full(inv(Ybus)); Zth Z(k, k); % k 是故障节点编号大系统中不要直接求逆而是对 Ybus 做 LU 分解再用单位向量回代求解 Z 的第 k 列。节点阻抗矩阵的对角元 Zkk 把网络拓扑信息压缩成了一个复数三相短路电流的幅值可以表示为I_f V_f / Zkk如果潮流计算给的 V_f 是标幺值I_f 也是标幺值。输出有名值时还需要乘以基准电流I_base S_base / (sqrt(3) * U_base)。很多计算程序跑出的短路电流大得离谱原因就是漏了这一步单位换算。3.2 三相短路电流计算的代码写法三相短路是对称故障只需要正序阻抗。下面是一个函数封装def three_phase_fault(V_pre, Z1): V_pre: 故障前节点电压(pu) Z1: 正序节点阻抗矩阵对角元(pu) 返回短路电流的幅值和相角 If V_pre / Z1 return abs(If), np.angle(If, degTrue)参数上要注意V_pre 是复数Z1 也是复数不能只取模相除。程序返回相角不是可选项保护整定时需要知道短路电流相对于参考节点的相位。比如 V_pre 1.02 0.05jZ1 0.01 0.2j短路电流会明显滞后电压这在单机无穷大系统里可以手算验证。3.3 不对称短路单相接地和两相短路的对称分量计算不对称短路必须引入负序 Z2 和零序 Z0。工程上通常认为负序阻抗约等于正序阻抗但零序阻抗差异很大它取决于变压器接线方式、中性点接地方式、线路是否带架空地线。计算程序里需要维护三个序网导纳矩阵不能只给一个 Ybus。核心公式如下三相短路: If_a E / Z1 单相接地: If_a 3E / (Z1 Z2 Z0) 两相短路: If_a sqrt(3) * E / (Z1 Z2)代码实现可以这样def asymmetric_fault(fault_type, E, Z1, Z2, Z0): if fault_type SLG: If 3 * E / (Z1 Z2 Z0) elif fault_type LL: If np.sqrt(3) * E / (Z1 Z2) else: raise ValueError(unsupported fault type) return abs(If), np.angle(If, degTrue)程序跑出来之后最先要怀疑的是零序参数。零序阻抗如果直接复用正序阻抗单相接地电流会小很多而且相位完全错误。我见过的短路电流计算程序里排第一的坑就是零序网络构建不完整。4. 把 PDF 里的计算步骤变成可运行脚本输入输出与数据导入4.1 “未在本地计算机上注册 Microsoft.ACE.OLEDB.12.0 提供程序”及更稳的数据导入方式计算程序写好后最耗时间的反而是把节点和线路参数送进去。很多人在 MATLAB 里用 xlsread 读 Excel报错信息是未在本地计算机上注册“microsoft.ace.oledb.12.0”提供程序。这个报错的原因通常是系统里没有安装 Microsoft Access Database Engine或者装了 64 位版本而 Excel 是 32 位版本位数不匹配。与其折腾 OLEDB我会直接用 CSV 文件或者让 Python 的 openpyxl 引擎去读 Excelimport pandas as pd branch pd.read_csv(branch.csv, sep,) node pd.read_excel(node.xlsx, engineopenpyxl)engine 参数指定为 openpyxl 后读取 xlsx 文件不依赖系统的 ADO/OLE DB 组件在 Linux 服务器上也能运行。如果项目要求必须在 MATLAB 里完成可以把 Excel 另存为 CSV再用 readtable 读取branch readtable(branch.csv); branch_array table2array(branch);这一小步改动能删掉一整个类别的环境依赖问题。数据文件里建议固定列顺序节点编号、线路首端、末端、R、X、B、基准电压。列名和单位在 PDF 计算书里写清楚程序才不会被脏数据污染。4.2 批量计算多个故障点的短路电流函数实际工程不是算一个节点而是要把目标母线列表全部扫一遍。可以用一个批量函数封装前面的单点计算def batch_fault_calc(fault_type, bus_list, Z1, Z2, Z0, V_pre_map): rows [] for k in bus_list: z1 Z1[k, k] z2 Z2[k, k] z0 Z0[k, k] v_pre V_pre_map[k] if fault_type 3ph: If v_pre / z1 elif fault_type SLG: If 3 * v_pre / (z1 z2 z0) elif fault_type LL: If np.sqrt(3) * v_pre / (z1 z2) else: continue rows.append({ bus: k, If_pu: abs(If), angle_deg: np.angle(If, degTrue), V_pre_real: v_pre.real, V_pre_imag: v_pre.imag }) return pd.DataFrame(rows)这个函数的好处是把潮流计算结果、序阻抗矩阵和故障类型组合成一个统一入口。后续如果想要输出短路电流有名值只需传入 S_base 和 U_base再乘以基准电流。故障类型用字符串参数区分扩展新故障类型时只需要在分支里增加一条公式。4.3 用已知算例校核你的程序程序写完不等于算得对。我会先用一个最小的单机-无穷大系统做基准测试设系统戴维南等值为 Z_th 0.1 0.3j故障前电压 V_f 1.0 pu那么三相短路电流理论幅值为 1.0 / 0.3162 3.16 pu。程序输出必须落到这个数值附近。潮流程序同样可以用残差法校核将输出的 V 和 theta 代回功率方程计算每个节点的注入功率与设定值比较差值应在 1e-6 量级。如果差值很大优先检查节点编号是否从 1 开始是否存在重复编号以及导纳矩阵是否漏掉了并联支路。校验对象理论值程序输出判断三相短路电流3.16 pu3.1601 pu通过节点2电压0.98 pu0.9803 pu通过零序电流0.85 pu0.93 pu失败复查零序参数零序这一项最容易翻车因为手册里的公式和线路零序参数常常不在同一个表里程序报错前先检查输入数据。5. 进阶技巧从潮流和短路电流结果里定位系统薄弱点5.1 用潮流雅可比矩阵的最小奇异值判断电压稳定性单独的电压幅值不足以说明母线强弱。我会把最后一次迭代的雅可比矩阵 J 做奇异值分解取最小奇异值 sigma_min 作为电压稳定性指标。sigma_min 越接近零说明系统越接近电压崩溃点。这个指标比看电压低于 0.95 pu 更灵敏因为它反映了网络结构对电压的支撑能力。程序中只需要对 J 调用np.linalg.svd把 sigma_min 与短路电流结果放在一起看。5.2 短路电流相角比幅值更早暴露数据错误短路电流幅值相除后受阻抗幅值影响大而相角主要由阻抗角决定。如果程序算出的相角偏离预期超过 10 度通常是序阻抗矩阵中的电阻和电抗比例输反了。单相接地故障时零序阻抗角错误会直接导致相角异常这类问题在幅值上不容易看出来。因此批量计算函数里一定要保留angle_deg字段校核时同时看幅值和相角不要只输出绝对值。5.3 把电压、短路电流和雅可比奇异值合并成一张诊断表实用技巧是把所有计算结果拼进同一个 Pandas DataFrame按照短路电流降序排列diag v_df.merge(i_df, onbus) diag[sigma_min] sigma_min diag.sort_values(If_pu, ascendingFalse, inplaceTrue)排序后排在前面的母线通常是系统短路容量最大的位置也是设备选型中最容易超标的点。把潮流电压低和短路电流大的母线找出来再对照 sigma_min会比单看任何一张表都高效。这个联动分析没有写在传统 PDF 计算书里但它是把两个程序模块变成真正决策工具的关键一步。本文还有配套的精品资源点击获取