交直流混合系统的潮流计算麻烦的从来不是交流那一半。交流部分的 PQ 分解、牛顿迭代、稀疏求解随便找个开源库都能跑。真正让人抓头的是直流侧那几个换流器它们既是交流系统的负荷又是电源还带着触发角、熄弧角、换相电抗这一堆相互纠缠的参数。你改一个直流电流交流侧两个节点的有功无功同时变直流电压也跟着变最后连控制方式都可能从一个切到另一个。我最初做这套东西是为了给一个含两端直流联络线的测试系统算稳态工作点。一开始图省事用了交替迭代交流算一遍、直流算一遍、来回倒腾。小扰动工况下能收敛一旦把直流功率推到接近额定值或者在逆变侧把熄弧角压到 16 度附近迭代就开始来回振荡残差卡在 1e-3 上下死活下不去。后来老老实实改成统一迭代法把交流节点电压、相角、直流电压、直流电流、控制角全部塞进同一个未知量向量用同一个牛顿步一起修正收敛性立刻不一样了——典型工况 4 到 6 次迭代就到 1e-8。这篇就把这套基于 Matlab 的统一迭代法交直流潮流计算从头到尾拆一遍模型怎么建、方程怎么配平、雅可比矩阵那四个块怎么拼、代码骨架长什么样、以及我在调试时踩过的那些坑。适合已经会写纯交流牛顿法潮流、想往交直流方向走一步的读者也适合手头有个混合系统要算稳态点、不想从零推公式的人直接抄结构。1. 为什么交直流潮流要单独搭一套求解器1.1 纯交流潮流模型在换流站面前会直接失灵纯交流潮流里节点只有三类PQ、PV、平衡节点。每个节点的注入功率是常数或者电压幅值给定雅可比矩阵的元素都是电压幅值和相角的函数形式干净利落。换流站一进来这个假设就崩了。换流站对交流系统来说有功和无功都不是定值。有功由直流侧传输功率决定无功由换流器的功率因数决定而功率因数又跟触发角、换相重叠角、直流电流挂钩。这意味着你在交流节点上写不出一个只含 V 和 θ 的功率方程——它里面必然夹着直流侧的未知量。更直接点说换流节点是一个特性由另一套方程描述的负荷你没法把它简化成恒功率或者恒阻抗。还有一层换流变压器变比、换相电抗、每极换流器个数、直流接线方式单极、双极、金属回线这些结构参数都会改变方程的形态。双极系统里正负极的直流电流在理想对称时相等但接地回路和金属回线的处理方式不同方程数量也不同。我见过有人直接把单极代码改个系数套到双极上结果雅可比维度都配不平跑起来直接报维度不匹配。所以结论很明确交直流潮流不是在交流潮流上加个直流负荷而是要重新建立一套包含两类变量、两类方程的联合模型。这一点想清楚了后面的工作才有方向。1.2 交替迭代和统一迭代的取舍逻辑工程上算交直流潮流有两条路。第一条是交替迭代也叫顺序求解先把直流侧等效成交流节点的注入功率解一遍交流潮流再把交流母线电压代回换流器方程解一遍直流侧两个结果来回传递直到前后两次的换流站功率差小于阈值。第二条是统一迭代把所有未知量放进一个向量所有方程放进一个残差向量用牛顿法一次性求解。交替迭代的优点是实现简单、代码复用率高交流部分可以直接拿现成的牛顿法程序。但它的收敛性依赖一个隐含条件交流侧和直流侧之间的耦合不能太强。直流功率越大、换流站电气距离越近、控制方式越接近切换边界这个耦合就越强交替迭代的收敛速度会急剧下降甚至出现振荡。我实测过一个直流功率 800 MW 的两端系统交替迭代要 20 多次才勉强收敛把功率提到 950 MW 就直接发散了。统一迭代的代价是要自己推导并组装扩展雅可比矩阵工作量大概多出两三天。换来的是牛顿法本身的二次收敛特性——只要初值不离谱残差平方级下降。而且因为耦合项显式出现在雅可比里迭代过程对控制方式切换、直流功率大范围变化的适应性明显更好。我现在这套代码直流功率从 0 扫到 1.05 倍额定值除了极端工况需要调整初值基本都能收敛。选型上我的建议是如果只是做个粗略估算、直流容量占比很小交替迭代够用但凡要做参数扫描、要做 N-1 分析、要研究控制方式切换直接上统一迭代省得后期反复返工。1.3 整套程序的结构划分代码我按四层来组织这个分层方式用下来最省心。第一层是数据层把交流网络参数支路、节点类型、负荷、直流网络参数换流站位置、换流变变比、换相电抗、控制目标值、直流电阻分开存成结构体。第二层是模型层负责根据当前状态变量算出所有残差和雅可比块的元素。第三层是求解层做稀疏线性方程求解和状态量更新。第四层是输出层把结果整理成表格、算换流站的功率因数、校验功率平衡。这么分层的好处是调试时能单独拎出某一层来验证。比如怀疑雅可比算错了就把模型层的数值雅可比和解析雅可比做一次差分对比误差应该在 1e-6 量级。这种自检手段在后期定位问题时非常关键比盯着不收敛的残差曲线瞎猜有效得多。整套代码的核心文件我控制在五个左右主脚本、数据文件、残差函数、雅可比函数、结果输出函数。文件太多反而增加跳转成本。2. 把物理系统拆成可求解的方程2.1 交流侧节点功率平衡方程怎么写交流侧跟纯交流潮流完全一样每个节点两组方程有功不平衡和无功不平衡。对 PQ 节点两个都写对 PV 节点只写有功、无功由电压幅值给定的条件替代平衡节点两个都不写。这部分我不赘述细节只强调一点换流节点必须按 PQ 节点处理因为它的有功无功都随直流状态变化不能指定电压幅值。节点注入功率写成标准形式P_i V_i * Σ_j V_j * (G_ij*cosθ_ij B_ij*sinθ_ij) Q_i V_i * Σ_j V_j * (G_ij*sinθ_ij - B_ij*cosθ_ij)换流节点上这个 P_i 和 Q_i 里要额外加上直流注入项。注意符号约定我统一取交流系统向换流站注入为正那么整流站功率从交流流向直流的注入是负的逆变站是正的。符号约定如果一开始没定死后面推导偏导时会非常痛苦很容易出现一整块雅可比符号全反的情况。2.2 换流器稳态方程与等值换相电阻这是整个模型的核心也是新手最容易栽跟头的地方。理想换流器的稳态方程忽略换相重叠时是U_d (3√2/π) * n_T * U_s * cosα考虑换相重叠后直流电压被吃掉一块用等值换相电阻表示U_d (3√2/π) * n_T * U_s * cosα - (3/π) * X_c * I_d这里 n_T 是换流变压器变比折算到交流侧电压基准X_c 是折算到同一基准下的换相电抗I_d 是直流电流。把 (3/π) * X_c 记作 R_c就是所谓的等值换相电阻它不代表真实的有功损耗只是描述换相过程对直流电压的压降效应。这个区分很重要——算功率平衡的时候R_c 上消耗的功率不是热损耗而是换流器无功消耗的一部分体现。逆变侧方程形式上要改符号因为直流电流的方向相对逆变器是流出的U_d (3√2/π) * n_T * U_s * cosγ R_c * I_dγ 是熄弧角跟触发角的关系是 γ π - α - μμ 是换相重叠角。实际工程里逆变侧一般直接拿 γ 当控制变量比用 α 更直观因为 γ 直接对应换相失败的风险。写代码时我把常数 K (3√2/π) * n_T 单独算出来存着每次迭代只更新 cosα 和 I_d能省一点重复计算。别小看这点优化几百次迭代、上千个节点的时候累加起来的开销是可观的。2.3 直流网络方程与控制方程的配对直流网络本身很简单。两端系统就是一条支路直流电阻 R_d两端电压 U_dr 和 U_di直流电流 I_dU_dr - U_di - R_d * I_d 0多端系统就写成节点电流平衡的形式跟交流网络的节点导纳矩阵完全一个思路只是直流侧只有电阻没有电抗导纳矩阵是实数对称阵。我建议直流网络也建成矩阵形式这样扩展到三端、四端的时候不用改代码结构只改数据。真正需要动脑子的是控制方程。直流系统必须有两个独立的控制量才能确定工作点缺一个方程就欠定多一个就超定。常见的组合有这么几种控制方式组合整流侧方程逆变侧方程适用场景定电流 定熄弧角I_d - I_set 0γ - γ_set 0传统直流输电最常见定功率 定直流电压U_dr*I_d - P_set 0U_di - U_set 0背靠背、柔性直流常见定电流 定直流电压I_d - I_set 0U_di - U_set 0一端定电流、一端稳压定触发角 定熄弧角α - α_set 0γ - γ_set 0特殊研究工况这里有个关键细节我要专门说未知量的选择直接决定方程形式。如果你把 I_d 当成一个全局变量单极单回线只有一条直流电流那么定电流控制就是把 I_d 钉在给定值上方程直接是 I_d - I_set 0非常干净。但如果两端各有独立的直流电流变量定电流方程就要分别写变量数也多出来容易配不平。我的做法是每条直流支路只有一个电流变量多端系统里用节点电流平衡代替支路电流这样变量数最少雅可比也最紧凑。2.4 未知量向量和方程向量的完整配平以双端单极系统、n 个交流节点为例其中一个平衡节点、npv 个 PV 节点未知量包括所有非平衡节点的电压相角n-1 个所有 PQ 节点的电压幅值n - 1 - npv 个直流侧变量U_dr、U_di、I_d、α或整流侧控制角、γ 共 5 个方程包括非平衡节点的有功不平衡方程n-1 个PQ 节点的无功不平衡方程n - 1 - npv 个直流回路方程 1 个整流侧换流器方程 1 个逆变侧换流器方程 1 个整流侧控制方程 1 个逆变侧控制方程 1 个数量刚好对上。我第一版代码就是在这上面栽了把逆变侧的换流器方程漏写了维度差 1Matlab 报的是矩阵接近奇异排查了大半天才发现是方程少了一个。建议在写代码前先在纸上把变量清单和方程清单各列一遍一一对应打勾。这一步花十分钟能省后面十个小时。3. 统一迭代法雅可比矩阵的拼接方法3.1 扩展雅可比的分块结构统一迭代的整体迭代格式还是牛顿法那一套J * Δx -F(x) x_new x Δx其中 F 是所有方程的残差向量J 是 F 对 x 的偏导矩阵。按交流变量和直流变量分块J 长这样J [ J_AA J_AD ] [ J_DA J_DD ]J_AA 就是传统交流潮流的雅可比这部分代码可以直接复用。J_AD 是交流方程对直流变量的偏导J_DA 是直流方程对交流变量的偏导J_DD 是直流方程对直流变量的偏导。这两个耦合块是统一迭代的灵魂也是交替迭代法里根本没有的东西。很多教程讲到这里就停了只说把耦合项加上但不告诉你具体加什么。下面我把每一块都摊开。3.2 交流侧对直流变量的耦合项如何推导换流节点上直流注入交流的有功是P_dc U_d * I_d注意这是有功的绝对值符号靠前面约定的方向来定。对整流侧注入是 -U_dI_d对逆变侧是 U_dI_d。所以∂P_dc/∂U_d ±I_d ∂P_dc/∂I_d ±U_d ∂P_dc/∂α 0有功部分很干净控制角不出现在有功表达式里。无功部分麻烦一些因为要过功率因数这一关。换流器的功率因数有个漂亮的关系式可以从换流器方程反推出来cosφ U_d / (K * U_s)推导过程很简单cosα (U_d R_cI_d)/(KU_s)cos(αμ) (U_d - R_cI_d)/(KU_s)两者相加除以二就是 cosφ。这个式子说明什么说明换流器的功率因数完全由直流电压和交流母线电压的比值决定跟触发角没有直接关系。那么无功功率Q_dc P_dc * tanφ U_d * I_d * sqrt(1 - (U_d/(K*U_s))²) / (U_d/(K*U_s)) I_d * sqrt((K*U_s)² - U_d²)这个形式很好用偏导直接能算∂Q_dc/∂U_d -I_d * U_d / sqrt((K*U_s)² - U_d²) ∂Q_dc/∂I_d sqrt((K*U_s)² - U_d²) ∂Q_dc/∂U_s I_d * K² * U_s / sqrt((K*U_s)² - U_d²)最后一项 ∂Q_dc/∂U_s 是容易被漏掉的关键项。它落的位置是 J_AA 块里换流节点无功方程对自身电压幅值的偏导。你如果只按传统公式写 J_AA这一项就丢了结果是雅可比不精确牛顿法退化成拟牛顿法收敛次数从 5 次涨到 15 次甚至不收敛。我第一版就是这个问题收敛慢得离谱后来做数值差分对比才发现漏了这块。3.3 直流侧对交流变量的耦合项怎么处理直流方程里显式含有交流母线电压的只有换流器方程。以整流侧为例残差写成f_conv U_d - K * U_s * cosα R_c * I_d 0各偏导∂f/∂U_d 1 ∂f/∂U_s -K * cosα ∂f/∂α K * U_s * sinα ∂f/∂I_d R_c逆变侧符号略有不同cosα 换成 cosγR_c 项符号反过来。直流回路方程对交流变量没有偏导因为它只涉及直流量。控制方程也是定电流、定熄弧角、定功率都不含交流电压。所以 J_DA 块其实只有换流器方程那两行是稀疏的非零其余全是零。这个稀疏性在组装的时候要利用起来不然矩阵会白白变稠。3.4 标幺化处理与稀疏矩阵组装标幺化这件事我必须单独强调。交直流混合系统里交流侧电压是几百千伏量级直流电流是千安量级功率是兆瓦量级如果不统一基准雅可比矩阵的条件数会非常糟糕线性求解器给出的解精度很差。我的处理方式交流侧沿用系统基准 S_base 和交流电压基准 U_base阻抗基准 Z_base U_base²/S_base。直流侧取直流电压基准 U_dB 等于换流变交流侧电压基准折算到直流侧的值直流电流基准 I_dB S_base/U_dB直流电阻基准 R_dB U_dB/I_dB这正好等于 Z_base。这样所有量都在 1 附近雅可比元素数量级一致求解稳定。代码里稀疏组装用三元组方式% 预分配按最大非零元数估算 max_nnz 4 * n 20; rows zeros(max_nnz, 1); cols zeros(max_nnz, 1); vals zeros(max_nnz, 1); idx 0; % 交流块填充示例实际按行遍历更省事 for i 1:n-1 % ... 逐行填充 J_AA end % 最后一次性组装 J sparse(rows(1:idx), cols(1:idx), vals(1:idx), nvar, nvar);用 sparse 而不是 zeros 建满矩阵10 个节点看不出区别100 个节点就能差出一个数量级的求解时间。而且 Matlab 的 \ 运算符对稀疏矩阵会自动选合适的求解器比手动调 lu 靠谱。4. Matlab 实现全流程与代码骨架4.1 数据结构怎么组织才不容易乱我用四个结构体bus交流节点、branch交流支路、dc直流系统、opt求解选项。bus struct(num, [], ... % 节点编号 type, [], ... % 1PQ, 2PV, 3平衡, 4换流站 V, [], ... % 电压幅值初值标幺 theta, [], ... % 相角初值弧度 Pg, [], Qg, [], Pl, [], Ql, []); dc struct(conv, [], ... % 换流站数组每站一个结构 line, []); % 直流线路参数 dc.conv(k).bus_ac 5; % 接入的交流节点号 dc.conv(k).nT 1.0; % 换流变变比 dc.conv(k).Xc 0.1; % 换相电抗标幺 dc.conv(k).mode rect; % 整流还是逆变 dc.conv(k).ctrl I; % 控制变量类型 dc.conv(k).ctrl_val 1.0; % 控制目标值用结构体数组而不是元胞数组好处是可以直接用点号访问代码读起来清楚。换流站数量变化时不用改任何逻辑程序里全部按 length(dc.conv) 循环。有个细节换流站的交流接入节点类型要标成 4这样在主循环里能快速定位哪些节点需要加直流注入项、哪些节点要跳过电压幅值更新。用 type 字段区分比单独维护一个索引数组更不容易出错。4.2 主迭代循环的完整骨架整个牛顿迭代循环结构如下max_iter 20; tol 1e-8; x init_state(bus, dc); % 初始状态向量 history zeros(max_iter, 1); for iter 1:max_iter [F, J] build_residual_jacobian(x, bus, branch, dc, Ybus, ...); err max(abs(F)); history(iter) err; if err tol fprintf(第 %d 次迭代收敛最大残差 %.2e\n, iter, err); break; end if iter max_iter warning(达到最大迭代次数仍未收敛当前残差 %.2e, err); end dx -J \ F; % 稀疏左除自动选求解器 x x dx; % 越限检查和钳位 x enforce_limits(x, dc); end几个实操要点第一dx -J \ F这一行是稀疏求解的核心Matlab 会自动判断用 LU 还是 Cholesky不要手写 inv(J)*F数值上差很多。第二enforce_limits这个函数判断控制角的合理性。触发角一般限制在 5 到 90 度熄弧角限制在 10 到 60 度。牛顿迭代偶尔会给出负的触发角或超过 90 度的值这时候如果直接钳位残差会跳变迭代可能被带偏我的做法是钳位的同时把对应的控制方程临时改成角度不变相当于人为把这个变量冻结一步让其他变量先收敛。第三迭代日志一定要存。我把每次的残差存进 history 数组收敛完成后画个半对数图能直接看出收敛类型——二次收敛是斜率越来越陡线性收敛是直线振荡是锯齿。这个图对调试的价值极高。4.3 换流器残差与雅可比的子函数实现把换流器的计算单独封装成函数输入是直流状态量、交流电压输出是残差和对应的雅可比元素function [f, Jrow] conv_equation(Ud, Id, Us, alpha, K, Rc, mode) if strcmp(mode, rect) f Ud - K * Us * cos(alpha) Rc * Id; Jrow struct(dUd, 1, ... dUs, -K * cos(alpha), ... dalpha, K * Us * sin(alpha), ... dId, Rc); else f Ud - K * Us * cos(alpha) - Rc * Id; Jrow struct(dUd, 1, ... dUs, -K * cos(alpha), ... dalpha, K * Us * sin(alpha), ... dId, -Rc); end end无功对电压幅值的偏导那个 sqrt((KU_s)² - U_d²) 分母有开方当 U_d 接近 KU_s 时cosφ 接近 1换流器运行在接近单位功率因数的极端状态这个偏导会趋于零反过来 ∂Q/∂U_d 会趋于无穷。实际系统里 U_d/(K*U_s) 一般在 0.85 到 0.95 之间离奇异点还有距离但代码里要加个保护ratio Ud / (K * Us); ratio min(max(ratio, 0.5), 0.999); % 防止开方为负或除零这个保护不是多余的。我遇到过一次逆变侧电压被迭代推到很低、ratio 超过 1sqrt 直接返回复数整个迭代崩掉。加上钳位之后至少不会出复数虽然结果不物理但能让你从残差曲线看出是哪一步跑偏的。4.4 初值选取与收敛控制牛顿法对初值敏感交直流系统尤其如此。我的初值策略是分层给交流侧全部给 1.0 幅值、0 相角这个标准做法基本够用。直流侧要稍微算一下直流电压初值按额定值给直流电流初值由目标功率反算 I_d P_set / U_d整流侧触发角给 15 度约 0.26 弧度逆变侧熄弧角给 18 度。为什么给 15 到 18 度这个范围因为这是工程上实际运行点的典型值牛顿法从这附近出发基本一步就能进收敛域。你要是给 0 度或者 80 度第一第二个迭代步会跑得很远有时候直接跑到非物理解那边去了。收敛判据我用最大残差而不是二范数。二范数在大系统里会被节点数量放大阈值不好定最大残差直观1e-8 就是每个方程都满足到 8 位有效数字物理意义明确。还有一个技巧前两次迭代可以用放宽的雅可比也就是用初值处的雅可比多迭代几步省下重新组装和分解的时间。但对交直流系统来说耦合项变化快这个技巧收益不大我最后没采用老老实实每步重算。5. 实际调试中踩过的坑5.1 不收敛的几种典型症状和定位路径残差曲线是最好的诊断工具。我总结了几种形态和对应的病因残差表现可能原因排查手段首次迭代后残差反而变大雅可比符号错误或初值严重偏离数值差分校验雅可比残差线性下降每次只降一个数量级耦合项漏写退化成拟牛顿检查 J_AD 和 J_DA 是否为零残差在 1e-3 附近来回振荡控制方式在切换边界反复跳打印每步触发角看是否触发限值残差突然变成 NaN 或复数开方项为负、除零、角度越界在各处加 assert 断言前三次下降很快之后停滞雅可比接近奇异算 cond(J)检查是否有零对角元数值差分校验雅可比是排查的第一步也是最有效的一步。做法很简单对每个状态变量给一个 1e-6 的扰动用前向差分算数值偏导跟解析偏导逐个比较相对误差超过 1e-4 的就是有问题。我第二次遇到收敛慢的问题时就是这么定位的查出来是 ∂Q/∂U_s 那一项少了个 K² 因子。5.2 收敛到非物理解的表现和规避牛顿法有个讨厌的特性它只管把残差推到零不管解的物理意义。交直流系统里很容易出现这种情况——迭代收敛了残差 1e-10但打印出来发现整流侧触发角是负的 20 度或者逆变侧熄弧角是 3 度这明显是换相失败工况不可能是稳态解。原因通常有两个。一是初值离物理解太远牛顿步跨过了物理可行域落到了数学解的另一支上。二是控制方式配错了比如整流侧明明应该是定电流你按定功率写了方程系统会找到一个数学上满足定功率但物理上不可行的点。规避办法每次迭代后检查所有控制角的物理范围超出范围就打印警告并记录收敛后如果仍有越界就说明这次求解无效需要换初值重来。我在主程序里加了自动重启机制如果首次求解出现越界就把触发角初值从 15 度调到 20 度重新算一次通常第二次就能落在合理区域。5.3 结果校验的三种手段算完之后必须校验不能收敛了就用。我固定做三件事第一全系统有功平衡校验。所有交流电源的有功加上直流注入的有功减去负荷有功和网络损耗应该接近零。误差在 1e-6 量级说明有功方程没问题。这里要注意等值换相电阻 R_c 上的功率不是损耗不能算进网损这是新手常犯的错误。第二换流站交直流功率对照。交流侧算出来的注入有功跟直流侧 U_d * I_d 应该相等。这两条路径算出来的值如果差得多说明换流器方程或者功率因数公式用错了。我一般要求相对误差小于 1e-6。第三跟交替迭代的结果做交叉验证。同一组数据用两套方法算结果应该一致到小数点后四五位。这个交叉验证做过一次之后后面改代码就有信心了。5.4 我的几条实操心得换相电抗的单位是最容易出错的地方。数据手册上给的往往是折算到换流变阀侧的欧姆值你必须先折算到交流母线侧再除以阻抗基准变成标幺值。少除一个变比平方结果能差出百分之几而且不容易发现因为照样收敛。直流电阻的处理也一样双极系统的直流线路电阻是单极的一半两极并联金属回线方式下回路电阻还要加上回线电阻。这些细节数据文件里写错了结果就是错的程序本身不会报错。还有一点换流变压器的变比定义有歧义。有的资料定义成阀侧额定电压/交流侧额定电压有的定义成交流侧电压/阀侧电压两者互为倒数。用错了直流电压直接差一个变比平方功率就差得更多。我的建议是在数据文件里明确注释清楚定义并且用一个已知工况反算验证一遍。6. 这套代码后续能怎么扩展现在这套代码基本框架是稳的往几个方向扩展都不难。多端直流就是把直流网络从单支路推广到节点导纳矩阵形式方程结构不变只是变量和方程数量增加雅可比矩阵的稀疏性要重新组织一下。加入直流网络损耗的精细模型也不难把线路电阻从常数改成随温度变化的函数或者加上线路的分布参数效应。如果要做含新能源接入的混合系统难点不在潮流本身而在于新能源的出力不确定性——那就变成概率潮流或者场景分析了需要外层套蒙特卡洛或者点估计法。这时候统一迭代法的优势更明显因为单次求解快几千次场景扫描才跑得动。另外一个我觉得值得做的方向是把换流器的详细模型引进来考虑换相重叠角的动态变化和阀的导通压降。这样算出来的稳态点更接近实际尤其是低功率运行的时候理想模型的误差会比较大。代价是方程更复杂雅可比的解析推导工作量翻倍可能需要借助符号计算工具辅助推导。最后分享一个我一直在用的小习惯每次改动代码之后先用一个 3 节点的小系统跑一遍冒烟测试确认能收敛、功率平衡满足再去跑 30 节点以上的实际算例。小系统的迭代过程可以直接在命令行里 print 出来一眼就能看出哪一步不对大系统一旦出问题定位成本高十倍。这个习惯帮我省下的时间比我优化代码省下的时间多得多。 基于Matlab统一迭代法的交直流混合系统潮流计算
交直流混合系统的潮流计算麻烦的从来不是交流那一半。交流部分的 PQ 分解、牛顿迭代、稀疏求解随便找个开源库都能跑。真正让人抓头的是直流侧那几个换流器它们既是交流系统的负荷又是电源还带着触发角、熄弧角、换相电抗这一堆相互纠缠的参数。你改一个直流电流交流侧两个节点的有功无功同时变直流电压也跟着变最后连控制方式都可能从一个切到另一个。我最初做这套东西是为了给一个含两端直流联络线的测试系统算稳态工作点。一开始图省事用了交替迭代交流算一遍、直流算一遍、来回倒腾。小扰动工况下能收敛一旦把直流功率推到接近额定值或者在逆变侧把熄弧角压到 16 度附近迭代就开始来回振荡残差卡在 1e-3 上下死活下不去。后来老老实实改成统一迭代法把交流节点电压、相角、直流电压、直流电流、控制角全部塞进同一个未知量向量用同一个牛顿步一起修正收敛性立刻不一样了——典型工况 4 到 6 次迭代就到 1e-8。这篇就把这套基于 Matlab 的统一迭代法交直流潮流计算从头到尾拆一遍模型怎么建、方程怎么配平、雅可比矩阵那四个块怎么拼、代码骨架长什么样、以及我在调试时踩过的那些坑。适合已经会写纯交流牛顿法潮流、想往交直流方向走一步的读者也适合手头有个混合系统要算稳态点、不想从零推公式的人直接抄结构。1. 为什么交直流潮流要单独搭一套求解器1.1 纯交流潮流模型在换流站面前会直接失灵纯交流潮流里节点只有三类PQ、PV、平衡节点。每个节点的注入功率是常数或者电压幅值给定雅可比矩阵的元素都是电压幅值和相角的函数形式干净利落。换流站一进来这个假设就崩了。换流站对交流系统来说有功和无功都不是定值。有功由直流侧传输功率决定无功由换流器的功率因数决定而功率因数又跟触发角、换相重叠角、直流电流挂钩。这意味着你在交流节点上写不出一个只含 V 和 θ 的功率方程——它里面必然夹着直流侧的未知量。更直接点说换流节点是一个特性由另一套方程描述的负荷你没法把它简化成恒功率或者恒阻抗。还有一层换流变压器变比、换相电抗、每极换流器个数、直流接线方式单极、双极、金属回线这些结构参数都会改变方程的形态。双极系统里正负极的直流电流在理想对称时相等但接地回路和金属回线的处理方式不同方程数量也不同。我见过有人直接把单极代码改个系数套到双极上结果雅可比维度都配不平跑起来直接报维度不匹配。所以结论很明确交直流潮流不是在交流潮流上加个直流负荷而是要重新建立一套包含两类变量、两类方程的联合模型。这一点想清楚了后面的工作才有方向。1.2 交替迭代和统一迭代的取舍逻辑工程上算交直流潮流有两条路。第一条是交替迭代也叫顺序求解先把直流侧等效成交流节点的注入功率解一遍交流潮流再把交流母线电压代回换流器方程解一遍直流侧两个结果来回传递直到前后两次的换流站功率差小于阈值。第二条是统一迭代把所有未知量放进一个向量所有方程放进一个残差向量用牛顿法一次性求解。交替迭代的优点是实现简单、代码复用率高交流部分可以直接拿现成的牛顿法程序。但它的收敛性依赖一个隐含条件交流侧和直流侧之间的耦合不能太强。直流功率越大、换流站电气距离越近、控制方式越接近切换边界这个耦合就越强交替迭代的收敛速度会急剧下降甚至出现振荡。我实测过一个直流功率 800 MW 的两端系统交替迭代要 20 多次才勉强收敛把功率提到 950 MW 就直接发散了。统一迭代的代价是要自己推导并组装扩展雅可比矩阵工作量大概多出两三天。换来的是牛顿法本身的二次收敛特性——只要初值不离谱残差平方级下降。而且因为耦合项显式出现在雅可比里迭代过程对控制方式切换、直流功率大范围变化的适应性明显更好。我现在这套代码直流功率从 0 扫到 1.05 倍额定值除了极端工况需要调整初值基本都能收敛。选型上我的建议是如果只是做个粗略估算、直流容量占比很小交替迭代够用但凡要做参数扫描、要做 N-1 分析、要研究控制方式切换直接上统一迭代省得后期反复返工。1.3 整套程序的结构划分代码我按四层来组织这个分层方式用下来最省心。第一层是数据层把交流网络参数支路、节点类型、负荷、直流网络参数换流站位置、换流变变比、换相电抗、控制目标值、直流电阻分开存成结构体。第二层是模型层负责根据当前状态变量算出所有残差和雅可比块的元素。第三层是求解层做稀疏线性方程求解和状态量更新。第四层是输出层把结果整理成表格、算换流站的功率因数、校验功率平衡。这么分层的好处是调试时能单独拎出某一层来验证。比如怀疑雅可比算错了就把模型层的数值雅可比和解析雅可比做一次差分对比误差应该在 1e-6 量级。这种自检手段在后期定位问题时非常关键比盯着不收敛的残差曲线瞎猜有效得多。整套代码的核心文件我控制在五个左右主脚本、数据文件、残差函数、雅可比函数、结果输出函数。文件太多反而增加跳转成本。2. 把物理系统拆成可求解的方程2.1 交流侧节点功率平衡方程怎么写交流侧跟纯交流潮流完全一样每个节点两组方程有功不平衡和无功不平衡。对 PQ 节点两个都写对 PV 节点只写有功、无功由电压幅值给定的条件替代平衡节点两个都不写。这部分我不赘述细节只强调一点换流节点必须按 PQ 节点处理因为它的有功无功都随直流状态变化不能指定电压幅值。节点注入功率写成标准形式P_i V_i * Σ_j V_j * (G_ij*cosθ_ij B_ij*sinθ_ij) Q_i V_i * Σ_j V_j * (G_ij*sinθ_ij - B_ij*cosθ_ij)换流节点上这个 P_i 和 Q_i 里要额外加上直流注入项。注意符号约定我统一取交流系统向换流站注入为正那么整流站功率从交流流向直流的注入是负的逆变站是正的。符号约定如果一开始没定死后面推导偏导时会非常痛苦很容易出现一整块雅可比符号全反的情况。2.2 换流器稳态方程与等值换相电阻这是整个模型的核心也是新手最容易栽跟头的地方。理想换流器的稳态方程忽略换相重叠时是U_d (3√2/π) * n_T * U_s * cosα考虑换相重叠后直流电压被吃掉一块用等值换相电阻表示U_d (3√2/π) * n_T * U_s * cosα - (3/π) * X_c * I_d这里 n_T 是换流变压器变比折算到交流侧电压基准X_c 是折算到同一基准下的换相电抗I_d 是直流电流。把 (3/π) * X_c 记作 R_c就是所谓的等值换相电阻它不代表真实的有功损耗只是描述换相过程对直流电压的压降效应。这个区分很重要——算功率平衡的时候R_c 上消耗的功率不是热损耗而是换流器无功消耗的一部分体现。逆变侧方程形式上要改符号因为直流电流的方向相对逆变器是流出的U_d (3√2/π) * n_T * U_s * cosγ R_c * I_dγ 是熄弧角跟触发角的关系是 γ π - α - μμ 是换相重叠角。实际工程里逆变侧一般直接拿 γ 当控制变量比用 α 更直观因为 γ 直接对应换相失败的风险。写代码时我把常数 K (3√2/π) * n_T 单独算出来存着每次迭代只更新 cosα 和 I_d能省一点重复计算。别小看这点优化几百次迭代、上千个节点的时候累加起来的开销是可观的。2.3 直流网络方程与控制方程的配对直流网络本身很简单。两端系统就是一条支路直流电阻 R_d两端电压 U_dr 和 U_di直流电流 I_dU_dr - U_di - R_d * I_d 0多端系统就写成节点电流平衡的形式跟交流网络的节点导纳矩阵完全一个思路只是直流侧只有电阻没有电抗导纳矩阵是实数对称阵。我建议直流网络也建成矩阵形式这样扩展到三端、四端的时候不用改代码结构只改数据。真正需要动脑子的是控制方程。直流系统必须有两个独立的控制量才能确定工作点缺一个方程就欠定多一个就超定。常见的组合有这么几种控制方式组合整流侧方程逆变侧方程适用场景定电流 定熄弧角I_d - I_set 0γ - γ_set 0传统直流输电最常见定功率 定直流电压U_dr*I_d - P_set 0U_di - U_set 0背靠背、柔性直流常见定电流 定直流电压I_d - I_set 0U_di - U_set 0一端定电流、一端稳压定触发角 定熄弧角α - α_set 0γ - γ_set 0特殊研究工况这里有个关键细节我要专门说未知量的选择直接决定方程形式。如果你把 I_d 当成一个全局变量单极单回线只有一条直流电流那么定电流控制就是把 I_d 钉在给定值上方程直接是 I_d - I_set 0非常干净。但如果两端各有独立的直流电流变量定电流方程就要分别写变量数也多出来容易配不平。我的做法是每条直流支路只有一个电流变量多端系统里用节点电流平衡代替支路电流这样变量数最少雅可比也最紧凑。2.4 未知量向量和方程向量的完整配平以双端单极系统、n 个交流节点为例其中一个平衡节点、npv 个 PV 节点未知量包括所有非平衡节点的电压相角n-1 个所有 PQ 节点的电压幅值n - 1 - npv 个直流侧变量U_dr、U_di、I_d、α或整流侧控制角、γ 共 5 个方程包括非平衡节点的有功不平衡方程n-1 个PQ 节点的无功不平衡方程n - 1 - npv 个直流回路方程 1 个整流侧换流器方程 1 个逆变侧换流器方程 1 个整流侧控制方程 1 个逆变侧控制方程 1 个数量刚好对上。我第一版代码就是在这上面栽了把逆变侧的换流器方程漏写了维度差 1Matlab 报的是矩阵接近奇异排查了大半天才发现是方程少了一个。建议在写代码前先在纸上把变量清单和方程清单各列一遍一一对应打勾。这一步花十分钟能省后面十个小时。3. 统一迭代法雅可比矩阵的拼接方法3.1 扩展雅可比的分块结构统一迭代的整体迭代格式还是牛顿法那一套J * Δx -F(x) x_new x Δx其中 F 是所有方程的残差向量J 是 F 对 x 的偏导矩阵。按交流变量和直流变量分块J 长这样J [ J_AA J_AD ] [ J_DA J_DD ]J_AA 就是传统交流潮流的雅可比这部分代码可以直接复用。J_AD 是交流方程对直流变量的偏导J_DA 是直流方程对交流变量的偏导J_DD 是直流方程对直流变量的偏导。这两个耦合块是统一迭代的灵魂也是交替迭代法里根本没有的东西。很多教程讲到这里就停了只说把耦合项加上但不告诉你具体加什么。下面我把每一块都摊开。3.2 交流侧对直流变量的耦合项如何推导换流节点上直流注入交流的有功是P_dc U_d * I_d注意这是有功的绝对值符号靠前面约定的方向来定。对整流侧注入是 -U_dI_d对逆变侧是 U_dI_d。所以∂P_dc/∂U_d ±I_d ∂P_dc/∂I_d ±U_d ∂P_dc/∂α 0有功部分很干净控制角不出现在有功表达式里。无功部分麻烦一些因为要过功率因数这一关。换流器的功率因数有个漂亮的关系式可以从换流器方程反推出来cosφ U_d / (K * U_s)推导过程很简单cosα (U_d R_cI_d)/(KU_s)cos(αμ) (U_d - R_cI_d)/(KU_s)两者相加除以二就是 cosφ。这个式子说明什么说明换流器的功率因数完全由直流电压和交流母线电压的比值决定跟触发角没有直接关系。那么无功功率Q_dc P_dc * tanφ U_d * I_d * sqrt(1 - (U_d/(K*U_s))²) / (U_d/(K*U_s)) I_d * sqrt((K*U_s)² - U_d²)这个形式很好用偏导直接能算∂Q_dc/∂U_d -I_d * U_d / sqrt((K*U_s)² - U_d²) ∂Q_dc/∂I_d sqrt((K*U_s)² - U_d²) ∂Q_dc/∂U_s I_d * K² * U_s / sqrt((K*U_s)² - U_d²)最后一项 ∂Q_dc/∂U_s 是容易被漏掉的关键项。它落的位置是 J_AA 块里换流节点无功方程对自身电压幅值的偏导。你如果只按传统公式写 J_AA这一项就丢了结果是雅可比不精确牛顿法退化成拟牛顿法收敛次数从 5 次涨到 15 次甚至不收敛。我第一版就是这个问题收敛慢得离谱后来做数值差分对比才发现漏了这块。3.3 直流侧对交流变量的耦合项怎么处理直流方程里显式含有交流母线电压的只有换流器方程。以整流侧为例残差写成f_conv U_d - K * U_s * cosα R_c * I_d 0各偏导∂f/∂U_d 1 ∂f/∂U_s -K * cosα ∂f/∂α K * U_s * sinα ∂f/∂I_d R_c逆变侧符号略有不同cosα 换成 cosγR_c 项符号反过来。直流回路方程对交流变量没有偏导因为它只涉及直流量。控制方程也是定电流、定熄弧角、定功率都不含交流电压。所以 J_DA 块其实只有换流器方程那两行是稀疏的非零其余全是零。这个稀疏性在组装的时候要利用起来不然矩阵会白白变稠。3.4 标幺化处理与稀疏矩阵组装标幺化这件事我必须单独强调。交直流混合系统里交流侧电压是几百千伏量级直流电流是千安量级功率是兆瓦量级如果不统一基准雅可比矩阵的条件数会非常糟糕线性求解器给出的解精度很差。我的处理方式交流侧沿用系统基准 S_base 和交流电压基准 U_base阻抗基准 Z_base U_base²/S_base。直流侧取直流电压基准 U_dB 等于换流变交流侧电压基准折算到直流侧的值直流电流基准 I_dB S_base/U_dB直流电阻基准 R_dB U_dB/I_dB这正好等于 Z_base。这样所有量都在 1 附近雅可比元素数量级一致求解稳定。代码里稀疏组装用三元组方式% 预分配按最大非零元数估算 max_nnz 4 * n 20; rows zeros(max_nnz, 1); cols zeros(max_nnz, 1); vals zeros(max_nnz, 1); idx 0; % 交流块填充示例实际按行遍历更省事 for i 1:n-1 % ... 逐行填充 J_AA end % 最后一次性组装 J sparse(rows(1:idx), cols(1:idx), vals(1:idx), nvar, nvar);用 sparse 而不是 zeros 建满矩阵10 个节点看不出区别100 个节点就能差出一个数量级的求解时间。而且 Matlab 的 \ 运算符对稀疏矩阵会自动选合适的求解器比手动调 lu 靠谱。4. Matlab 实现全流程与代码骨架4.1 数据结构怎么组织才不容易乱我用四个结构体bus交流节点、branch交流支路、dc直流系统、opt求解选项。bus struct(num, [], ... % 节点编号 type, [], ... % 1PQ, 2PV, 3平衡, 4换流站 V, [], ... % 电压幅值初值标幺 theta, [], ... % 相角初值弧度 Pg, [], Qg, [], Pl, [], Ql, []); dc struct(conv, [], ... % 换流站数组每站一个结构 line, []); % 直流线路参数 dc.conv(k).bus_ac 5; % 接入的交流节点号 dc.conv(k).nT 1.0; % 换流变变比 dc.conv(k).Xc 0.1; % 换相电抗标幺 dc.conv(k).mode rect; % 整流还是逆变 dc.conv(k).ctrl I; % 控制变量类型 dc.conv(k).ctrl_val 1.0; % 控制目标值用结构体数组而不是元胞数组好处是可以直接用点号访问代码读起来清楚。换流站数量变化时不用改任何逻辑程序里全部按 length(dc.conv) 循环。有个细节换流站的交流接入节点类型要标成 4这样在主循环里能快速定位哪些节点需要加直流注入项、哪些节点要跳过电压幅值更新。用 type 字段区分比单独维护一个索引数组更不容易出错。4.2 主迭代循环的完整骨架整个牛顿迭代循环结构如下max_iter 20; tol 1e-8; x init_state(bus, dc); % 初始状态向量 history zeros(max_iter, 1); for iter 1:max_iter [F, J] build_residual_jacobian(x, bus, branch, dc, Ybus, ...); err max(abs(F)); history(iter) err; if err tol fprintf(第 %d 次迭代收敛最大残差 %.2e\n, iter, err); break; end if iter max_iter warning(达到最大迭代次数仍未收敛当前残差 %.2e, err); end dx -J \ F; % 稀疏左除自动选求解器 x x dx; % 越限检查和钳位 x enforce_limits(x, dc); end几个实操要点第一dx -J \ F这一行是稀疏求解的核心Matlab 会自动判断用 LU 还是 Cholesky不要手写 inv(J)*F数值上差很多。第二enforce_limits这个函数判断控制角的合理性。触发角一般限制在 5 到 90 度熄弧角限制在 10 到 60 度。牛顿迭代偶尔会给出负的触发角或超过 90 度的值这时候如果直接钳位残差会跳变迭代可能被带偏我的做法是钳位的同时把对应的控制方程临时改成角度不变相当于人为把这个变量冻结一步让其他变量先收敛。第三迭代日志一定要存。我把每次的残差存进 history 数组收敛完成后画个半对数图能直接看出收敛类型——二次收敛是斜率越来越陡线性收敛是直线振荡是锯齿。这个图对调试的价值极高。4.3 换流器残差与雅可比的子函数实现把换流器的计算单独封装成函数输入是直流状态量、交流电压输出是残差和对应的雅可比元素function [f, Jrow] conv_equation(Ud, Id, Us, alpha, K, Rc, mode) if strcmp(mode, rect) f Ud - K * Us * cos(alpha) Rc * Id; Jrow struct(dUd, 1, ... dUs, -K * cos(alpha), ... dalpha, K * Us * sin(alpha), ... dId, Rc); else f Ud - K * Us * cos(alpha) - Rc * Id; Jrow struct(dUd, 1, ... dUs, -K * cos(alpha), ... dalpha, K * Us * sin(alpha), ... dId, -Rc); end end无功对电压幅值的偏导那个 sqrt((KU_s)² - U_d²) 分母有开方当 U_d 接近 KU_s 时cosφ 接近 1换流器运行在接近单位功率因数的极端状态这个偏导会趋于零反过来 ∂Q/∂U_d 会趋于无穷。实际系统里 U_d/(K*U_s) 一般在 0.85 到 0.95 之间离奇异点还有距离但代码里要加个保护ratio Ud / (K * Us); ratio min(max(ratio, 0.5), 0.999); % 防止开方为负或除零这个保护不是多余的。我遇到过一次逆变侧电压被迭代推到很低、ratio 超过 1sqrt 直接返回复数整个迭代崩掉。加上钳位之后至少不会出复数虽然结果不物理但能让你从残差曲线看出是哪一步跑偏的。4.4 初值选取与收敛控制牛顿法对初值敏感交直流系统尤其如此。我的初值策略是分层给交流侧全部给 1.0 幅值、0 相角这个标准做法基本够用。直流侧要稍微算一下直流电压初值按额定值给直流电流初值由目标功率反算 I_d P_set / U_d整流侧触发角给 15 度约 0.26 弧度逆变侧熄弧角给 18 度。为什么给 15 到 18 度这个范围因为这是工程上实际运行点的典型值牛顿法从这附近出发基本一步就能进收敛域。你要是给 0 度或者 80 度第一第二个迭代步会跑得很远有时候直接跑到非物理解那边去了。收敛判据我用最大残差而不是二范数。二范数在大系统里会被节点数量放大阈值不好定最大残差直观1e-8 就是每个方程都满足到 8 位有效数字物理意义明确。还有一个技巧前两次迭代可以用放宽的雅可比也就是用初值处的雅可比多迭代几步省下重新组装和分解的时间。但对交直流系统来说耦合项变化快这个技巧收益不大我最后没采用老老实实每步重算。5. 实际调试中踩过的坑5.1 不收敛的几种典型症状和定位路径残差曲线是最好的诊断工具。我总结了几种形态和对应的病因残差表现可能原因排查手段首次迭代后残差反而变大雅可比符号错误或初值严重偏离数值差分校验雅可比残差线性下降每次只降一个数量级耦合项漏写退化成拟牛顿检查 J_AD 和 J_DA 是否为零残差在 1e-3 附近来回振荡控制方式在切换边界反复跳打印每步触发角看是否触发限值残差突然变成 NaN 或复数开方项为负、除零、角度越界在各处加 assert 断言前三次下降很快之后停滞雅可比接近奇异算 cond(J)检查是否有零对角元数值差分校验雅可比是排查的第一步也是最有效的一步。做法很简单对每个状态变量给一个 1e-6 的扰动用前向差分算数值偏导跟解析偏导逐个比较相对误差超过 1e-4 的就是有问题。我第二次遇到收敛慢的问题时就是这么定位的查出来是 ∂Q/∂U_s 那一项少了个 K² 因子。5.2 收敛到非物理解的表现和规避牛顿法有个讨厌的特性它只管把残差推到零不管解的物理意义。交直流系统里很容易出现这种情况——迭代收敛了残差 1e-10但打印出来发现整流侧触发角是负的 20 度或者逆变侧熄弧角是 3 度这明显是换相失败工况不可能是稳态解。原因通常有两个。一是初值离物理解太远牛顿步跨过了物理可行域落到了数学解的另一支上。二是控制方式配错了比如整流侧明明应该是定电流你按定功率写了方程系统会找到一个数学上满足定功率但物理上不可行的点。规避办法每次迭代后检查所有控制角的物理范围超出范围就打印警告并记录收敛后如果仍有越界就说明这次求解无效需要换初值重来。我在主程序里加了自动重启机制如果首次求解出现越界就把触发角初值从 15 度调到 20 度重新算一次通常第二次就能落在合理区域。5.3 结果校验的三种手段算完之后必须校验不能收敛了就用。我固定做三件事第一全系统有功平衡校验。所有交流电源的有功加上直流注入的有功减去负荷有功和网络损耗应该接近零。误差在 1e-6 量级说明有功方程没问题。这里要注意等值换相电阻 R_c 上的功率不是损耗不能算进网损这是新手常犯的错误。第二换流站交直流功率对照。交流侧算出来的注入有功跟直流侧 U_d * I_d 应该相等。这两条路径算出来的值如果差得多说明换流器方程或者功率因数公式用错了。我一般要求相对误差小于 1e-6。第三跟交替迭代的结果做交叉验证。同一组数据用两套方法算结果应该一致到小数点后四五位。这个交叉验证做过一次之后后面改代码就有信心了。5.4 我的几条实操心得换相电抗的单位是最容易出错的地方。数据手册上给的往往是折算到换流变阀侧的欧姆值你必须先折算到交流母线侧再除以阻抗基准变成标幺值。少除一个变比平方结果能差出百分之几而且不容易发现因为照样收敛。直流电阻的处理也一样双极系统的直流线路电阻是单极的一半两极并联金属回线方式下回路电阻还要加上回线电阻。这些细节数据文件里写错了结果就是错的程序本身不会报错。还有一点换流变压器的变比定义有歧义。有的资料定义成阀侧额定电压/交流侧额定电压有的定义成交流侧电压/阀侧电压两者互为倒数。用错了直流电压直接差一个变比平方功率就差得更多。我的建议是在数据文件里明确注释清楚定义并且用一个已知工况反算验证一遍。6. 这套代码后续能怎么扩展现在这套代码基本框架是稳的往几个方向扩展都不难。多端直流就是把直流网络从单支路推广到节点导纳矩阵形式方程结构不变只是变量和方程数量增加雅可比矩阵的稀疏性要重新组织一下。加入直流网络损耗的精细模型也不难把线路电阻从常数改成随温度变化的函数或者加上线路的分布参数效应。如果要做含新能源接入的混合系统难点不在潮流本身而在于新能源的出力不确定性——那就变成概率潮流或者场景分析了需要外层套蒙特卡洛或者点估计法。这时候统一迭代法的优势更明显因为单次求解快几千次场景扫描才跑得动。另外一个我觉得值得做的方向是把换流器的详细模型引进来考虑换相重叠角的动态变化和阀的导通压降。这样算出来的稳态点更接近实际尤其是低功率运行的时候理想模型的误差会比较大。代价是方程更复杂雅可比的解析推导工作量翻倍可能需要借助符号计算工具辅助推导。最后分享一个我一直在用的小习惯每次改动代码之后先用一个 3 节点的小系统跑一遍冒烟测试确认能收敛、功率平衡满足再去跑 30 节点以上的实际算例。小系统的迭代过程可以直接在命令行里 print 出来一眼就能看出哪一步不对大系统一旦出问题定位成本高十倍。这个习惯帮我省下的时间比我优化代码省下的时间多得多。