简介：面向电力系统专业学生、研究人员与调度工程师，这份基于Matlab的潮流计算程序包围绕IEEE 33节点与69节点标准配电系统，解决电网电压、电流、功率分布的计算与分析问题，可作为配电网稳态分析、分布式电源接入评估及算法研究的基准工具。包内涵盖5个文件，其中4个m文件为可直接运行的Matlab源码，已按33节点、69节点分别实现潮流计算，并包含前推回代法程序；另附1个PDF文档，用于说明模型参数与程序调用方式，压缩包仅341KB，部署轻量、便于二次开发。目前已有459人学习下载，适用于本科及研究生教学实验、课程设计与科研验证。通过运行程序，可直观得到各节点电压幅值、相角及支路功率等关键参数，支持牛顿-拉夫森法、前推回代法等多种算法的对比实践，还可在同一框架下扩展分布式电源接入场景，帮助读者系统掌握IEEE标准配电系统建模思路与Matlab数值计算编程方法。

1. 为什么这个 IEEE33 节点潮流计算资源值得下载

做配电网规划时，上级只给一张单线图和两行负荷说明，却要求半小时内评估分布式光伏接入后的电压分布。这种情况下，我不会从零开始搭模型，而是直接打开 IEEE33 节点标准算例跑一次潮流计算。这个资源的价值就在这里：它把 IEEE33 节点和 IEEE69 节点的支路参数、负荷参数、拓扑连接整理成可直接使用的数据文件，配套前推回代潮流计算主程序，覆盖 loadflow 的全流程。适合正在写配电网论文的研究生，也适合刚接手配网规划的工程师。不需要重复录入数据，不需要自己踩一遍拓扑编号的坑，唯一要做的就是理解算法然后用起来。

2. 数据准备：把 IEEE33/69 节点参数变成可计算的数据

2.1 33 节点与 69 节点的标准拓扑差异

IEEE33 节点测试系统和 IEEE69 节点测试系统是配电网研究里引用率很高的两个基准模型。IEEE33 节点系统额定电压 12.66kV，包含 32 条常闭支路和 5 个联络开关；IEEE69 节点系统同样是 12.66kV，但支路数量明显增加，包含 68 条常闭支路和 7 个联络开关。两个系统在默认状态下都保持辐射状开环运行，联络开关处于断开状态，因此非常适合用前推回代法求解。

注意，不同文献给出的总负荷数据会有细微差异，例如 69 节点系统总有功负荷有写 3802.19 kW，也有写 3802.2 kW，差别不大但对照文献时要注意来源。这个资源里的 case33 和 case69 数据，总负荷值应该和上述数据基本一致。若差得太多，先检查是不是把联络开关支路也算进负荷了。

2.2 用 MATLAB 格式整理支路和负荷数据

拿到这份资源后，不要急着跑程序，先看数据结构。常见做法是把每个测试系统单独放一个数据文件，里面用矩阵保存三类信息：节点编号、支路阻抗、节点负荷。下面是我习惯用的 MATLAB 格式，这套格式也符合资源里多数脚本的输入约定。

% case33.m —— IEEE33 节点配电网数据 % 支路矩阵：行 = 支路，列 = [首端节点, 末端节点, R/Ohm, X/Ohm] % 这里列出前 3 条支路，完整 32 条在资源内 case33.m 中给出 branch = [ 1 2 0.0922 0.0470; 2 3 0.4930 0.2511; 3 4 0.3660 0.1864; % 4 5 ... 省略一大部分，下载后打开文件可看到全部数据 ]; % 负荷矩阵：行 = 节点，列 = [节点号, P/kW, Q/kvar] load = [ 1 0 0; 2 100 60; 3 90 40; % 4 5 ... 同样省略，完整数据在资源内 ]; % 根节点编号，用于后面前推回代排序 root = 1;

