做直流配电网的仿真时最头疼的往往不是模型本身而是找不到一套真正顺手、能看清每一步原理的潮流计算程序。交流配电网的潮流程序一抓一大把NR法、前推回代、PQ分解法全都有现成的可一旦把模型换成纯直流无功、相角、频率全都不存在了那些程序里一大半逻辑变成了摆设。实际上直流配电网的潮流计算比交流版本要干净得多——节点电压是唯一的状态变量牛顿拉夫逊法在这里只需要维护一张(n-1)阶的雅可比矩阵。这篇文章围绕经典33节点辐射状配电网算例的直流化改造把NR法的公式推演、MATLAB程序实现、收敛性分析和工程扩展完整拆开讲一遍适合做直流配电、直流微电网仿真或者课程设计的同学直接拿去改。1. 直流化到底改了啥从经典33节点算例到一张纯电阻网络1.1 交流潮流方程与直流潮流方程的本质差异先看交流系统里的节点功率方程。对任意节点i注入功率可以写成两部分有功部分P_i V_i × Σ_j [G_ij×cos(θ_i-θ_j) B_ij×sin(θ_i-θ_j)] × V_j无功部分Q_i V_i × Σ_j [G_ij×sin(θ_i-θ_j) - B_ij×cos(θ_i-θ_j)] × V_j这里有电压幅值V、相角θ、有功P、无功Q四个量所以交流NR法的雅可比矩阵是4个分块矩阵的组合分别处理P对θ、P对V、Q对θ、Q对V的偏导数计算规模直接翻倍。而且交流潮流还要考虑PV节点的无功越限、变压器分接头位置、线路充电电容一堆附加因素。直流配电网稳态运行时没有无功功率没有频率也没有相角差的概念。换流器控制直流电压和输送有功线路和负荷的数学模型只剩电阻。节点电流方程退化成I G × V节点注入功率则是P_i V_i × I_i V_i × Σ_j (G_ij × V_j)这里G就是纯电导矩阵没有B虚部。这意味着NR法的未知量只剩一个维度——节点电压幅值迭代方程的规模直接砍掉一半以上。这就是为什么直流配电网潮流程序能写得很精简也值得单独写一版而不是拿交流程序硬改。1.2 经典33节点算例的直流化改造思路标题里的“33节点”用的是配电领域非常经典的那个33节点辐射状测试系统原始参数包含线路电阻、电抗、有功负荷和无功负荷。做直流配电网研究时通常不会重新构建一套全新拓扑而是直接把这个成熟算例做直流化改造拓扑结构完全保留33个节点、32条支路外加5条可操作的联络开关支路线路参数只取电阻列电抗全部置零负荷只保留有功分量无功分量直接丢弃节点1作为平衡节点由电压源换流器维持直流母线电压恒定其余32个节点全部按恒功率节点处理。直流化之后的网络就是一张纯电阻网络电压等级我习惯按直流配电网常用的10kV设定基准功率取1MVA。这样每个节点的负荷功率数量级在0.1~0.4 p.u.之间潮流计算的结果分布比较典型也方便和后续工程数据做对比。1.3 节点类型大幅简化交流潮流里节点分三类平衡节点、PQ节点、PV节点。PV节点需要给定有功和电压幅值无功是待求量这给雅可比矩阵的构建和迭代过程中的无功越限处理增加了不少代码量。直流配电网里没有无功的概念节点类型只剩两种交流节点类型直流对应节点类型给定量待求量平衡节点直流平衡节点VSC恒压控制节点电压幅值注入功率PQ节点恒功率节点P节点注入有功功率节点电压幅值PV节点不存在除非有恒压控制变流器但那本质上是另一个平衡节点——这个简化听着不起眼但实际编码时省掉的事情非常多。不需要处理θ角排序、不需要分块雅可比、不需要做PV→PQ的越限切换迭代矩阵从(2n-2)×(2n-2)直接降到(n-1)×(n-1)逻辑链条短了一大截。顺便提一句如果某个DC/DC换流器或者储能变流器采用恒压控制它在模型里就是又一个“已知电压节点”本质上和平衡节点等价。在代码里处理多平衡节点时只需要把这些节点的电压也固定住把它们从迭代方程中剔除即可后面第5章会说到。2. 牛顿拉夫逊法的直流版推演一个雅可比矩阵说清全部逻辑2.1 功率不平衡方程与迭代修正的符号约定先把直流潮流的NR法完整推一遍。目标是找到一组节点电压V使得每个恒功率节点的注入功率等于给定值。定义功率不平衡量ΔP_i P_spec_i - P_calc_i其中P_spec_i是节点i给定的注入功率负荷节点为负值分布式电源节点为正值P_calc_i是根据当前电压计算得到的注入功率P_calc_i V_i × Σ_j (G_ij × V_j)牛顿法的核心思路是在当前电压点做一阶泰勒展开求一组电压修正量ΔV让不平衡量归零。写成矩阵形式J × ΔV ΔP这里J就是雅可比矩阵也就是P_calc对电压V的偏导数矩阵。迭代更新公式为V V ΔV很多人在这一步容易把符号搞反。常见的一种教材写法是定义F P_calc - P_spec然后解J×ΔV -F我这里定义的是ΔP P_spec - P_calc所以修正方程里没有负号。两种写法最终等价但一开始就必须固定下来否则程序里一会儿收敛一会儿发散最后还以为是算法问题。2.2 雅可比矩阵的解析表达式不对称才是正常状态对P_calc_i V_i × (G(i,:)·V)逐个求偏导。对角元∂P_calc_i / ∂V_i G(i,:)·V V_i × G_ii这里有项要特别注意G(i,:)·V本身含有V_i因为求和里包含G_ii×V_i这一项所以对V_i求导时会产生两项——一项来自当前电压点处的网络注入项另一项来自系数G_ii与V_i的乘积。非对角元∂P_calc_i / ∂V_j V_i × G_ij其中 i ≠ j非对角元的推导就简单得多G(i,:)·V对V_j求导得到系数G_ij再乘上前面的V_i因子即可。有一点需要单独强调直流网络的电导矩阵G本身是对称的G_ij G_ji但雅可比矩阵J并不对称。因为J(i,j) V_i×G_ij而J(j,i) V_j×G_ji V_j×G_ij两者只有在V_i V_j时才相等。迭代过程中各节点电压一般都不相等所以J基本上不会保持对称。如果用某些对称矩阵求解算法或者把J误写成对称矩阵结果一定不对。2.3 收敛判据与整体迭代流程收敛判据一般用功率不平衡量的无穷范数即所有节点ΔP_i的最大绝对值。我习惯设置为1e-8 p.u.对应到1MVA基准下就是0.01W量级对潮流计算来说精度足够且不会让迭代无谓地反复。完整迭代流程可以概括为初始化所有电压为1.0 p.u.平启动平衡节点电压置为给定值用当前V计算每个节点的P_calc计算ΔP向量取PQ节点的子集判断max(|ΔP|)是否小于收敛容差满足则结束构造雅可比矩阵J只取PQ行、PQ列求解修正方程ΔV J\ΔP更新V V ΔV回到第2步。交流NR法通常需要十几次迭代甚至几十次直流版由于没有相角变量问题规模小、非线性程度也低正常算例基本5~8次迭代就能收敛。如果超过15次还在动就要怀疑数据或者初值出问题了这一点后面专门讲。3. MATLAB程序分段拆解先跑通小算例再装进33节点拓扑3.1 数据组织表驱动而不是写死在公式里程序的第一件事是把网络数据和负荷数据组织成可替换的表。我最推荐的做法是用单独的数组维护拓扑和负荷这样换一个算例只需要改数据不需要动算法代码。负荷数据用一行数组元素顺序对应节点1到33% 节点有功负荷单位 kW顺序对应节点 1~33 % 直流化改造后按经典33节点算例常用参数整理 Pload_kW [0, 100, 90, 120, 60, 60, 200, 200, 60, 60, 45, 60, 60, 120, 60, 60, 60, ... 90, 90, 90, 90, 90, 90, 420, 420, 60, 60, 60, 120, 200, 150, 210, 60];支路数据用三列矩阵格式是[首端节点号末端节点号电阻值]。注意直流化之后没有电抗列这是和交流程序最大的数据差异% 支路数据[首端, 末端, 电阻欧姆] % 完整33节点支路电阻按经典算例录入这里展示前6条作为格式示例 branch [ 1 2 0.0922 2 3 0.4930 3 4 0.3660 4 5 0.3811 5 6 0.8190 6 7 0.1872 % ... 继续扩到全部支路和联络开关 ];基准值设置S_base 1e6; % 基准功率 1 MVA V_base 10e3; % 基准电压 10 kV Z_base V_base^2 / S_base; % 基准阻抗这里等于 100 欧姆这个基准值的好处是负荷功率都在0.04~0.42 p.u.的合理区间电压结果读起来直观。有特殊需求的话把V_base改成±10kV对称双极系统对应的极间电压即可其他逻辑不用动。3.2 电导矩阵形成与NR迭代主循环电导矩阵G的装配逻辑很简单每一条支路电阻r对应电导g 1/r自导纳加上g互导纳减去g。这才是真正理解“节点电压方程I G×V”的关键——G的每一行描述的是该节点与全网其他节点的电流电压关系。%% 形成电导矩阵标幺值 n_node 33; G zeros(n_node, n_node); for k 1:size(branch, 1) i branch(k, 1); j branch(k, 2); r_pu branch(k, 3) / Z_base; g 1 / r_pu; G(i, i) G(i, i) g; G(j, j) G(j, j) g; G(i, j) G(i, j) - g; G(j, i) G(j, i) - g; end G sparse(G); %% 生成节点注入功率标幺值 P_spec -Pload_kW / S_base; % 负荷为负注入 slack 1; P_spec(slack) 0; % 平衡节点功率由潮流解算后得到 %% 初始化与迭代参数 V ones(n_node, 1); V(slack) 1.0; % 平衡节点电压给定 PQ setdiff(1:n_node, slack); % 恒功率节点集合 tol 1e-8; max_iter 50; iter 0; %% 牛顿拉夫逊迭代主循环 while iter max_iter % 计算当前电压下的节点注入功率P_calc V .* (G * V) P_calc V .* (G * V); % 功率不平衡量 dP_all P_spec - P_calc; dP dP_all(PQ); % 收敛判断 if max(abs(dP)) tol break; end % 构造雅可比矩阵 J dP_calc / dV nPQ length(PQ); J zeros(nPQ, nPQ); for ii 1:nPQ i PQ(ii); % 对角元 J(ii, ii) G(i, :) * V V(i) * G(i, i); for jj 1:nPQ if jj ii continue; end j PQ(jj); % 非对角元 J(ii, jj) V(i) * G(i, j); end end % 求解修正方程并更新电压 dV J \ dP; V(PQ) V(PQ) dV; iter iter 1; end这一段代码是整个程序的骨架。P_calc V .* (G * V)这一步是MATLAB向量化运算的经典写法G*V先得到各节点的注入电流向量再逐点乘以电压得到注入功率简洁也能避免写内层循环。构造雅可比矩阵的两重循环对33节点来说完全够快因为矩阵只有32×32。如果你要扩展到几百上千节点可以再用向量化技巧或者稀疏矩阵存储优化但日常研究阶段完全没必要折腾。3.3 先跑通一个微型算例再放大到33节点直接上33节点数据调试有个问题一旦结果不对你很难分辨是算法错、数据错还是某个联络开关没处理对。我的习惯永远是先构造一个最小的可手算网络把算法里的每一个符号、每一个数值方向验证清楚再放大到目标算例。最小算例长这样平衡节点1电压1.0 p.u.节点2通过一条电阻R 100欧姆的线路接到节点1节点2带100kW恒功率负荷。基准阻抗Z_base 100欧姆所以R的标幺值正好是1.0100kW对应标幺值0.1。节点的功率方程是P_calc_2 V_2 × (G_22×V_2 G_21×V_1) V_2² - V_2因为这里G_22 1G_21 -1V_1 1。负荷100kW是吸收功率所以注入功率P_spec_2 -0.1。方程就成了V_2² - V_2 -0.1解出来有两个根V_2 0.8873 或 V_2 0.1127。潮流计算取高电压解也就是0.8873这是恒功率负荷下直流网络的正常运行点。用上面的主程序把n_node改成2branch填[1 2 100]Pload_kW填[0 100]从V_2 1.0平启动开始迭代前几步就会快速逼近0.8873。如果程序迭代轨迹和这个对不上说明符号约定或者雅可比构造有bug。这比盯着33节点网络找bug高效得多。3.4 结果输出与功率平衡校验迭代结束后的输出部分除了打印节点电压我还会额外做两项校验平衡节点注入功率和全网损耗这是判断结果是否可信的硬指标。%% 结果输出与校验 P_bal V(slack) * (G(slack, :) * V); % 平衡节点注入功率标幺 P_loss V * G * V; % 全网电阻损耗标幺 P_total_sp sum(P_spec); % 所有非平衡节点注入之和 fprintf(迭代次数%d\n, iter); fprintf(平衡节点注入功率%.4f p.u.%.2f kW\n, P_bal, P_bal * S_base / 1e3); fprintf(全网网损%.4f p.u.%.2f kW\n, P_loss, P_loss * S_base / 1e3); fprintf(功率平衡校验%.6f\n, P_bal P_total_sp - P_loss);功率平衡校验的物理意义是平衡节点注入功率加上其它所有节点注入功率的总和应该等于网络电阻损耗。因为所有节点注入功率之和等于V × G × V也就是支路电阻上的总损耗。校验值在1e-6以下基本说明结果自洽。节点电压分布可以直接画柱状图figure; bar(1:n_node, V, FaceColor, [0.3 0.5 0.7]); grid on; xlabel(节点编号); ylabel(电压 (p.u.)); title(33节点直流配电网潮流电压分布);33节点辐射状直流配电网的电压分布通常有一个明显的趋势离平衡节点越远电压越低线路越长、负荷越重的地方压降越明显。末端节点如果电压掉到0.93 p.u.以下就该考虑是不是负荷太重或者线径选细了。4. 收敛性分析恒功率负荷的临界负载与那些“发散得莫名其妙”的时刻4.1 恒功率负荷为什么会导致潮流无解直流配电网里恒功率负荷是一个很微妙的存在。它吸收的功率不随电压变化电压降低时电流反而增大这会给网络带来一种类似正反馈的压力。回到刚才的微型算例方程是V_2² - V_2 P_L 0其中P_L是节点2吸收的有功标幺值取正数。这是一个一元二次方程判别式是1 - 4×P_L。只有当P_L ≤ 0.25 p.u.时才有实数解也就是说这个单负荷直流网络最多只能带250kW基准1MVA。超过这个值潮流方程没有可行解NR法从任何初值出发都会发散。这不是程序bug而是网络的物理极限。恒功率负荷比重越高直流系统越容易出现这种电压稳定问题。实际工程中换流器或储能系统会提供电压支撑但纯潮流计算时如果负荷给得太重程序发散是正常现象得从网络参数和运行点去找原因。4.2 平启动失败时的三个排查方向NR法跑发散做完微型算例验证之后最可能的原因有三个按出现频率排序一是数据录入错误。支路首末端节点号写反、电阻数量级差10倍、负荷单位没从MW换算成kW这些错误会让雅可比矩阵方向没问题但解根本不存在。我见过最隐蔽的一次是联络开关支路的电阻填成了0.0001欧姆相当于把两个节点直接短接结果电压被拉得乱七八糟。二是初值离解太远。直流NR法平启动V1.0通常没问题但重负荷或者高阻抗网络下解的电压可能远低于1.0平启动迭代容易跑到低电压解或者直接发散。这时可以把初值改到0.9甚至0.8试试也可以先跑恒阻抗负荷模型得到一组相对稳定的电压把它当作恒功率模型的初值。三是负荷总量超过了网络承载能力。前面说的临界负载问题在网络规模变大后依然存在只是没有一个简单的判别式可以直接算出来。此时需要逐步增加负荷观察哪个点开始发散这实际上就是连续潮流法的雏形。4.3 每次迭代都记录最大不平衡量调程序时我会在循环里记录每一步的max(abs(dP))并打印出来这比只看到最终“不收敛”三个字有用得多。正常收敛时这个值是一条近似单调快速下降的曲线通常从10的-2次方量级一路掉到1e-8以下。如果看到曲线先下降几步又反弹或者在一个数值附近来回震荡多半是雅可比矩阵构造有误重点检查对角元是不是漏了G(i,:)*V这一项。如果曲线完全发散每一步都比上一步大先检查符号方向再看初值。另外求解修正方程时用J\b而不是inv(J)*b。MATLAB的反斜杠运算符自动选择合适的高斯消元求解方法数值稳定性比显式求逆好速度也更快。虽然32×32矩阵用哪种都能算出来但好的习惯可以避免以后扩展到大规模网络时踩坑。5. 三个趁手的扩展方向合环运行、ZIP负荷、分布式电源与下垂控制5.1 直接闭合联络开关NR法天然支持环网经典33节点系统带有5条联络开关支路平时断开保持辐射状运行故障或倒负荷时可以闭合形成弱环网。前推回代法处理环网很麻烦需要在断点处做补偿迭代NR法完全没这个问题。因为节点法建模的本质就是用节点电压方程描述整个网络不管支路怎么连接最终都归并到电导矩阵G里。闭合一条联络开关就是在branch表里增加一行[8 21 R]重新装配G迭代逻辑一个字都不用改。这也是为什么我在程序里坚持用“支路表驱动”的原因——拓扑变了数据变算法不动。如果想对比辐射状和合环两种运行方式的差异可以分两次计算分别输出末端节点电压和全网网损。通常合环后某些末端节点电压会抬升网损也可能因为潮流重新分配而变化但要注意合环带来的环流问题那是工程层面另外需要考虑的事。5.2 从恒功率到ZIP负荷模型恒功率负荷假设很常用但实际配网里大量负荷是恒阻抗和恒电流特性混合的。ZIP模型把负荷功率拆成三部分P_L(V) P_0 × [a b×(V/V_0) c×(V/V_0)²]其中a、b、c分别是恒功率、恒电流、恒阻抗占比三者之和为1。这个模型改动程序时P_spec不再是常数而是电压的函数。于是雅可比矩阵也要跟着改因为不平衡量本身对V有依赖J_total ∂P_calc/∂V - ∂P_spec/∂V具体来说恒电流部分的功率对电压导数是b×P_0/V_0恒阻抗部分是2×c×P_0×V/V_0²这些项需要加到雅可比矩阵的对角元上符号方向要看注入还是吸收。对于恒功率部分导数为0代码可以原样保留。这也是为什么我在程序里把P_spec单独放在一个向量里而不是写死在功率方程内部——后续扩展ZIP模型时只需要把“计算P_spec”和“修正雅可比对角元”两个地方改动即可其余迭代骨架完全不动。5.3 分布式电源接入与下垂控制节点光伏、储能接入直流配电网时在潮流模型里就是负的恒功率负荷直接在P_spec对应节点上加上电源注入功率即可。如果有多个换流器采用下垂控制V V_ref - k×P这些节点既不是严格的平衡节点也不是严格的恒功率节点它们的处理方式是在迭代方程里保留这些节点但在雅可比矩阵对角元中加入下垂系数带来的修正项。这样电压和功率会在迭代中自动权衡最终稳定在下垂曲线上。这类控制在直流微电网里非常常见也是直流配电网相比交流配电更有优势的场景——多条馈线通过不同换流器并联供电各自按下垂比例分担功率不需要像交流系统那样纠结无功电压控制。把这三个扩展方向都实现一遍这个程序基本就能应付大部分直流配电网潮流研究场景了。我在几个不同的直流算例上都跑过这套代码最大的体会是NR法本身非常简单真正决定程序质量的是数据组织方式和边界条件处理是否清晰。先把微型算例调到每一步都能看懂再往33节点、环网、ZIP模型逐步扩展比一开始就堆功能踏实得多。 
长河君 · 2026/10/13 11:07:52 · 阅读 1.2k 直流配电网潮流计算：牛顿拉夫逊法的MATLAB实现与33节点算例
做直流配电网的仿真时最头疼的往往不是模型本身而是找不到一套真正顺手、能看清每一步原理的潮流计算程序。交流配电网的潮流程序一抓一大把NR法、前推回代、PQ分解法全都有现成的可一旦把模型换成纯直流无功、相角、频率全都不存在了那些程序里一大半逻辑变成了摆设。实际上直流配电网的潮流计算比交流版本要干净得多——节点电压是唯一的状态变量牛顿拉夫逊法在这里只需要维护一张(n-1)阶的雅可比矩阵。这篇文章围绕经典33节点辐射状配电网算例的直流化改造把NR法的公式推演、MATLAB程序实现、收敛性分析和工程扩展完整拆开讲一遍适合做直流配电、直流微电网仿真或者课程设计的同学直接拿去改。1. 直流化到底改了啥从经典33节点算例到一张纯电阻网络1.1 交流潮流方程与直流潮流方程的本质差异先看交流系统里的节点功率方程。对任意节点i注入功率可以写成两部分有功部分P_i V_i × Σ_j [G_ij×cos(θ_i-θ_j) B_ij×sin(θ_i-θ_j)] × V_j无功部分Q_i V_i × Σ_j [G_ij×sin(θ_i-θ_j) - B_ij×cos(θ_i-θ_j)] × V_j这里有电压幅值V、相角θ、有功P、无功Q四个量所以交流NR法的雅可比矩阵是4个分块矩阵的组合分别处理P对θ、P对V、Q对θ、Q对V的偏导数计算规模直接翻倍。而且交流潮流还要考虑PV节点的无功越限、变压器分接头位置、线路充电电容一堆附加因素。直流配电网稳态运行时没有无功功率没有频率也没有相角差的概念。换流器控制直流电压和输送有功线路和负荷的数学模型只剩电阻。节点电流方程退化成I G × V节点注入功率则是P_i V_i × I_i V_i × Σ_j (G_ij × V_j)这里G就是纯电导矩阵没有B虚部。这意味着NR法的未知量只剩一个维度——节点电压幅值迭代方程的规模直接砍掉一半以上。这就是为什么直流配电网潮流程序能写得很精简也值得单独写一版而不是拿交流程序硬改。1.2 经典33节点算例的直流化改造思路标题里的“33节点”用的是配电领域非常经典的那个33节点辐射状测试系统原始参数包含线路电阻、电抗、有功负荷和无功负荷。做直流配电网研究时通常不会重新构建一套全新拓扑而是直接把这个成熟算例做直流化改造拓扑结构完全保留33个节点、32条支路外加5条可操作的联络开关支路线路参数只取电阻列电抗全部置零负荷只保留有功分量无功分量直接丢弃节点1作为平衡节点由电压源换流器维持直流母线电压恒定其余32个节点全部按恒功率节点处理。直流化之后的网络就是一张纯电阻网络电压等级我习惯按直流配电网常用的10kV设定基准功率取1MVA。这样每个节点的负荷功率数量级在0.1~0.4 p.u.之间潮流计算的结果分布比较典型也方便和后续工程数据做对比。1.3 节点类型大幅简化交流潮流里节点分三类平衡节点、PQ节点、PV节点。PV节点需要给定有功和电压幅值无功是待求量这给雅可比矩阵的构建和迭代过程中的无功越限处理增加了不少代码量。直流配电网里没有无功的概念节点类型只剩两种交流节点类型直流对应节点类型给定量待求量平衡节点直流平衡节点VSC恒压控制节点电压幅值注入功率PQ节点恒功率节点P节点注入有功功率节点电压幅值PV节点不存在除非有恒压控制变流器但那本质上是另一个平衡节点——这个简化听着不起眼但实际编码时省掉的事情非常多。不需要处理θ角排序、不需要分块雅可比、不需要做PV→PQ的越限切换迭代矩阵从(2n-2)×(2n-2)直接降到(n-1)×(n-1)逻辑链条短了一大截。顺便提一句如果某个DC/DC换流器或者储能变流器采用恒压控制它在模型里就是又一个“已知电压节点”本质上和平衡节点等价。在代码里处理多平衡节点时只需要把这些节点的电压也固定住把它们从迭代方程中剔除即可后面第5章会说到。2. 牛顿拉夫逊法的直流版推演一个雅可比矩阵说清全部逻辑2.1 功率不平衡方程与迭代修正的符号约定先把直流潮流的NR法完整推一遍。目标是找到一组节点电压V使得每个恒功率节点的注入功率等于给定值。定义功率不平衡量ΔP_i P_spec_i - P_calc_i其中P_spec_i是节点i给定的注入功率负荷节点为负值分布式电源节点为正值P_calc_i是根据当前电压计算得到的注入功率P_calc_i V_i × Σ_j (G_ij × V_j)牛顿法的核心思路是在当前电压点做一阶泰勒展开求一组电压修正量ΔV让不平衡量归零。写成矩阵形式J × ΔV ΔP这里J就是雅可比矩阵也就是P_calc对电压V的偏导数矩阵。迭代更新公式为V V ΔV很多人在这一步容易把符号搞反。常见的一种教材写法是定义F P_calc - P_spec然后解J×ΔV -F我这里定义的是ΔP P_spec - P_calc所以修正方程里没有负号。两种写法最终等价但一开始就必须固定下来否则程序里一会儿收敛一会儿发散最后还以为是算法问题。2.2 雅可比矩阵的解析表达式不对称才是正常状态对P_calc_i V_i × (G(i,:)·V)逐个求偏导。对角元∂P_calc_i / ∂V_i G(i,:)·V V_i × G_ii这里有项要特别注意G(i,:)·V本身含有V_i因为求和里包含G_ii×V_i这一项所以对V_i求导时会产生两项——一项来自当前电压点处的网络注入项另一项来自系数G_ii与V_i的乘积。非对角元∂P_calc_i / ∂V_j V_i × G_ij其中 i ≠ j非对角元的推导就简单得多G(i,:)·V对V_j求导得到系数G_ij再乘上前面的V_i因子即可。有一点需要单独强调直流网络的电导矩阵G本身是对称的G_ij G_ji但雅可比矩阵J并不对称。因为J(i,j) V_i×G_ij而J(j,i) V_j×G_ji V_j×G_ij两者只有在V_i V_j时才相等。迭代过程中各节点电压一般都不相等所以J基本上不会保持对称。如果用某些对称矩阵求解算法或者把J误写成对称矩阵结果一定不对。2.3 收敛判据与整体迭代流程收敛判据一般用功率不平衡量的无穷范数即所有节点ΔP_i的最大绝对值。我习惯设置为1e-8 p.u.对应到1MVA基准下就是0.01W量级对潮流计算来说精度足够且不会让迭代无谓地反复。完整迭代流程可以概括为初始化所有电压为1.0 p.u.平启动平衡节点电压置为给定值用当前V计算每个节点的P_calc计算ΔP向量取PQ节点的子集判断max(|ΔP|)是否小于收敛容差满足则结束构造雅可比矩阵J只取PQ行、PQ列求解修正方程ΔV J\ΔP更新V V ΔV回到第2步。交流NR法通常需要十几次迭代甚至几十次直流版由于没有相角变量问题规模小、非线性程度也低正常算例基本5~8次迭代就能收敛。如果超过15次还在动就要怀疑数据或者初值出问题了这一点后面专门讲。3. MATLAB程序分段拆解先跑通小算例再装进33节点拓扑3.1 数据组织表驱动而不是写死在公式里程序的第一件事是把网络数据和负荷数据组织成可替换的表。我最推荐的做法是用单独的数组维护拓扑和负荷这样换一个算例只需要改数据不需要动算法代码。负荷数据用一行数组元素顺序对应节点1到33% 节点有功负荷单位 kW顺序对应节点 1~33 % 直流化改造后按经典33节点算例常用参数整理 Pload_kW [0, 100, 90, 120, 60, 60, 200, 200, 60, 60, 45, 60, 60, 120, 60, 60, 60, ... 90, 90, 90, 90, 90, 90, 420, 420, 60, 60, 60, 120, 200, 150, 210, 60];支路数据用三列矩阵格式是[首端节点号末端节点号电阻值]。注意直流化之后没有电抗列这是和交流程序最大的数据差异% 支路数据[首端, 末端, 电阻欧姆] % 完整33节点支路电阻按经典算例录入这里展示前6条作为格式示例 branch [ 1 2 0.0922 2 3 0.4930 3 4 0.3660 4 5 0.3811 5 6 0.8190 6 7 0.1872 % ... 继续扩到全部支路和联络开关 ];基准值设置S_base 1e6; % 基准功率 1 MVA V_base 10e3; % 基准电压 10 kV Z_base V_base^2 / S_base; % 基准阻抗这里等于 100 欧姆这个基准值的好处是负荷功率都在0.04~0.42 p.u.的合理区间电压结果读起来直观。有特殊需求的话把V_base改成±10kV对称双极系统对应的极间电压即可其他逻辑不用动。3.2 电导矩阵形成与NR迭代主循环电导矩阵G的装配逻辑很简单每一条支路电阻r对应电导g 1/r自导纳加上g互导纳减去g。这才是真正理解“节点电压方程I G×V”的关键——G的每一行描述的是该节点与全网其他节点的电流电压关系。%% 形成电导矩阵标幺值 n_node 33; G zeros(n_node, n_node); for k 1:size(branch, 1) i branch(k, 1); j branch(k, 2); r_pu branch(k, 3) / Z_base; g 1 / r_pu; G(i, i) G(i, i) g; G(j, j) G(j, j) g; G(i, j) G(i, j) - g; G(j, i) G(j, i) - g; end G sparse(G); %% 生成节点注入功率标幺值 P_spec -Pload_kW / S_base; % 负荷为负注入 slack 1; P_spec(slack) 0; % 平衡节点功率由潮流解算后得到 %% 初始化与迭代参数 V ones(n_node, 1); V(slack) 1.0; % 平衡节点电压给定 PQ setdiff(1:n_node, slack); % 恒功率节点集合 tol 1e-8; max_iter 50; iter 0; %% 牛顿拉夫逊迭代主循环 while iter max_iter % 计算当前电压下的节点注入功率P_calc V .* (G * V) P_calc V .* (G * V); % 功率不平衡量 dP_all P_spec - P_calc; dP dP_all(PQ); % 收敛判断 if max(abs(dP)) tol break; end % 构造雅可比矩阵 J dP_calc / dV nPQ length(PQ); J zeros(nPQ, nPQ); for ii 1:nPQ i PQ(ii); % 对角元 J(ii, ii) G(i, :) * V V(i) * G(i, i); for jj 1:nPQ if jj ii continue; end j PQ(jj); % 非对角元 J(ii, jj) V(i) * G(i, j); end end % 求解修正方程并更新电压 dV J \ dP; V(PQ) V(PQ) dV; iter iter 1; end这一段代码是整个程序的骨架。P_calc V .* (G * V)这一步是MATLAB向量化运算的经典写法G*V先得到各节点的注入电流向量再逐点乘以电压得到注入功率简洁也能避免写内层循环。构造雅可比矩阵的两重循环对33节点来说完全够快因为矩阵只有32×32。如果你要扩展到几百上千节点可以再用向量化技巧或者稀疏矩阵存储优化但日常研究阶段完全没必要折腾。3.3 先跑通一个微型算例再放大到33节点直接上33节点数据调试有个问题一旦结果不对你很难分辨是算法错、数据错还是某个联络开关没处理对。我的习惯永远是先构造一个最小的可手算网络把算法里的每一个符号、每一个数值方向验证清楚再放大到目标算例。最小算例长这样平衡节点1电压1.0 p.u.节点2通过一条电阻R 100欧姆的线路接到节点1节点2带100kW恒功率负荷。基准阻抗Z_base 100欧姆所以R的标幺值正好是1.0100kW对应标幺值0.1。节点的功率方程是P_calc_2 V_2 × (G_22×V_2 G_21×V_1) V_2² - V_2因为这里G_22 1G_21 -1V_1 1。负荷100kW是吸收功率所以注入功率P_spec_2 -0.1。方程就成了V_2² - V_2 -0.1解出来有两个根V_2 0.8873 或 V_2 0.1127。潮流计算取高电压解也就是0.8873这是恒功率负荷下直流网络的正常运行点。用上面的主程序把n_node改成2branch填[1 2 100]Pload_kW填[0 100]从V_2 1.0平启动开始迭代前几步就会快速逼近0.8873。如果程序迭代轨迹和这个对不上说明符号约定或者雅可比构造有bug。这比盯着33节点网络找bug高效得多。3.4 结果输出与功率平衡校验迭代结束后的输出部分除了打印节点电压我还会额外做两项校验平衡节点注入功率和全网损耗这是判断结果是否可信的硬指标。%% 结果输出与校验 P_bal V(slack) * (G(slack, :) * V); % 平衡节点注入功率标幺 P_loss V * G * V; % 全网电阻损耗标幺 P_total_sp sum(P_spec); % 所有非平衡节点注入之和 fprintf(迭代次数%d\n, iter); fprintf(平衡节点注入功率%.4f p.u.%.2f kW\n, P_bal, P_bal * S_base / 1e3); fprintf(全网网损%.4f p.u.%.2f kW\n, P_loss, P_loss * S_base / 1e3); fprintf(功率平衡校验%.6f\n, P_bal P_total_sp - P_loss);功率平衡校验的物理意义是平衡节点注入功率加上其它所有节点注入功率的总和应该等于网络电阻损耗。因为所有节点注入功率之和等于V × G × V也就是支路电阻上的总损耗。校验值在1e-6以下基本说明结果自洽。节点电压分布可以直接画柱状图figure; bar(1:n_node, V, FaceColor, [0.3 0.5 0.7]); grid on; xlabel(节点编号); ylabel(电压 (p.u.)); title(33节点直流配电网潮流电压分布);33节点辐射状直流配电网的电压分布通常有一个明显的趋势离平衡节点越远电压越低线路越长、负荷越重的地方压降越明显。末端节点如果电压掉到0.93 p.u.以下就该考虑是不是负荷太重或者线径选细了。4. 收敛性分析恒功率负荷的临界负载与那些“发散得莫名其妙”的时刻4.1 恒功率负荷为什么会导致潮流无解直流配电网里恒功率负荷是一个很微妙的存在。它吸收的功率不随电压变化电压降低时电流反而增大这会给网络带来一种类似正反馈的压力。回到刚才的微型算例方程是V_2² - V_2 P_L 0其中P_L是节点2吸收的有功标幺值取正数。这是一个一元二次方程判别式是1 - 4×P_L。只有当P_L ≤ 0.25 p.u.时才有实数解也就是说这个单负荷直流网络最多只能带250kW基准1MVA。超过这个值潮流方程没有可行解NR法从任何初值出发都会发散。这不是程序bug而是网络的物理极限。恒功率负荷比重越高直流系统越容易出现这种电压稳定问题。实际工程中换流器或储能系统会提供电压支撑但纯潮流计算时如果负荷给得太重程序发散是正常现象得从网络参数和运行点去找原因。4.2 平启动失败时的三个排查方向NR法跑发散做完微型算例验证之后最可能的原因有三个按出现频率排序一是数据录入错误。支路首末端节点号写反、电阻数量级差10倍、负荷单位没从MW换算成kW这些错误会让雅可比矩阵方向没问题但解根本不存在。我见过最隐蔽的一次是联络开关支路的电阻填成了0.0001欧姆相当于把两个节点直接短接结果电压被拉得乱七八糟。二是初值离解太远。直流NR法平启动V1.0通常没问题但重负荷或者高阻抗网络下解的电压可能远低于1.0平启动迭代容易跑到低电压解或者直接发散。这时可以把初值改到0.9甚至0.8试试也可以先跑恒阻抗负荷模型得到一组相对稳定的电压把它当作恒功率模型的初值。三是负荷总量超过了网络承载能力。前面说的临界负载问题在网络规模变大后依然存在只是没有一个简单的判别式可以直接算出来。此时需要逐步增加负荷观察哪个点开始发散这实际上就是连续潮流法的雏形。4.3 每次迭代都记录最大不平衡量调程序时我会在循环里记录每一步的max(abs(dP))并打印出来这比只看到最终“不收敛”三个字有用得多。正常收敛时这个值是一条近似单调快速下降的曲线通常从10的-2次方量级一路掉到1e-8以下。如果看到曲线先下降几步又反弹或者在一个数值附近来回震荡多半是雅可比矩阵构造有误重点检查对角元是不是漏了G(i,:)*V这一项。如果曲线完全发散每一步都比上一步大先检查符号方向再看初值。另外求解修正方程时用J\b而不是inv(J)*b。MATLAB的反斜杠运算符自动选择合适的高斯消元求解方法数值稳定性比显式求逆好速度也更快。虽然32×32矩阵用哪种都能算出来但好的习惯可以避免以后扩展到大规模网络时踩坑。5. 三个趁手的扩展方向合环运行、ZIP负荷、分布式电源与下垂控制5.1 直接闭合联络开关NR法天然支持环网经典33节点系统带有5条联络开关支路平时断开保持辐射状运行故障或倒负荷时可以闭合形成弱环网。前推回代法处理环网很麻烦需要在断点处做补偿迭代NR法完全没这个问题。因为节点法建模的本质就是用节点电压方程描述整个网络不管支路怎么连接最终都归并到电导矩阵G里。闭合一条联络开关就是在branch表里增加一行[8 21 R]重新装配G迭代逻辑一个字都不用改。这也是为什么我在程序里坚持用“支路表驱动”的原因——拓扑变了数据变算法不动。如果想对比辐射状和合环两种运行方式的差异可以分两次计算分别输出末端节点电压和全网网损。通常合环后某些末端节点电压会抬升网损也可能因为潮流重新分配而变化但要注意合环带来的环流问题那是工程层面另外需要考虑的事。5.2 从恒功率到ZIP负荷模型恒功率负荷假设很常用但实际配网里大量负荷是恒阻抗和恒电流特性混合的。ZIP模型把负荷功率拆成三部分P_L(V) P_0 × [a b×(V/V_0) c×(V/V_0)²]其中a、b、c分别是恒功率、恒电流、恒阻抗占比三者之和为1。这个模型改动程序时P_spec不再是常数而是电压的函数。于是雅可比矩阵也要跟着改因为不平衡量本身对V有依赖J_total ∂P_calc/∂V - ∂P_spec/∂V具体来说恒电流部分的功率对电压导数是b×P_0/V_0恒阻抗部分是2×c×P_0×V/V_0²这些项需要加到雅可比矩阵的对角元上符号方向要看注入还是吸收。对于恒功率部分导数为0代码可以原样保留。这也是为什么我在程序里把P_spec单独放在一个向量里而不是写死在功率方程内部——后续扩展ZIP模型时只需要把“计算P_spec”和“修正雅可比对角元”两个地方改动即可其余迭代骨架完全不动。5.3 分布式电源接入与下垂控制节点光伏、储能接入直流配电网时在潮流模型里就是负的恒功率负荷直接在P_spec对应节点上加上电源注入功率即可。如果有多个换流器采用下垂控制V V_ref - k×P这些节点既不是严格的平衡节点也不是严格的恒功率节点它们的处理方式是在迭代方程里保留这些节点但在雅可比矩阵对角元中加入下垂系数带来的修正项。这样电压和功率会在迭代中自动权衡最终稳定在下垂曲线上。这类控制在直流微电网里非常常见也是直流配电网相比交流配电更有优势的场景——多条馈线通过不同换流器并联供电各自按下垂比例分担功率不需要像交流系统那样纠结无功电压控制。把这三个扩展方向都实现一遍这个程序基本就能应付大部分直流配电网潮流研究场景了。我在几个不同的直流算例上都跑过这套代码最大的体会是NR法本身非常简单真正决定程序质量的是数据组织方式和边界条件处理是否清晰。先把微型算例调到每一步都能看懂再往33节点、环网、ZIP模型逐步扩展比一开始就堆功能踏实得多。