MATLAB多算法潮流计算包：从IEEE9到IEEE300的对比实验

简介：面向电气工程本硕博及科研人员，这份资源聚焦IEEE 9节点至IEEE 300节点标准系统的潮流计算，覆盖直流法、牛顿拉夫逊法、快速PQ分解法和Gauss-Seidel法四类经典算法，并提供对应MATLAB实现脚本与操作录像，适合电力系统分析课程设计、算法对比研究及入门编程练习。压缩包内共44个文件，以41个MATLAB的.m脚本为主，涵盖潮流计算主函数、各算法子程序、节点与支路数据文件、结果输出与绘图辅助模块等；另有1个位图说明文件、1个AVI操作录像和1个文本说明，方便逐行调试与对照学习。包体仅472KB，轻量紧凑，便于下载后快速部署。已有1025人学习使用，资源通过Runme_.m一键运行，配合录像可迅速掌握从数据导入到结果分析的全流程，是理解不同潮流算法收敛特性与适用场景的实用工具。

1. 多算法潮流计算包：从 IEEE9 到 IEEE300 的一次完整对照实验

做电力系统研究的人几乎都绕不开潮流计算，但真正把直流法、牛顿拉夫逊法、快速 PQ 分解法、Gauss-Seidel 法放在同一套 IEEE 标准算例体系下逐一对齐比较的资源并不多。这个 MATLAB 程序包的价值在于：它把 IEEE9、14、30、39、57、118、300 等不同规模算例全部纳入了同一套计算框架，每种算法都有独立入口函数，横向对比收敛速度、迭代次数和计算结果变得非常直接。尤其值得注意的一点是，当系统规模从 IEEE9 放大到 IEEE300 时，牛顿拉夫逊法和快速 PQ 分解法在迭代次数上的差距并不大，但单次迭代耗时差异显著，这个观察只有在一口气跑完所有算例之后才容易形成。适合正在做课程设计或准备论文仿真部分的本科生、研究生，也适合刚接触 MATPOWER 二次开发、想理解内部迭代细节的工程师。

程序不是零散脚本的堆砌，而是按 MATPOWER 风格组织的完整工程：case9.m到case300.m定义了各规模算例的母线、支路和发电机参数，Runme_.m作为总入口统一调度四种算法，hhu_前缀文件封装了针对不同算法的求解器。接下来先厘清四种算法的数学边界，再进入代码层面逐段拆解。

2. 直流法、牛顿拉夫逊与快速 PQ 分解的数学边界

2.1 四种算法在电力系统方程组中的定位差异

潮流计算本质上是在给定母线注入功率和网络导纳矩阵的前提下，求解节点电压幅值与相角的非线性方程组：

P_i = V_i * Σ V_j * (G_ij * cos(θ_ij) + B_ij * sin(θ_ij)) Q_i = V_i * Σ V_j * (G_ij * sin(θ_ij) - B_ij * cos(θ_ij))

四种算法的差异不在于目标函数，而在于对雅可比矩阵的处理方式。直流法把问题彻底线性化：忽略无功、忽略电阻、假设所有电压幅值为 1.0 pu、相角差足够小使得 cos(θ) ≈ 1、sin(θ) ≈ θ，于是只剩下 P = B' * θ 这个线性方程；牛顿拉夫逊法保留完整雅可比矩阵并迭代求解修正量；快速 PQ 分解法利用输电网络高 X/R 比的特点，把有功和无功解耦，并使用恒定、对称的 B' 和 B'' 矩阵替代每次迭代都重新计算的雅可比；Gauss-Seidel 法则完全避开矩阵求逆，逐母线迭代更新电压。

从理论上看，直流法的适用范围是纯有功潮流分析，比如经济调度和输电能力评估；牛顿法是精度基准；PQ 分解法是牛顿法在内存和速度上的折中；Gauss-Seidel 更适合教学演示而非大规模系统。这个包把四种方法放在统一框架内跑同一批算例，实测结果的差异会直接印证理论边界。

2.2 四种算法的迭代机制对照

