三相不平衡潮流计算在配电网分析里是个绕不开的话题。最近我把自己一直在用的那套基于前推回代法与互阻抗模型的Matlab程序重新整理了一遍顺手把容易出错的细节和调试心得也记录了下来。这篇文章适合正在做配电网三相潮流仿真、或者打算自己写三相不平衡分析程序的同学它不只是一份代码说明更是一份从物理模型到代码实现的完整施工记录。为什么绕不开实际台区里单相负荷大量存在A相和C相可能挂几十户B相却几乎没有负载线路也因为不换位而参数不对称A相电流的变化会通过相间互感传递到其他两相。用简化单相模型去估算出来的末端电压和中性点偏移经常和现场实测对不上。我最早就是用单相前推回代后来换到三相互阻抗模型误差才真正降下来。这套程序的核心逻辑是把三相节点阻抗、前推回代迭代和收敛处理完整串起来最后能输出各相电压、支路电流和线损方便后续做低电压治理和三相不平衡度评估。文章里所有思路都基于Matlab环境实现代码块里的函数我尽量做了命名和注释可以按需裁剪。1. 三相不平衡问题为什么必须单独建模1.1 配电网的三相不平衡都从哪来配电系统的三相不平衡从来不是假设出来的而是实测得出来的。输电网讲究三相换位线路参数对称负荷也是三相对称所以经典潮流算法里可以直接用正序网络替代。配电网完全不是这回事低压台区从母线出来就是四线制或三线制馈线大量支路是单相两线用户接在哪一相基本看施工习惯。台区负荷越细单相比例越高三相不平衡度就越明显。除了负荷分布还有几个容易被忽略的来源。单相光伏逆变器、电动自行车充电桩这类设备本质上都是单相负荷接入相位如果没人管理很容易出现A相满负荷、C相接近空载的局面。我再举个实际现象晚上7点高峰时段某台区三相电流可能分别是180A、60A、150A中性线电流能到40A以上。这个电流在中性线上会产生电压降直接抬高中性点电位导致某一相电压偏低、另一相电压偏高。低电压用户投诉往往就是这么来的。1.2 单相等值模型在哪些场景下会失效单相模型通常的做法是把三相总负荷平均到一相上或者用总功率除以三只计算一相的电压降。它隐含了两个关键假设三相完全对称、相与相之间没有电磁耦合、中性点电位始终为零。实际配电线路这三条基本都不成立。拿一个典型场景说某条四线制馈线A相重载B相轻载。用单相模型看B相电压得到的结论是“很健康”但实际中A相重载电流流过线路时会在B相导线上感应出互感电压如果这个感应电压的方向与B相自身电流产生的电压降方向一致B相电压会被明显拉低甚至出现三相电压不平衡度超标。这种相间耦合效应单相模型完全看不到。中性线电流引起的中性点偏移也只有在三相模型中才会显现。单相模型天然认为中性点就是大地的零电位但实际中性线有阻抗、有电流中性点电位会漂移几伏甚至十几伏。研究单相接地、缺相运行、零序电流分布时单相模型给出的结果基本没有参考价值。1.3 哪些分析场景必须上三相模型我自己的经验是只要工作目标里出现下面任意一个关键词就必须换用三相模型电压不平衡度评估。国标要求电网公共连接点电压不平衡度不超过2%多数低压台区在高峰期早超了这个指标单相模型算不出来。中性线电流分析。校验中性线载流量、设计无功补偿和三相负荷调整方案时需要知道每一段中性线上的电流单相模型会直接无视它。单相分布式光伏接入评估。接了哪些相、每相多少容量都会影响三相电压只有三相潮流能给到每相的结果。低电压治理方案设计。如果只看单相结果会低估A相的电压跌落治理措施自然也会偏保守或选错相位。台区线损精细分摊。单相模型算出来的网损会漏掉中性线损耗换三相模型后往往发现实际损耗比原来估算高8%以上。如果只是粗略估算变压器低压侧总电流单相模型还能顶一阵一旦涉及逐节点、逐相的电压和损耗分析单相模型得出的结论几乎没法用于现场决策。这也是我坚持把互阻抗模型写进程序的原因。2. 前推回代法在配电网潮流中的定位与核心思路2.1 为什么用前推回代而不是牛顿法配电网潮流计算最常用的两类算法一类是牛顿类一类是前推回代。我优先选前推回代原因有三个。第一配电网辐射状结构为主支路多、深度大而且R/X比普遍偏高。传统牛顿法需要反复求解雅可比矩阵在R/X比大时矩阵病态严重迭代容易震荡或者直接发散。前推回代不需要构造雅可比矩阵也不做任何矩阵求逆只做两次网络遍历数值稳定性要好得多。第二前推回代的物理意义非常直白。从末端负荷倒推支路电流再用支路电流从首端往前推电压整个过程和实际功率流动方向一致。出现收敛问题时可以很直观地判断是哪个节点、哪条支路的数据出了问题。第三三相模型下每个节点有三个相电压雅可比矩阵规模直接翻三倍。前推回代虽然也需要处理3×1复电压向量和3×3阻抗矩阵但没有雅可比矩阵的构造和分解内存开销和计算量都小一个量级。不过要注意前推回代只适用于辐射状网络或者少量合环经过破环处理的网络。如果网络里有明显的环网就需要改用牛顿法或者在断开处叠加补偿电流不能硬套前推回代。2.2 前推回代的计算顺序与公式前推回代的基本流程可以拆成三块网络分层。把根节点作为第0层按照父节点关系给每个节点确定层号同时记录每个节点的父子关系和分支编号。前推。从最深层节点开始逐层向上先计算当前节点的三相注入电流再累加其所有下层子支路的电流得到离开本节点往根方向流动的支路电流。回代。从第0层根节点开始逐层向下用已经算好的支路电流和对应支路的三相阻抗矩阵计算电压降落更新子节点电压。核心公式不复杂。恒功率负荷在第k个节点的三相注入电流为I_inj_abc conj(S_abc ./ V_abc)其中S_abc是节点三相复功率向量V_abc是当前迭代的三相复电压向量两者都是3×1。注意这里必须按相做复数除法如果某一相没有负荷S对应为0电流自然就是0。支路电流的前推公式I_branch(k) I_inj_abc(k) sum(I_branch(children(k)))也就是说父节点往根方向流过的电流等于本节点自身负荷电流加上所有子支路电流之和。电压回代公式V_child V_parent - Z_abc_branch * I_branch这里的Z_abc_branch是父节点到子节点这一段线路的3×3相分量阻抗矩阵乘以3×1支路电流向量得到三相电压降落向量。由于互阻抗的存在A相的电压降落不只由A相电流决定还要叠加B相和C相电流通过互阻抗产生的贡献所以必须用矩阵乘法拆成三个单相独立计算就前功尽弃了。2.3 三相场景下几个容易被忽略的变量前推回代切到三相后有几个点特别容易踩空。第一是负荷的相别。写程序时我习惯把Sload存成n×3的复数矩阵每一行对应一个节点列对应A/B/C相。这样B相没有负荷时S0计算电流自动就是0逻辑非常清晰。千万不能把三相负荷压缩成一个总功率再按平均值分配到各相那就又退回单相模型了。第二是平衡节点。通常取配电变压器低压侧母线作为平衡节点三相电压设定为额定相电压的正序值也就是幅值相同、相位互差120度。但实际现场变压器低压出口三相电压往往已经存在一定偏移所以我的程序允许直接用实测电压作为初值。第三是初值问题。所有非电源节点初始电压都取额定电压三相各自带120度相位差。这个初值对常规辐射状网络是够用的但如果系统比较重载前几次迭代会出现比较明显的震荡后面我会专门讲怎么处理。3. 互阻抗模型如何落地从物理矩阵到程序数据3.1 三相线路阻抗矩阵的物理意义三相线路的参数不能只用一组正序、零序阻抗来代替。三相导线之间存在电磁耦合A相流过的电流不仅在本相产生电压降还会通过互感在B相、C相上产生电压降。要精确描述这种耦合最自然的方式就是用3×3复阻抗矩阵Zabc [ Zaa Zab Zac Zba Zbb Zbc Zca Zcb Zcc ]对角线上的Zaa、Zbb、Zcc是各相导线的自阻抗包括导线内阻抗和大地回路的返回阻抗。非对角线上的Zab、Zbc、Zca是相间互阻抗主要由导线之间的几何均距决定。由于配电线路一般不做换位三相导线布置不对称对角线元素并不相等互阻抗之间也有差异矩阵不能简化成完全对称的循环矩阵。这里面还有一层容易混淆的地方。互阻抗本质上是互易的所以物理上Zab应该等于Zba。但我们用Carson公式或者其他电磁参数计算工具得到的原始矩阵可能因为大地回路处理方式、分相段编号等细节出现微小数值不对称。遇到这种情况程序里最好取上下三角平均避免后面叠代时出现电压结果不对称。3.2 是否包含中性线Kron化简是捷径三相四线制低压台区实际上有四根导线严格的相分量模型应该是4×4矩阵包含三相导线和中性线。如果直接把中性线作为独立变量每个节点就有A/B/C/N四个电位程序结构和数据规模都会复杂不少而且中性线的接地方式、重复接地等问题处理起来特别容易出错。我的做法是优先采用Kron归减法把中性线消去只保留三相导线的等效阻抗矩阵。思路很简单先把4×4阻抗矩阵按三相导线和中性线分块Z_full [ Z_abc_self Z_an Z_na Z_nn ]其中Z_abc_self是3×3的相导线自互阻抗分块Z_an是3×1的中性线互阻抗列向量Z_na是它的转置Z_nn是中性线自阻抗。消去中性线后三相导线的等效阻抗矩阵为Zabc_eq Z_abc_self - Z_an * inv(Z_nn) * Z_na这个公式的物理解释很直观中性线电流并不是零但它流经中性线阻抗而后又通过互阻抗耦合回相导线最终效果可以折算成相导线之间的一个等效互阻抗。把中性线“吞”进3×3矩阵后程序就不再需要单独算中性线电压但中性线电流引起的电压降影响已经被包含进去了。不过要记住Kron化简成立的前提是中性线回路唯一。如果一条馈线有多个重复接地或多段中性线通过大地并联简单消去中性线就不够准确了。这种情况我建议显式建模中性线节点保留4×4矩阵直接把中性线也纳入前推回代的计算范围。3.3 程序里的数据组织与生成方式在Matlab里每个分支的三相阻抗矩阵我习惯用cell数组存储例如Zabc_branch{b}就是一个3×3复数矩阵对应第b条支路的等效相分量阻抗。手头没有精确的Carson计算工具时可以按典型参数近似。比如一段LGJ-120导线参数大致如下Zaa 0.27 j0.347 Ω/km Zbb 0.27 j0.347 Ω/km Zcc 0.27 j0.35 Ω/km Zab 0.08 j0.17 Ω/km Zbc 0.08 j0.18 Ω/km Zca 0.08 j0.17 Ω/km然后根据线路长度直接相乘function Zabc build_Zabc(r_self, x_self, r_mut, x_mut, length_km) Zabc diag(r_self 1j*x_self); Zabc(1,2) r_mut(1) 1j*x_mut(1); Zabc(2,1) Zabc(1,2); Zabc(2,3) r_mut(2) 1j*x_mut(2); Zabc(3,2) Zabc(2,3); Zabc(1,3) r_mut(3) 1j*x_mut(3); Zabc(3,1) Zabc(1,3); Zabc Zabc * length_km; end这里只放了一个简化版本实际工程里互阻抗最好用导线几何参数计算否则算出来的中性点偏移和实测偏差会很明显。但我先用典型参数把程序跑通再逐步替换精确数据这个路径是可行的。4. Matlab程序实现拓扑、迭代与收敛控制4.1 拓扑输入与分层处理写前推回代程序第一步不是写迭代公式而是把网络的父子关系理清楚。我用两个核心数组parent(k)节点k的父节点编号根节点父节点设为0。children{k}节点k的所有子节点编号的cell数组。有这两个结构就可以从根节点做广度优先遍历得到节点深度。示例函数如下function depth getDepth(parent, root) n length(parent); depth zeros(1, n); children cell(1, n); for k 1:n if parent(k) 0 children{parent(k)} [children{parent(k)}, k]; end end queue root; while ~isempty(queue) cur queue(1); queue(1) []; for ch children{cur} depth(ch) depth(cur) 1; queue [queue, ch]; end end end注意这里用队列做BFS不要用find(parentcur)逐次扫描整个数组节点一多会非常慢。我最早写第一版程序就是吃了这个亏30个节点还能忍到100个节点的台区就明显卡顿。在迭代过程中还需要一个按深度分层的节点列表。简单做法是用nodesByDepth{d}来存第d层的所有节点这样前推和回代都可以按层遍历逻辑清楚也不会漏节点。4.2 前推回代主循环核心迭代代码我精简成一个可以对照实现的结构V repmat(Vbase, n, 1); V V .* (ones(n,1) * [1 exp(-1j*2*pi/3) exp(1j*2*pi/3)]); Ibranch zeros(n, 3); for iter 1:maxIter % 前推从最深一层往上 for d maxDepth:-1:1 for k nodesByDepth{d} if parent(k) 0 continue; end Iinj conj(Sload(k, :) ./ V(k, :)); Ibranch(k, :) Iinj; for ch children{k} Ibranch(k, :) Ibranch(k, :) Ibranch(ch, :); end end end % 回代从根节点往下 V_old V; for d 1:maxDepth for k nodesByDepth{d} for ch children{k} dV (Zabc_branch{ch} * Ibranch(ch, :).).; V(ch, :) V_old(k, :) - dV; end end end % 收敛判断 if max(abs(V(:) - V_old(:))) tol break; end end里面有几个细节要特别说明。Sload(k,:)是节点k的三相复功率向量B相没有负荷就置0。conj(Sload ./ V)是恒功率负荷电流公式。电压初值我用了一个1×3的复数向量A相为1∠0°B相为1∠-120°C相为1∠120°再乘到每一行。这样比直接设成1∠0更接近真实工况迭代能少跑好几轮。回代的时候子节点电压更新用的是父节点上一轮的值V_old不是本轮已经更新的值。这一点很关键否则网络顺序会让某些节点的更新路径受同一轮内其他节点影响产生不必要的耦合误差。虽然最终都能收敛但会降低稳定性尤其是重载网络下容易出现振荡。4.3 收敛判据与迭代加速收敛判据我习惯用电压幅值偏差最大值小于1e-6标幺值。三相系统里三个相电压都要检查所以直接用max(abs(V(:)-V_old(:)))取全局最大值简单可靠。如果发现迭代很慢或震荡可以用阻尼因子V alpha * V (1 - alpha) * V_old;alpha可以设成0.5到0.9具体数值看网络重载程度。重载网络alpha取小一点反而稳。我一般先用0.8跑不收敛再降到0.5。另一个实用技巧是不要每次迭代都从初始平坦启动。可以先用上一轮收敛电压作为初值只做一次前推回代多迭代几轮。实际测试下来对同一网络这种热启动可以减少约30%的迭代次数。4.4 代码效率与Matlab实现细节三相前推回代看起来简单但如果不注意实现细节数据一多还是会慢。预先构建children和nodesByDepth避免在循环里用find。Ibranch和V都预先分配为zeros(n,3)不要在循环里动态扩展矩阵。矩阵乘法Zabc_branch{ch} * Ibranch(ch,:)会构造临时向量如果节点非常多可以考虑用循环嵌套计算三相电压降虽然代码长一点但内存访问更友好。用parfor并行化不同馈线时要注意Matlab的parfor里访问cell数组容易报错最好把线路阻抗矩阵改成普通三维数组或拆分变量。这些优化点在前推回代里属于“锦上添花”但如果程序要反复批量跑几个百节点的台区性能差距会非常明显。5. 算例验证典型三相不平衡馈线的结果对比5.1 测试算例参数我拿一个10节点低压台区做验证。变压器容量400kVA额定电压10/0.4kV低压侧线电压400V相电压230V。线路全部采用LGJ-120导线单位长度参数按上面3.3节给的数据。第1节点是电源节点后面9个节点带负荷。负荷分配刻意设计成严重不平衡A相重载B相基本空载C相中等负载。节点2到10的单相负荷数据如下节点A相(kVA)B相(kVA)C相(kVA)225j10020j8320j8015j6418j7012j5522j9014j6615j6010j4716j6012j5812j508j3910j406j2108j305j2各段线路长度依次为40m、50m、45m、55m、50m、40m、60m、45m、50m全部按低压四线制敷设中性线重复接地一次这个算例我按Kron化简后的三相四线模型计算。5.2 三相模型与单相模型的结果对比把同一个负荷数据折算成三相总功率后用单相前推回代算一遍得到的是统一的“平均电压”。三相互阻抗模型则会输出三个相电压。取几个代表节点作对比节点三相模型Va(pu)三相模型Vb(pu)三相模型Vc(pu)单相模型V(pu)20.9780.9950.9880.98750.9480.9820.9690.97380.9190.9720.9480.952100.8920.9650.9370.946很明显单相模型给出的末端电压0.946pu看起来“勉强可以接受”但三相模型显示A相末端电压只有0.892pu这已经是比较严重的低电压问题。如果治理方案只按单相模型来设计补偿容量和调整相别都会偏离实际需求。同时注意B相电压也比单相模型低虽然B相本身没有挂负荷。这就是互阻抗耦合的作用A相重载电流通过互阻抗在B相导线上产生了额外的电压降单相模型完全无视了这个物理过程。5.3 不平衡度和网损分析用对称分量法把末端节点三相电压分解得到负序分量和正序分量之比约为3.1%超过GB/T 15543规定的2%限值。也就是说这个台区在高峰期不仅低电压而且三相不平衡已经超标治理时必须做三相调补不能简单单相补无功。网损方面三相模型算出来的总损耗除了三根相线损耗外还包含中性线损耗。因为中性线电流比较大这部分损耗并不小。单相模型因为只计算一相且假设三相平衡算出来的总网损明显偏低。本例中三相模型总网损比单相模型高约9%这个差值主要就是中性线电流和互阻抗附加损耗贡献的。这个结果也说明一个道理只做单相潮流不只是电压结果偏差连线损分析都会低估。台区线损治理项目如果拿单相模型做基线很容易在改造后觉得“效果不明显”因为基线本身就不准确。6. 调试过程中踩过的坑与经验总结6.1 相别错位和负荷数据格式我最早调试程序时犯过最蠢的一个错是把负荷表里的A、B、C列顺序读反结果算出来某节点负序电压特别大。排查了半天才发现是列没对齐。这里给个建议程序入口处统一用n×3复数矩阵存储负荷并且把相别标清楚。如果输入数据来自Excel读进来后先打印前几行人工核对一遍。看似多花一分钟实际上能省掉后面好几个小时的排查。6.2 变压器低压侧电压初值不要教条前推回代并不强制要求根节点三相电压严格正序。实际台区里负序电流流过配电变压器时低压侧出口三相电压已经有不小的偏移。如果把根节点电压硬设成完全对称程序前几次迭代会强行把电压“拉回对称”这就跟现场实测对不上了。我的做法是增加一个开关允许手动输入根节点三相电压初值。如果现场有测量数据直接用实测值迭代很快且结果更符合实际。没有实测数据时再退回到标准正序初值。6.3 Kron化简在重复接地场景下会失真前面提到过Kron化简的前提是中性线电流只走过一条明确路径。但低压台区经常有多处重复接地比如变压器中性点接地、末端重复接地、用户侧保护接地。这种情况下中性线电流会分流到大地简单消去中性线的等效3×3矩阵会低估中性线回路阻抗进而导致中性点偏移算得偏小。碰到这种台区我建议别省事直接显式建模中性线节点。程序要额外增加一层电压变量前推回代中把中性线也当作一根“相”来处理虽然复杂但结果更可信。6.4 不收敛时的排查顺序三相前推回代不收敛原因往往比单相版本更多。我一般按这个顺序排查拓扑是否分层正确。重点看parent数组是否有环根节点是否唯一有没有孤立节点。线路长度单位是否一致。单位长度阻抗是Ω/km长度如果写成米电压降会被放大1000倍程序必然发散。负荷功率符号是否搞反。我的算例里所有负荷都为正但程序里加了光伏后光伏节点可能作为负负荷处理。如果光伏注入功率大于负荷支路电流方向会翻转再叠加互阻抗迭代曲线会来回跳。线路阻抗矩阵是否合理。曾遇到过互阻抗参数填得比自阻抗还大程序算出正的自损耗为负这种网络在物理上就不成立。还有一个经验把迭代过程中的最大电压偏差打印出来观察它是单调下降还是振荡。如果是振荡多半是阻尼因子太大或者初值太差如果是快速单调下降后突然反弹往往是某个节点负荷超过了线路传输极限这时候要从负荷数据找问题。6.5 第二层级的扩展建议前面所有讨论都集中在纯辐射状低压馈线。如果要把程序扩展到含配电变压器的中压系统或者接入多台分布式光伏还需要在三相潮流基础上增加变压器绕组模型和逆变器控制模型。变压器三相不平衡时Yyn0和Dyn11绕组的中性点电压处理方式不同不能只靠修改线路阻抗矩阵解决光伏逆变器通常按PQ节点建模仿但电压越限时可能切换成恒电压或恒无功模式这些都需要针对不同控制策略单独处理。我自己目前的做法是先保证基础三相前推回代程序足够可靠把所有线路和负荷数据管理好再逐步往后端功能扩展。主程序保持接口清晰方便把变压器和分布式电源模型作为独立模块加进去。最后再说说个人体会。三相不平衡潮流计算真正难的不是迭代公式而是模型细节——互阻抗归算是否准确、负荷相别是否对齐、变压器中性点怎么处理。程序本身几行代码就能搭起来但要让算例和现场数据对得上功夫全在线路参数和负荷获取上。如果你打算自己复现建议先从一个小台区开始核对每个分段的三相阻抗矩阵再逐步扩大到几十个节点在没验证之前尽量别直接用单相模型评估低电压问题。我还会继续在配电三相潮流这个方向上做一些整理后面有机会再单独写一篇关于配电变压器三相建模和含光伏接入的扩展内容。