这段代码里支路阻抗单位是欧姆，负荷单位是千瓦和千乏，不是标幺值。为什么这么设计？因为配网数据习惯直接给有名值，用有名值计算也便于和单线图上的设备参数对应。根节点就是变电站出口母线，IEEE33 节点系统一般把母线编号为 1，但部分文献把根节点写成 0，拿到数据后先用一个校核脚本确认。

% 校核数据连续性 N = max(max(branch(:,1:2)), max(load(:,1))); assert(min(branch(:)) >= 1, '节点编号必须大于等于 1'); P_total = sum(load(:,2)); Q_total = sum(load(:,3)); fprintf('总负荷 P=%.1f kW, Q=%.1f kvar\n', P_total, Q_total);

这段校核脚本主要做两件事：检查节点编号是否从 1 开始，统计总负荷是否落在合理区间。如果min(branch(:))返回 0，说明数据里的根节点是 0 号母线，后面所有节点编号都要整体加 1。别小看这一步，很多人在数据导入阶段就翻车，算出来的电压全是乱的。

2.3 69 节点数据入口

69 节点数据在资源里通常单独命名，比如 case69.m，结构和 case33 完全一致。替换方法很简单：加载case69.m后，把branch和load换成 69 节点系统的矩阵，根节点依然是 1，主程序不用改。

% 加载 33 或 69 节点数据 case33; % 或者 case69 % case69 中已经定义了 branch, load, root

换数据后要格外留意联络开关状态。IEEE69 节点在默认算例中 7 个联络开关全部断开，前推回代才能保持辐射状。如果某个联络开关被误当作闭合，程序会报不收敛或者电压结果异常。关于这一点，后面避坑章节会专门展开。

3. 算法选择：前推回代法的完整推导与代码实现

3.1 为什么配电网潮流计算优先用前推回代法

配电网和输电网最明显的区别是线路阻抗中电阻占比高，R/X 比值可能大于 1，而输电网里 X 通常远大于 R。牛顿-拉夫逊法在处理输电网 PQ 分解时依赖高 X/R 假设，直接搬到配电网容易出现雅可比矩阵病态，收敛困难。前推回代法不需要形成雅可比矩阵，只依赖树木层次结构，从末端往前推电流、从根节点往后推电压，物理概念清晰，对辐射状配网收敛性也稳定。

前推回代法有一个前提：网络必须是无环的辐射状网络。IEEE33 和 IEEE69 在联络开关断开时就满足这个条件。如果闭合联络开关，网络出现合环，就必须在算法里额外加补偿电流迭代，或者改用其他方法求解。这一点是选型时就要想清楚的。

3.2 先构建树的父节点和遍历顺序

前推回代第一步是确定网络拓扑中的父子关系。常见做法是用广度优先搜索从根节点开始遍历，记录每个节点的父节点和支路阻抗。下面这段代码可以独立保存成build_tree.m。

function [parent, Zparent, order] = build_tree(branch, root) % branch: [from, to, R, X] N = max(max(branch(:,1:2))); % 邻接表记录相连节点和支路阻抗 adj = cell(N,1); for k = 1:size(branch,1) f = branch(k,1); t = branch(k,2); Z = branch(k,3) + 1j*branch(k,4); adj{f}(end+1,:) = [t, Z]; adj{t}(end+1,:) = [f, Z]; end % 广度优先搜索 parent = zeros(N,1); Zparent = zeros(N,1); order = zeros(N,1); visited = false(N,1); queue = root; visited(root) = true; cnt = 0; while ~isempty(queue) u = queue(1); queue(1) = []; cnt = cnt + 1; order(cnt) = u; for i = 1:size(adj{u},1) v = adj{u}(i,1); if ~visited(v) visited(v) = true; parent(v) = u; Zparent(v) = adj{u}(i,2); queue(end+1) = v; end end end if any(~visited) error('存在孤立节点或环路，请检查联络开关是否全部断开'); end end

