1. 交直流潮流为什么最后落在了统一迭代上第一次把交直流混合系统往 Matlab 里塞的时候我的做法和大多数人一样交流侧用一套牛顿法算完把换流站交流侧的有功当已知量丢给直流网络直流网络解出来再回代更新换流站功率来回往复直到两次迭代的换流站功率差小于阈值。这套交替迭代法也有人叫它顺序求解法跑纯主从控制的简单算例时确实舒服代码短、调试直观、交流子程序和直流子程序各自独立出问题一眼就知道是哪半边。但只要换流站数量一多、控制方式一混搭、再加上换流器损耗随电压波动这种细枝末节交替迭代就开始露怯迭代次数从七八次涨到三四十次偶尔还会在两个解之间来回振荡你把收敛阈值放宽一点它就收敛了检查结果发现直流电压偏离参考值上千分之几功率平衡差得离谱。基于 Matlab 的交直流潮流计算统一迭代法建模与实现本质上就是在讲一件事把交流节点方程和直流节点方程塞进同一个牛顿迭代框架里用一个大雅可比矩阵同时求解用额外的矩阵维度换取收敛的鲁棒性和解的唯一性。它不算什么高深算法真正难的地方全在建模细节上——换流站的稳态方程怎么摆、标幺基准怎么统一、耦合偏导要不要留、直流松弛节点怎么处理。这些细节一旦搞错程序照样能跑出数来但结果是错的而且错得很隐蔽。这篇东西适合三类人看第一类是在做交直流混合系统相关课题、需要自己写潮流程序而不是只会点软件按钮的第二类是想搞清楚为什么我的交替迭代不收敛的第三类是有一定 Matlab 基础、能看懂矩阵运算但没接触过电力系统潮流内部结构的。我会把方程怎么来、雅可比怎么装、代码怎么写、坑在哪里踩一条线讲下来。默认你会基础的交流牛顿-拉夫逊潮流如果连 Ybus 怎么组、H 矩阵怎么来的都不清楚建议先把纯交流版本手推一遍再往下看不然读到雅可比那节会很痛苦。需要提前说明一点下面给出的算例参数和迭代记录是我为了讲清楚流程自己搭的一个小系统不是任何实际工程数据你在复现时数值对不上很正常重点看结构和逻辑。2. 方案选型交替迭代与统一迭代到底差在哪2.1 交替迭代法在工程里卡在什么地方交替迭代的思路非常工程师把大问题切成两个能各自求解的小问题轮流解用边界量互相传递。它的收敛性依赖一个隐含假设——交流侧和直流侧的耦合不能太强。主从控制下定直流电压的那个换流站几乎把两侧钉死了耦合确实弱交替两三轮就稳。但一旦上 droop 控制直流电压本身参与功率分配换流站有功不再由交流侧给定而是由直流电压决定此时交流方程里出现未知量 Vdc直流方程里出现未知量 Vm两个子系统互相牵着走交替迭代的收敛速度会肉眼可见地掉下来。我在一个五端直流的算例上做过对比全主从控制时交替迭代 6 次收敛阈值 1e-8把其中两端改成下垂控制后交替迭代需要 31 次而且下垂系数取大一点就发散。这不是算法实现有 bug而是交替迭代本质上是一个块高斯-赛德尔迭代耦合项的谱半径随耦合强度增大而逼近 1收敛自然就慢。统一迭代相当于把它换成了整体牛顿二次收敛的性质量级上不一样。还有一个更现实的问题交替迭代的收敛判据不好定。你按换流站功率差判断可能出现功率差已经很小但电压还在漂的情况你按电压差判断又可能在功率严重不平衡时误判收敛。统一迭代没有这个烦恼不平衡量向量统一范数一个判据管到底。2.2 统一迭代付出的代价值不值代价很直接雅可比维度变大装配逻辑变复杂调试难度上升。一个 nb 个交流节点、nd 个直流节点的系统纯交流雅可比是 (2nb-1) 阶左右加上直流后大概多出 nd 阶而且多出来的耦合块元素推导容易出错。更要命的是一旦某个耦合偏导符号写反程序不一定发散而是收敛到一个物理上不成立的解你拿万用表都量不出来。但从我自己的经验看这笔账是划算的。第一牛顿法对初值敏感但一旦进入收敛域就是二次收敛交直流系统初值可以很方便地用平启动 直流电压全置 1.0给出实际很少不收敛第二统一迭代天然支持各种控制方式混搭主从、下垂、定功率、定电压只需要改对应行的方程框架不动第三扩展性好以后要加换流器限流、交流侧无功限值、直流线路损耗都只是在已有方程上加项不像交替迭代那样要重新设计外层循环。代价还有一条容易被忽略统一迭代对模型的完整性要求高。交替迭代可以容忍某一侧用简化模型比如直流侧只当恒功率负载统一迭代不行你写进雅可比的每一项都必须和方程严格对应否则收敛性立刻崩。这其实是好事逼着我把换流器方程老老实实推了一遍。2.3 什么情况下我会劝你别用统一迭代不是所有场景都值得上统一迭代。如果系统里只有一个背靠背换流站两侧交流网络规模都很小交替迭代两三次就收敛代码量还只有统一迭代的三分之一那就没必要折腾。如果直流侧是完全理想的恒功率负载模型损耗忽略、控制理想交直流之间严格解耦统一迭代的雅可比是分块对角的等于把两个独立问题拼在一起解白白浪费计算量。真正值得上统一迭代的典型场景是多端直流、下垂控制、换流站数量大于等于三、需要严格校验功率平衡、后续要做 N-1 扫描或连续潮流。这种情况下交替迭代的迭代次数不稳定批量计算时总有几个算例卡住统一迭代的稳定收敛反而是省时间的。对比维度交替迭代法统一迭代法雅可比规模两个小矩阵分别求解一个大矩阵一次求解收敛阶数线性块高斯-赛德尔二次牛顿控制方式混搭需要改外层调度逻辑只改对应行方程下垂控制支持收敛慢易振荡与主从控制同框架调试难度低两侧独立高耦合项易错代码量我的实现约 300 行约 600 行扩展新设备改动大加行加列即可3. 模型拆解把换流站写成能进雅可比的方程3.1 节点分类与未知量编号统一迭代第一个要想清楚的问题不是方程而是谁是未知量。交流侧沿用经典分类PQ 节点未知量是电压幅值 Vm 和相角 VaPV 节点未知量只有相角幅值由控制给定平衡节点两个都给定。直流侧稍微特殊一点直流网络没有相角概念只有电压幅值 Vdc所以每个直流节点的未知量就是 Vdc 本身但定直流电压的那个节点例外——它的 Vdc 已知功率未知等价于交流侧的平衡节点。这个例外是很多人第一次写统一迭代时栽跟头的地方。直觉上你会觉得直流节点都应该有功率方程于是给定了电压的节点也列了一条 ΔPdc 方程结果雅可比对应行是零行或者接近零行矩阵直接奇异。正确做法是定电压节点的 Vdc 从未知量集合里剔除对应的 ΔPdc 方程也不列它的功率在收敛后由 Pdc Vdc * Σ G Vdc 反算出来。我用的是未知量索引映射的办法先扫一遍表把所有未知量编号再用映射数组把物理节点号翻译成未知量序号。这样后面装配雅可比时不用到处判断类型索引一查就行。具体做法在第四节给代码。3.2 换流站的稳态方程怎么摆换流站建模有两条路平均模型和开关模型。潮流计算用平均模型忽略开关过程和 PWM 谐波只保留基波功率平衡。我采用的是最常见的功率平衡写法交流侧注入电网的有功 Pc直流侧从直流网络取用的有功 Pdc。两者之差是换流器损耗Pdc Pc Ploss。交流侧无功 Qc 由控制方式决定定 Q 控制时直接给定。损耗模型取 Ploss a b·I c·I²其中 I 是换流站交流侧电流的标幺值I sqrt(Pc² Qc²) / Vm。选这个损耗模型不是因为它精确而是因为它同时包含了常数项、线性项和二次项能覆盖导通损耗和开关损耗的主要成分而且对 Vm 的偏导好推——这一点在装配耦合块时非常关键。有些文献用恒定效率模型Pdc Pc / η简单但那样 Ploss 对 Vm 的偏导恒为零交直流之间就没有耦合项了统一迭代退化成两个独立问题失去了演示意义。注意损耗模型的选择直接决定耦合块是否存在。如果你的算例里所有换流站都用恒定效率、所有控制都是定 P 定 Q那统一迭代的雅可比是分块对角的收敛性和交替迭代一样这时候没必要上统一迭代用它只是自己给自己找麻烦。3.3 标幺基准怎么统一交直流混合系统最容易出错的地方是基准值。交流侧一般取三相功率基准 Sb 和线电压基准 Vb直流侧如果直接沿用交流的 Sb直流电压基准要单独取——我习惯取 Vdb 使得直流侧基准阻抗与交流侧一致也就是 Vdb² / Sb Vb² / Sb这样 Vdb Vb这里的 Vb 是交流相电压基准若用线电压基准则要乘系数。如果不统一直流网络的电导矩阵 G 数值会偏离好几个数量级雅可比条件数爆炸。我吃过一次亏直流侧电导一开始用的是实际欧姆值没转标幺结果雅可比对角元量级从 1e1 跳到 1e4稀疏求解器虽然没报错但迭代过程出现了明显的数值抖动收敛判据在 1e-6 附近停住不动。改成标幺后一次就收敛了。所以我的建议是建模阶段就把所有量统一到同一套基准下代码里只出现标幺值物理量只在输入输出环节转换。量基准转换关系常见错误交流功率SbP_pu P / Sb单相/三相混淆交流电压Vb线电压Vm_pu V / Vb用相电压基准直流电压Vdb VbVdc_pu Vdc / Vdb直接用电网额定值直流电阻Zdb Vdb²/SbR_pu R / Zdb忘了转标幺换流器电流Ib Sb/(√3·Vb)I_pu I / Ib与功率混淆损耗系数按 Sb 归一a_pu a / Sb系数单位不统一3.4 完整不平衡量向量的构成把所有方程按顺序排好不平衡量向量 ΔF 分成三块第一块是交流有功不平衡 ΔP对所有非平衡交流节点列共 (nb - 1) 条。对纯负荷节点ΔP_i P_i^calc - P_i^spec对挂换流站的节点注入功率里要减去换流站交流侧取用的 Pc。第二块是交流无功不平衡 ΔQ对所有 PQ 节点列。挂换流站的节点如果换流站定 Q 控制那 Qc 作为给定值进 P_spec 的位置。第三块是直流功率不平衡 ΔPdc对所有非定电压直流节点列ΔPdc_k Vdc_k·Σ_m G_km·Vdc_m - Pdc_k^spec。这里有个符号约定的坑必须讲清楚。换流站整流交流到直流时 Pc 为正、Pdc 为负还是反过来不同文献写法不一样。我统一约定Pc 为换流站从交流电网吸收的有功正表示吸收Pdc 为换流站注入直流网络的有功正表示注入那么整流时 Pc 0、Pdc 0 的绝对值小于 Pc关系是 Pdc -(Pc - Ploss)。逆变时反过来。约定一旦定下来全套代码都要遵守中间改一次符号所有耦合项都得重推。4. 雅可比矩阵的结构与装配细节4.1 分块结构与稀疏索引统一迭代的雅可比是一个三行三列的分块矩阵J [H N Ndc; M L Ldc; Md Hd D ]H、N、M、L 是大家熟悉的交流子块分别是 ΔP 对 Va、ΔP 对 Vm、ΔQ 对 Va、ΔQ 对 Vm 的偏导。D 是直流子块ΔPdc 对 Vdc 的偏导。剩下四块是耦合块Ndc 是 ΔP 对 VdcLdc 是 ΔQ 对 VdcMd 是 ΔPdc 对 VaHd 是 ΔPdc 对 Vm。在我采用的换流站模型下Ndc 和 Ldc 恒为零——因为交流侧注入功率 Pc、Qc 是控制给定的不直接依赖直流电压。Md 也是零——直流功率不依赖交流相角。只有 Hd 非零它来自损耗项对交流电压幅值的依赖。这个结论很重要它说明统一迭代的耦合其实只有一个块在起作用如果你把这个块也简化掉损耗恒定矩阵就是分块下三角可以用两次回代求解本质上又退回了交替迭代。所以我的原则是耦合块哪怕数值很小也必须保留它是统一迭代区别于交替迭代的物理来源。实际算例里 Hd 的元素量级大概在 1e-3 到 1e-2确实小但正是这个小项在迭代初期提供了方向信息减少了来回振荡。4.2 交流子块的经典写法回顾交流部分沿用标准极坐标牛顿法公式我不重复推导只强调两个实现细节。第一Ybus 要用稀疏矩阵Matlab 里 sparse 存储装配功率时用 Ybus 乘以电压向量别用循环。第二对角元和非对角元的偏导公式不同装配时按行列分别处理不要图省事用一个循环套所有元素。H 的对角元H_ii -Q_i - B_ii·Vm_i² H 的非对角元H_ij Vm_i·Vm_j·(G_ij·sin θ_ij - B_ij·cos θ_ij)N 的对角元N_ii P_i/Vm_i G_ii·Vm_i N 的非对角元N_ij Vm_i·(G_ij·cos θ_ij B_ij·sin θ_ij)M、L 类似符号别搞反。这部分是教科书内容但我在实际 debug 中发现很多人栽在 θ_ij 的定义上——θ_ij θ_i - θ_j 还是 θ_j - θ_i两种写法都有人用混用会导致 H 的非对角元整体符号反了程序仍然收敛但结果完全错误。我的做法是全程用 θ_ij θ_i - θ_j并且写一个单元测试两节点系统手算一遍对拍。4.3 直流子块与耦合子块的偏导推导直流注入功率 Pdc_k Vdc_k · Σ_m G_km · Vdc_m。对 Vdc_m 求偏导m ≠ k 时∂Pdc_k/∂Vdc_m Vdc_k · G_kmm k 时∂Pdc_k/∂Vdc_k Σ_m G_km·Vdc_m Vdc_k·G_kk这个对角元很容易写错成 2·Vdc_k·G_kk或者漏掉求和项里的 k 项。正确形式是该节点注入电流 自身电导乘电压物理上就是电流对电压的灵敏度。耦合块 Hd 的推导ΔPdc_k 里包含 Pdc_k^spec而 Pdc_k^spec Pc_k Ploss_k按前面的符号约定取绝对值Ploss_k a b·I c·I²I sqrt(Pc² Qc²)/Vm_iVm_i 是换流站所接交流节点的电压幅值。∂Ploss/∂Vm (b 2c·I) · ∂I/∂Vm (b 2c·I) · (-I/Vm)所以 ∂ΔPdc_k/∂Vm_i -∂Ploss/∂Vm_i (b 2c·I)·I/Vm_i。符号为正物理意义是电压幅值升高、电流减小、损耗减小、直流侧所需注入功率增加合理。提示这个偏导的量级很小如果你在调试时发现它一直是零先检查是不是损耗系数 b、c 都设成了零或者 Vm 初值取 1.0 导致灵敏度看起来没变化。可以用数值差分验证一遍方法在下一节。4.4 用数值差分给雅可比做体检手推偏导最怕符号错。我的习惯是写完解析雅可比后立刻用中心差分对拍一次误差超过 1e-4 就说明有地方错了。代码不到二十行但能省掉几个晚上的调试function Jnum jac_check(x, Ffun, nx, h) % x : 当前状态向量 % Ffun: 不平衡量函数句柄 F Ffun(x) % nx : 未知量个数 % h : 差分步长取 1e-6 左右 F0 Ffun(x); Jnum zeros(numel(F0), nx); for k 1:nx xp x; xp(k) xp(k) h; xm x; xm(k) xm(k) - h; Jnum(:, k) (Ffun(xp) - Ffun(xm)) / (2*h); end end用法是先在第一次迭代时调用把解析雅可比 J 和 Jnum 相减取绝对值最大元素打印出来。如果最大值出现在耦合块对应的行说明耦合偏导有问题如果出现在交流子块说明 Ybus 或者 θ 约定有误。差分步长不宜太小1e-6 到 1e-8 之间试太大会有截断误差太小会撞上双精度舍入误差。5. Matlab 实现从数据结构到求解循环5.1 三张表的数据结构我不用面向对象也不用 struct 数组就用三个二维矩阵列定义写在注释里。好处是索引快、易于向量化、保存成 mat 文件方便。缺点是列含义靠记忆所以我把列定义固化成常量。% bus : 交流节点表每行一个节点 % 1 id | 2 type(1PQ,2PV,3平衡) | 3 Vm | 4 Va(rad) % 5 Pg | 6 Qg | 7 Pd | 8 Qd | 9 conv_id(0无换流站) % % dc : 直流节点表 % 1 id | 2 Vdc | 3 type(1定电压,2定功率,3下垂) | 4 Pdc_spec % 5 kdroop | 6 Vdc_ref | 7 conv_id % % conv : 换流站表 % 1 id | 2 ac_bus | 3 dc_node | 4 mode(1定P定Q,2定Vdc) % 5 Pset | 6 Qset | 7 a | 8 b | 9 ctype字段我坚持用整数而不是字符串因为后面要大量做逻辑判断和索引字符串比较慢且容易拼错。conv_id这个外键放在 bus 和 dc 表里是双向索引的关键装配耦合块时靠它从交流节点找到对应的直流节点反之亦然。5.2 未知量映射与主循环骨架主循环的结构固定算不平衡量 → 判断收敛 → 装雅可比 → 解修正方程 → 更新状态。真正需要小心的是状态更新的边界处理以及修正量对应的物理量。function res acdc_nr(bus, dc, conv, opt) nb size(bus,1); nd size(dc,1); % 交流未知量索引 iVa find(bus(:,2) ~ 3); % 非平衡节点的相角 iVm find(bus(:,2) 1); % PQ 节点的幅值 % 直流未知量索引剔除定电压节点 iVd find(dc(:,3) ~ 1); nx numel(iVa) numel(iVm) numel(iVd); % 状态向量 x 与物理表的映射 mapVa zeros(nb,1); mapVa(iVa) 1:numel(iVa); mapVm zeros(nb,1); mapVm(iVm) numel(iVa) (1:numel(iVm)); mapVd zeros(nd,1); mapVd(iVd) numel(iVa)numel(iVm) (1:numel(iVd)); Vm bus(:,3); Va bus(:,4); Vdc dc(:,2); Y make_ybus(bus); G make_gdc(dc); for it 1:opt.maxit [dP, dQ] ac_mismatch(Y, Vm, Va, bus, conv, Vdc); dPdc dc_mismatch(G, Vdc, dc, conv, Vm); F [dP(iVa); dQ(iVm); dPdc(iVd)]; if norm(F, inf) opt.tol break; end J build_jac(Y, Vm, Va, G, Vdc, bus, dc, conv, ... iVa, iVm, iVd, mapVa, mapVm, mapVd); dx J \ (-F); % 稀疏左除 % 回代更新 Va(iVa) Va(iVa) dx(mapVa(iVa)); Vm(iVm) Vm(iVm) dx(mapVm(iVm)); Vdc(iVd) Vdc(iVd) dx(mapVd(iVd)); end res pack_result(bus, dc, conv, Vm, Va, Vdc, it, norm(F,inf)); end这段骨架里有个细节值得说J \ (-F)用的是左除Matlab 会自动识别 J 是稀疏矩阵并调用合适的求解器。如果你写inv(J)*(-F)一是慢二是精度差三是在大系统上会直接内存溢出。另外norm(F, inf)取无穷范数比二范数更严格因为它盯的是单个最大的不平衡量不会被平均掉。5.3 交流不平衡量的计算function [dP, dQ] ac_mismatch(Y, Vm, Va, bus, conv, Vdc) V Vm .* exp(1j*Va); S V .* conj(Y * V); P real(S); Q imag(S); Ps (bus(:,5) - bus(:,7)) / 1.0; % 给定注入已标幺 Qs (bus(:,6) - bus(:,8)) / 1.0; % 换流站交流侧取用功率计入负荷 for k 1:size(conv,1) i conv(k,2); c k; [Pc, Qc, ~] conv_ac_power(conv(c,:), Vm(i), Vdc(conv(c,3))); Ps(i) Ps(i) - Pc; Qs(i) Qs(i) - Qc; end dP P - Ps; dQ Q - Qs; end这里conv_ac_power是核心它要根据控制方式返回换流站交流侧实际取用的有功和无功。定 P 定 Q 模式下 Pc Pset、Qc Qset直接返回定直流电压模式下Pc 是由直流侧功率反推的未知量不能在这里给死需要靠方程间的关系解决——这也是统一迭代相比交替迭代省事的地方未知量 Pc 不出现在状态向量里但它通过直流方程隐含求解。严格来说定电压站的 Pc 在每次迭代时应该由当前 Vdc 和损耗反算Pc Pdc - Ploss而 Pdc Vdc·ΣG·Vdc。这就是所谓的内层更新不影响雅可比结构但影响收敛速度。我一开始没做这个内层更新直接用上一次迭代的 Pc 值结果迭代次数多了快一倍。5.4 直流不平衡量与稀疏装配function dPdc dc_mismatch(G, Vdc, dc, conv, Vm) Idc G * Vdc; Pdc Vdc .* Idc; Pspec zeros(size(Vdc)); for k 1:size(conv,1) j conv(k,3); [~, ~, Pdc_k] conv_ac_power(conv(k,:), Vm(conv(k,2)), Vdc(j)); Pspec(j) Pdc_k; end dPdc Pdc - Pspec; end雅可比装配我采用三元组预分配的方式先用spalloc估一个非零元上界用三个长向量 irow、jcol、jval 累积最后一次性sparse。这种方式比逐元素赋值快一个数量级在节点数上百之后差别非常明显。nnz_est 4*(nb nd) 16*size(conv,1); irow zeros(nnz_est,1); jcol irow; jval irow; cnt 0; add (r,c,v) deal_helper(); % 这里用局部函数累加略实际写的时候我把累加封装成局部函数避免匿名函数带来的性能损失function push(r, c, v) cnt cnt 1; irow(cnt) r; jcol(cnt) c; jval(cnt) v; end这种写法在嵌套函数里能直接改外部变量比传结构体清爽。最后J sparse(irow(1:cnt), jcol(1:cnt), jval(1:cnt), nx, nx)。5.5 收敛判据与迭代控制判据我一般用两层主判据是无穷范数小于 1e-8标幺辅助判据是修正量范数小于 1e-6两个都满足才认为收敛。只判不平衡量的话在病态系统上可能出现不平衡量很小但电压还在缓慢漂移的情况。迭代次数上限设 30超过就认为不收敛打印最后一次的不平衡量分布看是哪个节点贡献最大。这个诊断信息非常有用——如果最大不平衡量集中在某个换流站节点八成是控制方式设置和方程不匹配如果分散在所有节点八成是初值或者标幺化的问题。还有一个实用技巧记录每次迭代的不平衡量范数画成半对数图。正常牛顿法应该是一条陡峭的直线下降如果出现平台期、抖动或者先降后升说明雅可比有问题或者系统本身无解。这个图我几乎每次调试新模型都会画比看数字直观得多。6. 一个可复现的小算例与结果校验6.1 算例配置我搭了一个五节点交流 三节点直流的系统规模小到能手算校验又能体现出交直流耦合的特征。交流侧节点 1 是平衡节点节点 2、3 是纯负荷节点 4、5 各接一个换流站。直流侧节点 1 是定电压站对应交流节点 4 的换流站节点 2、3 是定功率站对应交流节点 5 的换流站和一个纯直流负荷。参数数值说明交流基准功率100 MVA全系统统一交流基准电压10 kV线电压交流线路阻抗0.01 j0.08 pu环形结构负荷节点 20.6 j0.2 pu恒功率负荷节点 30.4 j0.15 pu恒功率换流站 1 控制定 Vdc 1.0主从控制的从站换流站 2 控制定 P -0.3 pu逆变运行直流线路电阻0.02 pu三节点三角形损耗系数 a0.002 pu常数项损耗系数 b0.01线性项损耗系数 c0.02二次项需要注意的是这里的定功率是站在直流侧说的换流站 2 从直流网络取 0.3 pu 的功率送到交流侧同时因为损耗它从直流网络实际取用的功率要略大于 0.3。6.2 迭代过程记录平启动所有电压置 1.0相角置 0。前三次迭代的不平衡量范数大致是迭代次数不平衡量范数说明06.2e-1初始不平衡主要是负荷节点18.7e-2大幅下降牛顿方向正确24.1e-3进入二次收敛区39.6e-6继续下降45.2e-11收敛这个下降模式是典型的牛顿法特征。如果你跑出来第一次迭代范数反而升了多半是雅可比符号错了或者修正量没有乘上正确的方向。我特意对比了一下忽略损耗b c 0的情况迭代次数变成 3 次收敛但收敛后的直流电压和考虑损耗时差了约 0.008 pu。这个差别在电压等级上就是 80 V 左右对于某些对电压精度要求高的场景是不能忽略的。6.3 结果校验的三种手段第一种是功率平衡校验。全系统的总有功注入应该等于总负荷加总损耗交直流两侧分别算误差应该在 1e-6 以内。这个方法能抓出大部分符号错误和单位错误。第二种是对拍。我用一个最简单的两节点交直流系统一个换流站、一条直流线路把方程手工展开成二元方程组用 fsolve 解一遍和自己的牛顿法结果对比。fsolve 不依赖雅可比结果可信度高。这个方法有点笨但对新模型非常有效。第三种是扰动测试。把某个负荷增加 1%看直流电压的变化方向是否符合物理直觉。逆变侧负荷增加直流电压应该下降因为要从直流网络多取功率线路压降增大。如果方向反了说明某个偏导符号错了哪怕收敛结果看起来挺合理。提示三种方法里我觉得扰动测试最有价值。功率平衡校验只能抓失衡抓不出符号错误对拍要写两套代码成本高扰动测试几行代码就能做而且直接检验物理正确性。7. 常见问题与排查实录7.1 不收敛的六类原因速查现象可能原因排查方法处理第一次迭代范数就变大雅可比符号反了数值差分对拍逐个块检查范数卡在 1e-3 不动耦合块漏了或为零打印耦合块元素补上损耗偏导迭代振荡不收敛下垂系数过大减小 kdroop 试或改用主从控制矩阵奇异报错定电压节点未剔除检查未知量索引剔除对应行列收敛但功率不平衡标幺基准不统一检查 G 矩阵量级统一转标幺收敛到明显错误解换流站符号约定混乱检查 Pc/Pdc 正方向全程统一约定这六类我基本都踩过。印象最深的是第三类在一个五端直流算例上下垂系数取 20 时怎么都不收敛取 5 就正常。当时的直觉是系数大响应快应该更好收敛实际上下垂系数大意味着直流电压对功率变化更敏感雅可比对角占优性变差牛顿法反而更容易跳出收敛域。这也是为什么实际工程里下垂系数都要留裕度。7.2 雅可比奇异与病态除了定电压节点不剔除导致的六异还有一种更隐蔽的情况两个换流站接在同一个交流节点上。这时交流侧注入功率是两者之和雅可比对应行的偏导也是两者叠加如果你的 conv_id 映射写错把两个站当成一个处理矩阵会出现重复行奇异但不报错Matlab 会给警告但很多人不看警告。病态则多半来自量级差异。交流侧的偏导量级在 1e0 到 1e1直流侧的电导偏导在 1e1 到 1e2耦合块在 1e-3。这种跨四个数量级的矩阵条件数很容易上 1e6虽然稀疏求解器能处理但精度会下降。我的做法是迭代前对雅可比做一次行列均衡Matlab 的equilibrate或者手动缩放实测能把条件数降一到两个数量级。7.3 结果看着对但不对的隐性错误最麻烦的不是不收敛而是收敛到一个错误解。我遇到过两次一次是 θ_ij 的定义搞反交流潮流的相角全部反号但因为系统对称功率分布看起来完全正常直到我检查某条线路的功率方向才发现不对。另一次是直流网络的 G 矩阵没有把对角元置为负的非对角元之和导致直流注入功率算出来偏大一倍。这个问题在纯直流测试时没暴露因为测试用例只有一条线路一接上交直流系统收敛后的直流电压全部偏低我才回头查 G 矩阵的组装。预防这类问题的办法只有一个每个模块写完就单独测别等整体跑通再测。Ybus 组完先跑个纯交流潮流G 矩阵组完先跑个纯直流网络都对上了再拼起来。拼起来之后再用手算小算例对拍一次。7.4 性能与工程化的几点建议节点数上到几百之后性能瓶颈会从算法转移到实现。几条经验Ybus 和 G 矩阵必须稀疏装配用三元组一次性构造别用循环逐元素赋值。不平衡量计算里避免 for 循环遍历换流站改成向量化索引换流站数量少的时候差别不大上百之后能差几倍。雅可比装配是热点能复用的部分比如交流子块的结构在迭代间不变可以只更新数值不重建稀疏结构。如果要批量扫算例比如 N-1 分析把外层循环改成 parfor注意每个 worker 内部别共享大矩阵内存会爆。收敛判据不要用相对误差标幺制下绝对误差更直观1e-8 对标幺量来说已经足够严格。8. 我个人踩过之后想说的几句这套代码我前后重构过三次。第一次是照着交替迭代的思路硬改成统一迭代结果雅可比装错了收敛是收敛了但结果和交替迭代对不上查了两天才发现耦合块的符号反了。第二次是把换流站方程重推了一遍用数值差分逐块验证这次对了但代码写得极其啰嗦每次加新控制方式都要改五六个地方。第三次才想明白控制方式的变化本质上是哪条方程换成哪个形式于是把方程抽象成可配置的模板主循环一行没动扩展新控制方式只加一个分支。如果让我给刚开始做这块的人一条建议那就是先把两节点交直流系统做到绝对正确再往上加节点。两节点系统里每一个偏导都能手算每一个符号错误都会立刻表现为不收敛或者明显偏离这时候调试成本最低。等系统大了错误会被网络的平均效应掩盖你看到的只是收敛慢或者结果差一点点根本定位不到问题。另外别迷信商用软件的结果。我用某商业仿真软件对拍过一个算例两边差了 0.3%查了半天发现是损耗模型不同——它默认用恒定效率我用的是二次损耗。这种差异不是 bug但如果你不知道就会以为是自己的程序错了然后去改一个本来正确的公式。模型对齐比结果对齐重要。最后说个小事Matlab 的稀疏左除在某些版本上对近奇异矩阵会给 warning 但不报错如果你的循环里没有检查lastwarn很可能跑完一整批算例都不知道其中有几个是错的。我现在的做法是每次求解后强制检查rank(full(J))小系统或者条件数估计超阈值就直接标记该算例为可疑不让它污染统计结果。这个习惯帮我省了很多返工。
