简介面向电力系统潮流计算学习者与工程人员这份MATLAB源码基于牛顿-拉夫森迭代法支持IEEE 6节点与9节点标准测试系统用于分析稳态下的电压幅值、相角以及线路有功无功潮流分布。压缩包中共有2个m文件整体大小仅2KB代码精简适合直接运行和修改。目前已有426人学习该资源。脚本覆盖了网络模型定义、初始值设置、迭代求解、收敛判定与结果输出等关键步骤清晰展示了PQ节点、PV节点及平衡节点的处理方式可帮助读者从代码层面理解潮流计算的数值解法本质。同时两个示例分别针对IEEE 6节点和9节点系统便于对比不同规模网络下的算法表现也可作为课程实验、毕业设计或科研算法验证的基础工具。整体来看这套MATLAB脚本结构清晰复用性强便于快速掌握核心流程。1. 从9节点还是6节点入手先搞清数据再谈潮流计算拿到 power_flow、IEEE 9节点、6节点这类任务的人通常不是不懂潮流方程而是卡在数据和算法的接缝上手里没有一份能直接喂给程序的母线表和支路表或者数据表跟算法对不上。IEEE 9节点WSCC 3机9节点和6节点测试系统是电力系统论文里最常用的两个小型基准算例前者参数全网高度统一、最容易跑通后者在不同教材和论文里存在多个版本、最容易踩坑。这篇笔记打算完整走一遍怎么把母线数据和支路数据整理成程序能读的格式怎么用牛顿-拉夫逊法自己写一个几十行的潮流求解器再用 MATLAB 下的 MATPOWER 做交叉验证。适合刚接手电力系统仿真、准备做配电网或微电网算例以及打算在论文里引用 IEEE 标准算例的工程师和研究生。2. IEEE 6节点与9节点测试系统数据字段、版本差异与选型2.1 两个系统的来头WSCC 3机9节点与教材经典6节点IEEE 9节点系统的正式名字叫 WSCC 3机9节点系统来自美国西部电网协调委员会的一套标准模型MATPOWER 里直接叫 case9。它最值钱的一点是数据收敛全世界的论文、课件、开源仓库用的几乎都是同一套母线阻抗和负荷参数你随便找一篇附录去核对数值都不会差太远。9条母线、9条支路、3台发电机拓扑简单清楚调压器、输电线路、负荷都覆盖到了做算法验证很合适。IEEE 6节点则要多留一个心眼。它最早出自 Wood 和 Wollenberg 合著的《Power Generation, Operation and Control》教材后来论文和课件里的 6 节点算例基本都从这套数据衍生出来但在负荷量、变压器抽头、并联电容上各有微调。你拿到手的 6 节点和另一个课题组拿到的 6 节点算出来的母线电压可能差 0.1% 甚至更多这不是程序错了是数据版本不同。为什么两个系统这么受潮流计算偏爱因为规模小到能做白盒验证。9 节点可以全部手工推导任何一步算错都能顺着矩阵找回来6 节点更小但里面带了两台调压变压器能把潮流计算里最阴险的非标变比问题暴露出来。很多 IEEE 电力系统会议论文在提新算法时都习惯先用这两个系统做自证再上大规模算例业界基本把这当成了约定俗成的流程。2.2 母线/支路数据怎么读一张字段表说清单位与基准不管从哪个渠道拿到数据电网潮流数据最终都可以归成两张表母线表bus和支路表branch。母线表每一行是一个节点支路表每一行是一条线路或一台变压器。理解这两张表的字段比理解算法本身更早一步因为字段单位不一致是最常见的翻车根源。表字段含义单位在算法里的角色母线Bus母线编号1 到 NYbus 的行列索引母线Type2slack1PV0PQ-决定哪些量已知、哪些要迭代母线Pg / Qg发电机有功 / 无功出力MW / MVar节点注入功率的正项母线Pd / Qd负荷有功 / 无功MW / MVar节点注入功率的负项母线Vmag电压幅值初值puPV 和 slack 固定PQ 待求母线Vang相角初值degslack 固定其余待求支路from / to首端 / 末端母线编号-Ybus 非对角元素位置支路R / X支路电阻 / 电抗pu计算串联导纳支路B线路总充电电纳pu组装对地半电容支路tap变压器变比pu非标变比时要修正 Ybus这里最关键的是基准值。通常这两个系统都以 100 MVA 为功率基准所有阻抗、导纳、电压都已经是标幺值只有功率是 MW/MVar 的绝对值。所以在把功率写进方程之前必须除以基准功率否则牛顿法算出来的结果会整体放大 100 倍这是新手必踩的坑。2.3 为什么先用9节点跑通再用6节点验泛化我的建议很明确如果你第一次写潮流程序先用 9 节点。原因是 9 节点数据没有歧义MATPOWER 自带 case9跑出来的结果全网可以互相核对而且它的结构非常清晰1 号母线是平衡节点2、3 号是 PV 节点4 到 9 号是负荷与联络节点一条发电—输电—负荷链路明明白白。用这套数据调通你的牛顿-拉夫逊代码遇到任何偏差都容易定位。6 节点适合作为第二关。它的规模更小一台普通笔记本做上万次潮流计算也就是秒级的事适合批量场景但它的变压器抽头和重负荷会让程序暴露出更多边界问题正好用来检验你的程序对新数据版本够不够稳。我一般把 9 节点当标准答案把 6 节点当压力测试两边都跑过才敢说这个程序可以拿去做别的算例。提示如果你拿到的是 MATPOWER 自带的 case6ww它就是教材那套 6 节点系统数据源头可靠。如果是论文附录里的 6 节点先看它是否写了变压器变比和充电电纳这两项最容易引入版本差异。3. 手写牛顿-拉夫逊潮流Ybus、Jacobian 与两个收敛细节3.1 从支路表到Ybus抽头变比k和充电电容怎么处理潮流计算的第一步永远是组装节点导纳矩阵 Ybus。Ybus 的对角线是各母线自导纳非对角线是互导纳所有支路的阻抗、充电电容、变压器变比都汇总在这里。组装时最容易出错的是两点线路充电电容 B 在数据里给的是总电纳实际并联在两端节点时要各取一半变压器非标变比 tap 不是 1.0 时首端和末端的自导纳修正公式不一样。下面这段代码以 IEEE 9 节点为例先把数据按统一格式定义好再组装 Ybus。注意我这里的支路数据里B 列填的是总充电电纳变比 tap 统一填 1.0因为 case9 的三台变压器都是标称变比。import numpy as np # ---------- IEEE 9节点数据基准 100MVA ---------- # 母线: [类型, Pg, Qg, Pd, Qd, V初值, θ初值(°)] # 类型: 0PQ, 1PV, 2slack bus9 np.array([ [2, 0, 0, 0, 0, 1.040, 0.0], # 1号 平衡母线 [1, 163, 0, 0, 0, 1.025, 0.0], # 2号 PV [1, 85, 0, 0, 0, 1.025, 0.0], # 3号 PV [0, 0, 0, 0, 0, 1.000, 0.0], # 4号 PQ [0, 0, 0, 90, 30, 1.000, 0.0], # 5号 PQ [0, 0, 0, 0, 0, 1.000, 0.0], # 6号 PQ [0, 0, 0, 100, 35, 1.000, 0.0], # 7号 PQ [0, 0, 0, 0, 0, 1.000, 0.0], # 8号 PQ [0, 0, 0, 125, 50, 1.000, 0.0], # 9号 PQ ]) # 支路: [首端, 末端, R(pu), X(pu), B总(pu), 变比tap] br9 np.array([ [0, 3, 0.000, 0.0576, 0.000, 1.0], # 1-4 变压器 [3, 4, 0.017, 0.092, 0.158, 1.0], # 4-5 线路 [4, 5, 0.039, 0.170, 0.358, 1.0], # 5-6 线路 [2, 5, 0.000, 0.0586, 0.000, 1.0], # 3-6 变压器 [5, 6, 0.012, 0.101, 0.209, 1.0], # 6-7 线路 [6, 7, 0.009, 0.072, 0.149, 1.0], # 7-8 线路 [7, 1, 0.000, 0.0625, 0.000, 1.0], # 8-2 变压器 [7, 8, 0.032, 0.161, 0.306, 1.0], # 8-9 线路 [8, 3, 0.010, 0.085, 0.176, 1.0], # 9-4 线路 ]) def build_ybus(bus, br): n len(bus) Y np.zeros((n, n), dtypecomplex) for f, t, R, X, B, tap in br: z complex(R, X) y 1.0 / z if z ! 0 else 1e-12 # 防止除以零 # 统一处理变比k1时退化为普通线路k!1时按变压器π模型 Y[f, f] y / (tap * tap) 0.5j * B Y[t, t] y 0.5j * B Y[f, t] - y / tap Y[t, f] - y / tap return Y逻辑说明支路阻抗 z 由 R 和 X 构成导纳 y 是它的倒数。充电电纳 B 在数据里是总量所以两端各挂一半。变比 tap 的修正只影响首端自导纳和两个互导纳首端自导纳除以 tap 的平方互导纳除以 tap。当 tap 等于 1.0 时这套公式默认退化成普通线路所以不需要单独写 if 分支。参数说明bus9 里每一行最后两个数是电压初值 1.0 pu 和相角初值 0 度这叫做平启动。对绝大多数传输系统平启动配合牛顿-拉夫逊法都能收敛不需要额外做复杂初始化。br9 里的 B 列注意是总充电电纳MATPOWER 的 case9 数据里 4-5 线路给的是 0.158对应每端 0.079我代码里直接对 B 取一半正好对上。3.2 Jacobian四块偏导怎么组装成线性方程组牛顿-拉夫逊潮流的核心思路是把节点注入功率方程在当前电压处做一阶泰勒展开得到失配量 ΔP、ΔQ 与修正量 Δθ、ΔV 之间的线性关系反复求解直到失配量小于阈值。Jacobian 矩阵由四块组成∂P/∂θ、∂P/∂V、∂Q/∂θ、∂Q/∂V。组装时要按状态量的取舍来决定保留哪些行列。平衡母线的 θ 和 V 都是固定的不参与迭代PV 母线的 V 固定、θ 待求无功失配不参与迭代PQ 母线的 θ 和 V 都待求。所以最终线性方程组的规模是角度方向包含所有 PV 和 PQ 母线电压方向只包含 PQ 母线。def newton_pf(bus, br, tol1e-10, max_iter40): n, base len(bus), 100.0 typ bus[:, 0].astype(int) pq np.where(typ 0)[0] # PQ 母线索引 pv np.where(typ 1)[0] # PV 母线索引 slack np.where(typ 2)[0][0] # 平衡母线索引 pvpq np.concatenate([pv, pq]) # 角度待求的母线不含 slack Vm bus[:, 5].astype(float).copy() # 电压幅值 Va np.deg2rad(bus[:, 6].astype(float)).copy() # 相角转弧度 Y build_ybus(bus, br) # 功率注入 发电 - 负荷转标幺 Sbase (bus[:, 1] - bus[:, 3]) / base 1j * (bus[:, 2] - bus[:, 4]) / base Pspec, Qspec Sbase.real, Sbase.imag for it in range(max_iter): V Vm * np.exp(1j * Va) S V * np.conj(Y V) # 当前电压下的注入功率 dP Pspec - S.real # 有功失配 dQ Qspec - S.imag # 无功失配 dF np.concatenate([dP[pvpq], dQ[pq]]) if np.max(np.abs(dF)) tol: break # ---- 组装 Jacobian 四块 ---- J11 np.zeros((n, n)); J12 np.zeros((n, n)) J21 np.zeros((n, n)); J22 np.zeros((n, n)) G, B Y.real, Y.imag for i in range(n): for j in range(n): dth Va[i] - Va[j] if i j: J11[i, i] -S.imag[i] - Vm[i]**2 * B[i, i] J12[i, i] S.real[i] / Vm[i] Vm[i] * G[i, i] J21[i, i] S.real[i] - Vm[i]**2 * G[i, i] J22[i, i] S.imag[i] / Vm[i] - Vm[i] * B[i, i] else: J11[i, j] Vm[i] * Vm[j] * (G[i,j]*np.sin(dth) - B[i,j]*np.cos(dth)) J21[i, j] -Vm[i] * Vm[j] * (G[i,j]*np.cos(dth) B[i,j]*np.sin(dth)) J12[i, j] Vm[i] * (G[i,j]*np.cos(dth) B[i,j]*np.sin(dth)) J22[i, j] Vm[i] * (G[i,j]*np.sin(dth) - B[i,j]*np.cos(dth)) # 按状态量取舍剔除 slack 的相角列和 PV 的电压列 J np.block([ [J11[np.ix_(pvpq, pvpq)], J12[np.ix_(pvpq, pq)]], [J21[np.ix_(pq, pvpq)], J22[np.ix_(pq, pq)]], ]) dx np.linalg.solve(J, dF) Va[pvpq] dx[:len(pvpq)] Vm[pq] dx[len(pvpq):] V Vm * np.exp(1j * Va) Qcal (V * np.conj(Y V)).imag # PV 和 slack 的无功出力 return V, Qcal逻辑说明每次迭代先算当前电压下的注入功率 S然后和给定值做差得到失配量 dF。对角项的四个偏导公式用的是标准形式非对角项通过 θ_i - θ_j 计算。组装时用 np.ix_ 做索引切片把 slack 的相角列和 PV 的电压列从矩阵里剔除这样线性方程组才是可解的方阵。参数说明tol 取 1e-10对应功率失配小于 0.00001 W 级别对 9 节点这种小系统已经足够严格。max_iter 设 40正常 4 到 6 次迭代就能收敛。PV 母线的无功不参与迭代但收敛后要从功率平衡方程里反推出来这就是返回值里 Qcal 的作用——用它可以核对发电机无功是否越限。3.3 收敛判据和初值平启动为什么够用很多人喜欢把收敛判据分成电压偏差和功率偏差两种实际工程里我只认功率失配量。因为电压差 0.0001 pu 在重负荷系统里可能对应很大的功率偏差而功率失配才是潮流方程真正要平衡的东西。我上面代码里用的是 dF 的最大绝对值也就是同时看 ΔP 和 ΔQ 的最大值这个判据简单直接不会出现电压看着收敛了、支路功率还差着一截的假象。初值方面平启动所有 PQ 母线 V1.0、θ0对 9 节点和 6 节点都够用。只有在系统重载到极限、或者支路参数特别悬殊时平启动才会出问题那时候可以把初值改成上一次收敛结果或者先用高斯-赛德尔法迭代几轮再切到牛顿法。对本文的两个系统平启动加牛顿法就够了不用额外做初值优化。4. 用MATPOWER做交叉验证最小算例与结果对拍4.1 loadcase runpf 跑通9节点的最小命令自己写的代码算出来必须有个可信的参照物。MATPOWER 是 MATLAB 环境下的开源潮流计算工具内置了 case9这是业界公认的参考实现。最小验证流程只有三步加载算例、跑潮流、打印结果。% 在 MATLAB 中运行需要先把 MATPOWER 添加到路径 mpc loadcase(case9); % 载入 IEEE 9节点自带算例 res runpf(mpc); % 默认牛顿-拉夫逊法 fprintf(收敛迭代次数: %d\n, res.iterations); for k 1:9 fprintf(母线 %d: V %.4f pu, theta %.3f deg\n, ... k, res.bus(k, 8), res.bus(k, 9)); end逻辑说明loadcase 把 case9.m 解析成 MATPOWER 内部的算例结构体 mpcrunpf 执行潮流计算返回的结果结构体 res 里bus 矩阵第 8 列是电压幅值第 9 列是相角。iterations 字段记录实际牛顿迭代次数正常 9 节点应该在 3 到 5 次之间。参数说明runpf 默认用的是牛顿法如果系统规模变大或者想对比算法可以传入第三个参数切换比如 runpf(mpc, mpoption(PF_ALG, 2)) 切到快速解耦法。但做交叉验证时我建议就用默认设置因为默认参数最接近你自己写的那套标准牛顿法。4.2 对拍三张表电压幅值、相角、支路有功交叉验证不是看一眼差不多就行要落到三张表母线电压幅值、母线相角、支路有功潮流。自己代码跑完把结果写进一个数组和 MATPOWER 的结果逐项做差看最大偏差落在了哪条母线、哪条支路上。# 自己代码的结果V_mine 是 newtown_pf 返回的复数电压向量 V_mine, Q_mine newton_pf(bus9, br9) # 从 MATPOWER res.bus 第8、9列抄回参考值 V_ref np.array([1.040, 1.025, 1.025, 1.026, 0.996, 1.013, 1.026, 1.016, 1.032]) ang_ref np.array([0.000, 9.280, 4.660, -2.220, -3.990, -3.690, 3.720, 0.730, -4.360]) err_v np.max(np.abs(np.abs(V_mine) - V_ref)) err_a np.max(np.abs(np.angle(V_mine, degTrue) - ang_ref)) print(最大电压幅值偏差(pu):, err_v) print(最大相角偏差(deg):, err_a)逻辑说明这段代码只是把 MATPOWER 的参考值手动抄进数组实际工作中更省事的做法是直接在 MATLAB 里把 res.bus 第 8、9 列导出成 CSV再用 Python 读进来避免手抄出错。对拍的标准是电压幅值偏差小于 1e-6 pu相角偏差小于 1e-4 度。如果偏差达到 1e-3 量级说明你的程序里还有模型细节没对上不能算通过。参数说明这里给的参考值是经典 case9 的典型收敛结果不同 MATPOWER 版本之间可能有最后一位小数的差异所以严格对拍时应以你自己机器上 runpf 的输出为准而不是拿我抄的这组数当标准。支路有功的对拍同样重要具体看 res.branch 表里的 PF、QF、PT、QT 四列分别对应首端和末端的有功、无功。4.3 6节点版本差异为什么你的结果和论文对不上6节点的交叉验证要更谨慎。如果你用 case6ww 跑完发现和某篇论文附录里的结果差了 0.1%先别急着怀疑程序和 MATPOWER先检查三件事一是负荷参数是否一致二是两台调压变压器的变比是否标注为 1.0三是线路充电电纳 B 是给了还是当成 0 处理。这三个变量任何一个不同都会让母线电压产生可以测量的偏差。行业里对 6 节点系统有一个共识它从来不是一个由 IEEE 标准化委员会钦定的唯一基准更多是靠教材和论文互相引用流传下来的。所以做算法对比时我习惯把 9 节点当锚点6 节点当辅助验证。如果你的论文需要引用算例建议同时给出 9 节点和 6 节点的结果并注明数据来源是 MATPOWER 自带版本还是教材附录省得审稿人拿另一套数据来质疑你。提示6节点数据里如果 tap 列不是 1.0但你用普通线路模型去算无功分布一定会出问题。这种偏差不会让程序报错只能靠对拍发现所以每次拿到新数据先打印一下支路表的前几行确认变比列的实际值。5. 潮流计算避坑6节点与9节点最常见的5个翻车现场5.1 节点编号从0开始导致Ybus行错位现象程序不报错但算出来的电压离谱比如某条母线电压高达 1.3 pu或者相角乱成一团。原因MATPOWER 和大多数论文数据里母线编号从 1 开始。你从数据文件读到程序里如果忘了减 1支路表里的1-4就被解释成 bus[1] 到 bus[4]实际应该指向 bus[0] 到 bus[3]。索引整体错一位Ybus 的拓扑关系全乱但矩阵本身还是方阵所以不会触发任何异常提示。解决数据定义阶段统一改成从 0 编号像我第 3 章那样手写数据时就写 0-based。如果是读外部文件构造支路表时强制减 1并在组装 Ybus 后打印对角线验证。# 排查技巧打印Ybus对角线检查母线i的自导纳是否远大于交叉导纳 Y build_ybus(bus9, br9) print(np.round(np.diag(Y).real, 4)) # 对角线应包含多条支路的贡献5.2 平衡节点放最后Jacobian奇异现象np.linalg.solve 报 Singular matrix直接中断。原因平衡节点的相角和电压都是固定值不参与迭代。组装 Jacobian 时如果忘了把 slack 对应行和列剔除矩阵会多出一个自由度变成奇异矩阵。特别是在把 slack 放在母线列表末位的数据里新手很容易因为没有删除最后一列而中招。解决用索引集合来管理状态量而不是靠位置。第 3 章代码里 pvpq 和 pq 都是从类型列里动态提取的slack 永远是 np.where(typ 2) 找出来的那个不管它在母线表里排第几都不会出问题。5.3 变压器支路忽略变比k无功怎么都对不上现象电压幅值和有功分布挺正常但无功功率和参考值偏差很大而且你调来调去都消不掉。原因6 节点系统里的调压变压器在不同数据版本里变比可能是 0.98 或者 1.05不是 1.0。你没处理 tapYbus 里就把变压器当成普通线路漏掉了非标变比带来的无功修正。9 节点恰好三台变压器都是 1.0所以这个坑在 9 节点上不暴露一换 6 节点就炸。解决build_ybus 里保留 tap 参数凡是 tap 不等于 1.0 的支路首端自导纳除以 tap 的平方互导纳除以 tap。注意 tap 的方向约定我一般默认变比位于首端母线侧如果数据文件用的是末端侧变比要取倒数再填进去。5.4 收敛判据设太粗PQ分解法结果不达标现象tol 设成 1e-4牛顿法正常收敛但对拍时发现支路有功差了 0.3 MW 左右怎么看都不像收敛到位。原因1e-4 的功率失配判据在标幺制下对应 0.01 MW 量级的偏差经过多条支路累加后线路潮流的误差会放大到可观测的水平。更重要的是快速解耦法因为忽略了电阻的影响本身存在模型近似误差判据如果再放宽偏差会进一步叠加。解决判据统一收紧到 1e-8 或者更小。牛顿法对这个规模的系统不差那几次迭代判据严一点结果才经得起和 MATPOWER 逐项对拍。如果用的是 PQ 分解法先和牛顿法结果对拍一次确认模型近似误差在可接受范围内再谈收敛判据。5.5 功率单位混用MW/MVar没除以基准值现象所有节点的计算结果整体放大了 100 倍发电机出力和线路潮流数值巨大但电压和相角看起来还算正常。原因系统阻抗和电压是标幺值但 Pg、Pd 这些功率还是 MW/MVar 的绝对值。把它们直接当成标幺功率塞进方程等于把整个系统的注入功率放大了基准值的倍数。解决进入迭代前统一除以基准功率。# 正确做法功率必须转标幺 base 100.0 P_pu (bus[:, 1] - bus[:, 3]) / base Q_pu (bus[:, 2] - bus[:, 4]) / base6. 进阶PQ分解法提速以及离线验证程序的三板斧如果你的算例规模从 9 节点涨到几百节点牛顿法每次迭代都要重新组装并分解 Jacobian计算量会明显上来。这时候可以换成快速解耦法也就是 PQ 分解法。它的核心假设是高压输电网络里支路电阻远小于电抗有功主要受相角影响无功主要受电压幅值影响于是把原来的四块 Jacobian 拆成两块独立的常数矩阵 B 和 B。B 对应有功-相角B 对应无功-电压两个矩阵离线组装一次迭代时反复用每次求解规模减半。# PQ分解法迭代核心B1、B2 是离线组装好的常数矩阵 # 迭代中只做两次小规模求解不再重新分解 for it in range(max_iter): V Vm * np.exp(1j * Va) S V * np.conj(Y V) dP (Pspec - S.real) / Vm # 有功失配除以电压幅值 dQ (Qspec - S.imag) / Vm # 无功失配除以电压幅值 Va[pvpq] np.linalg.solve(B1, dP[pvpq]) Vm[pq] np.linalg.solve(B2, dQ[pq])B 的组装一般只取支路电抗的倒数忽略电阻和充电电容B 则直接取 Ybus 的虚部再剔除 PV 母线和平衡母线的对应行列。换算法之后别忘了用同一套 9 节点数据做一次对拍PQ 分解法因为忽略了电阻结果和牛顿法会有微小差别知道差多少才不会误判程序故障。验证程序正确性我习惯用三板斧。第一板斧是标准算例对拍9 节点和 6 节点都必须过电压幅值偏差小于 1e-6 pu。第二板斧是功率平衡自检全网发电减去全网负荷再减去所有线路损耗应该约等于零这个检查能同时揪出支路功率计算和充电电容处理的问题。第三板斧是初值扰动测试把平启动的 1.0 改成 0.95 和 1.05重新跑一遍结果必须一致否则说明程序可能落到了一个错误的工作点。我自己早期吃过大亏拿 6 节点数据验证程序支路无功一直差 0.3%排查两天发现是变压器变比没处理。从那以后我养成了一个习惯任何潮流程序先跑 9 节点再跑 6 节点两边对上才敢说收敛。这个习惯一直留到现在希望帮到你。本文还有配套的精品资源点击获取