做直流配电网的仿真时，最头疼的往往不是模型本身，而是找不到一套真正顺手、能看清每一步原理的潮流计算程序。交流配电网的潮流程序一抓一大把，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.0；100kW对应标幺值0.1。节点的功率方程是：
P_calc_2 = V_2 × (G_22×V_2 + G_21×V_1) = V_2² - V_2
因为这里G_22 = 1，G_21 = -1，V_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改成2，branch填[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法平启动V=1.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模型逐步扩展，比一开始就堆功能踏实得多。