实践中的真实差异在 IEEE300 上体现得最充分：牛顿法通常 3~5 次迭代即可收敛到 1e-8 的功率偏差，PQ 分解法需要 8~15 次但每次迭代只做前代回代，整体时间反而可能更短；Gauss-Seidel 在 300 节点上经常要几百次迭代，且对初始电压很敏感。

3. Runme_.m 入口与 hhu_ 前缀封装：读懂这套代码的工程结构

3.1 总体工程布局

解压后第一件事不是打开某个求解函数，而是先看Runme_.m。这个入口脚本和hhu_enter.m、hhu_again.m、hhu_agains.m构成了整个程序的三级调度关系。Runme_.m负责加载算例数据、选择算法、调用对应的hhu_xxx函数并输出结果；hhu_enter.m是核心求解入口，按用户指定的算法进入不同分支；hhu_again.m和hhu_agains.m是断点续算或参数调整后的重跑入口。

Runme_.m % 唯一的用户入口 │ ├── hhu_enter.m % 算法分发中心 │ ├── hhu_newtonpf.m % 牛顿拉夫逊法实现 │ ├── hhu_pqpf.m % 快速PQ分解法实现 │ ├── hhu_dcpf.m % 直流法实现 │ ├── hhu_gspf.m % Gauss-Seidel法实现 │ └── hhu_runpf.m % 综合调度与结果输出 │ ├── case9.m / case14.m / case30.m / case39.m ├── case57.m / case118.m / case300.m % 各规模算例 │ ├── makeYbus.m % 构建节点导纳矩阵 ├── makeBdc.m % 构建直流法B'矩阵 ├── makeSbus.m % 构建注入功率向量 ├── dSbus_dV.m % 计算功率对电压的偏导 ├── dSbr_dV.m % 计算支路功率对电压的偏导 └── bustypes.m / pfsoln.m / printpf.m % 节点类型识别、结果回代、打印

hhu_runpf.m是贯穿所有算法的主调度，它调用了makeYbus.m构建导纳矩阵，调用bustypes.m识别 PQ、PV、平衡节点，再根据算法名分发给newtonpf.m、fdpf.m、dcpf.m、gausspf.m或它们对应的hhu_版本。

3.2 运行入口的配置逻辑

Runme_.m中常见的配置块如下：

%% 选择算例与算法 mpc = loadcase('case300.m'); % 载入IEEE300节点系统 mpopt = mpoption('PF_ALG', 1, 'VERBOSE', 2, 'OUT_ALL', 1); % PF_ALG: 1-牛顿法, 2-快速PQ分解, 3-直流法, 4-Gauss-Seidel results = runpf(mpc, mpopt);

参数含义拆开来说：PF_ALG=1走牛顿法、2走fdpf、3走dcpf，而4对应gausspf；VERBOSE=2控制输出详细程度，设成0可以关闭所有中间打印；OUT_ALL=1表示输出完整的results结构体，包括bus、branch、gen字段。值得留意的是，文件列表里同时存在runpf.m、newtonpf.m、fdpf.m、dcpf.m和hhu_runpf.m、hhu_newtonpf.m等成对文件，形式上看像是 MATPOWER 原生函数与本封装函数并存。实际运行时Runme_.m通过hhu_前缀版本走完整个流程，原因是原生函数对输入参数格式的要求更严苛，hhu_版本放宽了部分检查、增加了断点续算支持。想理解算法本质，优先看hhu_版本；想核对标准结果，对比原生版本。

4. hhu_newtonpf 核心迭代：Jacobian 组装、稀疏求解与收敛判定

4.1 从功率偏差到修正方程的推导

牛顿拉夫逊法的每一步迭代核心是求解：

[ H N ] [Δθ] [ΔP] [ M L ] * [ΔV] = [ΔQ]

其中H = ∂P/∂θ、N = ∂P/∂V、M = ∂Q/∂θ、L = ∂Q/∂V。在hhu_newtonpf.m中，这个分块矩阵由dSbus_dV.m计算，它输出的是复功率对电压相量和幅值的偏导。代码实现的关键在于把复导纳矩阵 Ybus 的实虚部拆开利用：