这段代码返回三个关键变量：parent表示每个节点的父节点，Zparent表示该节点与父节点之间的线路阻抗，order是 BFS 遍历顺序。回代时要用order的逆序，因为逆序保证先处理末端节点，再向根节点累加电流。代码最后对孤立节点做了检查，如果visited里还有 false，说明网络不是一棵完整树，要么有断线，要么存在无法从根节点到达的环路。

3.3 前推回代主循环

构建完树之后，主迭代就变得很直接。每次迭代分三步：计算节点注入电流，回代累加支路电流，前推更新节点电压。完整代码放在run_pf.m中。

function [V, iter] = run_pf(branch, load, root, tol, max_iter) [parent, Zparent, order] = build_tree(branch, root); N = length(order); % 把负荷整理成复数功率矩阵 S = zeros(N,1); for i = 1:size(load,1) b = load(i,1); S(b) = load(i,2) + 1j*load(i,3); end V = ones(N,1); I_line = zeros(N,1); iter = 0; mismatch = 1e10; while mismatch > tol && iter < max_iter V_old = V; % 第 1 步：计算节点注入电流 I_node = zeros(N,1); for n = 2:N I_node(n) = conj(S(n) / V(n)); end % 第 2 步：回代，从末端向根节点累加支路电流 I_line = zeros(N,1); for i = N:-1:2 n = order(i); p = parent(n); I_line(n) = I_node(n); for c = 1:N if parent(c) == n I_line(n) = I_line(n) + I_line(c); end end end % 第 3 步：前推，从根向末端更新电压 for i = 2:N n = order(i); p = parent(n); V(n) = V(p) - Zparent(n) * I_line(n); end mismatch = max(abs(V - V_old)); iter = iter + 1; end if mismatch > tol warning('迭代未收敛，最终偏差 %g', mismatch); end end

这段代码里，S(n)是节点 n 的负荷复功率，conj(S(n)/V(n))是注入电流共轭。回代部分用两层循环寻找子节点，虽然复杂度是 O(N^2)，但 33 和 69 节点规模下运行时间可忽略。前推部分直接利用Zparent(n)计算电压降落，根节点电压始终保持 1.0 pu。tol一般设1e-6，max_iter设 100 到 200。配电网正常算例通常在 3 到 6 次迭代内收敛。

3.4 调用方式与结果输出

主程序封装好之后，调用非常简单。下面这段代码展示完整调用流程，同时把电压标幺值曲线画出来。

% 加载数据 case33; % 运行潮流 [V, iter] = run_pf(branch, load, root, 1e-6, 100); fprintf('迭代次数：%d\n', iter); % 输出最低电压节点和幅值 [Vmin, idx] = min(abs(V)); fprintf('最低电压：%.4f pu，节点：%d\n', Vmin, idx); % 画电压曲线 figure; plot(abs(V), 'o-'); xlabel('节点'); ylabel('电压/pu'); grid on;

这段代码里的abs(V)是电压幅值标幺值。IEEE33 节点标准算例的电压最低点通常在 18 号节点附近，幅值约 0.9038 pu。如果你跑出来的结果明显偏离，优先检查负荷矩阵和支路阻抗有没有录入错误，而不是怀疑代码。

4. 从 33 节点扩到 69 节点：改数据、改拓扑的相关代码改动

4.1 数据文件替换与初始检查

从 33 节点扩展到 69 节点，算法层面不需要改动，难的是数据替换后的校验。case69.m里的branch是 68 行矩阵，load是 68 行，根节点同样是 1。替换后先运行数据校核脚本，确认总负荷和节点连续性。

% 切换到 69 节点数据 case69; % 校验数据 N = max(max(branch(:,1:2))); assert(N == 69, '节点数量不是 69，请检查数据'); P_total = sum(load(:,2)); Q_total = sum(load(:,3)); fprintf('69 节点系统总负荷：P=%.1f kW, Q=%.1f kvar\n', P_total, Q_total);

