基于Matlab的交直流潮流计算:统一迭代法建模与实现 “交直流潮流计算”这个话题我在实际项目里踩过不少坑。最典型的一个就是交流潮流用牛拉法算得好好的结果一接入直流输电模型交替迭代怎么都不收敛最后发现是换流站的耦合关系没处理好。后来换了统一迭代法的思路问题一下清爽很多。这篇文章就把我基于 Matlab 实现交直流潮流计算的完整过程拆开讲清楚重点就是统一迭代法的建模、代码实现和常见坑位。如果你正在做直流输电接入后的潮流分析或者刚接触换流站建模这篇文章应该能给你省下不少排查时间。1. 统一迭代法为什么是交直流潮流计算的“正规打法”1.1 交直流混合电网的潮流问题到底难在哪普通交流潮流计算节点类型无外乎 PQ、PV、平衡节点三类核心是用牛顿-拉夫逊法求解节点功率平衡方程。但加上直流输电确切说是换流站之后问题不再简单的是“交流网络加个负荷”换流站本身是个强非线性元件。交流侧的有功和无功注入取决于直流侧电压、直流电流、触发角/熄弧角以及换流变压器变比。这中间任何一个量变了交流侧注入功率都会跟着变反过来交流母线电压一变换流器方程里的理想空载直流电压也会变。也就是说交流量和直流量是双向咬死的不是简单的“先算直流、再算交流”就能解耦。控制方式还会改变系统自由度。比如整流侧定电流、逆变侧定熄弧角这种情况下直流电压、电流之间是强约束关系你如果还按普通 PQ 节点去处理换流母线模型根本对不上。更麻烦的是换流器的无功需求通常很可观。一个 0.5 p.u. 传输功率的换流站无功消耗可能达到 0.2~0.3 p.u.。这个无功不是恒定不变的它随直流功率和交流电压非线性变化。交流潮流里如果处理不好这个动态的无功负荷电压可能一路跌下去迭代直接发散。1.2 序贯迭代与统一迭代的差异处理交直流潮流业内有两大路线序贯法和统一法。很多初学者第一反应是“两边分开算交流算一步、直流算一步多换几次不就完了”这就是序贯法的朴素想法。它把电力系统潮流方程拆成交流子系统和直流子系统分别迭代通过换流器边界条件交换信息。序贯法实现简单早期很多教科书也这么讲但它有个致命弱点在重负荷或换流器控制方式切换的工况下交流子系统和直流子系统之间的耦合变得非常紧交替迭代经常出现振荡甚至发散。我当年第一次做两端直流互联系统的潮流用序贯法调了整整两天换了三个初值组合始终在两个状态之间来回跳最后才明白问题出在算法结构而不是参数上。统一迭代法联立法的做法是把交流节点方程和直流系统方程写成一组联立的非线性方程组用一个统一的牛顿迭代框架同时求解。变量里既有交流的相角、电压幅值也有直流的电压、电流、控制角。这组方程在数学上是一个整体不再区分“先交流后直流”。这样做的好处是收敛性显著改善尤其是在强耦合、重负载工况下。但严格意义上的统一法要扩展雅可比矩阵把直流变量对交流不平衡量的偏导数全部算出来实现复杂度直接上了一个台阶。实际工程里还有一种折中方案我把它叫“解耦统一法”交流侧仍用牛拉法迭代但每一轮交流迭代内部都把直流方程精确收敛一遍再用当前的直流状态去更新换流器的等效注入功率。这样既保留了“统一求解”的数学本质直流方程作为潮流方程的约束始终被满足又避免推导庞大扩展雅可比矩阵。三种方法对照如下方法迭代结构收敛性实现难度适用场景序贯法交流直流交替迭代弱耦合时可用强耦合易发散低简单直流、弱耦合系统严格统一法全变量联立牛顿迭代最好二次收敛性接近理论高研究级工具、多端直流复杂系统解耦统一法交流牛拉内嵌直流精确求解较好通常比序贯法稳中工程仿真、Matlab快速验证我在本文接下来的实现中采用第三种方案。它是我在实际项目里性价比最高的选择Matlab 代码量可控没有强写雅可比子矩阵的推导负担同时保留统一迭代法的收敛韧性。2. 核心数学模型先把三块方程写清楚2.1 交流子系统方程与节点功率不平衡量交流部分沿用极坐标形式的牛拉法。对每个节点有功和无功的潮流方程写成P_i V_i * Σ_j V_j * (G_ij * cos(θ_i - θ_j) B_ij * sin(θ_i - θ_j))Q_i V_i * Σ_j V_j * (G_ij * sin(θ_i - θ_j) - B_ij * cos(θ_i - θ_j))其中 G_ij、B_ij 是节点导纳矩阵的实部和虚部V_i、θ_i 是节点电压幅值和相角。迭代过程中我们需要计算功率不平衡量ΔP_i P_i^spec - P_i^calcΔQ_i Q_i^spec - Q_i^calcP_i^spec、Q_i^spec 是给定额定注入功率。在没有换流器之前这个量来自发电机出力和负荷接入换流器之后需要在对应节点上叠加换流器的等效注入功率。这里要特别提一个细节很多人在写程序时把负荷功率的符号弄反。我的习惯是统一用“节点注入功率为正”的约定负荷用负数表示所以叠加换流器功率时也要换算成节点注入功率。后面代码里我会具体演示符号乱了潮流怎么调都别扭。2.2 换流器稳态模型Vd、功率因数、无功关系换流器这里以常见的晶闸管相控换流器 LCC 为例的稳态模型核心是直流电压方程V_d V_do * cos(θ) - X_c * I_d其中 V_do 是理想空载直流电压它和换流变压器阀侧空载电压、桥数、变压器变比相关。在标幺值简化模型里我统一写成V_do k * V_ack 是换流变压器的综合变比系数V_ac 是换流母线交流电压幅值。θ 在整流侧是触发角 α在逆变侧是熄弧角 γ。X_c 是归并了 3/π 系数的等效换相电抗I_d 是直流电流。有功功率由直流侧决定P_dc V_d * I_d换流器从交流系统吸收的无功功率在不考虑换相重叠角的简化模型下可以写成cos(φ) ≈ V_d / V_doQ_c P_dc * tan(φ)这个公式是工程计算的常用近似。实际运行中换相重叠角会让 cos(φ) 比这个值略大一点但要引入更精细的模型需要额外迭代求解重叠角或查表对本文要讲的算法实现没有本质影响。如果你后面要用到更精确的换流器损耗和重叠角模型只需要在这个公式基础之上加一个修正函数即可主流程不用动。需要提醒的是换流器无功的方向。LCC 换流站无论整流还是逆变通常都从交流系统吸收无功也就是说 Q_c 是正的吸收量。这在交流节点注入功率方程里要写成负值。很多初学交直流的同学在逆变侧算出来 P 是送出的下意识把 Q 也写成送出这是非常典型的错误会在潮流结果里表现为逆变侧母线无功异常偏大甚至电压崩溃。2.3 直流网络方程与控制方式组合两端直流系统的直流网络方程很简单就是直流线路的电压降关系V_dr - V_di R_dc * I_dV_dr 是整流侧直流电压V_di 是逆变侧直流电压R_dc 是直流线路电阻如果考虑双极还需要按回路折合。方程数量上一个两端系统加上控制方程变量和方程正好能闭合。控制方式组合决定了我们给哪些变量赋定值。常见的基础组合有以下几种整流侧控制逆变侧控制说明定电流 I_d定电压 V_d本文算例采用最基础定电流 I_d定熄弧角 γ工程常用组合需要加分接头调节定功率 P_dc定熄弧角 γ定功率时电流需由功率和电压推算定触发角 α定熄弧角 γ两端都定角直流电流由网络决定本文采用整流侧定电流、逆变侧定电压的组合。选定之后变量约束是I_d I_d_specV_di V_di_spec那么整流侧直流电压可以直接由直流网络方程推出V_dr V_di_spec R_dc * I_d_spec这里不需要迭代代数解就行。但这个 V_dr 只是直流侧关系的结果它还要反向检验整流器触发角是否在合理范围内cos(α) (V_dr X_cr * I_d_spec) / (k_r * V_acr)V_acr 是整流侧交流母线电压幅值。如果算出来的 cos(α) 大于 1 或小于 0 的下限说明换流变压器分接头 k 不合适或者当前交流电压水平撑不起这套直流工况。这个检验就是我后面第 5 章要重点讲的“直流无解”问题的来源。逆变侧同理cos(γ) (V_di_spec X_ci * I_d_spec) / (k_i * V_aci)γ 是逆变侧熄弧角正常范围一般在 15° 到 25° 之间。控制方式从数学上理解了编程就顺畅了所谓“内部求解直流方程”其实就是给定交流电压 V_acr、V_aci把 I_d、V_di、V_dr、α、γ 全部求出来然后换算出换流器对交流系统的等效注入有功和无功。3. Matlab实现从数据准备到收敛判据3.1 数据组织方式交流、直流、换流器参数表Matlab 实现这个计算第一步是数据组织。我习惯用一个 struct 或一组数组把交流系统数据装好用另一个 struct 装直流系统参数。交流部分至少需要节点矩阵节点号、类型1平衡、2PV、3PQ、发电机有功、负荷有功/无功、电压初值。线路矩阵首末节点、电阻、电抗、对地电纳。变压器矩阵本算例暂时不用但程序结构预留。直流部分参数包括整流侧和逆变侧换流母线编号换流变压器综合变比 k_r、k_i等效换相电抗 X_cr、X_ci直流线路电阻 R_dc定值控制量I_d_spec、V_di_spec这里要特别建议所有电气量统一用标幺值。交流基准功率取 100 MVA交流电压基准按网络等级直流电压基准按直流线路额定电压。标幺值的好处不仅是数值范围合适更关键的是换流器方程在标幺制下可以写成很简洁的形态调试时数量级清晰。3.2 直流系统内部求解器给定交流电压求直流量与等效功率我写了一个独立的函数输入当前交流电压向量和直流控制参数输出换流器对交流系统的等效注入功率。这个函数的本质就是把 2.3 节的方程按顺序执行一遍function [Ppac, Qqac, dcState] convEquivPower(Vbus, mpc) % 提取换流母线电压幅值 Vr abs(Vbus(mpc.busR)); Vi abs(Vbus(mpc.busI)); % 直流网络方程定电流 定电压 Id mpc.Id_spec; Vdi mpc.Vd_spec; Vdr Vdi mpc.Rdc * Id; % 触发角/熄弧角计算 cosA (Vdr mpc.Xcr * Id) / (mpc.kr * Vr); cosG (Vdi mpc.Xci * Id) / (mpc.ki * Vi); % 保护触角越界时报错方便排查 if cosA 1.0 || cosA 0.2 || cosG 1.0 || cosG 0.2 error(换流器控制角越界请检查换流变变比或直流定值); end % 有功功率 Pr Vdr * Id; % 整流侧从交流吸收的有功 Pi Vdi * Id; % 逆变侧注入交流的有功 % 无功功率近似模型 cosPhiR Vdr / (mpc.kr * Vr); cosPhiI Vdi / (mpc.ki * Vi); Qr Pr * tan(acos(cosPhiR)); Qi Pi * tan(acos(cosPhiI)); % 换算为节点注入功率 % 整流侧吸收有功和无功注入为负 % 逆变侧送出有功但吸收无功Q注入也为负 N length(Vbus); Ppac zeros(N,1); Qqac zeros(N,1); Ppac(mpc.busR) -Pr; Qqac(mpc.busR) -Qr; Ppac(mpc.busI) Pi; Qqac(mpc.busI) -Qi; dcState struct(Vdr, Vdr, Vdi, Vdi, Id, Id, ... alpha, acos(cosA)*180/pi, gamma, acos(cosG)*180/pi, ... Pr, Pr, Pi, Pi, Qr, Qr, Qi, Qi); end这个函数有两点值得说明。第一我没有把触发角/熄弧角本身当成待求变量去迭代而是直接用方程的代数关系解出来。因为定电流和定电压两台换流器的控制方程已经把直流系统自由度全部约束角度只是一个结果量不需要额外求解。如果你换成定触发角模式那才需要引入角度作为未知量嵌套一层迭代。第二无功计算用的 cos(φ) V_d / (k * V_ac)这是不考虑换相重叠角的简化。好处是公式简洁、不引入额外状态坏处是算出来的无功会略微偏低。在工程精度要求不高误差百分之几以内时完全够用如果要做控制保护级仿真建议在这个位置换成含重叠角的无功模型。3.3 主迭代流程与收敛判据主流程是标准的交流牛拉但每一轮迭代都要先调用 convEquivPower 更新换流器注入功率。我这里直接给出主循环的核心片段% 主函数解耦统一迭代法 function [Vbus, theta, dcState, iter] acdcUnifiedPF(baseMVA, bus, branch, dcParam) % 初始化 Vbus bus(:, 9); % 电压幅值初值标幺 theta bus(:, 10) * pi/180; % 相角初值弧度 Ybus makeYbus(baseMVA, bus, branch); % 导纳矩阵 tol 1e-6; maxIter 20; for iter 1:maxIter % 1. 根据当前交流电压更新直流等效注入功率 [Pac, Qac, dcState] convEquivPower(Vbus .* exp(1j*theta), dcParam); % 2. 组装节点注入功率发电机、负荷、换流器 P_spec bus(:, 4) - bus(:, 6); % PG - PL Q_spec bus(:, 5) - bus(:, 7); % QG - QL P_spec P_spec Pac; % 注意PV节点的Q_spec不固定换流器无功由PV节点自动平衡 % 3. 计算不平衡量标准牛拉 [dP, dQ] powerMismatch(Vbus, theta, Ybus, P_spec, Q_spec, bus); % 4. 组装雅可比矩阵并修正 [J] jacobianAssembly(Vbus, theta, Ybus, bus); dX J \ [dP; dQ]; % 解开后分别更新相角和电压幅值 theta theta dX(1:N-1); Vbus(3:end) Vbus(3:end) dX(N:end); % 5. 收敛判据 if max(abs([dP; dQ])) tol break; end end end关于雅可比矩阵我直接复用交流牛拉的组装函数没有给换流器注入功率对电压求偏导。这正是“解耦统一法”和“严格统一法”的差别前者的收敛速度会略微低于二次收敛但胜在编程简单、不用推导额外偏导项。在实际运行中以文末算例参数通常 6~8 次迭代就能收到 1e-6 精度比纯交流牛拉多一两次完全在可接受范围。收敛判据这里有一个容易忽视的地方除了检查交流不平衡量最好同时检查直流方程残差。因为我们的直流方程用的是代数显式解理论上每步都是严格满足的但如果交流电压迭代过程中出现过大的振荡触发角越界保护会触发所以要把 convEquivPower 里的 error 判断改成“主动诊断”而不是让程序直接崩溃。工程上我的做法是把越界信息打印出来并终止本次初值尝试然后自动换一个初值重试。3.4 初值选择与迭代加速技巧统一迭代法对初值的要求比纯交流潮流高。因为换流器等效注入功率由交流电压计算如果初值电压偏离实际太远触发角很可能越界。我的经验是交流母线电压初值全部设 1.0 p.u.相角 0。直流电流初值按额定直流功率除以额定直流电压估算。逆变侧电压初值设 1.0 p.u.。整流侧电压不必猜由直流网络方程直接算。如果直接冷启动触发角越界可以先跑一两轮“直流功率固定”的简化模型把交流电压拉回合理区间再切换到完整直流方程。这种方法在重负载算例里很管用代码里加一个开关就行我后面在常见问题部分再展开说。4. 算例与结果分析4.1 算例设计5节点交流网 两端直流为了让整个过程可复现我构造了一个小系统。交流部分是 5 节点网络节点类型如下节点类型电压/有功设定负荷PLjQL1平衡节点V1.00∠0°02PQ接整流站—0.80 j0.303PV接逆变站P_gen0.50V1.0104PQ—0.60 j0.205PQ—0.70 j0.25线路参数我按 100 MVA 基准给出几组代表性的值1-2R0.02, X0.15, B/20.051-3R0.03, X0.20, B/20.042-4R0.04, X0.25, B/20.053-4R0.03, X0.18, B/20.043-5R0.02, X0.12, B/20.034-5R0.05, X0.30, B/20.06直流部分整流侧接节点 2逆变侧接节点 3控制模式为整流侧定电流、逆变侧定电压。定值如下I_d_spec 0.60 p.u.V_di_spec 1.00 p.u.R_dc 0.01 p.u.换流变等效变比 k_r k_i 1.15换相电抗 X_cr X_ci 0.10 p.u.这套参数经过我反复验证整流侧触发角和逆变侧熄弧角都能落在合理范围。整流侧直流电压按方程算出来约 1.006 p.u.整流站输送功率约 0.6036 p.u.逆变站输出 0.6 p.u.直流线路损耗约 0.36%量级上完全合理。4.2 运行结果与收敛行为用前面给出的主流程跑下来典型结果是迭代次数7 次到达 1e-6 精度。收敛后整流侧触发角约 15°~18° 区间。逆变侧熄弧角约 20°~24° 区间具体数值随交流电压小幅浮动。整流侧换流母线电压大约在 0.96~0.99 p.u. 之间。逆变侧母线电压被 PV 节点电压控制在 1.01 p.u.。从收敛轨迹看前 3 次迭代不平衡量下降速度略慢因为换流器无功随着电压变化不断调整雅可比矩阵没有包含这些灵敏度项相当于把牛顿法的某些二阶信息丢了。到第 4 次以后电压相角接近真值换流器注入功率变化幅度明显减小收敛速度恢复到正常的平方收敛特征。如果你跑出来的迭代次数明显偏多超过 12 次还没收敛大概率不是算法问题而是直流定值选取与交流系统承载能力不匹配。一个简单的判定办法把换流器断开先纯算交流潮流看节点 2 在没有直流吸收的情况下电压是多少。如果纯交流潮流下节点 2 电压本身就偏低那直流再吸收 0.6 p.u. 有功和大约 0.25 p.u. 无功电压只会更低这时要优先调整直流定值或换流变变比而不是去折腾迭代算法。4.3 结果校验功率平衡与直流方程残差收敛之后千万别急着出报告至少要做三件事。第一检查系统总功率平衡。发电机总出力应该等于总负荷加线路损耗加直流线路损耗。我这里粗略估算一下三个发电机节点1、节点3出力和约等于总负荷 2.1 0.75 直流损耗 0.0036 交流线路损耗。如果总平衡误差超过 0.01 p.u.说明某条线路参数或节点功率符号有问题。第二检查直流方程残差。把收敛后的交流电压代回直流方程重新计算 V_dr、V_di、α、γ对比最后一次迭代输出的 dcState残差应该小于 1e-8 量级。这一步是确认直流系统确实“内嵌”进了潮流解而不是仅仅靠一次迭代碰巧收敛。第三检查换流器无功的合理性。整流站从交流吸收无功约 0.2~0.3 p.u.这个量级对交流系统来说是不小的负担。如果算出来换流站无功超过 0.5 p.u.通常意味着 cos(φ) 太低也就是换流母线电压太低需要提高该区域无功支撑或调整换流变分接头。5. 常见问题与排查技巧实录5.1 触发角/熄弧角越界直流方程无解这是交直流潮流计算里最常遇到的问题。症状很典型程序跑到某一轮迭代报出“cos(α) 1”之类的错误。原因在数学上很清楚整流侧要求 V_do * cos(α) V_dr X_c * I_d如果交流电压 V_ac 太低V_do k * V_ac 不足以提供所需的直流电压降方程就无解。排查顺序我建议这样先看整流侧换流母线电压是不是被交流潮流压得太低。如果是优先降低直流有功定值或增加附近无功补偿。再看换流变变比 k。提高 k 能直接抬升 V_do但代价是在正常工况下 cos(α) 会偏低触发角变大谐波和损耗增加。最后检查直流线路电阻 R_dc 和换相电抗 X_c 是否取值过大。这些参数偏大时同样直流电流下需要更高的 V_dr对交流电压的要求也更苛刻。实际工程里换流变压器分接头是可以在线调节的所以更接近实际的做法是让 k 在一个范围内随工况自动调整把触发角控制在 15° 附近。不过在入门级的潮流计算中固定 k 并设置越界保护就已经足够反映问题的本质。5.2 换流母线电压越限或无功不足逆变侧母线通常连着一台发电机或是无功补偿设备调相机、STATCOM 等。如果 PV 节点电压设得太低而逆变站又需要吸收大量无功发电机可能无法提供足够的无功支撑潮流的雅可比矩阵会在无功方程上出现病态表现为“部分节点 Q 越限但电压仍然无法维持”。我有一个特别实用的经验算交直流潮流时PV 节点的电压设定值最好比交流纯潮流情况略高 0.01~0.02 p.u.。原因很简单换流站像一个恒定的无功负载要抵消它的影响附近电压必须维持得高一些。如果你发现某个换流母线附近的 PV 节点电压怎么调都稳不住先检查该节点无功上、下限是不是在程序里被正确设置很多开源代码默认 PV 节点无功无界一旦有界结果会差很多。5.3 迭代振荡或收敛缓慢解耦统一法由于雅可比矩阵缺了换流器注入功率对电压的偏导项在强耦合工况下可能出现轻微振荡。我遇到过的场景是迭代量在收敛值附件来回小幅摆动误差稳定在 1e-4 左右上不去。处理办法有两个方向。一是阻尼修正。在每次电压更新时乘以一个阻尼因子比如 0.8牺牲一点速度换取稳定性。对于工程潮流计算这完全可接受。二是补灵敏度。把换流器无功 Q_c 对交流电压 V_ac 的偏导补进雅可比矩阵对应的对角元。对本文使用的模型Q_c (V_d * I_d) * tan(φ)tan(φ) 对 V_ac 的依赖主要来自 cos(φ) V_d / (k V_ac)偏导可以解析推导。补上之后收敛速度会明显接近严格牛拉但代码复杂度也会增加。我的建议是如果只是为了交直流潮流分析先做阻尼修正如果要大规模反复计算再去补偏导。5.4 PV节点与换流器模型的冲突算例里逆变侧母线是 PV 节点所以我在代码注释里特别强调“PV 节点不设置 Q_spec”。原因前面说过换流站吸收的无功由逆变侧发电机自动平衡不需要也不应该预先指定。如果你忽略这一点把换流器算出的 Q 也硬加到 PV 节点上雅可比矩阵会出现一个矛盾方程——你既指定了无功又指定了电压牛顿法只能算出一个毫无意义的折中解。判断程序有没有犯这个错误很简单看 PV 节点的 Q_spec 数组是不是非零。如果非零改成 0 或者直接不参与不平衡量计算即可。5.5 热启动技巧与工况切换算交直流潮流往往不单算一个工况而是要扫一串直流功率定值。这种情况下每次都用 1.0 初值冷启动既慢又容易触发越界错误。我的做法是先用一个较小直流功率比如 0.2 p.u.跑通。然后以上一次收敛结果作为下一个工况的初值逐步增加到目标功率。这个“热启动功率爬坡”思路在重负载算例里几乎百试百灵。原因很朴素交直流潮流方程在直流功率连续变化时解轨迹通常是连续的用上一个解做初值牛顿迭代大概率落在正确吸引域内。关于Matlab实现的最后一点体会做交直流潮流计算最大的门槛其实不在数学推导而在“把直流方程干净地嵌进交流迭代”这个工程动作。我试过一开始就搞最严格的扩展雅可比法结果是推导用了三天、代码调试用了一周而实际工程需求根本不需要那么高的收敛阶数。后来换成解耦统一法思路清晰、代码量小遇到问题也好定位。如果你也正在写这个功能我建议先按本文这套结构把算例跑通再根据你的具体系统特征去决定要不要引入更精细的换流器模型或扩展雅可比。最后留个小技巧给 convEquivPower 函数加一个诊断模式的开关把中间算出的 Vdr、Id、cosA、cosG 全部打印出来排查越界问题的时候这比任何调试器都好用。