function [dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V) % 计算功率偏差对电压幅值和相角的偏导 Ibus = Ybus * V; % 节点注入电流 diagV = sparse(1:length(V), 1:length(V), V); diagIbus = sparse(1:length(Ibus), 1:length(Ibus), Ibus); % 对相角求偏导: dS/dVa = j * diag(V) * conj(diag(Ibus) - Ybus*diag(V)) dSbus_dVa = 1j * diagV * conj(diagIbus - Ybus * diagV); % 对幅值求偏导: dS/dVm = diag(V) * conj(Ybus * diag(V)) + conj(diagIbus) * diagV dSbus_dVm = diagV * conj(Ybus * diagV) + conj(diagIbus) * diagV; end

这段代码的逻辑是：先从Ybus * V算出节点注入电流，再构造对角矩阵diagV和diagIbus，最后按复功率对电压相量求偏导的链式法则，把结果拆成对相角和对幅值两个偏导矩阵。注意这里用了sparse稀疏矩阵存储，因为 IEEE300 的 Ybus 稠密度不到 1%，用全矩阵会导致内存膨胀数倍。dSbus_dVa对应雅可比矩阵的 H 和 M 分块，dSbus_dVm对应 N 和 L 分块。

4.2 迭代循环与收敛判据

完整迭代从pfsoln.m的输出回溯可见：先初始化 V 为平启动值（幅值 1.0、相角 0），随后进入迭代，每轮先算功率偏差并判断是否小于容差，不满足则组装雅可比求修正量。

%% 牛顿法主迭代循环 V = V0; % V0为平启动初始电压 tol = 1e-8; % 收敛容差 iter = 0; max_iter = 30; while iter < max_iter %% 计算功率偏差 [Pcal, Qcal] = hhu_power_balance(V, Ybus); % 按当前V计算注入功率 dP = Pspec - Pcal; % 有功偏差向量 dQ = Qspec - Qcal; % 无功偏差向量 %% 收缩到自由节点 dVa = dP(pv); % PV和PQ节点的有功偏差 dVm = dQ(pq); % 仅PQ节点的无功偏差 if max(abs([dVa; dVm])) < tol break; % 满足收敛条件 end %% 组装雅可比并求解修正方程 [dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V); J = [ real(dSbus_dVa(:, pv)) real(dSbus_dVm(:, pq)); imag(dSbus_dVa(:, pv)) imag(dSbus_dVm(:, pq)) ]; dVa = J \ [dVa; dVm]; % 稀疏LU分解求解 end

J \ [dVa; dVm]是 MATLAB 内置的稀疏线性求解，内部自动选择 LU 分解策略。pv和pq索引向量来自bustypes.m，它把母线分为三类：平衡母线只给电压初值、PV 母线给定有功和电压幅值、PQ 母线给定有功和无功。实际运行时不同版本会直接在文档注释或hhu_runpf.m中给出鼓励的推荐入口——优先走hhu_enter.m进入不同算法的总入口，不要跳过调度层直接调用迭代函数。

注意：直接运行newtonpf.m或fdpf.m这类被调用文件会报“未定义变量”错误，因为它们的输入参数（Ybus、V0、ref、pv、pq等）由上层脚本传入。遇到报错先检查当前文件夹是否在工程根目录。

5. fdpf 与 gausspf 的实现差异：近似的代价和 Gauss-Seidel 的收敛短板

5.1 快速 PQ 分解法的 B' 与 B'' 矩阵构造

makeBdc.m负责构建直流法的 B' 矩阵，而fdpf.m中 PQ 分解法的 B' 和 B'' 矩阵由它派生。B' 取导纳矩阵虚部的有功功率相关部分，B'' 取无功功率部分。差异在于：B' 不含并联支路和变压器非标准变比的影响，B'' 则计入这些分量。makeBdc.m的关键代码逻辑如下：