69 节点系统比 33 节点多出一批分支，尤其是几条长线路末端节点，电压偏低现象更明显。跑完潮流后不要只看迭代次数，要画电压曲线确认最低电压点的位置。常见结果是最低电压出现在某条重载分支末端，如果画出来的曲线在某个中段节点附近出现异常凹陷，往往说明该节点负荷数据多录了一个数量级。

4.2 多分支节点的子节点表优化

前文提到回代部分用两层循环找子节点，33/69 节点时没问题，但如果你想扩展到几百节点的实际馈线，建议改成预先建立子节点表。改动很小，但能明显减少内层循环的无效遍历。

% 在 build_tree 函数中补充 children 表 children = cell(N,1); for n = 2:N p = parent(n); children{p}(end+1) = n; end

有了children表后，回代部分可以改成下面这样：

for i = N:-1:2 n = order(i); p = parent(n); I_line(n) = I_node(n); for c = children{n} I_line(n) = I_line(n) + I_line(c); end end

这个优化没有改变算法本质，只是把查找孩子节点的操作提前到迭代前完成。资源里如果直接用了for c = 1:N的写法，性能上没问题，但逻辑上不够友好。拆资源时我习惯顺手改成children表，后续加分布式电源或扩展网络时更好维护。

4.3 接入分布式光伏：把光伏当负负荷处理

配电网潮流最常见的扩展场景是接入分布式光伏。光伏逆变器在没有参与无功调节时，可以简单建模为恒功率负负荷，也就是在对应节点上把S(n)从负荷功率减去光伏输出功率。

% 在 run_pf 主循环前加入光伏修正 PV = [ 17 500 0; % 17 号节点光伏 500 kW，无功 0 18 300 0; % 18 号节点光伏 300 kW ]; for i = 1:size(PV,1) b = PV(i,1); S(b) = S(b) - (PV(i,2) + 1j*PV(i,3)); end

这段代码放在主循环开始之前，把S矩阵修正后再进入前推回代迭代。光伏注入用减法，因为S的方向是负荷方向。这里有一个容易搞错的地方：光伏发出的有功功率，在负荷复功率表达式中是负值。如果你不小心写成了加法，算出来的结果是光伏变成了负载，电压反而下降，这就完全反了。

4.4 闭合联络开关时的弱环处理

前推回代法无法直接处理合环网络。若想闭合 IEEE33 的某个联络开关，比如 8-21 开关，网络就出现一个闭环，此时需要采用弱环补偿电流法。常见做法是：先按纯辐射状网络计算一次前推回代，然后求解一个环路补偿电流，再把它叠加到相关支路电流上，重复迭代直到收敛。

% 伪代码示意，实际实现需要先构建环路增广矩阵 % 在每次前推回代后计算环路电压残余量和补偿电流 V_loop = 0; % 环路开关两端的电压差 for it = 1:max_iter V = run_pf_once(branch_radial, load, root); deltaV = V(open_from) - V(open_to); I_comp = deltaV / Z_loop_total; % 将补偿电流叠加到对应支路 end

这段只是思路示意，实际问题里需要根据开关两端节点建立环路阻抗矩阵。如果你做的是闭环运行模式下的潮流，还是建议在资源基础上增加配网特有的电压控制策略，或者直接换用牛顿法求解器。前推回代法在闭环场景下盲目扩写容易出数值问题。

5. 避坑：潮流计算里的五个常见翻车点

5.1 迭代始终不收敛，偏差在 1e-3 附近震荡

现象：run_pf跑完后报警告，mismatch停在 1e-3 级别上下跳动。

原因：最常见的是联络开关没有全部分断。IEEE33 有 5 个联络开关，如果在branch里把这些开关对应的支路留下，网络就会出现环，前推回代算法无法收敛。另外，最大迭代次数设太小时，程序在没稳定时就退出了。

解决：检查branch矩阵的行数。IEEE33 标准算例只有 32 条闭合支路，如果看到 37 条，说明 5 个联络开关被错误地当成闭合支路。把联络开关对应的支路从branch中删掉，再重新运行。

5.2 总负荷和文献值对不上

