简介基于Matlab的牛顿拉夫逊法基波潮流计算程序面向电力系统专业学生、科研人员及工程技术人员以何仰赞第四版《电力系统分析》为理论支撑适用于任意规模纯交流电网的潮流计算。程序支持节点与支路灵活增删可接入多个风电、光伏等分布式电源子函数模块涵盖节点导纳矩阵、雅克比矩阵n-1m维等核心计算环节输出格式对标Matpower的runpf()函数误差低于10^-3并附带步骤讲解教程与详细注释便于快速理解算法推导、掌握编程实现要点也可直接用于课程设计或科研验证。压缩包共14个文件以13个.m脚本文件为主包含主程序、核心子函数及case14、case33bw等典型算例另附1个xlsx格式的计算结果表整体仅43KB。目前已有2553人学习下载是学习潮流计算和牛拉法编程的实用参考工具尤其适合电力系统专业的初学者与进阶开发者。1. 牛拉法基波潮流计算在MATLAB里到底要解什么做电力系统潮流计算牛顿-拉夫逊法下面直接叫牛拉法是绕不开的基准算法可一旦matlab里把基波潮流跑通最耗时间的往往不是迭代框架而是雅可比矩阵的符号、PV节点和平衡节点的取舍以及电压初值怎么给。基波潮流说的是工频稳态下的潮流所有电气量都用相量和标幺值表示最终归结成一组实数非线性方程组——每个PQ节点有两个方程每个PV节点有一个有功方程平衡节点只提供电压参考。这篇文章用一个3节点测试系统把牛拉法基波潮流的MATLAB实现从Y矩阵到主迭代按步骤写出来每个关键循环都带注释最后给出收敛判断和调试手段。2. 牛拉法基波潮流计算的数学模型节点功率方程与雅可比矩阵2.1 从复功率到节点注入方程基波潮流在解哪些量对节点 i设节点电压为V_i V_i∠θ_i节点导纳矩阵元素为Y_ij G_ij jB_ij。注入复功率可以写成S_i V_i * conj(Σ Y_ij * V_j)展开实部和虚部后得到节点有功、无功注入方程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)其中θ_ij θ_i - θ_j。这是牛拉法需要求解的非线性方程组。基波潮流不考虑谐波和不对称所有公式都在正序相量域里成立。节点类型决定了哪些方程需要参与迭代。表中列出常见三类节点节点类型给定量未知量方程组PQ节点P, QV, θ三相平衡时两个方程PV节点P, VQ, θ仅一个有功方程平衡节点V, θP, Q不建方程只提供相位参考平衡节点之所以必要是因为潮流方程只有差分量相角需要有一条绝对基准系统有功无功不平衡量最终由平衡节点吸收。2.2 修正方程的分块雅可比矩阵H、N、J、L 的定义与符号牛拉法每步迭代求解[ΔP; ΔQ] J * [Δθ; ΔV/V]注意这里ΔP P_spec - P_calcΔQ Q_spec - Q_calcJ 是运行点对状态量的偏导数矩阵。采用ΔV/V而不是ΔV电压迭代更新为V_new V V .* dV好处是雅可比对角线更简单且数值条件数更稳。雅可比矩阵按功能分成四个子块分块含义i≠j 非对角项ij 对角项H∂P/∂θV_i V_j(G_ij sinθ_ij - B_ij cosθ_ij)-Q_i - B_ii V_i^2NV∂P/∂VV_i V_j(G_ij cosθ_ij B_ij sinθ_ij)P_i G_ii V_i^2M∂Q/∂θ-V_i V_j(G_ij cosθ_ij B_ij sinθ_ij)P_i - G_ii V_i^2LV∂Q/∂VV_i V_j(G_ij sinθ_ij - B_ij cosθ_ij)Q_i - B_ii V_i^2这套公式是后面MATLAB代码的直接来源。容易出错的是符号尤其是 M 的非对角项带负号若抄错迭代会在第一次修正后直接发散。2.3 用 3 节点测试系统让公式落到数据上为了验证代码这里构造一个3节点测试系统所有数值都是标幺值基准容量 100 MVA。节点类型V初值θ初值(rad)P_specQ_spec1平衡1.050002PQ1.000-0.40-0.153PV1.0000.300支路参数按[始端, 末端, r, x, b]组织b 为线路全线的对地电纳标幺值支路rxb1-20.0100.0500.0201-30.0150.0600.0252-30.0200.0800.030这个系统的设计保证潮流解存在PV节点电压保持 1.0无功不越限。后面所有代码都围绕这套数据展开。3. MATLAB实现牛拉法基波潮流从节点支路数据到导纳矩阵3.1 数据布局与节点导纳矩阵组装在matlab里写潮流程序数据布局决定了后面所有索引的复杂度。这里使用两个矩阵bus放节点信息line放支路信息。% bus: [type, V0, theta0, Pspec, Qspec] % type: 3平衡节点, 2PV节点, 1PQ节点 bus [ 3, 1.05, 0, 0, 0; 1, 1.00, 0, -0.4, -0.15; 2, 1.00, 0, 0.3, 0; ]; % line: [i, j, r, x, b] % i,j 为节点编号r/x 为串联阻抗标幺值b 为线路对地电纳标幺值 line [ 1, 2, 0.01, 0.05, 0.02; 1, 3, 0.015, 0.06, 0.025; 2, 3, 0.02, 0.08, 0.03; ]; n size(bus,1); Y zeros(n,n); for k 1:size(line,1) i line(k,1); j line(k,2); z line(k,3) 1j*line(k,4); y 1/z; b line(k,5); % 自导纳加上本支路导纳和对地电纳的一半 Y(i,i) Y(i,i) y 1j*b/2; Y(j,j) Y(j,j) y 1j*b/2; % 互导纳为负支路导纳 Y(i,j) Y(i,j) - y; Y(j,i) Y(j,i) - y; end G real(Y); B imag(Y); disp(Y);这段代码把每条支路的串联阻抗转成导纳再加入并联对地电纳。注意b是全线总电纳分配到两端各取一半这是高压线路模型的标准做法。Y(i,j)取负号因为节点导纳矩阵互导纳定义为负的支路导纳。3.2 初值设定与待求变量索引牛拉法看起来不需要太好的初值但初值太差会让雅可比矩阵接近奇异。通用做法是 flat startPQ节点电压幅值取 1.0所有节点相角取 0PV节点使用给定的电压幅值。V bus(:,2); % 电压幅值初值 theta bus(:,3); % 相角初值单位弧度 type bus(:,1); Pspec bus(:,4); Qspec bus(:,5); slack find(type 3); % 平衡节点 non_slack find(type ~ 3); % 需要theta方程的节点 PQ find(type 1); % 需要电压幅值方程的节点 PV find(type 2); % 有功方程固定电压幅值索引数组后面反复使用。non_slack排除平衡节点的相角PQ排除PV节点的无功方程PV仅用于更新后强制恢复电压幅值。3.3 导纳矩阵构建容易踩的坑实际写这段代码时最容易出错的不是循环而是忘记把线路电容分成两端。b/2只影响自导纳若只给一侧潮流结果看起来也像样但对地充电功率会少一半在长线路、高电压等级系统中误差立刻暴露。另一个坑是标幺值。r、x必须基于同一基准电压和基准功率换算不能把工程单位欧姆直接填进去。matlab不会告警但 Y 矩阵数量级不对时雅可比矩阵也会失控。第三个坑是互导纳符号。Y(i,j)取-y而不是y否则自导纳和互导纳正负关系全反迭代几乎不可能收敛。4. 牛顿-拉夫逊主迭代循环的完整MATLAB实现带详细注释4.1 计算不平衡量 ΔP / ΔQ每次迭代的第一步用当前电压和相角计算节点注入功率。P zeros(n,1); Q zeros(n,1); for i 1:n % theta(i)-theta 表示 θ_i - θ_j P(i) V(i) * sum( V .* ( G(i,:).*cos(theta(i)-theta) B(i,:).*sin(theta(i)-theta) ) ); Q(i) V(i) * sum( V .* ( G(i,:).*sin(theta(i)-theta) - B(i,:).*cos(theta(i)-theta) ) ); end % 给定值减去计算值 dP Pspec - P; dQ Qspec - Q; % 只保留参与求解的方程非平衡节点P方程 PQ节点Q方程 dF [dP(non_slack); dQ(PQ)];平衡节点的Pspec和Qspec是占位 0不参与方程。PV节点的Qspec也是占位 0因为无功不是给定值。dF是最终用于收敛判断的残差向量。4.2 解析组装雅可比矩阵雅可比矩阵的组装是牛拉法核心。这里先用稠密矩阵完整建一个2n × 2n矩阵再裁剪成实际需要的子矩阵方便对照 2.2 节的公式。% 完整雅可比矩阵前n行为P方程后n行为Q方程 % 列前n列为θ后n列为ΔV/V对应的偏导 J zeros(2*n, 2*n); for i 1:n for j 1:n if i j % 对角项 J(i,i) -Q(i) - B(i,i)*V(i)^2; J(i,ni) P(i) G(i,i)*V(i)^2; J(ni,i) P(i) - G(i,i)*V(i)^2; J(ni,ni) Q(i) - B(i,i)*V(i)^2; else Vij V(i)*V(j); th theta(i) - theta(j); % H块dP_i/dθ_j J(i,j) Vij * (G(i,j)*sin(th) - B(i,j)*cos(th)); % N块V_j * dP_i/dV_j J(i,nj) Vij * (G(i,j)*cos(th) B(i,j)*sin(th)); % M块dQ_i/dθ_j J(ni,j) -Vij * (G(i,j)*cos(th) B(i,j)*sin(th)); % L块V_j * dQ_i/dV_j J(ni,nj) Vij * (G(i,j)*sin(th) - B(i,j)*cos(th)); end end end % 裁剪行和列 row_idx [non_slack, n PQ]; col_idx [non_slack, n PQ]; J_red J(row_idx, col_idx);为什么能直接裁剪因为平衡节点的θ固定PV节点的V固定对应列不会出现在修正方程里平衡节点的P/Q方程、PV节点的Q方程也不参与残差。这样构造出的J_red是方阵维度等于length(non_slack) length(PQ)。4.3 求解修正量与更新把 4.1 和 4.2 的内容组装进主循环形成完整可运行的牛拉法基波潮流程序。tol 1e-8; maxiter 50; for iter 1:maxiter % ---------- 计算P,Q和残差 ---------- P zeros(n,1); Q zeros(n,1); for i 1:n P(i) V(i) * sum( V .* ( G(i,:).*cos(theta(i)-theta) B(i,:).*sin(theta(i)-theta) ) ); Q(i) V(i) * sum( V .* ( G(i,:).*sin(theta(i)-theta) - B(i,:).*cos(theta(i)-theta) ) ); end dP Pspec - P; dQ Qspec - Q; dF [dP(non_slack); dQ(PQ)]; % ---------- 收敛判断 ---------- if max(abs(dF)) tol fprintf(第 %d 次迭代收敛\n, iter); break; end % ---------- 组装并裁剪雅可比矩阵 ---------- J zeros(2*n, 2*n); for i 1:n for j 1:n if i j J(i,i) -Q(i) - B(i,i)*V(i)^2; J(i,ni) P(i) G(i,i)*V(i)^2; J(ni,i) P(i) - G(i,i)*V(i)^2; J(ni,ni) Q(i) - B(i,i)*V(i)^2; else Vij V(i)*V(j); th theta(i) - theta(j); J(i,j) Vij * (G(i,j)*sin(th) - B(i,j)*cos(th)); J(i,nj) Vij * (G(i,j)*cos(th) B(i,j)*sin(th)); J(ni,j) -Vij * (G(i,j)*cos(th) B(i,j)*sin(th)); J(ni,nj) Vij * (G(i,j)*sin(th) - B(i,j)*cos(th)); end end end row_idx [non_slack, n PQ]; col_idx [non_slack, n PQ]; J_red J(row_idx, col_idx); % ---------- 解线性方程组 ---------- dx J_red \ dF; % ---------- 拆解修正量 ---------- dtheta zeros(n,1); dV_scaled zeros(n,1); dtheta(non_slack) dx(1:length(non_slack)); dV_scaled(PQ) dx(length(non_slack)1:end); % ---------- 更新状态量 ---------- theta theta dtheta; V V V .* dV_scaled; % PV节点电压幅值强制保持设定值 V(PV) bus(PV,2); end if iter maxiter warning(达到最大迭代次数未收敛); end这段代码使用matlab里的反斜杠算子求解线性方程组。对于稠密小矩阵反斜杠会自动选LU分解节点数较小时性能足够。dV_scaled对应的是ΔV/V因此更新电压时要乘回V这是最容易和直接V V dV混淆的地方。5. 结果校验、稀疏化与牛顿-拉夫逊调试技巧5.1 功率平衡校验与PV节点无功输出迭代收敛不代表结果正确。我一般会先用复数形式从解出的电压反推一次注入功率和节点功率方程交叉验证。Vc V .* exp(1j*theta); I Y * Vc; S_check Vc .* conj(I); fprintf(平衡节点注入功率: %.6f j%.6f\n, real(S_check(slack)), imag(S_check(slack))); fprintf(PV节点无功输出: %.6f\n, imag(S_check(3)));S_check是基于Y矩阵和最终电压直接计算得到的注入复功率如果P、Q计算和雅可比公式里没有符号错误它应基本等于收敛后的P 1j*Q。另外可以改变初始电压再做一次潮流看是否收敛到同一解。牛拉法可能收敛到不同的运行点但在这个小系统里重复试验能暴露初值相关的错误。5.2 用稀疏矩阵重写雅可比去掉双层循环的提速方式当系统规模超过几百个节点时稠密J_red的内存和求解时间会迅速膨胀。常见做法是在组装时直接用稀疏三元组。rows []; cols []; vals []; for i 1:n for j 1:n if Y(i,j) ~ 0 % 这里 val 按4.2中的公式计算 rows [rows; i]; cols [cols; j]; vals [vals; val]; end end end J_sp sparse(rows, cols, vals, 2*n, 2*n);实际工程可以在循环里用preallocate预分配数组避免反复扩容。matlab的\对稀疏矩阵会自动选择稀疏LU分解不需要手动指定算法。还有一个很实用的调试技巧当第一次迭代没有收敛而残差序列出现震荡时多半是初值太差或雅可比符号有问题。在确认正负号无误后可以给修正量加阻尼。alpha 0.8; theta theta alpha * dtheta; V V V .* (alpha * dV_scaled);alpha从 1 开始残差增大时改成 0.5 甚至 0.2往往能把发散拖回来。这个阻尼思路和计算机视觉里的Levenberg-Marquardt加正则的思路类似对牛拉法潮流初值敏感的问题很有效。把小系统跑顺之后再换IEEE标准节点数据前也值得多看一眼收敛过程中的norm(dF,inf)是否单调下降。本文还有配套的精品资源点击获取 MATLAB实现牛拉法基波潮流计算：从节点导纳矩阵到雅可比矩阵
简介基于Matlab的牛顿拉夫逊法基波潮流计算程序面向电力系统专业学生、科研人员及工程技术人员以何仰赞第四版《电力系统分析》为理论支撑适用于任意规模纯交流电网的潮流计算。程序支持节点与支路灵活增删可接入多个风电、光伏等分布式电源子函数模块涵盖节点导纳矩阵、雅克比矩阵n-1m维等核心计算环节输出格式对标Matpower的runpf()函数误差低于10^-3并附带步骤讲解教程与详细注释便于快速理解算法推导、掌握编程实现要点也可直接用于课程设计或科研验证。压缩包共14个文件以13个.m脚本文件为主包含主程序、核心子函数及case14、case33bw等典型算例另附1个xlsx格式的计算结果表整体仅43KB。目前已有2553人学习下载是学习潮流计算和牛拉法编程的实用参考工具尤其适合电力系统专业的初学者与进阶开发者。1. 牛拉法基波潮流计算在MATLAB里到底要解什么做电力系统潮流计算牛顿-拉夫逊法下面直接叫牛拉法是绕不开的基准算法可一旦matlab里把基波潮流跑通最耗时间的往往不是迭代框架而是雅可比矩阵的符号、PV节点和平衡节点的取舍以及电压初值怎么给。基波潮流说的是工频稳态下的潮流所有电气量都用相量和标幺值表示最终归结成一组实数非线性方程组——每个PQ节点有两个方程每个PV节点有一个有功方程平衡节点只提供电压参考。这篇文章用一个3节点测试系统把牛拉法基波潮流的MATLAB实现从Y矩阵到主迭代按步骤写出来每个关键循环都带注释最后给出收敛判断和调试手段。2. 牛拉法基波潮流计算的数学模型节点功率方程与雅可比矩阵2.1 从复功率到节点注入方程基波潮流在解哪些量对节点 i设节点电压为V_i V_i∠θ_i节点导纳矩阵元素为Y_ij G_ij jB_ij。注入复功率可以写成S_i V_i * conj(Σ Y_ij * V_j)展开实部和虚部后得到节点有功、无功注入方程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)其中θ_ij θ_i - θ_j。这是牛拉法需要求解的非线性方程组。基波潮流不考虑谐波和不对称所有公式都在正序相量域里成立。节点类型决定了哪些方程需要参与迭代。表中列出常见三类节点节点类型给定量未知量方程组PQ节点P, QV, θ三相平衡时两个方程PV节点P, VQ, θ仅一个有功方程平衡节点V, θP, Q不建方程只提供相位参考平衡节点之所以必要是因为潮流方程只有差分量相角需要有一条绝对基准系统有功无功不平衡量最终由平衡节点吸收。2.2 修正方程的分块雅可比矩阵H、N、J、L 的定义与符号牛拉法每步迭代求解[ΔP; ΔQ] J * [Δθ; ΔV/V]注意这里ΔP P_spec - P_calcΔQ Q_spec - Q_calcJ 是运行点对状态量的偏导数矩阵。采用ΔV/V而不是ΔV电压迭代更新为V_new V V .* dV好处是雅可比对角线更简单且数值条件数更稳。雅可比矩阵按功能分成四个子块分块含义i≠j 非对角项ij 对角项H∂P/∂θV_i V_j(G_ij sinθ_ij - B_ij cosθ_ij)-Q_i - B_ii V_i^2NV∂P/∂VV_i V_j(G_ij cosθ_ij B_ij sinθ_ij)P_i G_ii V_i^2M∂Q/∂θ-V_i V_j(G_ij cosθ_ij B_ij sinθ_ij)P_i - G_ii V_i^2LV∂Q/∂VV_i V_j(G_ij sinθ_ij - B_ij cosθ_ij)Q_i - B_ii V_i^2这套公式是后面MATLAB代码的直接来源。容易出错的是符号尤其是 M 的非对角项带负号若抄错迭代会在第一次修正后直接发散。2.3 用 3 节点测试系统让公式落到数据上为了验证代码这里构造一个3节点测试系统所有数值都是标幺值基准容量 100 MVA。节点类型V初值θ初值(rad)P_specQ_spec1平衡1.050002PQ1.000-0.40-0.153PV1.0000.300支路参数按[始端, 末端, r, x, b]组织b 为线路全线的对地电纳标幺值支路rxb1-20.0100.0500.0201-30.0150.0600.0252-30.0200.0800.030这个系统的设计保证潮流解存在PV节点电压保持 1.0无功不越限。后面所有代码都围绕这套数据展开。3. MATLAB实现牛拉法基波潮流从节点支路数据到导纳矩阵3.1 数据布局与节点导纳矩阵组装在matlab里写潮流程序数据布局决定了后面所有索引的复杂度。这里使用两个矩阵bus放节点信息line放支路信息。% bus: [type, V0, theta0, Pspec, Qspec] % type: 3平衡节点, 2PV节点, 1PQ节点 bus [ 3, 1.05, 0, 0, 0; 1, 1.00, 0, -0.4, -0.15; 2, 1.00, 0, 0.3, 0; ]; % line: [i, j, r, x, b] % i,j 为节点编号r/x 为串联阻抗标幺值b 为线路对地电纳标幺值 line [ 1, 2, 0.01, 0.05, 0.02; 1, 3, 0.015, 0.06, 0.025; 2, 3, 0.02, 0.08, 0.03; ]; n size(bus,1); Y zeros(n,n); for k 1:size(line,1) i line(k,1); j line(k,2); z line(k,3) 1j*line(k,4); y 1/z; b line(k,5); % 自导纳加上本支路导纳和对地电纳的一半 Y(i,i) Y(i,i) y 1j*b/2; Y(j,j) Y(j,j) y 1j*b/2; % 互导纳为负支路导纳 Y(i,j) Y(i,j) - y; Y(j,i) Y(j,i) - y; end G real(Y); B imag(Y); disp(Y);这段代码把每条支路的串联阻抗转成导纳再加入并联对地电纳。注意b是全线总电纳分配到两端各取一半这是高压线路模型的标准做法。Y(i,j)取负号因为节点导纳矩阵互导纳定义为负的支路导纳。3.2 初值设定与待求变量索引牛拉法看起来不需要太好的初值但初值太差会让雅可比矩阵接近奇异。通用做法是 flat startPQ节点电压幅值取 1.0所有节点相角取 0PV节点使用给定的电压幅值。V bus(:,2); % 电压幅值初值 theta bus(:,3); % 相角初值单位弧度 type bus(:,1); Pspec bus(:,4); Qspec bus(:,5); slack find(type 3); % 平衡节点 non_slack find(type ~ 3); % 需要theta方程的节点 PQ find(type 1); % 需要电压幅值方程的节点 PV find(type 2); % 有功方程固定电压幅值索引数组后面反复使用。non_slack排除平衡节点的相角PQ排除PV节点的无功方程PV仅用于更新后强制恢复电压幅值。3.3 导纳矩阵构建容易踩的坑实际写这段代码时最容易出错的不是循环而是忘记把线路电容分成两端。b/2只影响自导纳若只给一侧潮流结果看起来也像样但对地充电功率会少一半在长线路、高电压等级系统中误差立刻暴露。另一个坑是标幺值。r、x必须基于同一基准电压和基准功率换算不能把工程单位欧姆直接填进去。matlab不会告警但 Y 矩阵数量级不对时雅可比矩阵也会失控。第三个坑是互导纳符号。Y(i,j)取-y而不是y否则自导纳和互导纳正负关系全反迭代几乎不可能收敛。4. 牛顿-拉夫逊主迭代循环的完整MATLAB实现带详细注释4.1 计算不平衡量 ΔP / ΔQ每次迭代的第一步用当前电压和相角计算节点注入功率。P zeros(n,1); Q zeros(n,1); for i 1:n % theta(i)-theta 表示 θ_i - θ_j P(i) V(i) * sum( V .* ( G(i,:).*cos(theta(i)-theta) B(i,:).*sin(theta(i)-theta) ) ); Q(i) V(i) * sum( V .* ( G(i,:).*sin(theta(i)-theta) - B(i,:).*cos(theta(i)-theta) ) ); end % 给定值减去计算值 dP Pspec - P; dQ Qspec - Q; % 只保留参与求解的方程非平衡节点P方程 PQ节点Q方程 dF [dP(non_slack); dQ(PQ)];平衡节点的Pspec和Qspec是占位 0不参与方程。PV节点的Qspec也是占位 0因为无功不是给定值。dF是最终用于收敛判断的残差向量。4.2 解析组装雅可比矩阵雅可比矩阵的组装是牛拉法核心。这里先用稠密矩阵完整建一个2n × 2n矩阵再裁剪成实际需要的子矩阵方便对照 2.2 节的公式。% 完整雅可比矩阵前n行为P方程后n行为Q方程 % 列前n列为θ后n列为ΔV/V对应的偏导 J zeros(2*n, 2*n); for i 1:n for j 1:n if i j % 对角项 J(i,i) -Q(i) - B(i,i)*V(i)^2; J(i,ni) P(i) G(i,i)*V(i)^2; J(ni,i) P(i) - G(i,i)*V(i)^2; J(ni,ni) Q(i) - B(i,i)*V(i)^2; else Vij V(i)*V(j); th theta(i) - theta(j); % H块dP_i/dθ_j J(i,j) Vij * (G(i,j)*sin(th) - B(i,j)*cos(th)); % N块V_j * dP_i/dV_j J(i,nj) Vij * (G(i,j)*cos(th) B(i,j)*sin(th)); % M块dQ_i/dθ_j J(ni,j) -Vij * (G(i,j)*cos(th) B(i,j)*sin(th)); % L块V_j * dQ_i/dV_j J(ni,nj) Vij * (G(i,j)*sin(th) - B(i,j)*cos(th)); end end end % 裁剪行和列 row_idx [non_slack, n PQ]; col_idx [non_slack, n PQ]; J_red J(row_idx, col_idx);为什么能直接裁剪因为平衡节点的θ固定PV节点的V固定对应列不会出现在修正方程里平衡节点的P/Q方程、PV节点的Q方程也不参与残差。这样构造出的J_red是方阵维度等于length(non_slack) length(PQ)。4.3 求解修正量与更新把 4.1 和 4.2 的内容组装进主循环形成完整可运行的牛拉法基波潮流程序。tol 1e-8; maxiter 50; for iter 1:maxiter % ---------- 计算P,Q和残差 ---------- P zeros(n,1); Q zeros(n,1); for i 1:n P(i) V(i) * sum( V .* ( G(i,:).*cos(theta(i)-theta) B(i,:).*sin(theta(i)-theta) ) ); Q(i) V(i) * sum( V .* ( G(i,:).*sin(theta(i)-theta) - B(i,:).*cos(theta(i)-theta) ) ); end dP Pspec - P; dQ Qspec - Q; dF [dP(non_slack); dQ(PQ)]; % ---------- 收敛判断 ---------- if max(abs(dF)) tol fprintf(第 %d 次迭代收敛\n, iter); break; end % ---------- 组装并裁剪雅可比矩阵 ---------- J zeros(2*n, 2*n); for i 1:n for j 1:n if i j J(i,i) -Q(i) - B(i,i)*V(i)^2; J(i,ni) P(i) G(i,i)*V(i)^2; J(ni,i) P(i) - G(i,i)*V(i)^2; J(ni,ni) Q(i) - B(i,i)*V(i)^2; else Vij V(i)*V(j); th theta(i) - theta(j); J(i,j) Vij * (G(i,j)*sin(th) - B(i,j)*cos(th)); J(i,nj) Vij * (G(i,j)*cos(th) B(i,j)*sin(th)); J(ni,j) -Vij * (G(i,j)*cos(th) B(i,j)*sin(th)); J(ni,nj) Vij * (G(i,j)*sin(th) - B(i,j)*cos(th)); end end end row_idx [non_slack, n PQ]; col_idx [non_slack, n PQ]; J_red J(row_idx, col_idx); % ---------- 解线性方程组 ---------- dx J_red \ dF; % ---------- 拆解修正量 ---------- dtheta zeros(n,1); dV_scaled zeros(n,1); dtheta(non_slack) dx(1:length(non_slack)); dV_scaled(PQ) dx(length(non_slack)1:end); % ---------- 更新状态量 ---------- theta theta dtheta; V V V .* dV_scaled; % PV节点电压幅值强制保持设定值 V(PV) bus(PV,2); end if iter maxiter warning(达到最大迭代次数未收敛); end这段代码使用matlab里的反斜杠算子求解线性方程组。对于稠密小矩阵反斜杠会自动选LU分解节点数较小时性能足够。dV_scaled对应的是ΔV/V因此更新电压时要乘回V这是最容易和直接V V dV混淆的地方。5. 结果校验、稀疏化与牛顿-拉夫逊调试技巧5.1 功率平衡校验与PV节点无功输出迭代收敛不代表结果正确。我一般会先用复数形式从解出的电压反推一次注入功率和节点功率方程交叉验证。Vc V .* exp(1j*theta); I Y * Vc; S_check Vc .* conj(I); fprintf(平衡节点注入功率: %.6f j%.6f\n, real(S_check(slack)), imag(S_check(slack))); fprintf(PV节点无功输出: %.6f\n, imag(S_check(3)));S_check是基于Y矩阵和最终电压直接计算得到的注入复功率如果P、Q计算和雅可比公式里没有符号错误它应基本等于收敛后的P 1j*Q。另外可以改变初始电压再做一次潮流看是否收敛到同一解。牛拉法可能收敛到不同的运行点但在这个小系统里重复试验能暴露初值相关的错误。5.2 用稀疏矩阵重写雅可比去掉双层循环的提速方式当系统规模超过几百个节点时稠密J_red的内存和求解时间会迅速膨胀。常见做法是在组装时直接用稀疏三元组。rows []; cols []; vals []; for i 1:n for j 1:n if Y(i,j) ~ 0 % 这里 val 按4.2中的公式计算 rows [rows; i]; cols [cols; j]; vals [vals; val]; end end end J_sp sparse(rows, cols, vals, 2*n, 2*n);实际工程可以在循环里用preallocate预分配数组避免反复扩容。matlab的\对稀疏矩阵会自动选择稀疏LU分解不需要手动指定算法。还有一个很实用的调试技巧当第一次迭代没有收敛而残差序列出现震荡时多半是初值太差或雅可比符号有问题。在确认正负号无误后可以给修正量加阻尼。alpha 0.8; theta theta alpha * dtheta; V V V .* (alpha * dV_scaled);alpha从 1 开始残差增大时改成 0.5 甚至 0.2往往能把发散拖回来。这个阻尼思路和计算机视觉里的Levenberg-Marquardt加正则的思路类似对牛拉法潮流初值敏感的问题很有效。把小系统跑顺之后再换IEEE标准节点数据前也值得多看一眼收敛过程中的norm(dF,inf)是否单调下降。本文还有配套的精品资源点击获取