简介：《Matlab牛拉法计算潮流》是一份面向电力系统专业学生与工程技术人员的牛拉法潮流计算源码包，利用牛顿-拉弗森迭代法求解电力网络稳态潮流，解决节点电压、功率分布及线路潮流的计算问题。压缩包共13个文件，含3个m格式源代码、9个txt格式数据与结果文件，以及1个docx说明文档，整体大小约529KB，按代码、输入数据、输出结果和文档分类，层次清晰。目前已有2077人学习下载。三个m程序对应不同电力网络配置的实例，完整覆盖模型建立、初始猜测、Jacobian矩阵构造、迭代更新与收敛判定等核心环节；txt文件提供线路导纳、节点负荷等输入参数及各节点电压、支路功率等结果数据，docx文档补充了使用说明与理论背景。读者可对照运行、逐步拆解，既能深入掌握牛拉法原理，也能提升Matlab编程与电力系统分析实操能力，适合作为课堂配套或自学参考资料。
1. 牛拉法潮流计算：拿到zip包前先想清楚的事
一份命名规整的“Matlab牛拉法计算潮流.zip”解压之后，里面通常是一个或几个.m文件、一张节点支路数据表，运气好还有一份IEEE标准算例。但真正要解决的不是“跑出一个P和Q”，而是理解牛顿-拉夫逊法在潮流方程里到底在迭代什么：把非线性功率方程在初值处泰勒展开，保留一阶项，反复求解线性修正方程，直到节点功率不平衡量小于阈值。这个思路和电气工程本科教材里的经典算法一致，但落地到Matlab代码时，雅可比矩阵的排列、PV节点无功越限处理、收敛判据的选取都直接影响结果。
这篇文章适合正在做课程设计、毕业论文或电力系统入门仿真的人，也适合想从“调用MATPOWER”转向“自己写一遍”的工程师。它能帮你落地一个可复现的极坐标牛拉法潮流程序，并理解为什么初值差一点就发散、为什么平衡节点必须存在、为什么PV节点要检查无功越限。这些细节在MATPOWER里是黑盒，在手写代码里是必须面对的取舍。
2. 牛拉法潮流计算的数学模型与极坐标雅可比矩阵
2.1 节点分类与功率方程
潮流计算的第一步是给节点贴标签。电力系统里的节点分成三类：PQ节点是有功、无功都给定的负荷节点；PV节点是有功和电压幅值给定、无功待求的发电机节点；平衡节点（也叫Vθ节点）则承担系统功率差额，电压幅值和相角都是已知量。绝大多数算例里，平衡节点只有1个，其余节点根据运行方式在PQ和PV之间切换。
极坐标形式的节点功率方程为：
- P_i = U_i * Σ U_j (G_ij cosθ_ij + B_ij sinθ_ij)
- Q_i = U_i * Σ U_j (G_ij sinθ_ij - B_ij cosθ_ij)
其中θ_ij = θ_i - θ_j。对PQ节点，已知P_i、Q_i，未知θ_i、U_i；对PV节点，已知P_i、U_i，未知θ_i、Q_i；平衡节点已知θ_i、U_i，不需要参加迭代。因此方程组中，PQ节点有2个方程，PV节点有1个有功方程，待求变量就是所有非平衡节点的θ以及所有PQ节点的U。
这段方程的理解直接决定代码数组的编号方式。我一般会把节点重新编号：先把PQ节点排在一起，再排PV节点，最后是平衡节点。这样做的好处是雅可比矩阵的分块结构非常清晰，PowerWorld和MATPOWER的内部数据格式也类似，方便后续对比。
2.2 雅可比矩阵的构造
令待求量x = [θ_PQ; U_PQ; θ_PV]，方程F = [ΔP; ΔQ]，其中ΔP对PQ和PV节点的有功不平衡量，ΔQ只对PQ节点。牛拉法的修正方程为：
[J] * [Δx] = -[F]
矩阵J是2×2分块形式：
- J1 = ∂ΔP / ∂θ
- J2 = ∂ΔP / ∂U，此时右侧乘U修正量，即ΔU/U
- J3 = ∂ΔQ / ∂θ
- J4 = ∂ΔQ / ∂U
手写代码时，最容易出错的是对角元和非对角元公式。以极坐标形式为例，雅可比矩阵的非对角元（i ≠ j）：
- H_ij = ∂P_i / ∂θ_j = -U_i U_j (G_ij sinθ_ij - B_ij cosθ_ij)
- N_ij = U_j * ∂P_i / ∂U_j = -U_i U_j (G_ij cosθ_ij + B_ij sinθ_ij)
- J_ij = ∂Q_i / ∂θ_j = U_i U_j (G_ij cosθ_ij + B_ij sinθ_ij)
- L_ij = U_j * ∂Q_i / ∂U_j = -U_i U_j (G_ij sinθ_ij - B_ij cosθ_ij)
对角元则要从原始求和公式对θ_i或U_i求偏导，会额外多出一项与i节点自导纳和注入功率相关的项。具体来说：
- H_ii = -Q_i - B_ii * U_i^2
- N_ii = P_i + G_ii * U_i^2
- J_ii = P_i - G_ii * U_i^2
- L_ii = Q_i - B_ii * U_i^2
这里Q_i和P_i是该节点当前迭代点的注入功率。很多初学者会直接套用非对角元公式到对角元位置，结果第一次迭代就偏差很大。
2.3 修正方程与收敛判据
修正方程是一个实线性方程组，维度等于2n_PQ + n_PV。Matlab里直接用反斜杠运算符求解：
dx = -J \ F;解出的dx里前半部分是Δθ，后半部分是ΔU/U。注意对于PV节点，ΔQ不需要解，所以F里对应位置只有ΔP。实际迭代时，为了改善收敛性，可以引入阻尼因子，比如把修正量乘以0.5或1.0，看系统是否朝功率不平衡量减小的方向走。
收敛判据常用两种：一种是所有节点功率不平衡量的绝对值最大值小于给定阈值，例如1e-6；另一种是Δx的范数小于阈值。前者物理意义更直接，因为潮流最终要求的是节点注入功率和线路损耗平衡。我一般同时监控两个值，并且把阈值设在1e-8以匹配双精度浮点极限。
3. 用Matlab从零实现牛拉法潮流的核心函数
3.1 数据准备：节点与支路参数数组
手写潮流代码的第一步是把算例数据组织成Matlab友好的数组。节点数组的每一行代表一个节点，列含义依次为：节点编号、类型（1表示PQ，2表示PV，3表示平衡）、有功注入P、无功注入Q、电压幅值初值、电压相角初值、无功下限Qmin、无功上限Qmax。支路数组每一行代表一条支路，列为：起始节点、终止节点、电阻R、电抗X、对地电纳B（半导通纳，实际单侧为B/2）、变压器变比k（非变比则为1）。
下面是一组典型的IEEE 14节点部分数据格式，完整14节点可以在MATPOWER的case14.m里找到并转成这个表格。
% 节点 [编号, 类型, P, Q, Vmag, Vang(rad), Qmin, Qmax] bus = [ 1 3 0 0 1.06 0 -999 999; 2 2 0.217 0.127 1.045 0 -0.4 0.5; 3 2 0.942 0.190 1.010 0 -0.4 0.4; 4 1 0.478 -0.039 1.0 0 0 0; 5 1 0.076 0.016 1.0 0 0 0; ]; % 支路 [首端, 末端, R, X, B/2, 变比] branch = [ 1 2 0.01938 0.05917 0.0264 1; 1 5 0.05403 0.22304 0.0246 1; 2 3 0.04699 0.19797 0.0219 1; 2 4 0.05811 0.17632 0.0187 1; 2 5 0.05695 0.17388 0.0170 1; ];参数说明：
- 类型字段不会直接参与功率方程，但决定了待求变量和方程数。
- 电压幅值初值：PQ节点通常给1.0（标幺值），PV节点给实际给定值。
- 相角初值：全部给0，牛拉法对平启动的鲁棒性算是不错的。
3.2 主迭代流程代码
完整的潮流求解函数如下。函数接收bus和branch矩阵，返回节点电压幅值、相角以及迭代信息。
function [V, theta, iter, error] = nr_power_flow(bus, branch, tol, max_iter) % 牛拉法极坐标潮流计算 % 输入bus, branch; tol为收敛阈值; max_iter为最大迭代次数 n = size(bus, 1); % 节点重新编号：PQ在前，PV次之，平衡最后 pq = find(bus(:,2) == 1); pv = find(bus(:,2) == 2); sl = find(bus(:,2) == 3); n_pq = length(pq); n_pv = length(pv); % 提取参数 Y = build_ybus(bus, branch); % 构建节点导纳矩阵，稍后实现 G = real(Y); B = imag(Y); V = bus(:,5); % 电压幅值向量 theta = bus(:,6); % 相角向量 % 已知注入功率（标幺值） P_spec = bus(:,3); Q_spec = bus(:,4); % 记录迭代 error = zeros(max_iter, 1); for iter = 1:max_iter % 计算当前注入功率（极坐标公式） U = V; TH = theta; P_calc = zeros(n,1); Q_calc = zeros(n,1); for i = 1:n for j = 1:n dth = TH(i) - TH(j); P_calc(i) = P_calc(i) + U(i)*U(j)*(G(i,j)*cos(dth) + B(i,j)*sin(dth)); Q_calc(i) = Q_calc(i) + U(i)*U(j)*(G(i,j)*sin(dth) - B(i,j)*cos(dth)); end end % 功率不平衡量：PQ节点算P和Q，PV节点只算P dP = P_spec - P_calc; dQ = Q_spec - Q_calc; % 组装F向量：先PQ的dP和dQ，再PV的dP F = [dP(pq); dQ(pq); dP(pv)]; % 计算雅可比矩阵 J = jacobian_calculation(U, TH, G, B, pq, pv); % 求解修正方程 dx = -J \ F; % 拆分修正量：前n_pq个是dtheta_pq，接着n_pq个是dU_pq/U_pq，最后n_pv个是dtheta_pv idx_theta = [pq; pv]; % 所有非平衡节点的相角顺序 dtheta = zeros(n,1); dtheta(idx_theta) = dx(1:n_pq+n_pv); dU = zeros(n,1); dU(pq) = dx(n_pq+n_pv+1:end); % 更新状态 theta = theta + dtheta; V(pq) = V(pq) .* (1 + dU(pq)); % 检查收敛 err = max(abs(F)); error(iter) = err; if err < tol V_ret = V; theta_ret = theta; return; end % PV节点无功越限检查（常放在迭代收敛后，但此处可提前处理） % 实际工程中会在每轮迭代后计算PV节点无功，若越限则转为PQ节点 end warning('未收敛，最大迭代次数 %d 达到', max_iter); V = []; theta = []; end逻辑说明：
- 雅可比矩阵是用独立函数
jacobian_calculation计算的，这样主循环更容易读。 - 更新电压的方式是
V(pq) * (1 + dU)，因为牛顿法的修正量约等于ΔU/U，这样做相当于把修正量映射成相对变化，数值上更稳定。 - 如果模型包含PV节点无功越限，需要在每次迭代收敛后计算Q_pv，若超出限值，则将该节点转为PQ节点并重新迭代。这段代码里先预留了处理空间，避免一次引入过多逻辑干扰主线。
3.3 雅可比矩阵与功率不平衡量的计算函数
雅可比矩阵计算子函数如下：
function J = jacobian_calculation(U, TH, G, B, pq, pv) % 计算极坐标牛拉法雅可比矩阵 non_slack = [pq; pv]; n1 = length(non_slack); % theta维度 n2 = length(pq); % U维度 n = n1 + n2; % 初始化分块 H = zeros(n1, n1); % dP/dtheta N = zeros(n1, n2); % dP/dU，乘U M = zeros(n2, n1); % dQ/dtheta L = zeros(n2, n2); % dQ/dU，乘U nodes = [pq; pv; find(~ismember(1:length(U), [pq; pv]))]; % 这个不是最终顺序 % 更稳妥做法：直接遍历所有节点并映射到位置 idx_theta = [pq; pv]; % 补齐，之后填充 for ii = 1:length(idx_theta) i = idx_theta(ii); for jj = 1:length(idx_theta) j = idx_theta(jj); if i == j % 对角元：需要依赖当前注入功率 % 先计算该节点当前注入功率 Pi = 0; Qi = 0; for k = 1:length(U) dth = TH(i) - TH(k); Pi = Pi + U(i)*U(k)*(G(i,k)*cos(dth) + B(i,k)*sin(dth)); Qi = Qi + U(i)*U(k)*(G(i,k)*sin(dth) - B(i,k)*cos(dth)); end H(ii,jj) = -Qi - B(i,i)*U(i)^2; M(ii,jj) = Pi - G(i,i)*U(i)^2; % 注意此处的行索引对齐后面会重新整理 else dth = TH(i) - TH(j); H(ii,jj) = U(i)*U(j)*(G(i,j)*sin(dth) - B(i,j)*cos(dth)); N(ii,jj) = -U(i)*U(j)*(G(i,j)*cos(dth) + B(i,j)*sin(dth)); M(ii,jj) = -U(i)*U(j)*(G(i,j)*cos(dth) + B(i,j)*sin(dth)); L(ii,jj) = -U(i)*U(j)*(G(i,j)*sin(dth) - B(i,j)*cos(dth)); end end end % 重新按PQ/PV顺序编制成的H是完整n1×n1，但对PQ节点的ΔP和ΔQ分块不同 % 这里由于前面用同一个循环填充，会导致M错位。需要更仔细的行列映射。 % 完整实现下面单独给一个按公式逐块构建的版本。 end上面的代码为了展示易错点，故意留了一个映射不严谨的版本。实际项目里，我会按更清晰的分块方式写：
function J = jacobian_calculation(U, TH, G, B, pq, pv) % 雅可比矩阵，按[PQ dtheta, PV dtheta, PQ dU]的顺序排列 n_pq = length(pq); n_pv = length(pv); n_theta = n_pq + n_pv; H = zeros(n_theta, n_theta); N = zeros(n_theta, n_pq); M = zeros(n_pq, n_theta); L = zeros(n_pq, n_pq); % 计算公式：对每个非平衡节点i（包括pq和pv），每个节点j for i_idx = 1:n_theta i = [pq; pv](i_idx); for j_idx = 1:n_theta j = [pq; pv](j_idx); if i == j Pi = 0; Qi = 0; for k = 1:length(U) dth = TH(i) - TH(k); Pi = Pi + U(i)*U(k)*(G(i,k)*cos(dth) + B(i,k)*sin(dth)); Qi = Qi + U(i)*U(k)*(G(i,k)*sin(dth) - B(i,k)*cos(dth)); end H(i_idx, j_idx) = -Qi - B(i,i)*U(i)^2; else dth = TH(i) - TH(j); H(i_idx, j_idx) = U(i)*U(j)*(G(i,j)*sin(dth) - B(i,j)*cos(dth)); end % N的列只对应PQ节点 if ismember(j, pq) j_col = find(pq == j); if i == j Pi = 0; for k = 1:length(U) dth = TH(i) - TH(k); Pi = Pi + U(i)*U(k)*(G(i,k)*cos(dth) + B(i,k)*sin(dth)); end N(i_idx, j_col) = Pi + G(i,i)*U(i)^2; else dth = TH(i) - TH(j); N(i_idx, j_col) = -U(i)*U(j)*(G(i,j)*cos(dth) + B(i,j)*sin(dth)); end end end end % M和L行对应PQ节点 for i_idx = 1:n_pq i = pq(i_idx); for j_idx = 1:n_theta j = [pq; pv](j_idx); if i == j Pi = 0; for k = 1:length(U) dth = TH(i) - TH(k); Pi = Pi + U(i)*U(k)*(G(i,k)*cos(dth) + B(i,k)*sin(dth)); end M(i_idx, j_idx) = Pi - G(i,i)*U(i)^2; else dth = TH(i) - TH(j); M(i_idx, j_idx) = -U(i)*U(j)*(G(i,j)*cos(dth) + B(i,j)*sin(dth)); end if ismember(j, pq) j_col = find(pq == j); if i == j Qi = 0; for k = 1:length(U) dth = TH(i) - TH(k); Qi = Qi + U(i)*U(k)*(G(i,k)*sin(dth) - B(i,k)*cos(dth)); end L(i_idx, j_col) = Qi - B(i,i)*U(i)^2; else dth = TH(i) - TH(j); L(i_idx, j_col) = -U(i)*U(j)*(G(i,j)*sin(dth) - B(i,j)*cos(dth)); end end end end J = [H N; M L]; end逻辑说明：
- N和L矩阵的列索引需要从节点编号映射回雅可比矩阵中的位置，这里用
find(pq==j)实现。 - 对角元里使用当前迭代点的注入功率Pi、Qi，所以每次迭代都需要重新计算。
- 这个版本的公式适用于标幺制，电压单位不是kV而是相对值，所以G和B也来自标幺导纳。
4. 在IEEE 14节点和实际算例上验证与排错
4.1 用IEEE 14节点数据跑通程序
完整IEEE 14节点算例包含14个节点、20条支路，其中平衡节点1个（1号）、PV节点4个（2、3、6、8号）、PQ节点9个。把前面写的nr_power_flow函数和build_ybus函数放在同一个工作目录，然后运行：
[V, theta, iter, err] = nr_power_flow(bus14, branch14, 1e-8, 30); % 打印结果 for i = 1:14 fprintf('节点%d 电压幅值%.6f 相角%.6f度\n', i, V(i), theta(i)*180/pi); end如果数据正确，2~3次迭代就能看到不平衡量降到1e-6以下，5次以内收敛到1e-10。典型的前几次误差序列大致是1e-2、1e-4、1e-8，呈二次收敛特征。如果第二次迭代误差反而变大，说明雅可比矩阵有笔误，重点查对角元符号。
build_ybus函数的实现是标准做法，需要处理变压器变比和线路对地导纳折算到哪一侧，这里不展开完整代码，但有一个关键点：变压器支路的等效导纳矩阵非对称，变比在首端时，Y12 = -y/k，Y21 = -y/k，Y11 = (y + y/2)/k^2，Y22 = y + y/2。
4.2 常见收敛失败原因与参数调整
牛拉法发散几乎都逃不开几个原因。初值问题：电压初值给0或者给负值，会导致导纳矩阵计算出现奇异值，必须先给定合理平启动值。数据单位问题：电阻、电抗以标幺值给，但功率和电压用有名值，算出来的注入功率可能差10倍，收敛曲线会剧烈波动。PV节点处理缺失：系统里PV节点的无功功率没有参与迭代，若某台发电机无功越限却仍保持PV类型，雅可比矩阵会把电压硬性固定在给定值，最终结果与实际运行点不一致，有时会导致相邻PQ节点电压崩掉。
处理PV节点无功越限的通用做法是：
% 在每轮迭代收敛后，计算PV节点无功 Q_pv = compute_q(bus, V, theta, Y); % 用导纳矩阵计算 for i = pv' if Q_pv(i) < bus(i,7) || Q_pv(i) > bus(i,8) % 转为PQ节点，将Q设为限值 bus(i,2) = 1; bus(i,4) = max(bus(i,7), min(bus(i,8), Q_pv(i))); disp(['节点', num2str(i), ' 无功越限，转为PQ']); end end这个逻辑放在迭代收敛判断之后、返回结果之前。另一个常见问题是平衡节点有功、无功输出超出合理范围，需要检查系统总负荷和总出力是否平衡。
4.3 与MATPOWER结果对比验证
自己写的程序结果靠不靠谱，和MATPOWER对比是最快的方式。MATPOWER的runpf函数默认用牛顿法，也是可以指定的。执行：
mpc = case14; res = runpf(mpc);然后对比各节点电压幅值和相角。最大误差在1e-5以内基本说明算法正确。如果差异较大，优先检查导纳矩阵是否一致。可以把自建Ybus和MATPOWER里makeYbus的结果对比：
Y_mp = makeYbus(mpc); Y_self = build_ybus(bus14, branch14); disp(max(abs(Y_mp - Y_self), [], 'all'));数值不为0时，最常见的是变压器变比归算方向搞反了，或者对地电纳没有除以2。这类问题单看一个算例很难发现，和数据规模无关。
5. 提高收敛性和计算效率的三个技巧
5.1 用平启动和逐步升压方式加速收敛
平启动是最经典的初值选择：所有非平衡节点电压幅值设1.0，相角设0。但病态重负荷系统在这种初值下可能不收敛。我习惯先把所有节点当成PQ节点跑一次，再把符合条件的节点改回PV节点重新迭代，相当于两阶段启动。另一种做法是把发电机节点电压设到1.05或1.06，因为实际运行中PV节点电压通常略高于额定，这比1.0初值更接近真实运行点。
5.2 稀疏化雅可比矩阵求解
中小算例用完整矩阵的J \ F没问题，节点数量到500以上时，内存占用和计算量会明显上涨。J本身是高度稀疏的，尤其是非对角元素只存在于网络拓扑相连的节点对之间。可以使用Matlab的稀疏矩阵存储：
J = sparse(J); dx = -J \ F;在构造雅可比矩阵时一开始就分配稀疏矩阵，例如J = sparse(n, n)，然后再填非对角元。注意J \ F对稀疏矩阵会自动选用LU分解，求解速度比满载矩阵快一个数量级。这个方法在处理2000节点以上的输电网时收益明显。
5.3 用容差自动调整和可视化验证
收敛容差从1e-6改成1e-8通常不会增加太多迭代次数，但能明显提升结果用于后续最优潮流计算的稳定性。迭代完成后，把电压分布画出来，能立刻看到是否有节点电压落在合理区间之外：
figure; bar(1:14, V); xlabel('节点编号'); ylabel('电压幅值（标幺）'); title('牛拉法潮流计算结果电压分布');如果出现某节点电压低于0.9或高于1.1，要回到运行方式数据检查无功平衡，而不是盲目调算法参数。潮流结果的可信度，最终还是由网络参数和运行方式决定的，牛拉法只是把这些数据翻译成收敛的电压断面。
本文还有配套的精品资源，点击获取
