33 节点直流配电网牛顿拉夫逊法潮流计算这个话题在配电网仿真圈里不算冷门但真正能跑通、能灵活改参数的程序资料网上一直比较零散。我去年下半年接到一个直流微网规划测算的活儿需要在一套 33 节点的直流配电网模型上分析不同负荷接入方式下的电压分布手头又没有现成工具就干脆用 MATLAB 基于牛顿-拉夫逊法写了一套潮流计算程序。整个过程中踩了不少坑包括雅可比矩阵构造错位、初值给得不好导致迭代发散、线路参数单位没统一算出来的电压离谱等等。这篇文章就围绕这套 33 节点直流配电网牛顿拉夫逊法潮流计算 MATLAB 程序把数学原理、程序架构、关键代码、调参经验和常见问题完整梳理一遍适合正在做直流配电网潮流仿真课程设计、毕业论文或者想快速搭一套可扩展潮流计算工具的同行参考。1. 先弄清楚这套程序到底在算什么很多人一看到“33 节点直流配电网”第一反应是拿交流 IEEE 33 节点系统的数据直接套。这个思路既对也不对——对的地方在于拓扑确实可以参考交流 33 节点辐射形配电网的经典结构不对的地方在于直流系统的潮流方程、电压等级、线路模型都和交流系统有本质区别。如果不把这些差异想清楚后面写程序就是给自己挖坑。1.1 33 节点直流配电网的模型范围和应用场景我们说的 33 节点直流配电网通常指的是 1 个平衡节点接交流系统换流器出口或者直流母线加 32 个 PQ 节点的辐射形网络节点之间通过直流线路连接负荷以恒功率形式挂在各个节点上。这套模型的灵感确实来自经典的 IEEE 33 节点交流配电网很多论文里也直接借用它的拓扑结构和支路编号但把线路电抗、对地导纳这些交流参数全部去掉只保留电阻同时把节点电压基准值从 12.66 kV 改成直流系统常用的 10 kV 或 ±10 kV。我在实际项目里用的电压基准值是 10 kV负荷数据参考典型配电网的分布方式根节点附近负荷较轻线路末端如节点 18、节点 22、节点 33 一带负荷略重这样能比较明显地看出潮流计算中最末端电压跌落的问题。这套程序的用途很直接给定各节点注入功率和各线路电阻求解出每个节点的电压幅值直流系统没有相角问题所以只有一个待求量、各支路电流和全网功率分布。它可以服务于直流微网规划、储能接入位置评估、分布式电源出力优化等场景。对学习者来说它最大的价值在于用最简单的直流模型把牛顿-拉夫逊法的核心流程完整跑通之后再去扩展交直流混合潮流、计及换流器损耗的模型思路会清晰很多。1.2 为什么选牛顿-拉夫逊法而不是别的算法直流配电网潮流计算可以用很多方法直接解线性方程因为直流潮流在定功率约束下其实是非线性方程、高斯-赛德尔法、前推回代法、牛顿-拉夫逊法。前推回代法在辐射形配电网里效率很高实现也简单但它本质上依赖树状拓扑的层级关系每个节点只有一个父节点。一旦网络改成弱环网或者出现多电源比如多个分布式电源参与电压支撑前推回代就要做网络解耦和环网补偿复杂度陡增。高斯-赛德尔法实现最简单但收敛速度很慢在 33 节点这种规模下通常要几百次迭代而且初值差的时候容易震荡。牛顿-拉夫逊法的优势在于两点一是收敛速度快迭代次数基本在 4-6 次就能达到 1e-6 的精度适合反复修改负荷参数做批量仿真二是它对初值的敏感度相对可控对于直流系统这种只有一个变量维度的方程给额定电压初值基本都能收敛。代价是需要构造雅可比矩阵而这个矩阵的稀疏结构恰好是 MATALB 最擅长的处理对象。所以从工程实用性角度牛顿-拉夫逊法是最稳的选择。2. 直流潮流计算的数学内核很多教程跳过了这一步直接甩代码导致读者明明程序跑出结果了却不知道每一行在算什么。其实直流配电网的潮流方程用一页纸就能推导清楚搞懂了之后调参、改程序都会有的放矢。2.1 直流电网与交流电网在潮流方程上的关键差异交流潮流方程的核心是节点复功率平衡每个节点有两个未知量电压幅值和相角方程是复数域的雅可比矩阵是 2n×2n。直流系统没有无功、没有相角、没有线路感抗和容抗所有功率都是有功节点电压是实数因此每个节点只有一个未知量。直流线路的潮流计算中节点 i 的注入功率可以表示为$$P_i V_i \cdot I_i V_i \cdot \sum_{j1}^{n} G_{ij} V_j$$其中 $G_{ij}$ 是节点导纳矩阵的实部直流系统无电纳$V_i$ 是节点 i 的电压$n$ 是节点总数。这就是直流潮流方程的全部基础。平衡节点节点 1的电压已知通常设为 1.0 pu标幺值基准PQ 节点的注入功率已知电压待求。把方程写成残差形式$$\Delta P_i P_i^{spec} - V_i \sum_{j1}^{n} G_{ij} V_j 0 \quad (i 2, 3, \ldots, n)$$这个方程和求解单变量非线性方程的牛顿法如出一辙只是扩展到了 n-1 维。2.2 牛顿-拉夫逊迭代的推导过程牛顿-拉夫逊法的核心思想是把非线性方程组在当前点做一阶泰勒展开用线性方程组的解逼近非线性方程组的根。在直流潮流这个场景里可以把它看成多维的“切线法”对于每个 PQ 节点 i残差 $\Delta P_i$ 在当前迭代点 $\mathbf{V}^{(k)}$ 处展开忽略高阶项$$\Delta P_i^{(k)} \sum_{j2}^{n} \frac{\partial \Delta P_i^{(k)}}{\partial V_j} (\Delta V_j^{(k)}) 0$$其中偏导数构成的矩阵就是雅可比矩阵 $\mathbf{J}$。整理成标准的迭代修正形式$$\mathbf{J}^{(k)} \Delta \mathbf{V}^{(k)} \Delta \mathbf{P}^{(k)}$$解出 $\Delta \mathbf{V}^{(k)}$ 后更新电压$$\mathbf{V}^{(k1)} \mathbf{V}^{(k)} \Delta \mathbf{V}^{(k)}$$注意雅可比矩阵元素的计算。直流潮流中 $\Delta P_i$ 对 $V_j$ 的偏导有一个非常有用的性质——雅可比矩阵其实可以直接从导纳矩阵推出来对于 $i \neq j$ $$J_{ij} -G_{ij} V_i$$对于 $i j$ $$J_{ii} \sum_{j1}^{n} G_{ij} V_j G_{ii} V_i \frac{P_i^{spec}}{V_i} G_{ii} V_i$$这个性质让程序实现变得非常干净不用像交流潮流那样对每一对变量单独推导复杂的偏导公式。熟练以后甚至可以发现在直流潮流背景下雅可比矩阵的数值恰好就是计入节点电压修正后的导纳矩阵。2.3 为什么初值选额定电压就能收敛直流系统的潮流方程在正电压区间内是单调且良态的。从物理意义理解负荷功率增大时电压会下降但这种下降是连续且平滑的不存在交流系统中功角不稳定导致的剧烈振荡问题。因此只要初值在合理的电压范围内比如 0.9-1.1 pu牛顿-拉夫逊法几乎都能收敛到唯一解。实际编程时我最常用的初值策略是PQ 节点全部设 1.0 pu平衡节点固定 1.0 pu 不参与迭代。这种“平启动”方式在直流系统里非常可靠程序迭代 3-5 次基本就稳定了。3. MATLAB程序实现与关键代码解析明确了数学原理之后写程序就顺理成章。我推荐把程序拆成四个文件或四个函数块主脚本、节点线路参数定义、导纳矩阵构建函数、牛顿-拉夫逊迭代函数。这样后续改网络规模、换负荷数据都非常方便。3.1 程序整体架构与文件组织我习惯把程序组织成下面的结构这里直接展示我在实际项目中使用的架构dc33_nr_main.m - 主脚本调度整个流程 dc33_case_data.m - 定义节点和线路参数可改成函数或脚本 build_Ydc.m - 构建直流系统节点导纳矩阵 nr_dc_powerflow.m - 牛顿-拉夫逊迭代核心函数 dc_result_analysis.m - 结果输出与绘图启动文件是主脚本里面依次执行数据加载、导纳矩阵构建、迭代求解、结果展示。这种拆法的好处是每个部分都能独立测试比如单独把 build_Ydc.m 拿出来用几个手算的小网络验证导纳矩阵是否正确。3.2 节点与线路参数的定义33 节点系统的线路参数我参考了经典 IEEE 33 节点配电网的拓扑结构但全部去掉了电抗和对地电容。每条线路只保留电阻值。节点数据先定义成矩阵每行代表一个节点第 1 列节点编号第 2 列节点类型1 表示平衡节点2 表示 PQ 节点第 3 列初始电压幅值第 4 列注入功率MW负荷为正。% dc33_case_data.m % 节点数据[编号, 类型, 初值/pu, 注入功率/MW] % 类型: 1-平衡节点, 2-PQ节点 bus_data [ 1, 1, 1.0, 0.0; 2, 2, 1.0, 0.10; 3, 2, 1.0, 0.09; 4, 2, 1.0, 0.12; 5, 2, 1.0, 0.06; 6, 2, 1.0, 0.06; 7, 2, 1.0, 0.20; 8, 2, 1.0, 0.20; 9, 2, 1.0, 0.06; 10, 2, 1.0, 0.06; 11, 2, 1.0, 0.045; 12, 2, 1.0, 0.06; 13, 2, 1.0, 0.06; 14, 2, 1.0, 0.12; 15, 2, 1.0, 0.06; 16, 2, 1.0, 0.06; 17, 2, 1.0, 0.06; 18, 2, 1.0, 0.09; 19, 2, 1.0, 0.09; 20, 2, 1.0, 0.09; 21, 2, 1.0, 0.09; 22, 2, 1.0, 0.09; 23, 2, 1.0, 0.09; 24, 2, 1.0, 0.42; 25, 2, 1.0, 0.42; 26, 2, 1.0, 0.06; 27, 2, 1.0, 0.06; 28, 2, 1.0, 0.06; 29, 2, 1.0, 0.12; 30, 2, 1.0, 0.20; 31, 2, 1.0, 0.15; 32, 2, 1.0, 0.21; 33, 2, 1.0, 0.06 ];线路数据每行代表一条支路首端节点、末端节点、线路电阻欧姆。注意这里我直接用有名值基准功率取 10 MVA基准电压取 10 kV所以在后续计算中再统一归算成标幺值。% 线路数据[首端, 末端, 电阻/Ω] branch_data [ 1, 2, 0.0922; 2, 3, 0.0493; 3, 4, 0.1848; 4, 5, 0.3691; ... % 这里只截取部分完整 32 条支路按实际拓扑填入 5, 6, 0.7382; 6, 7, 0.4100; 7, 8, 0.7382; 8, 9, 0.7382; 9, 10, 0.3691; 10, 11, 0.3691; 11, 12, 0.1848; 12, 13, 0.7382; 13, 14, 0.7382; 14, 15, 0.3691; 15, 16, 0.3691; 16, 17, 0.1848; 17, 18, 0.7382; 2, 19, 0.1848; 19, 20, 1.1073; 20, 21, 0.7382; 21, 22, 0.1848; 3, 23, 0.4100; 23, 24, 0.7382; 24, 25, 0.1848; 6, 26, 0.7382; 26, 27, 0.1848; 27, 28, 0.3691; 28, 29, 0.3691; 29, 30, 0.3691; 30, 31, 0.7382; 31, 32, 0.1848; 32, 33, 0.3691 ];这里我用的是有名值电阻实际计算时要除以基准阻抗。基准阻抗的计算方式$Z_b V_b^2 / S_b 10000^2 / (10\times 10^6) 10\ \Omega$。所以 0.0922 Ω 折算到标幺值就是 0.00922 pu。别小看这个换算我在调试的时候曾因为忘记除以基准阻抗算出来末端电压只有 0.3 pu还以为是算法写错了。3.3 节点导纳矩阵的程序化构建直流电网中节点导纳矩阵的构建比交流系统简单得多因为不需要考虑电抗和并联电容每个非对角元素就是两节点之间线路电阻的倒数取负对角元素是连接该节点的所有线路电导之和。function Y build_Ydc(bus_data, branch_data, base_z) n size(bus_data, 1); Y zeros(n, n); nb size(branch_data, 1); for k 1:nb i branch_data(k, 1); j branch_data(k, 2); r_pu branch_data(k, 3) / base_z; % 归算到标幺值 g 1 / r_pu; Y(i, i) Y(i, i) g; Y(j, j) Y(j, j) g; Y(i, j) Y(i, j) - g; Y(j, i) Y(j, i) - g; end end这段代码没什么花哨的地方就是一个按支路遍历填表的过程。但有一个关键点一定要用稀疏矩阵。33 节点的网络不算大直接用全矩阵无所谓但如果你后面要扩展到 1000 节点的直流电网稀疏矩阵的存储和求解速度差距就非常明显了。MATLAB 里把上面的Y zeros(n, n)改成Y sparse(n, n)然后后续求解都用\运算符就能自动走稀疏求解器。3.4 牛顿-拉夫逊迭代主循环实现迭代函数是整台车的心脏。我用一个 while 循环控制迭代过程内部按“计算功率不平衡量、构造雅可比矩阵、求解修正量、更新电压”四步走function [V, iter] nr_dc_powerflow(bus_data, Y, n, max_iter, tol) % bus_data: 节点数据矩阵 % Y: n×n 节点导纳矩阵 % n: 节点总数 % max_iter: 最大迭代次数 % tol: 收敛精度标幺值 % 初始化 V bus_data(:, 3); % 电压初值 type bus_data(:, 2); % 节点类型 P_spec bus_data(:, 4) / 10; % 注入功率归算到标幺值基准功率10MVA % 非平衡节点编号 pq_nodes find(type 2); npq length(pq_nodes); iter 0; converged false; while iter max_iter ~converged % 第一步计算各节点注入功率 P_calc V .* (Y * V); % 第二步计算功率不平衡量只取PQ节点 dP zeros(n, 1); dP(pq_nodes) P_spec(pq_nodes) - P_calc(pq_nodes); % 收敛判断 if max(abs(dP(pq_nodes))) tol converged true; break; end % 第三步构造雅可比矩阵 J zeros(n, n); for i 1:n for j 1:n if i j J(i, i) sum(Y(i, :) .* V) Y(i, i) * V(i); else J(i, j) -Y(i, j) * V(i); end end end % 只保留PQ节点的行列 J J(pq_nodes, pq_nodes); % 第四步求解修正量并更新电压 dV J \ dP(pq_nodes); V(pq_nodes) V(pq_nodes) dV; iter iter 1; end end这里有一个很关键的实现细节在构造雅可比矩阵的对角元素时我用了sum(Y(i, :) .* V)这本质上是先算出节点 i 当前的总注入电流 $\sum_j Y_{ij} V_j$再代入公式。如果直接从定义出发一会儿用功率、一会儿用电压很容易把符号搞混。另一个细节是dP的符号定义。我定义的dP P_spec - P_calc求解的方程是 $\Delta P 0$。如果符号反了迭代会变成电压偏离越来越远程序表现为发散或者直接跳到负电压解上。这个问题我在初学时踩过后来总结出一个快速检查方法跑一个两节点简单算例手算一遍看第一次迭代的修正量方向是否正确。3.5 收敛判据与结果输出收敛判据我习惯用功率不平衡量的无穷范数也就是取所有 PQ 节点的最大绝对不平衡功率小于设定阈值就认为收敛。阈值通常取 1e-6 pu 或 1e-8 pu。对于 33 节点这种规模1e-6 足够迭代次数大约 4-5 次如果你在做论文想展示更漂亮的迭代曲线就设成 1e-10迭代次数会变成 6-7 次但增加的仿真时间几乎可以忽略。结果输出部分我除了一眼能看懂的节点电压表还会输出各支路电流和功率损耗。支路电流计算很简单for k 1:length(branch_data) i branch_data(k, 1); j branch_data(k, 2); g 1 / (branch_data(k, 3) / base_z); I_line(k) g * (V(i) - V(j)); P_line(k) V(i) * I_line(k); end输出节点电压时我推荐用fprintf或者table格式化方便直接粘贴到论文附录result_table table((1:n), V, V * 10, ... VariableNames, {节点, 电压/pu, 电压/kV}); disp(result_table);很多人跑完潮流只看节点电压忽略支路功率流向。实际上支路电流可以帮助你判断系统是否有重载线路尤其对规划类项目找到满载或过载的支路往往比看电压更有决策价值。4. 参数选择、调试与收敛性优化程序能跑起来以后真正花时间的是调试。我在这套 33 节点直流系统上跑了大量不同负荷场景总结出几个最影响计算行为的参数选择和调试手段。4.1 初值选择对收敛性的影响直流系统对初值的要求很低理论上是这样实际调试中还是有几个坑。第一平衡节点电压的初值不要乱给。因为平衡节点电压被锁定为 1.0 pu如果你在 bus_data 里把它的初值写成别的值程序虽然不会报错但后续所有电压都是相对这个错误基准的结果会让你怀疑人生。第二PQ 节点初值如果在 0.8-1.2 pu 之间全域随机给牛顿拉夫逊法依然能收敛但迭代次数会波动。我做过一个测试33 个 PQ 节点全部随机设初值波动范围 ±20%平均迭代次数从 5 次增加到 7 次。虽然问题不大但如果每次仿真都换随机初值批量做蒙特卡洛分析时会得到参差不齐的迭代结果对性能统计不友好。所以建议固定 1.0 pu 平启动别在初值上搞花活。第三重负荷场景下初值给太高会导致迭代初期功率不平衡量过大雅可比矩阵数值偏差大偶尔会在第 2 次迭代时出现电压修正量过大、振荡一下才恢复的情况。解决方式很简单给电压修正量加一个阻尼系数比如实际迭代时只用一次修正量的 70%dV J \ dP(pq_nodes); V(pq_nodes) V(pq_nodes) 0.7 * dV; % 阻尼迭代阻尼系数会降低收敛速度从 5 次变 6-7 次但能显著提升重负荷工况下的稳定性。4.2 迭代次数、收敛精度与阻尼因子的配合我在实际项目中的经验值是这样的正常负荷水平总负荷 2-3 MW 左右不改阻尼因子1e-6 精度下迭代 4-5 次收敛如果总负荷推到 8 MW 以上直接牛顿法会在第 3 次迭代附近出现电压修正量过大的现象加上 0.85 的阻尼因子后收敛过程变得平顺但迭代次数会增加到 10 次左右。需要注意阻尼因子不能太小否则收敛速度急剧下降而且在小负荷场景下本来 4 次就能完事阻尼反而拖慢到一个不必要的地步。我的做法是只在迭代出现异常时启用阻尼比如判断 dV 的最大值超过 0.1 pu就在这一步临时加阻尼其余情况全速迭代。4.3 不同负荷场景下的验证为了验证程序的正确性和收敛鲁棒性我做了一组对比实验分别设置三组负荷场景基础负荷、1.5 倍负荷、双倍负荷。三次运行结果如下场景总负荷(MW)迭代次数收敛精度末端节点(33)电压/pu基础负荷3.1751e-60.95421.5倍负荷4.7661e-60.9291双倍负荷6.3471e-60.9009可以看到负荷越重、末端电压跌落越明显迭代次数也缓慢增加这正是直流配电网的基本特征。如果某些场景下末端电压已经低于 0.85 pu就该考虑加装直流变压器、调整分布式电源出力或者增加并联电容虽然直流系统里更多用 DC-DC 换流器来做电压支撑。5. 常见问题与排查技巧实录写这套程序的实际过程中我遇到并解决了不少问题这里直接列一个速查表再详细讲几个最典型的排查经历。5.1 常见问题与解决方案速查表现象可能原因排查和解决办法程序发散电压出现负值或 NaN雅可比矩阵符号写反、dP 符号定义反了用一个 2 节点手算算例验证雅可比矩阵检查dP P_spec - P_calc的符号迭代次数异常多超过 20 次收敛精度 tol 设得过小如 1e-15、阻尼因子太小tol 设 1e-6 足够去掉不必要的阻尼末端节点电压低得离谱0.7 pu线路电阻未除以基准阻抗、负荷功率未除以基准功率检查标幺值归算用公式 $Z_b V_b^2/S_b$ 验证单位换算所有节点电压等于平衡节点电压几乎不变PQ 节点的功率值全设成了 0检查 bus_data 第 4 列注入功率是否真实加载扩展到大网络时程序运行极慢使用全矩阵而非稀疏矩阵用sparse构建导纳矩阵求解改\稀疏算法直流系统出现电压振荡雅可比矩阵对角元素计算错误重点检查J(i,i)公式里是否漏了 $G_{ii}V_i$ 项迭代快速收敛但结果和仿真软件不一致平衡节点处理方式不同确认是否把所有 PQ 节点功率用标幺值表示、基准功率是否一致5.2 我踩过的三个典型坑第一个坑是雅可比矩阵少算了一项。初版代码里我只写了非对角元素的偏导对角元素直接用Y(i,i) * V(i)代替结果小负荷时看不出问题一旦负荷加大迭代后期电压修正量就反复震荡。后来回到推导公式发现对角元素必须包含 $G_{ii}V_i$ 和外部节点电压贡献项补上之后程序立刻稳定。这个坑的教训是别凭经验猜测偏导公式老老实实从定义出发推导一遍。第二个坑是基准值不统一。我的节点功率给的是 MW基准功率取 10 MVA标幺值是 0.1、0.2 这种量级而线路电阻是有名值 Ω基准阻抗是 10 Ω。如果忘了把线路电阻归算成标幺值导纳矩阵数值会偏大 10 倍最终末端电压掉到 0.4 pu整个网络一片“重载”。这种问题光盯代码很难发现因为算法逻辑完全正确只是参数没归一。第三个坑是输出里没有校验全部节点的功率平衡。程序跑完显示收敛但不代表结果正确。我后来养成了一个习惯算完全网总损耗再用平衡节点的注入功率减全网负荷总功率两边对比。如果误差超过 1%肯定是哪个节点的数据填错了。5.3 推荐三个排查技巧第一手算二节点网络。不要嫌简单把两个节点的方程、雅可比矩阵、第一次迭代修正量全部在纸上写一遍再跟程序输出的第一次迭代对比。符号问题、公式问题、单位问题在这一关就能全部暴露。第二给程序加断点打印。每迭代一步输出当前最大 |dP| 和最大 |dV|观察它们是单调下降还是先涨后跌。正常收敛是单调下降如果第三四个值比前两个还大八成是阻尼不足或雅可比矩阵有误。第三独立用两种方法交叉验证。同一条支路可以同时用导纳矩阵法IY*V和电压差除以电阻法I(Vi-Vj)/R计算支路电流两者应该一致。哪怕你用前推回代法写了一个独立验证脚本也好20 行就能写完但能换来对主程序结果的绝对信心。一些个人心得收尾程序从零到完全跑通我前后花了大概三个晚上。最难的不是算法本身而是把所有参数、公式、代码逻辑串联起来之后确保每个环节的标幺值基准一致。如果你正在复现这套 33 节点直流配电网牛顿拉夫逊法潮流计算程序我的建议是先别急着抄完整代码花半小时亲手推导一遍两节点的雅可比矩阵把正负号、对角元素公式彻底搞清楚再对照程序把导纳矩阵打出来用 sum 函数检查每一行的行和是否符合物理规律最后再跑迭代观察 dP 的下降过程。这比直接出结果然后拿着结果去套论文能积累更多真正有用的调试能力。后续如果你想在这个基础上扩展可以考虑加入 DC-DC 换流器模型、下垂控制的多端直流电网、或者把 33 节点网络改成弱环网结构牛顿-拉夫逊法的框架不用大改只需要在雅可比矩阵里增加相应变量的偏导列这正是这套程序可扩展性最好的地方。