做交直流混合配电网潮流计算时我最开始用的也是教科书里常见的交替迭代法AC网和DC网各算各的换流器功率在两个网络之间来回传递。但碰到直流母线多、交流侧分布式电源出力波动大的场景迭代经常会走到死胡同要么震荡要么直接发散。后来换成统一求解法再把牛顿-拉夫逊的全流程用Matlab实现才算把这套系统稳稳跑起来。这篇文章就记录这次实现的完整过程包括数学模型、编码细节和几个我花了不少时间才填平的坑。交直流混合配电网不是把交流潮流和直流潮流两套程序拼起来就完事。交流部分的节点电压相角、直流节点的电压幅值、换流器的控制变量彼此之间存在强耦合。这篇文章适合正在做配电网潮流、分布式光伏/储能并网、柔性直流配网研究的同学尤其是想自己写代码而不是只调工具箱的人。看完之后你能理解统一求解法的方程组织方式也能照着把Matlab主程序、雅可比矩阵组装、算例验证这条链路搭起来。1. 交直流混合配电网潮流计算的难点与统一求解法的切入点1.1 交流系统和直流系统的耦合关系到底耦合在哪配电网里的直流部分通常通过电压源型换流器VSC接入交流母线直流母线上再挂光伏、储能和直流负荷。交流侧和直流侧在换流器这个位置并不是简单地“功率相加”而是有一套完整的状态约束交流侧注入换流器的有功功率决定了直流侧母线功率的输入直流侧的电压水平又反过来影响换流器交流侧的无功支撑能力。换句话说交流节点的电压相角、直流节点的电压幅值、换流器内部控制量这三个层面的变量必须放在同一个方程体系里才能互相制约。我举个例子直流母线电压因为光伏出力突然上升而升高VSC为了维持直流电压会调整有功输出方向这个调整会立刻改变交流节点注入功率交流节点电压和相角重新分布而交流侧电压的变化又会让换流器内部电流、损耗以及无功上限全部跟着变。这种双向反馈如果拆开迭代A网络算完得到边界条件传给B网络B网络算完再传回A网络每一步都只能“凭上一次的结果猜下一次”强耦合时很容易来回震荡。这也是我后来坚持用统一求解法的根本原因。1.2 为什么我最终选了统一求解法而不是交替求解法交替求解法的程序结构确实诱人交流潮流用现成的Matpower或牛顿法算完把换流器有功、无功当作直流侧的电源注入直流网络用线性方程G_dc*U_dcI解一遍然后把直流功率反馈到交流侧循环直到偏差小于阈值。问题在于这个方法只适合“直流部分简单、耦合弱”的情况。一旦直流侧有多个VSC、多个控制模式或者交流侧重载、R/X比较高交替迭代的步长很容易绕着一个极限环打转怎么调松弛因子都压不下来。统一求解法的思路是放弃“两个网络分开算”的幻觉直接把交流节点功率平衡方程、直流节点功率平衡方程、换流器控制方程全部写成一组非线性方程组F(x)0然后用牛顿-拉夫逊法同步更新所有状态变量。这个方法代码量确实大一些但收敛性和解的一致性都明显更好。下面是两套方案在我实际测试中的对比对比维度交替求解法统一求解法数学模型AC、DC各自成组边界交换全部联立成一组非线性方程收敛特性弱耦合时稳定强耦合时易震荡牛顿法局部二阶收敛整体稳健编程复杂度较低模块可复用需要组装AC-DC耦合雅可比矩阵适用场景单端直流、连接点少、负荷轻多端直流、VSC控制复杂、DG高渗透对初值敏感度每轮子迭代相对宽容对状态初值更敏感但可用平启动缓解我并不是说交替求解法一无是处。直流系统规模很小、只想快速看一个趋势的时候交替法半小时就能写完。但如果你要为“交直流混合配电网”这个题目做一套能反复使用的分析工具统一求解法能省掉后面无穷无尽的调参时间。2. 统一求解法的数学模型构建2.1 交流节点功率平衡方程交流部分仍然沿用经典极坐标形式的功率平衡方程。对交流节点i[ P_i^{sp} - V_i\sum_{j\in i}V_j(G_{ij}\cos\theta_{ij}B_{ij}\sin\theta_{ij})0 ][ Q_i^{sp} - V_i\sum_{j\in i}V_j(G_{ij}\sin\theta_{ij}-B_{ij}\cos\theta_{ij})0 ]其中 (P_i^{sp}P_{G,i}-P_{L,i})(Q_i^{sp}Q_{G,i}-Q_{L,i})。配电网里负荷如果直接用恒功率模型低压场景下容易让雅可比矩阵对角项偏弱所以我一般会在算例里预留ZIP负荷的转化接口把一部分恒功率负荷按当前电压转换成恒阻抗和恒电流再叠加到(P_i^{sp})上。这样做的好处是收敛半径更大也更贴近实际台区负荷特征。交流节点类型有三种松弛节点平衡节点、PV节点和PQ节点。交直流混合潮流里松弛节点一般选在变电站二次侧母线PV节点对应分布式光伏逆变器或其他可控无功源。换流器的交流侧母线的类型不完全取决于换流器本身还要看换流器采用什么控制模式。比如定交流电压控制的VSC它所在的交流节点的无功注入就不是“给定值”而是方程解出来的这一点在分组方程时必须想清楚。2.2 直流网络方程与换流器模型直流网络在稳态下本质上是一个线性电阻网络[ I_{dc,k} \sum_{m\in k} G_{dc,km}U_{dc,m} ]对于直流节点k注入功率平衡可以写成[ P_{dc,k}^{sp} - U_{dc,k}\sum_{m\in k}G_{dc,km}U_{dc,m}0 ](P_{dc,k}^{sp})是直流侧电源注入减去直流负荷后的净功率。光伏和储能如果通过DC/DC变换器接入直流母线在潮流计算里可以当成直流节点上的注入功率源如果接在交流侧那就按交流节点上的PQ或PV节点处理。换流器模型是统一求解法的灵魂。我使用的是经过简化的VSC稳态模型只关心基波正序忽略谐波和换流站内部电磁暂态把损耗统一归算成有功电流的函数[ P_{loss}a bI_{conv} cI_{conv}^2 ]其中(I_{conv})是流过VSC的直流侧电流(a)为空载损耗(b)为线性损耗系数(c)为二次损耗系数。交流侧有功、直流侧有功和损耗之间满足[ P_{ac,conv}P_{dc,conv}P_{loss} ]无功方面VSC可以向交流母线提供或吸收无功具体数值由控制目标决定。常见的控制模式有定有功功率、定无功功率、定直流电压、定交流电压四种。在数学上每种控制模式对应一条额外方程比如定直流电压模式对应(U_{dc,conv}-U_{dc}^{ref}0)定无功模式对应(Q_{conv}-Q^{ref}0)。2.3 联立方程组的牛顿-拉夫逊求解框架统一求解法的关键是把所有方程放在同一个失配量向量F里。我习惯把状态向量排列为[ x [\theta_{ac};\ V_{ac};\ U_{dc};\ I_{conv}] ]其中(\theta_{ac})是交流节点相角(V_{ac})是交流节点电压幅值(U_{dc})是直流节点电压(I_{conv})是换流器内部电流也可以根据需要换成调制比或触发角。失配量向量对应[ F [\Delta P_{ac};\ \Delta Q_{ac};\ \Delta P_{dc};\ \Delta F_{conv}] ]交流节点和直流节点的失配量就是上面讲的功率平衡方程的右边减左边换流器失配量包括有功平衡方程、无功控制方程、直流电压控制方程等。牛顿法的迭代格式是[ F(x^k)J(x^k)\Delta x^k0 ][ x^{k1}x^k-\Delta x^k ]其中雅可比矩阵(J\partial F/\partial x)。这个矩阵天然是分块结构[ J \begin{bmatrix} J_{ac,ac} J_{ac,dc} \ J_{dc,ac} J_{dc,dc} \end{bmatrix} ]交流系统和直流系统的交叉耦合全部体现在(J_{ac,dc})和(J_{dc,ac})这两个非对角块里。交替法之所以在强耦合时拉胯就是因为它没有这两个交叉块统一法把交叉偏导显式写出来等于把“两个系统如何互相影响”这件事一次算清。3. Matlab实现关键细节3.1 数据结构设计节点、支路、换流器我最初想写一套“面向对象”风格的程序后来发现配电网节点规模通常只有几十到几百个用Matlab的struct数组反而更直观、调试更方便。我建议至少建立三张表交流节点表、直流节点表、换流器表。交流节点表每个元素包含节点编号、节点类型、负荷有功无功、电源有功无功、电压初值、相角初值。交流支路表每个元素包含首末端节点、电阻、电抗、对地导纳和变压器变比。直流节点表包含节点编号、直流负荷/电源注入功率、直流电压初值。换流器表最重要的字段是控制模式我用两个状态枚举量表示“有功/直流电压控制模式”和“无功/交流电压控制模式”再存对应的参考值。另外必须建立起三个索引映射交流节点编号到交流雅可比行列位置的映射、直流节点编号到直流雅可比行列位置的映射、换流器变量到全局状态向量位置的映射。这个映射关系如果建不好后面组装雅可比矩阵时全是坑。我的经验是提前画一张“变量排布图”先排交流相角再排交流电压幅值再排直流电压最后排换流器内变量写代码时照着表格填就行。3.2 雅可比矩阵的组装与修正组装雅可比矩阵是最容易出错的一步但也是最值得花时间做的一件事。交流潮流部分可以直接参照传统牛顿法的写法对(P_i)关于(\theta_j)、(V_j)求偏导对(Q_i)关于(\theta_j)、(V_j)求偏导。直流部分对直流节点功率方程求偏导结果直接就是直流电导矩阵(G_{dc})的某个元素。真正需要手推的是换流器相关的交叉偏导。比如换流器有功平衡方程(P_{ac,conv}P_{dc,conv}P_{loss})失配量对交流相角(\theta_i)的偏导实际上就是交流潮流方程中换流器注入有功对(\theta_i)的偏导失配量对直流电压(U_{dc,k})的偏导会涉及(U_{dc}I_{conv})项以及损耗(P_{loss})对(I_{conv})的导数。这些交叉偏导必须一项一项写清楚不能偷懒用差分近似代替。解析雅可比在收敛速度和稳定性上的优势非常明显尤其矩阵维数几百的时候差分法既慢又容易引入数值噪声。我常用的一个技巧是把雅可比矩阵按“分块矩阵”方式维护先算AC块再算DC块再算换流器块最后把交叉偏导数填充进去。所有矩阵都用(sparse)存储求解线性方程时直接写“dx J \ F”千万不要写成“inv(J)*F”因为Matlab对稀疏矩阵的除法运算在数值稳定性和计算速度上都远好于显式求逆。下表是我整理的常见偏导关系方便对照失配量类型对交流相角的偏导对交流电压的偏导对直流电压的偏导对换流器电流的偏导(\Delta P_{ac})来自交流潮流方程来自交流潮流方程换流器有功注入项换流器方程电流项(\Delta Q_{ac})来自交流潮流方程来自交流潮流方程通常为0取决于换流器无功方程(\Delta P_{dc})换流器直流侧功率耦合项通常为0来自直流潮流方程(U_{dc})项(\Delta F_{conv})换流器交流功率方程项换流器交流功率方程项直流电压控制/功率方程项损耗及电压方程项3.3 收敛判据与初值选取统一求解法的收敛判据我建议同时看两个指标失配量范数和状态修正量范数。只盯着失配量有可能出现“F很小但dx还很大”的临界情况只盯dx又可能在迭代初期误判收敛。我实际使用的判据是[ \max(|F|) 1e^{-6} \quad \text{且} \quad \max(|\Delta x|) 1e^{-8} ]初值方面交流PQ节点电压平启动通常设为(1.0\angle 0^\circ)直流节点电压设为1.0标幺值换流器电流初值根据参考功率和直流电压估算为(I_{conv,0}P^{ref}/U_{dc,0})。这个初值策略对大多数配电网算例都够用但如果系统重载或者控制方向比较激进我会在牛顿迭代前多做一次“纯交流潮流”预迭代先把换流器当作PQ节点交流网单独收敛一次再把直流电压初值设为换流器交流侧功率对应的直流电压这样整个统一迭代的初值质量会好很多。4. 代码实现步骤与核心函数4.1 主程序流程我把整个主程序拆成六个步骤每一步都对应一个独立函数方便单独调试定义算例数据交流节点、交流支路、直流节点、换流器参数构建交流导纳矩阵(Y_{ac})和直流电导矩阵(G_{dc})初始化全局状态向量(x)进入牛顿迭代循环计算失配量F组装雅可比J求解线性方程更新状态检查收敛收敛后计算支路功率、换流器损耗、直流网络损耗输出并可视化节点电压和功率分布。这里特别强调第4步的循环上限要设置合理我一般设50次。如果50次还不收敛与其继续迭代不如停下来分析问题硬跑300次只会浪费时间。4.2 关键代码片段示例构建交流导纳矩阵的代码比较基础我直接给出一个不带变压器抽头的简化版function Y buildYac(branch) % branch: [from, to, r, x, b]单位标幺值 % 输出稀疏交流导纳矩阵 n max(max(branch(:,1:2))); Y sparse(n, n); for k 1:size(branch,1) f branch(k,1); t branch(k,2); z branch(k,3) 1j*branch(k,4); y 1 / z; b_shunt branch(k,5); Y(f,f) Y(f,f) y 1j*b_shunt/2; Y(t,t) Y(t,t) y 1j*b_shunt/2; Y(f,t) Y(f,t) - y; Y(t,f) Y(t,f) - y; end end失配量计算函数的骨架如下这里只展示最关键的结构function F calcF(x, idx, data) % 状态向量x依次包含交流相角、交流电压、直流电压、换流器电流 % 先把x拆开 theta x(idx.theta); V_ac x(idx.V_ac); V_dc x(idx.V_dc); I_conv x(idx.I_conv); F zeros(idx.F_total, 1); % 计算交流节点失配量 [Pcal, Qcal] power_ac(V_ac, theta, data.Y_ac); F(idx.F_P) data.P_inj - Pcal; F(idx.F_Q) data.Q_inj - Qcal; % 计算直流节点失配量 F(idx.F_DC) data.P_dc_inj - V_dc .* (data.G_dc * V_dc); % 计算换流器失配量 % 换流器有功平衡P_ac_conv P_dc_conv P_loss F(idx.F_conv) data.P_conv_ref - (V_dc(data.vsc_dc) .* I_conv ... data.loss_a data.loss_b .* I_conv data.loss_c .* I_conv.^2); end主循环代码我习惯写成for iter 1:max_iter F calcF(x, idx, data); J buildJac(x, idx, data); dx -J \ F; x x dx; if max(abs(F)) tol_F max(abs(dx)) tol_dx break; end end每轮迭代我还会把(max(abs(F)))打印出来这样能实时判断到底是交流方程卡住、直流方程卡住还是换流器方程卡住。这个习惯让我后面排查问题节省了大量时间。4.3 结果输出与算例验证验证程序我用的算例是改造后的IEEE 33节点系统在节点18处接了一个三端直流网络两个VSC分别连接交流节点18和25其中一个定直流电压另一个定有功功率直流网络里接两台光伏逆变器当作直流电源。整体算例规模不大但交流重载和直流多源并存足够暴露问题。我对比了三种结果纯交流潮流结果、统一求解法带直流系统后的结果、以及把换流器有功功率设为0后的统一求解法结果。后两者在换流器有功设为0时应该完全一致这是最快的自检方法。实际算例中接入直流系统后交流节点18的电压从0.9738下降到0.9652节点25的电压从0.9811下降到0.9776这说明直流系统从交流侧吸收了有功功率压降方向符合物理直觉。直流网络内部定直流电压节点电压维持在1.002定功率节点电压略低换流器损耗占传输功率的2.3%左右这个量级对VSC是合理的。如果你要复用我的验证思路我建议先跑一个最小算例3个交流节点加2个直流节点手工用计算器算一组近似解再和Matlab程序对比。这样能确保不是“阴差阳错跑出了合理结果”而是每个方程、每个偏导都确实写对了。5. 常见问题与排查技巧实录5.1 换流器控制模式组合导致雅可比矩阵奇异性这是我在项目里踩得最狠的一个坑。直流网络必须有一个“直流电压锚点”就像交流系统必须有平衡节点一样。如果两个换流器同时定直流电压直流电压变量被重复约束雅可比矩阵会出现冗余行如果一个换流器都没有定直流电压直流网络电压就会整体漂移雅可比矩阵又会出现零特征值。两种情况的结果都是(J)奇异线性方程解出一堆奇怪的大数。排查方法很简单迭代之前先统计各换流器的有功/直流电压控制模式确保有且仅有一个定直流电压控制其余换流器要么定有功要么定下垂控制。交直流混合配电网实际工程中也基本会留一个主换流器维持直流母线电压。如果确实需要多端支撑可以把其中一个定电压模式改成下垂控制方程里加入一条下垂系数约束这样既能维持电压又不会造成冗余。5.2 交直流两侧标幺值不统一带来的“隐形错误”我刚写完代码时怎么检查结果都不对交流侧节点电压正常直流电压也正常但换流器两侧的有功功率差了20%。后来发现是直流侧的电压基准选错了。交流侧我用线电压标幺值直流侧却直接用额定直流电压的有名值代入换流器公式导致功率在交直流转换时悄悄出现了一个比例系数误差。解决这个问题的原则是全系统统一采用同一个功率基准(S_{base})交流电压基准(V_{ac,base})取线电压或相电压都可以但导纳折算必须用对应的(Z_{base}V_{base}^2/S_{base})直流电压基准(V_{dc,base})取额定直流电压。所有有名值在进入方程组之前一次性转换成标幺值不要每行方程里零零散散地转换。我甚至专门写了一个“baseValue”结构体里面有S_base、V_ac_base、V_dc_base、Z_ac_base、Z_dc_base换流器方程中所有功率和电压都从标幺值再回到有名值检查一遍确保一致。5.3 迭代发散时的三阶段定位法统一求解法发散时不要急着改初值或调阻尼。我自己的排查顺序分三步第一步把换流器有功功率全部设成0直流负荷也设成0这时换流器只剩下无功方程。如果程序能收敛到纯交流潮流解说明交流部分的雅可比矩阵没有问题。第二步把直流网络单独拿出来用线性直流潮流解一次检查直流电压是否正确这一步可以定位直流电导矩阵是否填错。第三步再把换流器有功功率逐渐加上去每加一步都判断迭代是否收敛。哪一步开始发散问题就在哪一步涉及的方程和雅可比块里。除了这个方法我还会在每轮迭代里打印下面这类信息最大失配量属于哪个方程、对应哪个变量、最近两轮的修正方向是否来回跳。如果发现某个失配量始终降不下去十有八九是雅可比矩阵里对应的偏导公式写错了一个正负号。6. 一点扩展思考与个人经验6.1 统一求解法向三相不平衡、动态仿真的扩展配电网三相不平衡是常态统一求解法完全可以在三相模型上继续用。状态变量从单相电压相角变成三相电压相角交流支路要用三相阻抗矩阵换流器模型也要换成abc坐标或正负零序分量形式。雅可比矩阵规模会明显变大、非零元更多但统一求解法“联立全部方程”的框架不用变只是J矩阵里的交叉项更丰富了。关于动态仿真统一求解法的思想也很有用。把储能变换器、同步机转子运动方程这些微分方程用隐式梯形积分离散后和潮流代数方程一起组成差分代数方程组再统一做牛顿迭代。本质上就是把“代数潮流”和“微分动态”当成一个更大的方程系统来联立求解。我在实际工程中体会最深的一点是只要把F(x)和J(x)写成不依赖全局变量的纯函数后面从静态潮流扩展到动态仿真是水到渠成的事根本不用重写核心迭代器。6.2 与智能算法结合的松紧耦合策略交直流混合配电网里换流器控制参数、储能出力曲线通常有大量可调自由度很多人喜欢用粒子群、遗传算法这类智能算法去寻优。我建议把统一求解法包成一个“潮流计算内核”负责快速计算每一组控制参数下的运行状态把节点电压和换流器损耗作为适应度值返回。这个内层求解器必须足够快、足够稳定否则外层智能算法跑几百次全卡在潮流发散上。我个人的体会是统一求解法配合Matlab的稀疏矩阵对一个几十节点的交直流配电网单次潮流计算通常在毫秒级以内完全能作为智能算法的内层适应度函数。后续如果你要我给一条最实在的建议那就是所有雅可比矩阵偏导都写成解析式对所有状态变量都做一次中心差分校验一旦偏导正确整个程序的稳定性会骤然提升。如果你打算自己复现我建议先把3节点交流加2节点直流的最小模型手算跑通再扩展到33节点系统。那时候你才能真正体会到统一求解法的坑和乐趣全都在雅可比矩阵的交叉项里。