简介这套基于MATLAB实现的改进牛顿拉夫逊法程序面向电力系统研究人员、电气相关专业学生及配电网仿真工程师专门解决三相不平衡配电网在不同节点规模下的潮流计算问题。算法在传统牛顿拉夫逊法基础上做了改进收敛性与适应性更好适合各类节点情况可直接替换数据进行迁移使用。压缩包共22个文件以12个m脚本文件为主包含主函数、节点支路数据文件及核心迭代求解模块另附9个自动备份文件与1份使用说明文档整体仅28KB短小精悍便于快速加载运行。目前已有124人学习下载适合搭建三相不平衡潮流模型、对比算法效果或开展配电网仿真研究的读者。配套说明文档对代码结构和运行流程作了梳理遇到问题还可联系作者获得运行咨询、期刊复现、算法定制等后续支持可直接用于课程设计或课题预研。1. 三相不平衡配电网的潮流计算为什么牛顿拉夫逊法需要“改进”用单相潮流程序去算含单相光伏、单相负荷的三相配电网节点电压往往能落在合理区间但支路电流和功率损耗会偏得离谱。这不是牛顿法本身失效而是模型把三相系统简化成了一相相间互感、单相接入、中性点位移这些因素被直接抹掉了。改进牛顿拉夫逊法在原框架上把节点导纳矩阵扩成 3N×3N并对每一相单独写功率失配方程和雅可比偏导才能在 MATLAB 里直接处理 6、12、36 这类多节点算例。实际使用中这套代码适合两类人一是做配电网三相潮流、分布式电源接入方向的研究生拿 6 节点算例跑通流程再去换自己的网络二是做电网仿真或系统规划的工程师需要一个能嵌入主程序的快速三相潮流求解器。关键不在迭代公式多抽象而在如何把三相网络原始数据整理成可计算的矩阵结构以及在收敛失败时知道往哪个方向调。2. 三相潮流建模与节点导纳矩阵构造从单相到 3N 阶的扩展2.1 为什么三相导纳矩阵不是 N 阶而是 3N 阶配电网与输电网的显著区别是线路短、R/X 比高而且大量用户是单相接入。即使变压器出口电压对称由于三相负荷不均衡中性点电位会偏移各相电流和功率不再独立。单相潮流程序只保留一个节点电压相量相当于把三相电压折叠到一起自然丢失了相间耦合信息。因此三相潮流要做的第一件事就是把网络拓扑扩展成每个节点对应 A、B、C 三根相线节点导纳矩阵 Y 的维度变成 3N×3N其中 N 是拓扑节点数。在这个维度下每个节点对应一个 3×3 子块。对角线子块是该节点三相的自导纳非对角线子块是节点之间的相间互导纳。即使某节点只带了单相负荷其他两相的自导纳和互导纳也必须保留因为它们会通过线路和变压器耦合影响该节点其余相的电压。这就是原包中originalData6.m到originalData36.m都要把支路按相别展开的原因数据量不是 N 条支路而是按相展开后的 3 倍支路记录。从工程实现角度看3N 阶矩阵不要直接写成zeros(3*N)尤其在 36 节点时是 108×108 的稠密矩阵后续雅可比矩阵还会进一步膨胀。更合理的方式是用sparse预分配再按支路循环累加这样存储和求解速度都会好很多。2.2 从 originalData6.m 读取支路参数并生成三相导纳矩阵拿到originalData6.m后第一件事不是看牛顿迭代而是看它的矩阵格式。这类文件通常保存的是支路起始节点、终止节点、相别标志、电阻、电抗、对地导纳。常见列布局如下列号含义示例1起始节点编号12终止节点编号23相别标志0 表示三相1 表示 A 相2 表示 B 相3 表示 C 相4每相电阻 R / Ω0.0025每相电抗 X / Ω0.0046对地导纳 B / S0originalData12.m和originalData36.m只是节点数不同列格式应当保持一致。开始实现前先用size(originalData6)确认行数再打印前几行这比直接跑主函数更能避免索引越界。基于这种格式构造三相导纳矩阵的典型代码如下function Y build_Y3(br_data) N max(max(br_data(:,1:2))); Y sparse(3*N, 3*N); nbr size(br_data, 1); for k 1:nbr n1 br_data(k,1); n2 br_data(k,2); phase_flag br_data(k,3); R br_data(k,4); X br_data(k,5); B br_data(k,6); z R 1i*X; y 1/z; ysh 1i*B/2; [idx1, idx2] get_phase_idx(n1, n2, phase_flag, N); for m 1:length(idx1) for n 1:length(idx2) Y(idx1(m), idx2(n)) Y(idx1(m), idx2(n)) - y; end end for m 1:length(idx1) Y(idx1(m), idx1(m)) Y(idx1(m), idx1(m)) y ysh; end for m 1:length(idx2) Y(idx2(m), idx2(m)) Y(idx2(m), idx2(m)) y ysh; end end end function [idxA, idxB] get_phase_idx(n1, n2, flag, N) switch flag case 0 base1 (n1-1)*3; base2 (n2-1)*3; idxA base1 (1:3); idxB base2 (1:3); otherwise idxA (n1-1)*3 flag; idxB (n2-1)*3 flag; end end这段代码里sparse初始化后通过索引累加会自动把零元素变成非零最终仍保持稀疏存储。get_phase_idx把节点编号和相别映射到全局索引例如节点 i 的 A 相索引是(i-1)*31B 相是(i-1)*32C 相是(i-1)*33。循环中先减去互导纳再加到自导纳这个顺序不能调换。若有多条支路并联后累加的支路会在同一位置叠加最终得到正确的等效导纳。2.3 改进牛顿法对雅可比矩阵提出了什么新要求传统单相牛顿拉夫逊法在输电网中表现很好因为输电线路 R 远小于 X有功和无功耦合较弱可以采用 P-Q 分解简化。但配电网 R/X 比高三相负荷不平衡又引入了额外的相间耦合快速分解法的近似条件不再满足。所以原包中的New_NR_method.m走的是极坐标牛顿法路线但每个节点上的三相功率失配量必须同时参与迭代。改进的关键点之一是计算雅可比矩阵时把同节点不同相之间的偏导也纳入。比如节点 i 的 A 相有功失配量对节点 j 的 B 相电压相角的偏导不能近似忽略而要完整考虑互导纳实部和虚部组合出来的三角项。这个改动带来的收益是当某一相电压偏移很大时迭代仍能保持接近二阶收敛的特性。后续代码中出现的qiu_PQ.m就是功率失配量计算模块输入三相电压相量、三相导纳矩阵和负荷功率输出每一相的 P、Q 不平衡量作为牛顿迭代的右端项。3. 改进牛顿拉夫逊法的迭代核心失配量、雅可比与阻尼控制3.1 三相功率失配方程与 qiu_PQ.m 设计令节点 i 的 A、B、C 三相电压相量分别为 (V_i^a)、(V_i^b)、(V_i^c)写成复数形式。三相网络的节点注入功率在复平面内可以表示为[ S_i^p V_i^p \cdot \mathrm{conj}\left(\sum_{q \in {a,b,c}} \sum_{j1}^N Y_{ij}^{pq} V_j^q\right) ]其中 p 是本侧相别q 是对侧相别(Y_{ij}^{pq}) 表示节点 i 的 p 相与节点 j 的 q 相之间的互导纳。把计算功率与给定负荷相减就得到失配量[ \Delta P_i^p P_{i,\mathrm{given}}^p - \mathrm{Re}(S_i^p) ] [ \Delta Q_i^p Q_{i,\mathrm{given}}^p - \mathrm{Im}(S_i^p) ]qiu_PQ.m的任务就是批量完成上述计算。我更喜欢用向量化写法而不是三层循环因为 36 节点三相网络有 108 个复数电压量循环容易拖慢速度function [dP, dQ] qiu_PQ(U, Y, S_load) N length(S_load) / 3; U reshape(U, 3, N).; I Y * U(:); I reshape(I, 3, N).; S_calc conj(I) .* U; S_load reshape(S_load, 3, N).; dS S_load - S_calc; dP real(dS(:)); dQ imag(dS(:)); end这里Y * U(:)借助稀疏矩阵一次完成三相全网络的电流累加随后按相别 reshape 回 N×3得到各相计算功率。S_load的顺序必须与节点编号、相别顺序严格对应否则某相负荷会被算到另一相上最终导致潮流怎么迭代都不收敛。使用这段代码时建议先在 6 节点算例上验证 dP、dQ 的量级正常情况应在 1e-3 到 1e-2 起步然后随迭代逐步下降而不是一开始就出现 1e6 这种异常值。3.2 雅可比矩阵组装同一支路上三相耦合项的处理把失配量对电压幅值和相角求偏导得到的雅可比矩阵 J 可以分成四个子块H 对应有功对相角N 对应有功对电压幅值M 对应无功对相角L 对应无功对电压幅值。在三相网络中每个子块还要按相别展开因此 H、N、M、L 都是 3N×3N 的稀疏矩阵。实际组装时常见做法是先在一个两层循环内对节点对 i、j 处理再在内部对相别 p、q 做 9 次遍历。由于同一对节点在三相网络中可能出现 9 组偏导数我一般把 H、N、M、L 定义为稀疏矩阵并提前分配结构再逐个填入function J build_jacobian(U, Y, N) Ud abs(U); Ua angle(U); H sparse(3*N, 3*N); Nmat sparse(3*N, 3*N); Mmat sparse(3*N, 3*N); Lmat sparse(3*N, 3*N); for i 1:N for j 1:N for p 1:3 for q 1:3 [H, Nmat, Mmat, Lmat] fill_block_ij(... H, Nmat, Mmat, Lmat, i, j, p, q, Ud, Ua, Y); end end end end J [H Nmat; Mmat Lmat]; end四重循环在 36 节点下大约有 36×36×9≈11664 次块填充单次迭代在 MATLAB 2020b 下可接受。若网络到上百节点建议把节点按连接关系分组只对非零导纳块填充。fill_block_ij内部需要判断 i 和 j 是否相同再决定使用对角公式还是非对角公式。对角元素中存在“本相自导纳”和“他相互导纳”的混合公式推导时最容易漏掉的就是交叉项。这里有一个容易踩的坑电压幅值在牛顿法中不取摸长绝对值而是直接使用归一化后的标幺值参与除法。如果某相电压跌到接近 0.1 pu雅可比矩阵中的某些元素会变得很大可能导致迭代步长异常。遇到这种情况先检查负荷是否单相过重而不是急着调收敛阈值。3.3 迭代循环与阻尼因子结合主迭代流程在New_NR_method.m中循环结构大致如下% 平启动初值 U0 ones(3*N, 1); U U0; tol 1e-6; max_iter 50; for iter 1:max_iter [dP, dQ] qiu_PQ(U, Y, S_load); F [dP; dQ]; if norm(F, inf) tol break; end J build_jacobian(U, Y, N); dx J \ (-F); dU dx(1:3*N); dTheta dx(3*N1:end); % 幅值和相角分别更新 Ud_new abs(U) dU; Ua_new angle(U) dTheta; % 阻尼修正先试探再缩小步长 alpha 1.0; while alpha 0.05 U_test Ud_new .* exp(1i*Ua_new); [dP2, dQ2] qiu_PQ(U_test, Y, S_load); if norm([dP2; dQ2], inf) norm(F, inf) break; end alpha alpha * 0.5; Ua_new angle(U) alpha * dTheta; Ud_new abs(U) alpha * dU; end U Ud_new .* exp(1i*Ua_new); end阻尼牛顿法在三相不平衡严重时非常有效。假如某相电压相角跨过 180°计算功率会产生跳变直接更新会导致下一步失配量不减反增。这里的 while 循环保证每次迭代后的失配量单调下降alpha从 1 开始搜索范围到 0.05 为止。收敛判据不能只看有功失配量norm(F, inf)同时包含了有功和无功比单独检查 P 或 Q 更稳妥。配电网工程上 1e-4 已够用但做算法对比时我通常固定为 1e-6。4. 从 6/12/36 节点算例到换入新网络数据文件结构与替换方法4.1 三个算例数据文件的实际含义与常见误用压缩包里的originalData6.m、originalData12.m、originalData36.m分别对应 6 节点、12 节点、36 节点三相配电网。文件名中的数字就是拓扑节点个数。以 6 节点为例脚本中定义的是支路参数矩阵、负荷矩阵和平衡节点编号。运行main.m之前主程序会通过run指令将这些变量加载到工作区。.asv文件是 MATLAB 自动保存的历史版本例如originalData12.asv只是编辑过程中留下的备份运行主程序时优先选择.m而不是同名.asv否则数据不同步会导致结果对不上。这个细节在交付代码时很容易被忽略但不影响实际使用只需注意别把.asv当成另一组算例。替换数据前应核对以下内容核对项检查内容常见错误节点编号是否从 1 开始连续编号支路中出现孤立节点相别标志三相用 0 还是 3把 0 当成故障相负荷单位是三相总负荷还是单相负荷把三相功率直接填到单相上R/X 单位欧姆还是标幺值混用导致导纳数量级偏差平衡节点数是否只有一个多个平衡节点使雅可比奇异originalData36.m中的负荷矩阵行序一般与节点编号一致即第 i 行代表节点 i 的 A/B/C 三相负荷。如果从 IEEE 33 节点改造为三相网络原单相负荷数据需要拆成三相对称负荷再加入单相光伏或单相负荷的不平衡分量。这个处理方式在使用说明文档中有提醒实际执行时建议先做一次功率之和校验确认三相总功率等于所有负荷相加再进入潮流计算。4.2 修改 main.m 的启动方式避免反复切换算例使用说明文档描述的操作步骤是将所有文件放到 MATLAB 当前文件夹打开main.m点击运行。脚本形式的数据文件好处是直观缺点是多次运行后工作区会残留旧变量。我一般会在main.m开头做成可选配置clear; clc; % 选择算例originalData6 / originalData12 / originalData36 data_name originalData36; run(data_name); % 统一变量名后面都用 branch、load_pq、slack_bus branch eval(data_name); % 若原文件变量名就是文件名可直接取上述代码中eval会把字符串转成变量访问不够优雅但足够直接。更安全的做法是让每个数据文件都统一赋值给branch这样main.m只保留一行run(originalData36.m)即可修改连锁最少。原包里的12jiedianzhiluguanlianjuzhen.m用来生成 12 节点的支路关联矩阵作用是在有环网或双回线时检查拓扑连接关系。配电网大多开环运行但如果算例包含联络开关闭合后的闭环场景这个函数可以帮助检查是否出现孤立节点或环网冲突。4.3 换入自己的网络时参数怎么改才能一次跑通换数据不只是替换 R/X 矩阵负荷节点类型和平衡节点设置也要对齐。三相配电网中平衡节点通常设在变压器低压侧母线作为电压参考点其余节点多为 PQ 节点给定三相有功和无功负荷。若算例包含分布式电源可能需要把部分节点改为 PV 节点但三相 PV 节点的处理比单相复杂因为它既要维持正序电压又要分配各相无功出力。运行参数建议按下表调整参数推荐值修改位置收敛精度 tol1e-6main.m 或 New_NR_method.m最大迭代次数 max_iter50New_NR_method.m电压初值1.0 0i三相平启动main.m 中 U0 赋值负荷基准容量1e6 W计算标幺值时统一除以基准值从 6 节点换到 36 节点时最常见的失败原因是迭代次数不足。36 节点网络规模增大潮流方程非线性明显增强如果电压初值给成 0 或 1.2j0.5前三轮失配量会非常大达到 1e10 也不奇怪。更好的方法是用平启动电压即所有节点取 1.0∠0°让牛顿法自己去修正。若初始相角随机而不是从 0 开始某些支路会进入相位模糊区收敛速度显著下降。5. 收敛性判断与数值陷阱初值、阻尼与 MATLAB 版本差异5.1 电压初值不能随意给平启动是最稳的选择三相配电网中平衡节点电压通常设为 1.0∠0°其余节点从 1.0∠0° 开始迭代。有些刚接触牛顿法的人会把初值设成 0.95j0.1 来模拟重负荷这其实会破坏雅可比矩阵的初始条件。原因是幅值修正量dU在迭代中要加到幅值数组上如果初值偏离 1.0 太多第一步的功率失配量可能超过 100%阻尼因子必须反复缩小迭代次数成倍增加。我通常只在两个位置允许非平启动一是已经知道某个节点带有大型感应电动机需要给定 0.98∠-2° 作为热启动初值二是做连续潮流或时域仿真时把上一时刻的解直接作为下一时刻初值。静态潮流计算中平启动永远是第一选择。5.2 失配量不下降时先查三相导纳矩阵再查负荷如果发现norm(F, inf)从第一轮开始就振荡优先检查sparse矩阵的结构而不是调阻尼。可以用spy(Y)查看非零元分布正常三相网络应该有三个对角带且每个支路对应 9 个耦合元素。若出现整行全零说明某个节点编号在支路数据中不存在导致孤立节点出现。另一种常见情况是负荷数据中混入了 NaN 或复数虚部。qiu_PQ.m里S_calc一旦包含 NaNnorm(F, inf)也会变成 NaN。在main.m中加一行断言即可快速定位assert(~any(isnan(S_load(:))), S_load contains NaN);运行到这一行时MATLAB 2020b 会直接中断并提示错误位置省去翻数据文件的时间。5.3 阻尼调节的边界什么时候该放弃二分搜索阻尼 while 循环从 alpha1 开始每次减半搜索最多约 5 轮。如果 alpha 已经到 0.05 失配量仍然不下降继续缩小意义不大因为此时步长过小迭代基本停滞。这种情况通常是模型本身出了问题比如有两个平衡节点、某相负荷为负值或者变压器变比没有折算进导纳矩阵。可以把max_iter暂时调大到 100并去掉阻尼观察失配量在无阻尼时的变化趋势能更快暴露问题。另外不同 MATLAB 版本对稀疏矩阵的\求解算法选择不同。R2020b 中J \ (-F)自动使用 UMFPACK但 J 是复数稀疏矩阵时内存占用会比实数矩阵高。36 节点三相网络规模不大通常没有性能瓶颈若换到上百节点建议把雅可比矩阵改为实部虚部分离的增广形式避免复数稀疏求解产生额外开销。这一步不用改算法逻辑只改数据装配但能明显降低内存峰值。5.4 一个直接可用的收敛性检查模板rho zeros(max_iter, 1); for iter 1:max_iter [dP, dQ] qiu_PQ(U, Y, S_load); rho(iter) norm([dP; dQ], inf); if rho(iter) tol break; end J build_jacobian(U, Y, N); dx J \ -[dP; dQ]; U update_with_damping(U, dx, Y, S_load); end % 输出每次失配量便于观察二阶收敛 disp(rho(1:iter));正常情况rho会按 1e-2 → 1e-4 → 1e-8 的节奏下降末段呈现接近二阶收敛的陡降。如果看到 1e-2 → 1e-3 → 1e-2 的波动优先回到 5.2 节检查矩阵和负荷。这组判断不仅在算例调试时有用当你想把三相牛顿法接入更大规模仿真时也能作为算法可靠性指标。本文还有配套的精品资源点击获取
三相不平衡配电网潮流计算的改进牛顿拉夫逊法及MATLAB实现
简介这套基于MATLAB实现的改进牛顿拉夫逊法程序面向电力系统研究人员、电气相关专业学生及配电网仿真工程师专门解决三相不平衡配电网在不同节点规模下的潮流计算问题。算法在传统牛顿拉夫逊法基础上做了改进收敛性与适应性更好适合各类节点情况可直接替换数据进行迁移使用。压缩包共22个文件以12个m脚本文件为主包含主函数、节点支路数据文件及核心迭代求解模块另附9个自动备份文件与1份使用说明文档整体仅28KB短小精悍便于快速加载运行。目前已有124人学习下载适合搭建三相不平衡潮流模型、对比算法效果或开展配电网仿真研究的读者。配套说明文档对代码结构和运行流程作了梳理遇到问题还可联系作者获得运行咨询、期刊复现、算法定制等后续支持可直接用于课程设计或课题预研。1. 三相不平衡配电网的潮流计算为什么牛顿拉夫逊法需要“改进”用单相潮流程序去算含单相光伏、单相负荷的三相配电网节点电压往往能落在合理区间但支路电流和功率损耗会偏得离谱。这不是牛顿法本身失效而是模型把三相系统简化成了一相相间互感、单相接入、中性点位移这些因素被直接抹掉了。改进牛顿拉夫逊法在原框架上把节点导纳矩阵扩成 3N×3N并对每一相单独写功率失配方程和雅可比偏导才能在 MATLAB 里直接处理 6、12、36 这类多节点算例。实际使用中这套代码适合两类人一是做配电网三相潮流、分布式电源接入方向的研究生拿 6 节点算例跑通流程再去换自己的网络二是做电网仿真或系统规划的工程师需要一个能嵌入主程序的快速三相潮流求解器。关键不在迭代公式多抽象而在如何把三相网络原始数据整理成可计算的矩阵结构以及在收敛失败时知道往哪个方向调。2. 三相潮流建模与节点导纳矩阵构造从单相到 3N 阶的扩展2.1 为什么三相导纳矩阵不是 N 阶而是 3N 阶配电网与输电网的显著区别是线路短、R/X 比高而且大量用户是单相接入。即使变压器出口电压对称由于三相负荷不均衡中性点电位会偏移各相电流和功率不再独立。单相潮流程序只保留一个节点电压相量相当于把三相电压折叠到一起自然丢失了相间耦合信息。因此三相潮流要做的第一件事就是把网络拓扑扩展成每个节点对应 A、B、C 三根相线节点导纳矩阵 Y 的维度变成 3N×3N其中 N 是拓扑节点数。在这个维度下每个节点对应一个 3×3 子块。对角线子块是该节点三相的自导纳非对角线子块是节点之间的相间互导纳。即使某节点只带了单相负荷其他两相的自导纳和互导纳也必须保留因为它们会通过线路和变压器耦合影响该节点其余相的电压。这就是原包中originalData6.m到originalData36.m都要把支路按相别展开的原因数据量不是 N 条支路而是按相展开后的 3 倍支路记录。从工程实现角度看3N 阶矩阵不要直接写成zeros(3*N)尤其在 36 节点时是 108×108 的稠密矩阵后续雅可比矩阵还会进一步膨胀。更合理的方式是用sparse预分配再按支路循环累加这样存储和求解速度都会好很多。2.2 从 originalData6.m 读取支路参数并生成三相导纳矩阵拿到originalData6.m后第一件事不是看牛顿迭代而是看它的矩阵格式。这类文件通常保存的是支路起始节点、终止节点、相别标志、电阻、电抗、对地导纳。常见列布局如下列号含义示例1起始节点编号12终止节点编号23相别标志0 表示三相1 表示 A 相2 表示 B 相3 表示 C 相4每相电阻 R / Ω0.0025每相电抗 X / Ω0.0046对地导纳 B / S0originalData12.m和originalData36.m只是节点数不同列格式应当保持一致。开始实现前先用size(originalData6)确认行数再打印前几行这比直接跑主函数更能避免索引越界。基于这种格式构造三相导纳矩阵的典型代码如下function Y build_Y3(br_data) N max(max(br_data(:,1:2))); Y sparse(3*N, 3*N); nbr size(br_data, 1); for k 1:nbr n1 br_data(k,1); n2 br_data(k,2); phase_flag br_data(k,3); R br_data(k,4); X br_data(k,5); B br_data(k,6); z R 1i*X; y 1/z; ysh 1i*B/2; [idx1, idx2] get_phase_idx(n1, n2, phase_flag, N); for m 1:length(idx1) for n 1:length(idx2) Y(idx1(m), idx2(n)) Y(idx1(m), idx2(n)) - y; end end for m 1:length(idx1) Y(idx1(m), idx1(m)) Y(idx1(m), idx1(m)) y ysh; end for m 1:length(idx2) Y(idx2(m), idx2(m)) Y(idx2(m), idx2(m)) y ysh; end end end function [idxA, idxB] get_phase_idx(n1, n2, flag, N) switch flag case 0 base1 (n1-1)*3; base2 (n2-1)*3; idxA base1 (1:3); idxB base2 (1:3); otherwise idxA (n1-1)*3 flag; idxB (n2-1)*3 flag; end end这段代码里sparse初始化后通过索引累加会自动把零元素变成非零最终仍保持稀疏存储。get_phase_idx把节点编号和相别映射到全局索引例如节点 i 的 A 相索引是(i-1)*31B 相是(i-1)*32C 相是(i-1)*33。循环中先减去互导纳再加到自导纳这个顺序不能调换。若有多条支路并联后累加的支路会在同一位置叠加最终得到正确的等效导纳。2.3 改进牛顿法对雅可比矩阵提出了什么新要求传统单相牛顿拉夫逊法在输电网中表现很好因为输电线路 R 远小于 X有功和无功耦合较弱可以采用 P-Q 分解简化。但配电网 R/X 比高三相负荷不平衡又引入了额外的相间耦合快速分解法的近似条件不再满足。所以原包中的New_NR_method.m走的是极坐标牛顿法路线但每个节点上的三相功率失配量必须同时参与迭代。改进的关键点之一是计算雅可比矩阵时把同节点不同相之间的偏导也纳入。比如节点 i 的 A 相有功失配量对节点 j 的 B 相电压相角的偏导不能近似忽略而要完整考虑互导纳实部和虚部组合出来的三角项。这个改动带来的收益是当某一相电压偏移很大时迭代仍能保持接近二阶收敛的特性。后续代码中出现的qiu_PQ.m就是功率失配量计算模块输入三相电压相量、三相导纳矩阵和负荷功率输出每一相的 P、Q 不平衡量作为牛顿迭代的右端项。3. 改进牛顿拉夫逊法的迭代核心失配量、雅可比与阻尼控制3.1 三相功率失配方程与 qiu_PQ.m 设计令节点 i 的 A、B、C 三相电压相量分别为 (V_i^a)、(V_i^b)、(V_i^c)写成复数形式。三相网络的节点注入功率在复平面内可以表示为[ S_i^p V_i^p \cdot \mathrm{conj}\left(\sum_{q \in {a,b,c}} \sum_{j1}^N Y_{ij}^{pq} V_j^q\right) ]其中 p 是本侧相别q 是对侧相别(Y_{ij}^{pq}) 表示节点 i 的 p 相与节点 j 的 q 相之间的互导纳。把计算功率与给定负荷相减就得到失配量[ \Delta P_i^p P_{i,\mathrm{given}}^p - \mathrm{Re}(S_i^p) ] [ \Delta Q_i^p Q_{i,\mathrm{given}}^p - \mathrm{Im}(S_i^p) ]qiu_PQ.m的任务就是批量完成上述计算。我更喜欢用向量化写法而不是三层循环因为 36 节点三相网络有 108 个复数电压量循环容易拖慢速度function [dP, dQ] qiu_PQ(U, Y, S_load) N length(S_load) / 3; U reshape(U, 3, N).; I Y * U(:); I reshape(I, 3, N).; S_calc conj(I) .* U; S_load reshape(S_load, 3, N).; dS S_load - S_calc; dP real(dS(:)); dQ imag(dS(:)); end这里Y * U(:)借助稀疏矩阵一次完成三相全网络的电流累加随后按相别 reshape 回 N×3得到各相计算功率。S_load的顺序必须与节点编号、相别顺序严格对应否则某相负荷会被算到另一相上最终导致潮流怎么迭代都不收敛。使用这段代码时建议先在 6 节点算例上验证 dP、dQ 的量级正常情况应在 1e-3 到 1e-2 起步然后随迭代逐步下降而不是一开始就出现 1e6 这种异常值。3.2 雅可比矩阵组装同一支路上三相耦合项的处理把失配量对电压幅值和相角求偏导得到的雅可比矩阵 J 可以分成四个子块H 对应有功对相角N 对应有功对电压幅值M 对应无功对相角L 对应无功对电压幅值。在三相网络中每个子块还要按相别展开因此 H、N、M、L 都是 3N×3N 的稀疏矩阵。实际组装时常见做法是先在一个两层循环内对节点对 i、j 处理再在内部对相别 p、q 做 9 次遍历。由于同一对节点在三相网络中可能出现 9 组偏导数我一般把 H、N、M、L 定义为稀疏矩阵并提前分配结构再逐个填入function J build_jacobian(U, Y, N) Ud abs(U); Ua angle(U); H sparse(3*N, 3*N); Nmat sparse(3*N, 3*N); Mmat sparse(3*N, 3*N); Lmat sparse(3*N, 3*N); for i 1:N for j 1:N for p 1:3 for q 1:3 [H, Nmat, Mmat, Lmat] fill_block_ij(... H, Nmat, Mmat, Lmat, i, j, p, q, Ud, Ua, Y); end end end end J [H Nmat; Mmat Lmat]; end四重循环在 36 节点下大约有 36×36×9≈11664 次块填充单次迭代在 MATLAB 2020b 下可接受。若网络到上百节点建议把节点按连接关系分组只对非零导纳块填充。fill_block_ij内部需要判断 i 和 j 是否相同再决定使用对角公式还是非对角公式。对角元素中存在“本相自导纳”和“他相互导纳”的混合公式推导时最容易漏掉的就是交叉项。这里有一个容易踩的坑电压幅值在牛顿法中不取摸长绝对值而是直接使用归一化后的标幺值参与除法。如果某相电压跌到接近 0.1 pu雅可比矩阵中的某些元素会变得很大可能导致迭代步长异常。遇到这种情况先检查负荷是否单相过重而不是急着调收敛阈值。3.3 迭代循环与阻尼因子结合主迭代流程在New_NR_method.m中循环结构大致如下% 平启动初值 U0 ones(3*N, 1); U U0; tol 1e-6; max_iter 50; for iter 1:max_iter [dP, dQ] qiu_PQ(U, Y, S_load); F [dP; dQ]; if norm(F, inf) tol break; end J build_jacobian(U, Y, N); dx J \ (-F); dU dx(1:3*N); dTheta dx(3*N1:end); % 幅值和相角分别更新 Ud_new abs(U) dU; Ua_new angle(U) dTheta; % 阻尼修正先试探再缩小步长 alpha 1.0; while alpha 0.05 U_test Ud_new .* exp(1i*Ua_new); [dP2, dQ2] qiu_PQ(U_test, Y, S_load); if norm([dP2; dQ2], inf) norm(F, inf) break; end alpha alpha * 0.5; Ua_new angle(U) alpha * dTheta; Ud_new abs(U) alpha * dU; end U Ud_new .* exp(1i*Ua_new); end阻尼牛顿法在三相不平衡严重时非常有效。假如某相电压相角跨过 180°计算功率会产生跳变直接更新会导致下一步失配量不减反增。这里的 while 循环保证每次迭代后的失配量单调下降alpha从 1 开始搜索范围到 0.05 为止。收敛判据不能只看有功失配量norm(F, inf)同时包含了有功和无功比单独检查 P 或 Q 更稳妥。配电网工程上 1e-4 已够用但做算法对比时我通常固定为 1e-6。4. 从 6/12/36 节点算例到换入新网络数据文件结构与替换方法4.1 三个算例数据文件的实际含义与常见误用压缩包里的originalData6.m、originalData12.m、originalData36.m分别对应 6 节点、12 节点、36 节点三相配电网。文件名中的数字就是拓扑节点个数。以 6 节点为例脚本中定义的是支路参数矩阵、负荷矩阵和平衡节点编号。运行main.m之前主程序会通过run指令将这些变量加载到工作区。.asv文件是 MATLAB 自动保存的历史版本例如originalData12.asv只是编辑过程中留下的备份运行主程序时优先选择.m而不是同名.asv否则数据不同步会导致结果对不上。这个细节在交付代码时很容易被忽略但不影响实际使用只需注意别把.asv当成另一组算例。替换数据前应核对以下内容核对项检查内容常见错误节点编号是否从 1 开始连续编号支路中出现孤立节点相别标志三相用 0 还是 3把 0 当成故障相负荷单位是三相总负荷还是单相负荷把三相功率直接填到单相上R/X 单位欧姆还是标幺值混用导致导纳数量级偏差平衡节点数是否只有一个多个平衡节点使雅可比奇异originalData36.m中的负荷矩阵行序一般与节点编号一致即第 i 行代表节点 i 的 A/B/C 三相负荷。如果从 IEEE 33 节点改造为三相网络原单相负荷数据需要拆成三相对称负荷再加入单相光伏或单相负荷的不平衡分量。这个处理方式在使用说明文档中有提醒实际执行时建议先做一次功率之和校验确认三相总功率等于所有负荷相加再进入潮流计算。4.2 修改 main.m 的启动方式避免反复切换算例使用说明文档描述的操作步骤是将所有文件放到 MATLAB 当前文件夹打开main.m点击运行。脚本形式的数据文件好处是直观缺点是多次运行后工作区会残留旧变量。我一般会在main.m开头做成可选配置clear; clc; % 选择算例originalData6 / originalData12 / originalData36 data_name originalData36; run(data_name); % 统一变量名后面都用 branch、load_pq、slack_bus branch eval(data_name); % 若原文件变量名就是文件名可直接取上述代码中eval会把字符串转成变量访问不够优雅但足够直接。更安全的做法是让每个数据文件都统一赋值给branch这样main.m只保留一行run(originalData36.m)即可修改连锁最少。原包里的12jiedianzhiluguanlianjuzhen.m用来生成 12 节点的支路关联矩阵作用是在有环网或双回线时检查拓扑连接关系。配电网大多开环运行但如果算例包含联络开关闭合后的闭环场景这个函数可以帮助检查是否出现孤立节点或环网冲突。4.3 换入自己的网络时参数怎么改才能一次跑通换数据不只是替换 R/X 矩阵负荷节点类型和平衡节点设置也要对齐。三相配电网中平衡节点通常设在变压器低压侧母线作为电压参考点其余节点多为 PQ 节点给定三相有功和无功负荷。若算例包含分布式电源可能需要把部分节点改为 PV 节点但三相 PV 节点的处理比单相复杂因为它既要维持正序电压又要分配各相无功出力。运行参数建议按下表调整参数推荐值修改位置收敛精度 tol1e-6main.m 或 New_NR_method.m最大迭代次数 max_iter50New_NR_method.m电压初值1.0 0i三相平启动main.m 中 U0 赋值负荷基准容量1e6 W计算标幺值时统一除以基准值从 6 节点换到 36 节点时最常见的失败原因是迭代次数不足。36 节点网络规模增大潮流方程非线性明显增强如果电压初值给成 0 或 1.2j0.5前三轮失配量会非常大达到 1e10 也不奇怪。更好的方法是用平启动电压即所有节点取 1.0∠0°让牛顿法自己去修正。若初始相角随机而不是从 0 开始某些支路会进入相位模糊区收敛速度显著下降。5. 收敛性判断与数值陷阱初值、阻尼与 MATLAB 版本差异5.1 电压初值不能随意给平启动是最稳的选择三相配电网中平衡节点电压通常设为 1.0∠0°其余节点从 1.0∠0° 开始迭代。有些刚接触牛顿法的人会把初值设成 0.95j0.1 来模拟重负荷这其实会破坏雅可比矩阵的初始条件。原因是幅值修正量dU在迭代中要加到幅值数组上如果初值偏离 1.0 太多第一步的功率失配量可能超过 100%阻尼因子必须反复缩小迭代次数成倍增加。我通常只在两个位置允许非平启动一是已经知道某个节点带有大型感应电动机需要给定 0.98∠-2° 作为热启动初值二是做连续潮流或时域仿真时把上一时刻的解直接作为下一时刻初值。静态潮流计算中平启动永远是第一选择。5.2 失配量不下降时先查三相导纳矩阵再查负荷如果发现norm(F, inf)从第一轮开始就振荡优先检查sparse矩阵的结构而不是调阻尼。可以用spy(Y)查看非零元分布正常三相网络应该有三个对角带且每个支路对应 9 个耦合元素。若出现整行全零说明某个节点编号在支路数据中不存在导致孤立节点出现。另一种常见情况是负荷数据中混入了 NaN 或复数虚部。qiu_PQ.m里S_calc一旦包含 NaNnorm(F, inf)也会变成 NaN。在main.m中加一行断言即可快速定位assert(~any(isnan(S_load(:))), S_load contains NaN);运行到这一行时MATLAB 2020b 会直接中断并提示错误位置省去翻数据文件的时间。5.3 阻尼调节的边界什么时候该放弃二分搜索阻尼 while 循环从 alpha1 开始每次减半搜索最多约 5 轮。如果 alpha 已经到 0.05 失配量仍然不下降继续缩小意义不大因为此时步长过小迭代基本停滞。这种情况通常是模型本身出了问题比如有两个平衡节点、某相负荷为负值或者变压器变比没有折算进导纳矩阵。可以把max_iter暂时调大到 100并去掉阻尼观察失配量在无阻尼时的变化趋势能更快暴露问题。另外不同 MATLAB 版本对稀疏矩阵的\求解算法选择不同。R2020b 中J \ (-F)自动使用 UMFPACK但 J 是复数稀疏矩阵时内存占用会比实数矩阵高。36 节点三相网络规模不大通常没有性能瓶颈若换到上百节点建议把雅可比矩阵改为实部虚部分离的增广形式避免复数稀疏求解产生额外开销。这一步不用改算法逻辑只改数据装配但能明显降低内存峰值。5.4 一个直接可用的收敛性检查模板rho zeros(max_iter, 1); for iter 1:max_iter [dP, dQ] qiu_PQ(U, Y, S_load); rho(iter) norm([dP; dQ], inf); if rho(iter) tol break; end J build_jacobian(U, Y, N); dx J \ -[dP; dQ]; U update_with_damping(U, dx, Y, S_load); end % 输出每次失配量便于观察二阶收敛 disp(rho(1:iter));正常情况rho会按 1e-2 → 1e-4 → 1e-8 的节奏下降末段呈现接近二阶收敛的陡降。如果看到 1e-2 → 1e-3 → 1e-2 的波动优先回到 5.2 节检查矩阵和负荷。这组判断不仅在算例调试时有用当你想把三相牛顿法接入更大规模仿真时也能作为算法可靠性指标。本文还有配套的精品资源点击获取