function [Bdc, Bdcf] = makeBdc(Ybus) % 直流法B'矩阵构造: 取导纳矩阵虚部并排除平衡节点 B = imag(Ybus); % 导纳矩阵虚部就是电纳部分 Bdc = B; % 先复制全矩阵 Bdc(:, ref) = []; % 删除平衡节点对应列 Bdc(ref, :) = []; % 删除平衡节点对应行 Bdc = -Bdc; % 取负号得到 B' 矩阵 end

这里容易忽略的是符号处理：潮流方程中 P = -B' * θ，所以代码里做了取负。同时imag(Ybus)直接取 Ybus 虚部而不是对每支路单独处理，这是直流法和快速 PQ 分解法在工程实现上的共同捷径。hhu_pqpf.m交替求解两个低阶线性方程组：先固定 V 求 Δθ，再固定 θ 求 ΔV，交替迭代直到收敛。由于 B' 和 B'' 在迭代中保持不变，只需做一次 LU 分解，后续迭代都是前代回代，速度优势在这里体现。

%% PQ分解法交替迭代 [Bp, Bpp] = hhu_bprime(Ybus, ref, pv, pq); % 常数矩阵，只算一次 [Lp, Up] = lu(Bp); % 对B'做LU分解 [Lpp, Upp] = lu(Bpp); % 对B''做LU分解 while iter < max_iter dP = Pspec - Pcal; % 有功偏差 dVa = U \ (L \ dP(pv)); % 回代求解相角修正 Va(pv) = Va(pv) - dVa; % 更新相角 dQ = Qspec - Qcal; % 无功偏差 dVm = Upp \ (Lpp \ dQ(pq)); % 回代求解幅值修正 Vm(pq) = Vm(pq) - dVm; % 更新幅值 end

与牛顿法每次迭代重新组装雅可比矩阵相比，hhu_pqpf.m只在进入迭代前做两次lu分解，循环内仅执行U \ (L \ x)两次前代回代。代码注释里通常写的是这是工程常用套路：在较高 R/X 比的输电网中，P 主要受 θ 影响、Q 主要受 V 影响，因此解耦造成的误差可接受。

5.2 Gauss-Seidel 的电压更新机制与发散风险

gausspf.m走的完全是另一条路线：

function V = gausspf(Ybus, Sbus, V0, ref, pv, pq, max_iter, tol) V = V0; for iter = 1:max_iter Vprev = V; % 保存上一次迭代结果 for i = 1:n % 逐母线更新 % 从第2个母线开始，计算第i个母线的注入功率 sum_j = Ybus(i,:) * V; % 全网络电压对i的贡献 V(i) = (conj(Sbus(i)/V(i)) - sum_j + Ybus(i,i)*V(i)) / Ybus(i,i); end % PV节点电压幅值修正回设定值 V(pv) = abs(V0(pv)) .* (V(pv) ./ abs(V(pv))); if max(abs(abs(V) - abs(Vprev))) < tol break; end end end

最大软肋有三处：一是逐母线串行更新，没有矩阵层面的并行加速；二是 PV 节点处理靠迭代结束后强制拉回电压幅值这一后处理操作完成，但每次强制修正都会引入新的无功不匹配，导致收敛曲线振荡；三是对重负荷系统或病态网络，Gauss-Seidel 经常出现相邻两次迭代电压差值无法下降到容差以下的现象。实际跑 IEEE118 和 IEEE300 时，gausspf在重负荷工况下经常需要 500 次以上迭代才能收敛，甚至比牛顿法慢一个数量级。这是正常的，不是程序 bug。

6. 从 IEEE9 到 IEEE300 的准入参数：负荷、电压初值与收敛行为映射

6.1 各算例的规模与运行参数

代码包中case9.m到case300.m的算例数据组织方式完全一致，都包含bus、branch、gen三个基础矩阵。bus矩阵前四列分别是母线编号、类型（1=PQ，2=PV，3=平衡）、有功负荷、无功负荷；branch矩阵前四列是从端、至端、电阻、电抗；gen矩阵包含发电机母线编号、有功出力、电压幅值设定值等。既然所有算例共享同一套数据结构，准入参数的核心就是管理好电压初值和负荷缩放方式。