现象：用校核脚本算出总负荷不是 3715 kW 和 2300 kvar，而是偏差很大。

原因：负荷单位不统一。有些数据来源把负荷写成了标幺值，但程序里按 kW 处理；或者根节点负荷没有清零。

解决：先确认资源里的数据说明，看负荷是标幺值还是有名值。如果是标幺值，需要先乘以系统基准功率再放入load矩阵。根节点一般不挂负荷，若load第 1 行不是 0，要检查原因。

5.3 节点编号从 0 开始导致索引报错

现象：运行run_pf时 MATLAB 报V(0)索引错误，或者数组越界。

原因：不同文献的 IEEE33 系统根节点编号不一致，有的用 0，有的用 1。数据文件里的母线编号整体偏移了一位。

解决：在调用主程序前做一次整体平移：

if min(min(branch(:,1:2))) == 0 || min(load(:,1)) == 0 branch(:,1:2) = branch(:,1:2) + 1; load(:,1) = load(:,1) + 1; end

平移后再跑，索引问题就能消除。

5.4 接入光伏后迭代次数增加甚至震荡

现象：在某个节点加入 PV 后，收敛速度从原来 4 次变成 20 次，有时还会卡在容差边缘。

原因：把光伏当成恒定电压源节点处理了，前推回代法本身不支持 PV 节点。PV 节点需要额外的无功功率迭代修正，标准前推回代没有这个机制。

解决：把光伏建模成 PQ 负负荷。如果光伏并网点要求电压控制，则需要在每轮迭代后用电压偏差修正无功功率，再重新计算电流。不要指望单纯增加迭代次数能解决。

5.5 总损耗和文献基准对不上

现象：算出来的 33 节点系统总有功损耗是 180 kW，但文献结果是 202 kW 左右，怎么调都差一截。

原因：收敛容差设得太宽松，例如tol=1e-3，电压还没有稳定就停止迭代；或者根节点电压被错误地设置成了 0.98 pu。

解决：先把tol设为1e-8，max_iter设为 200，重新计算。确认根节点电压标幺值为 1.0。如果还是对不上，检查支路阻抗数据来源，某些文献给出的 33 节点阻抗是标幺值，需要还原成欧姆。

6. 验证与进阶：用 IEEE33 基准结果校验你的潮流程序

拿到一份潮流资源后，最先要做的事情是验证，而不是急着改参数。IEEE33 节点系统的公开基准结果比较统一，我把它当作验收标准。

你的程序跑出来的结果如果和这张表对不上，先别怀疑代码，按下面步骤排查一遍：第一步，确认load总负荷是 3715+j2300；第二步，确认branch行数是 32；第三步，用tol=1e-8重跑；第四步，检查根节点电压是否为 1.0。这套排查流程能解决九成的不一致问题。

验证通过之后，可以做两个进阶练习。第一个是绘制电压剖面，观察 DG 接入前后电压曲线变化：

case33; V_base = run_pf(branch, load, root, 1e-8, 200); % 在 18 号节点接入 300 kW 光伏 load_modified = load; load_modified(18, 2) = load_modified(18, 2) - 300; V_pv = run_pf(branch, load_modified, root, 1e-8, 200); figure; plot(abs(V_base), 'o-', abs(V_pv), 's-'); legend('原始系统', '18号节点接入300kW光伏'); xlabel('节点'); ylabel('电压/pu'); grid on;

这段代码用两种负荷场景对比电压曲线，可以直接看出光伏接入对末端电压的抬升作用。第二种进阶是把 33 节点数据换成 69 节点，用同样的主程序跑一遍，并将最低电压点和总损耗与文献对比，验证程序的可移植性。

我自己习惯的做法是，在run_pf外面再包一层检查函数，每次运行前自动校验节点编号连续性、联络开关状态、总负荷数值，然后才进主循环。这个习惯救过我很多次，曾经有一次就是因为没检查开关状态，盲调了一个小时也没找出原因。希望这篇笔记能让你避开同样的坑。

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