简介：面向电力系统研究人员与学生的MATPOWER实用代码包，聚焦IEEE300节点系统在直角坐标系下的牛顿拉夫逊法潮流计算；MATPOWER作为MATLAB电力系统分析工具箱，在潮流计算、稳定性研究与优化问题求解中应用广泛，特别适合处理大规模网络拓扑。资源仅含一个M文件，压缩包大小仅2KB，结构精简，便于直接阅读、运行和二次修改。该文件基于MATPOWER框架搭建完整的IEEE300节点模型，模型包含发电机、负荷、变压器等多种设备，并以直角坐标形式表示各节点电压，通过牛顿拉夫逊迭代逐步逼近非线性功率平衡方程的解，每一步计算并更新雅可比矩阵，最终得到各节点电压幅值、相角及支路潮流结果。直角坐标形式直观易读，适合初学者理解牛顿拉夫逊算法的数学推导与编程实现；三百节点规模则能反映大型实际电网的拓扑复杂性和计算需求，可作为标准算例用于教学演示、算法对比、科研验证以及后续的极坐标扩展、最优潮流等研究基础。已有六百零五人学习下载，对正在学习电力系统潮流计算、复现标准测试系统或希望基于该案例开展更深层研究的读者，都具有直接的参考价值。
1. 拆开ieee300.rar：直角坐标牛顿法在MATPOWER里的正确打开方式
解压ieee300.rar，里面只有一个ieee300.m文件，连case300.m都没带。第一次拿到的人可能觉得不值，但这份脚本的价值不在数据，而在它展示了一条清晰的实现路径：如何利用MATPOWER的IEEE300节点系统数据，在直角坐标系下从零写出牛顿-拉夫逊（NR）潮流求解器。MATPOWER自带的runpf在极坐标下求解，收敛快、代码稳，但内部细节被封装得严严实实。而直角坐标版NR的雅可比矩阵结构、PV节点的约束方式、稀疏填充技巧都与极坐标有明显差异，这些恰恰是研究最优潮流、状态估计、含FACTS设备潮流调整时绕不开的底层能力。这篇文章从ieee300.m出发，适合已经能跑通MATPOWER标准算例、想深入理解NR内核的人，也适合需要把潮流算法改造成特殊模型的研究生和工程师。
2. IEEE300数据模型与MATPOWER的加载边界
2.1 case300.m和ieee300.m的关系
MATPOWER自带的case300.m是IEEE300节点系统的数据文件，定义了母线、支路、发电机、负荷以及变压器参数。而ieee300.m通常是一个脚本，它负责把这份数据读入内存，然后调用自定义的直角坐标NR代码来求解潮流。很多下载包里没有附带case300.m，因为MATPOWER安装目录的data文件夹下已经有了，你只需要保证MATLAB搜索路径包含该目录即可。
我一般会在ieee300.m开头加上两行：
define_constants; mpc = loadcase('case300');define_constants会导入BUS_TYPE、PD、QD、VM、VA等列索引常量，避免用魔法数字。loadcase返回结构体mpc，字段包含bus、branch、gen、baseMVA等。注意这里不要用MATPOWER老版本习惯的load('case300')，因为新版数据格式可能不兼容，loadcase会自动进行格式转换和校验，这是官方推荐的入口。
2.2 bus、branch、gen矩阵中与NR直接相关的列
直角坐标NR需要从mpc中提取三类信息：网络拓扑与参数、节点注入功率、发电机约束。以mpc.bus为例，每一行代表一个节点，常用列的用途如下：
| 列索引 | 字段名 | 在直角坐标NR中的作用 |
|---|---|---|
| BUS_I | 节点编号 | 建立内部连续编号映射 |
| BUS_TYPE | 节点类型 | 1=PQ，2=PV，3=REF，决定未知量组合 |
| PD、QD | 有功/无功负荷 | 计算节点注入功率指定值 |
| GS、BS | 对地电导/电纳 | 叠加到导纳矩阵对角元 |
| VM、VA | 电压幅值/相角初值 | 计算初始e、f |
| BASE_KV | 基准电压 | 标幺值转换辅助 |
branch矩阵的F_BUS、T_BUS、BR_R、BR_X、BR_B、TAP等字段会传给makeYbus，由它生成复导纳矩阵Ybus。gen矩阵中GEN_BUS、PG、QG、QMAX、QMIN、VG是NR迭代中的关键：PG和QG决定注入功率，QMAX/QMIN用于无功越限后PV/PQ节点类型转换，VG则是PV节点电压幅值设定值。
在写NR之前，我会先做一个快速检查，确认数据维度正确：
fprintf('bus: %dx%d, branch: %dx%d, gen: %dx%d\n', ... size(mpc.bus,1), size(mpc.bus,2), ... size(mpc.branch,1), size(mpc.branch,2), ... size(mpc.gen,1), size(mpc.gen,2));这一步能提前发现数据文件被截断或字段错位的问题，尤其从第三方仓库下载的case文件经常缺少列数据。
2.3 makeYbus与自定义求解器的接口
MATPOWER把网络拓扑到导纳矩阵的转换封装在makeYbus中，这是自定义NR可以直接复用的部分。调用方式如下：
Ybus = makeYbus(mpc.baseMVA, mpc.bus, mpc.branch);mpc.baseMVA是系统基准容量，典型的IEEE300系统取100MVA。makeYbus返回一个稀疏复矩阵，维度为节点数×节点数。它已经考虑了线路充电电容、变压器变比、对地导纳等，不需要你再手动修正。拿到Ybus之后，就可以用find提取非零元素的坐标和数值，为组装雅可比矩阵做准备。
至于runpf，它和自定义NR的分工不要混淆。runpf是一套完整的极坐标求解器，直接调用它只会得到结果，无法看到中间过程的雅可比矩阵。而ieee300.m中实现的是“自己控制每一步”的求解流程：数据加载用MATPOWER，Ybus用MATPOWER，但迭代核心是自己写的。这样既保留了数据模型的可靠性，又能把算法内核握在手里，方便加自定义输出或改造成其他迭代格式。
3. 直角坐标下的NR实现：功率方程、雅可比矩阵与迭代控制
3.1 直角坐标功率方程与不平衡量构造
设节点i的电压实部和虚部分别为e_i、f_i，即V_i = e_i + j f_i。注入功率P_i、Q_i的表达式为：
P_i = e_i * Σ_j (G_ij * e_j - B_ij * f_j) + f_i * Σ_j (G_ij * f_j + B_ij * e_j)
Q_i = f_i * Σ_j (G_ij * e_j - B_ij * f_j) - e_i * Σ_j (G_ij * f_j + B_ij * e_j)
定义有功不平衡量ΔP_i为指定注入有功减去当前计算有功，ΔQ_i同理。对于PQ节点，未知量是e_i和f_i，方程是ΔP_i = 0和ΔQ_i = 0。对于PV节点，电压幅值给定，因此ΔQ_i方程替换为幅值约束方程：
e_i² + f_i² = V_i_spec²
在ieee300.m中，第一步通常是建立节点索引分类：
pq = find(mpc.bus(:, BUS_TYPE) == PQ); pv = find(mpc.bus(:, BUS_TYPE) == PV); ref = find(mpc.bus(:, BUS_TYPE) == REF);注意，ref节点在直角坐标NR中不参与迭代，其e、f固定不变。pq和pv列表的长度决定了未知量的总数：PQ节点有两个未知数，PV节点也有两个未知数，但PV节点少一个无功方程，多一个幅值方程，所以总方程数仍然等于总未知数。
构造节点注入功率时，不能直接用mpc.bus中的PD和QD，因为发电机输出功率在mpc.gen中。常规做法是先分配发电机功率到对应节点上：
nb = size(mpc.bus, 1); Pgen = zeros(nb, 1); Qgen = zeros(nb, 1); Pgen(mpc.gen(:, GEN_BUS)) = mpc.gen(:, PG); Qgen(mpc.gen(:, GEN_BUS)) = mpc.gen(:, QG); P_spec = Pgen - mpc.bus(:, PD); Q_spec = Qgen - mpc.bus(:, QD);这里P_spec和Q_spec是给定注入功率，后续在迭代中与计算值求差。注意如果存在多个发电机连接到同一节点，mpc.gen(:, GEN_BUS)的赋值会覆盖前面发电机的值，此时应该改用accumarray进行累加。IEEE300系统每个节点最多挂一台发电机，覆盖写法足够，但换到IEEE57等系统时就要小心多个发电机并联的情况。
3.2 雅可比矩阵的解析组装与稀疏化
直角坐标NR的修正方程是：
J * Δx = - [ΔP; ΔQ]
J矩阵是分块稀疏结构。每个非对角块对应一对节点之间的偏导数，对角块额外叠加本节点相关项。以PQ节点i和节点j为例，对e_j求偏导时，需要区分j是否等于i。
直接投公式容易漏项，更稳妥的方式是用数值差分验证一次解析雅可比。但为了在300节点系统上获得可接受的性能，最终还是要用稀疏矩阵组装。这里给出一个常见做法：先提取Ybus的非零元素坐标，用accumarray快速生成J。
Ybus = makeYbus(mpc.baseMVA, mpc.bus, mpc.branch); G = real(Ybus); B = imag(Ybus); [Yi, Yj, ~] = find(Ybus);之后在一个循环里处理所有非零元素，计算四个2x2分块并填充到稀疏矩阵。需要对e和f分别求导。为了减少重复计算，可以预先计算Ge = G * e、Bf = B * f等向量。
关键的一点：J必须使用sparse函数创建，例如：
J_pf = sparse(nb, nb); J_pq = sparse(nb, nb); J_qf = sparse(nb, nb); J_qe = sparse(nb, nb);这4个子矩阵分别对应ΔP对e、ΔP对f、ΔQ对e、ΔQ对f的偏导。最终J按节点顺序装配成2nb × 2nb的大稀疏矩阵。
如果不想在稀疏组装的细节上花太多时间，也可以先写出稠密雅可比验证正确性，再改成稀疏版本。IEEE300节点的稠密雅可比是600×600，不算太大，但每次迭代重建一次稠密矩阵在速度上明显慢于稀疏版本，所以作为第二步优化很有必要。
3.3 迭代循环、收敛判据与初值处理
完整的直角坐标NR迭代骨架如下：
V = mpc.bus(:, VM) .* exp(1j * mpc.bus(:, VA) * pi / 180); e = real(V); f = imag(V); tol = 1e-8; max_iter = 20; for iter = 1:max_iter V = e + 1j * f; S_calc = V .* conj(Ybus * V); P_calc = real(S_calc); Q_calc = imag(S_calc); dP = P_spec - P_calc; dQ = Q_spec - Q_calc; % 对PV节点，用电压幅值方程代替无功方程 dQ(pv) = - (e(pv).^2 + f(pv).^2 - mpc.bus(pv, VM).^2); % 组装J J = assemble_jacobian(e, f, G, B, pq, pv); % 约束向量 dS = [dP(pq); dP(pv); dQ(pq)]; % 注意顺序要对应 dx = -J \ dS; % 更新e、f idx = 1:length(pq); e(pq) = e(pq) + dx(idx); f(pq) = f(pq) + dx(idx+length(pq)); idx = idx + 2*length(pq); e(pv) = e(pv) + dx(idx); f(pv) = f(pv) + dx(idx+length(pv)); if max(abs(dP)) < tol && max(abs(dQ)) < tol break; end end这里的dS排列顺序必须和雅可比矩阵的行列排列一致。一种常见约定是：先所有PQ节点的ΔP，再所有PV节点的ΔP，再所有PQ节点的ΔQ，最后是PV节点的电压幅值偏差。注意幅值方程的右边是“给定值减当前平方和”，所以dQ(pv)写成负号形式。
迭代初值方面，MATPOWER默认使用平启动，即PQ节点的电压幅值为1、相角为0。直角坐标NR对初值不敏感，但对于电压等级跨度大的系统，如果一开始就用极坐标下得到的初值，收敛速度会明显更快。通常平启动已经能让IEEE300在8次左右收敛到1e-8。
3.4 直角坐标与极坐标NR的取舍
MATPOWER的runpf默认使用极坐标NR，变量是电压幅值和相角。极坐标的好处是所有PV节点的无功方程天然被幅值约束替代，方程规模略小，且雅可比矩阵对角占优更明显，收敛性通常更好。但直角坐标NR也有不可替代的使用场景：处理含移相器、储能变流器、FACTS装置时，这些元件的功率注入方程在直角坐标下呈现多项式形式，求偏导更容易；另外在最优潮流和状态估计中，直角坐标下的等式约束是二次函数，更容易被内点法和序列二次规划利用。ieee300.m展示的正是这条“非默认但重要”的路线。
4. IEEE300系统上的收敛调试与大型系统扩展
4.1 用IEEE300观察NR的典型收敛曲线
把迭代过程中的最大不平衡量记录下来，画成曲线是调试的第一步。在ieee300.m中加入：
mismatch_history(iter) = max([norm(dP, inf), norm(dQ, inf)]);IEEE300在正确的实现下，前两次迭代不平衡量下降两个数量级，之后每轮下降约一个数量级，最终收敛到1e-10以内。如果曲线出现振荡、收敛停滞，常见原因有三种：雅可比矩阵某一行符号反了；PV节点幅值方程的dQ符号写反；更新解时把PQ和PV节点的未知量顺序搞混。
还有一种不太明显的错误：节点编号不是从1开始的连续整数时，直接使用mpc.bus(:, BUS_I)作为索引会导致Ybus维度错位。makeYbus默认按内部连续编号构建，你需要用MAX_BUS或自动重编号方式建立映射。IEEE300的节点编号恰好是1～300连续，所以不会有问题，但如果把脚本用到IEEE118（节点编号不连续）时，就必须先做重编号。
4.2 PV/QV节点转换与无功越限处理
大型系统中发电机无功越限是导致NR发散的最常见原因。迭代过程中如果发电机的无功功率超出QMAX或QMIN，该节点必须从PV节点变成PQ节点，无功用界值固定。MATPOWER的runpf内部有这一处理，而自己的实现里很容易漏掉。
在每轮迭代后，估算发电机无功：
Qgen_est = Q_calc + mpc.bus(:, QD);然后对每个发电机节点，检查是否越限：
genbus = mpc.gen(:, GEN_BUS); Qmax = mpc.gen(:, QMAX); Qmin = mpc.gen(:, QMIN); violated = (Qgen_est(genbus) > Qmax) | (Qgen_est(genbus) < Qmin);如果存在越限节点，就把对应节点的BUS_TYPE从PV改成PQ，同时更新该节点指定的无功值：
mpc.bus(genbus(violated), BUS_TYPE) = PQ; mpc.bus(genbus(violated), QD) = clamp(Qgen_est(genbus(violated)));同时要更新P_spec和Q_spec以及pq、pv索引列表，重新组装雅可比矩阵。这个处理相当于在迭代中动态改变方程结构，实现起来要特别小心，不能只在原来的矩阵上打补丁。
4.3 稀疏LU分解与混合迭代策略
对于300节点系统，直接用J \ dx已经足够快。如果想把算法扩展到数千节点级，需要关注两个性能点：第一，雅可比矩阵必须保持稀疏存储，不要用full()转换；第二，每轮迭代的线性求解开销很大，可以用lu分解并复用分解结果。
[L, U, P, Q] = lu(J); dx = Q * (U \ (L \ (P * dS)));注意这里P和Q是置换矩阵，如果MATLAB版本较新，还可以用decomposition(J, 'lu')获得一个可复用的分解对象。不过当节点类型发生转换时，J矩阵已改变，分解必须重新做。
初值优化也是提升收敛性的常用手段。平启动对IEEE300完全可行，但很不适合重负荷系统。一种混合策略是先用高斯-赛德尔法迭代5次得到一个近似解，再作为NR的初值。高斯-赛德尔法对初值不敏感，但收敛慢，正好适合“拉”一把。实现时只需要额外写一个简单的GS迭代器，在NR主循环前调用即可。
5. 把ieee300.m改造成可复用模块的验证与排错技巧
5.1 封装成nr_cartesian函数
不要满足于一份只能算IEEE300的脚本。把核心迭代逻辑封装成函数，输入是MATPOWER数据结构和控制参数，输出是电压向量、收敛标志和迭代次数，这样就能直接复用到IEEE118、IEEE230等标准算例。
function [V, success, iter] = nr_cartesian(mpc, tol, max_iter)函数内部先调用makeYbus，再执行迭代。调用示例：
mpc = loadcase('case300'); [V, ok, it] = nr_cartesian(mpc, 1e-10, 30);封装时要注意把define_constants放在函数外部或使用mpc.bus(:)列索引时避免全局常量依赖。我倾向于在函数开头调用define_constants，但不影响调用方。
5.2 与runpf极坐标结果进行误差对照
验证正确性的最直接方法，是和MATPOWER的runpf结果比较。runpf返回的结果数据也存放在bus结构中：
mpc2 = runpf(mpc); err_Vm = max(abs(abs(V) - mpc2.bus(:, VM))); err_Va = max(abs(angle(V) - mpc2.bus(:, VA) * pi / 180)); fprintf('Vm误差: %e, Va误差: %e\n', err_Vm, err_Va);对于IEEE300，误差在1e-6以下说明雅可比矩阵和迭代逻辑基本正确。有一点需要注意：极坐标下PV节点电压幅值严格取设定值，而直角坐标下由于收敛精度限制可能在小数点后8位有偏差，所以误差检查使用最大绝对值比均方根更合适。
5.3 功率不平衡量自检方法
即使没有runpf可对照，也可以用功率平衡方程自检。收敛后重新计算注入功率：
S = V .* conj(Ybus * V); P_res = P_spec - real(S); Q_res = Q_spec - imag(S);打印最大有功和无功不平衡量，如果都在1e-8量级，说明解正确。对于PV节点，还要额外检查幅值误差：
err_V = max(abs(abs(V(pv)) - mpc.bus(pv, VM)));这个小技巧在新接入自定义设备模型时特别有用，因为修改后的雅可比矩阵往往有细节错误，功率不平衡量能快速暴露问题所在。
5.4 别让“NR”搜索指向5G新空口
在电力系统语境中，NR是Newton-Raphson的简写，这一点在搜索资料时很容易踩坑。当你搜索“NR 流程”“NR 参数”时，结果会被5G新空口（New Radio）的大量内容淹没，包括PSS/SSS同步信号、载波聚合流程等。这些和潮流计算完全不相关。建议在ieee300.m的代码注释里统一写成“Newton-Raphson”，变量名用nr而不是NR，可以减少自己和协作者的混淆。当你需要查找资料时，用“Newton-Raphson power flow”或“直角坐标牛顿法”作为关键词，能直接命中电力系统技术文档。
本文还有配套的精品资源，点击获取