6.2 算例规模与算法适用性的实践映射

这些数字是一般趋势，实际取值受容差设置和负荷水平影响会有所浮动。对 NR 与 PQ 分解法在部分国产教材中被称为牛拉法与快速解耦法，实际使用时在hhu_pqpf.m中可看到两个矩阵的构建逻辑：B' 剔除平衡节点行，B'' 同时剔除平衡和 PV 节点行，这是造成两者矩阵维度不同的原因，稳定推论。

6.3 参数调整与避坑技巧

运行Runme_.m之前，最常遇到的两个改动需求一是换算例规模，二是改收敛精度：

%% 在 Runme_.m 中切换算例 mpc = loadcase('case118.m'); % 把 case300 换成 case118 即可 %% 调整收敛容差 mpopt = mpoption('PF_TOL', 1e-10); % 默认1e-8，调小提高精度

换成大算例时，makeYbus.m自动构建对应维度的 Ybus，无需优化代码。但有个前提：MATLAB 当前工作目录必须切换到工程根目录，否则loadcase会因找不到case300.m而报错。操作录像中演示的就是这个环节，录像文件操作录像0023.avi中有完整的文件夹切换和运行演示。

三个高频报错按经验排序：

• 直接点运行了newtonpf.m等子函数，弹出“未定义函数或变量 Ybus”——回到Runme_.m运行；
• 当前文件夹不对导致loadcase失败——用cd切到工程目录确认case300.m文件可见；
• 车载负荷等设备需在某些父层目录下找资源，版本差异导致mpoption语法报错——MATLAB 2021a 以下版本没有PF_ALG选项，用mpoption('PF_ALG', 2)会直接报错，需要改用runpf(mpc, 'dc')或升级到 2021a 以上版本。

6.4 一键跑完多算例的对比脚本

换算例与算法的手动操作本质上完全可自动化。推荐一个比手动换参数更有效率的做法：写一段脚本循环跑四种算法在所有算例上的结果，集中比较迭代次数和耗时：

%% 批量对比脚本 algs = {'NR', 'PQ', 'DC', 'GS'}; % 四种算法标签 cases = {'case9.m', 'case14.m', 'case30.m', ... 'case39.m', 'case57.m', 'case118.m', 'case300.m'}; result_table = cell(length(cases), length(algs)); for i = 1:length(cases) mpc = loadcase(cases{i}); for j = 1:length(algs) mpopt = mpoption('PF_ALG', j, 'VERBOSE', 0, 'OUT_ALL', 0); tic; r = runpf(mpc, mpopt); elapsed = toc; result_table{i, j} = [r.iterations, elapsed]; end end

r.iterations是 MATPOWER 求解器输出的迭代次数，toc得到的是包括矩阵组装和求解在内的完整耗时。跑完case300.m后你会发现：直流法的计算结果和牛顿法在电压幅值上差异低于 1% 左右，但相角在重负荷母线上可能差到 3~5 度。这正好为“输电能力评估用直流法、详细运行方式分析用牛顿法”提供了一个直观的数据证据。

6.5 收敛效果的最后验证方法

判断一个算法是否真正收敛到物理上有意义的解位，最直接的方法是验证功率平衡方程是否满足：

%% 结果对比验证 V = results.bus(:, 8) .* exp(1j * results.bus(:, 9) * pi / 180); Sbus = makeSbus(results.baseMVA, results.bus, results.gen); Scalc = V .* conj(results.bus(:, 1) / results.baseMVA); % 对比注入功率 max_residual = max(abs(Sbus - V .* conj(results.bus(:, 1) / 100)));

其中makeSbus.m把bus和gen中的负荷与发电数据合并成净注入功率向量，residual的计算检查当前电压下功率方程两侧的最大偏差。如果不同算法在同一算例上的 max_residual 都是 1e-8 量级，说明它们收敛到了同一个解（或至少是数值上不可区分的解）；如果某算法 residual 很大但迭代次数显示已收敛，说明该算法内部使用的收敛判据和本地评估不一致。

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