1. 从物理意义到方程组为什么交直流潮流必须用统一迭代法交直流潮流计算光听名字就知道比纯交流潮流麻烦。纯交流系统里节点只有电压幅值、相角两类状态量求解核心是PQ分解或牛顿法但一旦加入换流站、直流线路问题立刻变脏电力电子器件引入了新的控制变量直流侧的电压、电流和功率与交流侧的电压、相角相互耦合而且换流站的运行模式定功率、定电压、定熄弧角等还会改变方程的结构。如果你想“先把交流算完再把直流结果代进去”迭代过程中大概率会来回震荡甚至直接发散。我最早接触这个课题时用的就是交替求解法交流系统算一轮把换流站交流母线功率提取出来折算成直流侧的注入量再去解直流网络方程算完再反馈回交流侧。理论上这个思路很直观但实际跑起来问题一大堆。最主要的矛盾是换流站两侧的方程本来就是联立的交直流接口处的功率、电压、控制量必须同时满足两侧关系你硬要拆开算就相当于把一个二次耦合问题降级成了两个循环嵌套的一阶逼近收敛性完全看初值脸色。统一迭代法的思路恰好相反——把换流站方程、直流网络方程、交流网络方程全部塞进同一个牛顿-拉夫逊迭代框架里联立求解。代价是雅可比矩阵规模变大但只要初始值给得靠谱收敛速度和稳定性远优于交替法。另一个让我当时忽略、后来吃了亏的细节是统一迭代法并不是简单地把直流节点“附加”到交流节点矩阵里就完事它要求你重新思考每个节点的功率平衡方程。换流站交流母线不再是普通PQ节点它的注入功率由直流侧传输功率和换流器损耗共同决定而这个关系本身就是一个隐式方程。你必须在牛顿迭代的每一步同时更新交流状态量和直流状态量否则“统一”就名不副实。说到底选统一迭代法不是因为它新潮而是因为从物理本质上讲交流系统和直流系统之间不存在时间尺度上的解耦它们的稳态解必须同时满足。下面我按实际建模顺序把这套方法从方程搭建到Matlab实现整个拆开讲。2. 模型搭建第一部分换流站稳态方程与直流网络方程2.1 换流站的基本方程交直流潮流里最核心的元件就是电压源换流器VSC或者电网换相换流器LCC。目前学术界和工程界讨论最多的其实是VSC-HVDC因为它能独立控制有功和无功而且不需要额外的换相电压支撑。VSC换流站的稳态模型通常用以下三个方程描述交流侧与直流侧的功率平衡关系 (P_{ac}P_{dc}P_{loss})换流器输出电压与直流电压的关系 (U_c \frac{\sqrt{3}}{2\sqrt{2}} M U_{dc}) 调制比为M换流变压器与电抗器上的电压降落关系 (U_c U_s - (R_c jX_c)I_c)这三个方程是基础但真正用的时候需要根据控制模式作出调整。比如定有功功率控制时(P_{dc})是已知量定直流电压控制时(U_{dc})是已知量定无功功率控制时交流侧无功注入是已知量。每种控制模式会改变雅可比矩阵中对应行的结构这也是统一迭代法实现时最容易出错的地方。2.2 直流网络方程直流网络相对简单没有频率、没有相角只有节点电压和注入电流的关系。用节点导纳矩阵 (Y_{dc}) 描述 (I_{dc}Y_{dc}U_{dc})每个直流节点的注入功率 (P_{dc}U_{dc}I_{dc})。这些方程在统一迭代法里会作为独立的功率失配方程进入牛顿迭代。注意直流侧没有无功功率的概念所以每个直流节点只提供一个功率方程。这意味着如果你把整个系统放在同一个牛顿法框架里每个直流节点贡献一个状态量 (U_{dc}) 和一个失配方程而换流站则根据控制模式贡献不同数量的状态量和方程。2.3 为什么要引入归一化和标幺值我在第一次写代码时直接用有名值结果交流侧电压是110kV量级直流侧电压是±200kV量级雅可比矩阵的条件数大得吓人迭代几步就出现数值溢出。后来老老实实全部转成标幺值问题立刻缓解。潮流计算里标幺值不只是一个“单位换算”问题它直接决定了牛顿法中海森矩阵的数值稳定性。建议交流侧以100MVA为基准直流侧同样以100MVA为基准电压基准取各自的额定电压。这样处理后所有变量都在0.8到1.2这个量级附近牛顿法的收敛行为会稳定得多。3. 模型搭建第二部分统一迭代法的失配方程与雅可比矩阵构造3.1 失配方程组的整体结构统一迭代法的本质就是构造一个整体的非线性方程组 (F(x)0)然后用牛顿-拉夫逊法迭代求解。这个方程组的未知量分为三块交流节点状态量交流母线电压幅值 (V) 和相角 (\theta)直流节点状态量直流母线电压 (U_{dc})换流站内部状态量换流器输出电压的相角 (\delta) 和幅值 (U_c)或者等效为调制比 (M) 和移相角对应地失配方程也分三块交流节点的有功、无功功率失配方程直流节点的有功功率失配方程换流站交直流接口处的功率平衡方程和控制方程这里有个容易被忽视的细节换流站内部状态量不是凭空多出来的。VSC换流站的交流侧电压 (U_c) 的幅值和相角是可调的它必须满足“换流变压器和电抗器上的电压降落”方程同时还要满足控制目标方程。这些方程和状态量必须在雅可比矩阵里一一对应否则矩阵就是奇异的。3.2 雅可比矩阵的分块结构统一迭代法的雅可比矩阵可以写成如下分块形式[ J \begin{bmatrix} J_{AC-AC} J_{AC-DC} \ J_{DC-AC} J_{DC-DC} \end{bmatrix} ]其中 (J_{AC-AC}) 是传统交流潮流的雅可比矩阵加上换流站功率对交流状态量的偏导数(J_{AC-DC}) 是交流失配方程对直流状态量的偏导数(J_{DC-AC}) 是直流失配方程对交流状态量的偏导数(J_{DC-DC}) 是直流失配方程对直流状态量的偏导数。很多人第一次写代码时会偷懒把交流雅可比矩阵当成常量只在迭代过程中更新右端项。这个做法在纯交流潮流里可能勉强能用在统一迭代法里绝对不行。因为换流站功率与直流电压直接相关而直流电压本身是迭代变量你必须把交流失配方程对直流电压的偏导数也计算出来。我踩过的坑就在这里一开始没算交叉偏导数矩阵接近奇异迭代始终在10的-2次方精度附近打转怎么都下不去。3.3 偏导数的推导示例以VSC换流站为例设交流母线电压为 (U_s \angle \theta_s)换流器输出电压为 (U_c \angle \theta_c)两者之间的等效阻抗为 (Z R jX)。从交流母线注入换流站的有功功率为[ P_s \frac{U_s U_c}{Z} \sin(\theta_s - \theta_c \alpha) - \frac{U_s^2}{Z} \sin \alpha ]这里 (\alpha \arctan(R/X))。这个功率既与 (U_s)、(\theta_s) 有关也与 (U_c)、(\theta_c) 有关。其中 (U_c) 又与直流电压 (U_{dc}) 满足 (U_c \frac{\sqrt{3}}{2\sqrt{2}} M U_{dc})。因此失配方程对 (U_{dc}) 的偏导数必须通过链式法则求得。这部分推导虽然繁琐但谁绕过去谁就会在矩阵奇异和收敛失败上栽跟头。3.4 控制模式的切换问题换流站是定有功功率还是定直流电压会直接影响雅可比矩阵里对应行的结构。定功率模式下换流站的有功功率失配方程是硬约束而直流电压是自由变量定电压模式下直流电压是已知量不再作为未知量而换流站的有功功率失配方程要替换成电压偏差方程。这个切换在编程时最容易产生索引错位。我的习惯是把换流站的控制模式定义成一个枚举变量比如1代表定有功功率2代表定直流电压3代表定无功功率在构造未知量索引和雅可比矩阵计算时都通过这个枚举变量来判断。这样虽然代码里分支判断多了一些但至少结构清晰调试时也容易定位问题。4. Matlab代码实现从数据输入到统一迭代求解4.1 数据结构设计写交直流潮流程序第一件事不是写牛顿迭代而是先设计好数据结构。我建议用Matlab的struct来组织整个系统的数据不要用零散的全局变量。下面是我常用的数据结构示例% 交流系统数据 system.ac.bus [ 1 1.06 0.0 0 0 0 0 1; 2 1.00 0.0 0 0 50 30 2; 3 1.00 0.0 0 0 60 40 2; ]; % bus数据格式: [节点编号, 电压幅值初值, 相角初值, 发电机有功, 发电机无功, 负荷有功, 负荷无功, 节点类型] % 节点类型: 1平衡节点, 2PQ节点, 3PV节点 % 换流站数据 system.vsc struct(... bus_ac, 3, ... % 交流侧接入节点 bus_dc, 4, ... % 直流侧节点编号 mode_p, 1, ... % 有功控制模式: 1定功率, 2定直流电压 p_ref, 20, ... % 有功功率参考值(MW) u_dc_ref, 1.0, ... % 直流电压参考值(标幺值) x_t, 0.15, ... % 换流变压器电抗(标幺值) r_c, 0.005, ... % 换流器等效电阻(标幺值) mode_q, 1); % 无功控制模式: 1定无功, 2定交流电压 % 直流网络数据 system.dc.bus [ 4 1.0 0; 5 1.0 0 ]; % 直流节点数据格式: [节点编号, 电压初值, 节点类型(0无控制)] system.dc.branch [ 4 5 0.01 0.1; ]; % 直流支路数据格式: [起始节点, 终止节点, 电阻, 电抗(直流无电抗, 这里填0)]4.2 初始化与初值设定统一迭代法对初值的要求比纯交流潮流要高。我的经验是交流母线电压初值取1.0∠0°直流电压初值取1.0换流器输出电压初值取交流母线电压的0.95倍左右相角初值取比交流母线相角滞后10°左右。这个初值策略在绝大多数VSC-HVDC算例里都能保证迭代收敛。如果初值给得太离谱比如直流电压初值给0.5雅可比矩阵可能会在早期迭代中出现奇异。% 初始化状态向量x % x的顺序: [交流节点相角(除平衡节点), 交流节点电压幅值(PQ节点), 直流节点电压, 换流站附加状态量] n_ac size(system.ac.bus, 1); n_pv_pq_ac ...; % 根据节点类型统计 n_dc size(system.dc.bus, 1); n_vsc length(system.vsc); x zeros(n_pv_pq_ac n_dc n_vsc * 2, 1); % 按顺序填充初值 % 交流相角初值 x(1:n_ac-1) 0.0; % 交流电压幅值初值(PQ节点) if n_pq 0 x(n_ac:n_acn_pq-1) 1.0; end % 直流电压初值 x(idx_dc_start:idx_dc_startn_dc-1) 1.0; % 换流站附加状态量初值 for k 1:n_vsc x(idx_vsc_start (k-1)*2) 0.95; % 换流器输出电压幅值 x(idx_vsc_start (k-1)*2 1) -0.1; % 换流器输出电压相角 end4.3 牛-拉夫逊迭代主循环统一迭代法的核心就是一个标准的牛顿法迭代循环计算失配量F计算雅可比矩阵J求解修正方程J·Δx -F更新x重复直到收敛。下面给出主循环框架max_iter 30; tol 1e-8; for iter 1:max_iter % 1. 计算失配方程 [F, ~] compute_mismatch(system, x); % 2. 计算雅可比矩阵 J compute_jacobian(system, x); % 3. 求解修正方程 dx -J \ F; % 4. 更新状态量 x x dx; % 5. 检查收敛 if norm(F, inf) tol fprintf(迭代收敛于第%d次\n, iter); break; end if iter max_iter error(迭代未收敛); end end这里我强烈建议用norm(F, inf)而不是norm(F, 2)来判断收敛。最大范数能直接反映最差节点的功率失配量二范数会把误差“平均”掉可能掩盖某个节点失配很大的问题。工程上通常要求功率失配量小于 (10^{-6}) 标幺值如果只是学习用途(10^{-6}) 到 (10^{-8}) 都可以接受。4.4 失配方程计算的详细实现失配方程是整个程序的核心也是容易出错的地方。下面给出compute_mismatch函数的简化实现思路function F compute_mismatch(system, x) % 从x中提取交流状态量 % 计算交流节点注入功率P_calc, Q_calc % 计算直流网络节点注入功率P_dc_calc % 计算换流站接口功率 % 交流节点功率失配 F_ac_p P_spec - P_calc; % 有功失配 F_ac_q Q_spec - Q_calc; % 无功失配 % 直流节点功率失配 F_dc_p P_dc_spec - P_dc_calc; % 换流站方程失配 % 这里需要根据控制模式组合方程 F_vsc compute_vsc_mismatch(system, x); % 组装 F [F_ac_p; F_ac_q; F_dc_p; F_vsc]; end这里有个关键点交流节点功率 (P_{calc})、(Q_{calc}) 的计算要考虑换流站的注入功率。换流站接入的交流母线其注入功率不再只是发电机和负荷的净值还要叠加换流站从交流侧吸收的功率。我在实现时是把换流站当成一个“可变功率注入源”每次迭代用当前状态量计算它的注入功率然后叠加到对应节点的功率失配里。这个处理方式比直接把换流站当成特殊节点更简洁而且不需要修改原始的交流潮流计算函数。4.5 换流站方程的详细实现VSC换流站的失配方程通常包括以下三类有功控制方程(P_{dc} - P_{ref} 0) 定功率模式 或 (U_{dc} - U_{dc,ref} 0) 定直流电压模式无功控制方程(Q_{ac} - Q_{ref} 0)定无功模式 或 (V_{ac} - V_{ac,ref} 0)定交流电压模式物理约束方程换流器交流侧电压、变压器阻抗压降、直流电压三者的关系第三类方程是统一迭代法特有的。它把换流器内部的状态量与交直流两侧的状态量联系起来是保证方程组完备性的关键。如果不写这个方程换流器输出电压幅值和相角就成了无法确定的自由变量雅可比矩阵必然奇异。function F_vsc compute_vsc_mismatch(system, x) F_vsc []; for k 1:length(system.vsc) vsc system.vsc(k); % 提取交流母线状态量 V_ac ...; theta_ac ...; % 提取直流节点电压 U_dc ...; % 提取换流器状态量 U_c ...; theta_c ...; % 计算交流注入功率 [P_ac, Q_ac] calc_vsc_ac_power(V_ac, theta_ac, U_c, theta_c, vsc); % 计算直流功率 P_dc U_dc * ...; % 由直流网络方程给出 % 有功平衡方程 F_vsc [F_vsc; P_ac - P_dc]; % 控制方程 if vsc.mode_p 1 F_vsc [F_vsc; P_dc - vsc.p_ref]; elseif vsc.mode_p 2 F_vsc [F_vsc; U_dc - vsc.u_dc_ref]; end % 无功控制方程 if vsc.mode_q 1 F_vsc [F_vsc; Q_ac - vsc.q_ref]; end end end4.6 雅可比矩阵的数值计算很多教材会花大量篇幅推导解析雅可比矩阵但实际用Matlab开发时我推荐先用数值差分验证解析结果或者直接采用数值雅可比矩阵来完成第一版功能。所谓数值雅可比就是对每个状态量加一个小扰动δ重新计算失配量然后用差分近似偏导数。function J compute_jacobian_numerical(system, x) n length(x); F0 compute_mismatch(system, x); J zeros(length(F0), n); delta 1e-7; for j 1:n x_pert x; x_pert(j) x_pert(j) delta; F_pert compute_mismatch(system, x_pert); J(:, j) (F_pert - F0) / delta; end end这种做法的优点是代码简洁、不易出错缺点是计算量大。对于小规模系统几十个节点以内完全够用但如果是大型系统建议在数值雅可比验证通过后再改写成解析雅可比。我个人的开发路径是先用数值雅可比跑通整体流程然后逐步把关键偏导数解析化每替换一块就用数值结果对比验证误差在1e-6以内就认为正确。这个习惯帮我省去了大量调试时间。5. 调试心得常见问题与排查技巧实录5.1 问题1迭代发散残差越来越大最常见的发散原因有两个初值不合适或者雅可比矩阵奇异。排查方法很简单在每次迭代后打印雅可比矩阵的条件数如果条件数在不断增大说明方程组本身有问题可能是指标缺失或者方程冗余。如果条件数一直正常但残差反复震荡大概率是初值距离真实解太远。我的处理顺序是先用扁平化系统把换流站替换成恒定功率注入跑一遍纯交流潮流验证交流网络本身没有建模错误。然后在纯交流潮流结果的基础上逐渐引入直流网络把换流站功率从恒定值改成迭代值。这样分步调试比一上来就上完整模型高效得多。5.2 问题2雅可比矩阵奇异雅可比矩阵奇异通常意味着未知量个数与方程个数不匹配或者存在冗余方程。检查方法计算秩或者行列式再看看哪些方程是线性相关的。最常见的情况是VSC换流站的无功控制方程与交流母线的电压幅值状态量没有形成正确的偏导关系。比如定交流电压模式下控制方程是 (V_{ac} - V_{ac,ref}0)但交流母线的电压幅值本身可能不是状态量如果该母线是平衡节点这时就必须引入换流站输出电压作为控制量否则方程无解。5.3 问题3收敛精度不够始终停在1e-4这种情况往往是失配方程中某个变量没参与迭代或者迭代过程中某处用了旧值没有更新。我调试时会在每次迭代后打印关键节点的电压幅值和功率观察它们是否在按照合理的趋势变化。如果发现某个变量一直不变检查它是否被错误地从状态向量中剔除或者雅可比矩阵中对应行全为零。5.4 问题4算例结果与文献对不上交直流潮流的结果验证我通常用两个方式一是把换流站的损耗设为零把直流网络电阻设为零这样直流侧总有功功率应该等于交流侧注入功率之差可以验证功率平衡二是把直流电压固定为1.0换流站定功率模式交流潮流结果应该与把换流站等效成恒功率负荷的结果一致。这两个“退化测试”能快速定位是模型错误还是数值错误。6. 一个完整算例三节点交流系统接两端直流网络为了让你能直接照猫画虎我给出一个最简单的算例。交流系统为三节点系统节点1是平衡节点节点2带有负荷节点3通过VSC换流站接入直流网络。直流网络是两端结构节点4通过直流线路接到节点5节点5连接另一个VSC换流站该换流站定直流电压控制支撑直流网络电压。系统参数如下交流基准容量100MVA交流节点1平衡节点电压1.06∠0°交流节点2PQ节点负荷50MWj30Mvar交流节点3PQ节点负荷60MWj40Mvar同时接入VSC换流站VSC1接入交流节点3定有功功率控制P_ref20MW定无功功率控制Q_ref0MvarVSC2接入交流节点5直流侧定直流电压控制U_dc_ref1.0直流线路节点4到节点5电阻R0.01标幺值运行统一迭代法正确结果应该是交流节点3的电压略低于1.0因为换流站从该节点吸取了20MW有功直流节点4的电压会略低于1.0因为有线路压降VSC1的换流器输出电压幅值在0.93到0.97之间相角比交流母线滞后几度。如果这些量级出现异常说明代码里有模型错误。这个算例虽然简单五脏俱全包含了交流网络、直流网络、两种控制模式定功率和定直流电压、换流站接口方程。跑通这个算例之后扩展到多端直流、多换流站系统就只是增加节点和分支的问题了。7. 从会写到会用几个提高效率的实用技巧7.1 把失配方程函数写成向量化形式Matlab的循环效率远不如向量化计算。如果你的系统规模不大循环无所谓但一旦节点数量上百建议把交流导纳矩阵计算、功率失配计算都写成矩阵运算形式。比如用稀疏矩阵存储导纳矩阵用矢量化公式一次性计算所有节点的注入功率可以提速一个数量级以上。7.2 使用匿名函数简化控制模式切换控制模式切换可以用匿名函数或者函数句柄数组来组织。比如定义两个函数句柄calc_p_error (P_dc, U_dc, ref) ...然后在主循环里根据模式调用不同的句柄。代码会更清爽调试时也更容易单步跟踪。7.3 用disp和fprintf记录迭代轨迹统一迭代法的调试比纯交流潮流更依赖迭代轨迹信息。我习惯在每次迭代后打印迭代次数、最大失配量、直流节点电压、换流站输出电压幅值。这样一旦发散从打印数据里能快速看出是交流侧先发散还是直流侧先发散。很多问题看一眼轨迹曲线就能定位不用逐行查代码。7.4 学会用退化测试验证代码退化测试是代码正确性的试金石。最简单的一个把直流线路电阻设为零两个换流站都设定有功功率控制其他什么都不变。这时直流网络相当于一个无损耗功率传输通道两侧换流站的有功功率应该完全一致交流系统的功率平衡也应该是精确的。如果算出来两侧功率有微小偏差说明直流方程的符号或者方向定义有问题。8. 一些基于经验的小忠告交直流潮流计算的难点不在数学而在建模。很多人一上来就盯着雅可比矩阵的公式推导其实真正花时间的是理解每个控制模式对应的方程结构。控制模式一变状态量的索引要变失配方程要变雅可比矩阵的稀疏结构也要变。如果一开始就设计成“控制模式可配置”的程序结构后面扩展多端直流、混合直流都会省力很多。另一个容易被低估的问题是单位制的统一。我见过有人在Matlab里用有名值又混入标幺值参数结果迭代出来一些诡异的数值查了一晚上才发现是单位混乱。建议整个程序从头到尾统一用标幺值只在输入输出边界做转换。直流电压基准和交流电压基准分别取各自的额定电压功率基准统一取系统基准容量这样最稳妥。最后再分享一个小技巧调试时不要同时验证多个新功能。我每次只改一个点比如先只验证“换流站定有功功率模式”跑通了再加“定直流电压模式”再加“无功控制模式”。每加一个模式都用退化测试验证结果没有回归。这个习惯看上去慢但实际总耗时最少。交直流潮流是一个经典但很容易写错的地方希望这整套流程能帮你少走弯路。
交直流潮流统一迭代法的建模与Matlab实现
1. 从物理意义到方程组为什么交直流潮流必须用统一迭代法交直流潮流计算光听名字就知道比纯交流潮流麻烦。纯交流系统里节点只有电压幅值、相角两类状态量求解核心是PQ分解或牛顿法但一旦加入换流站、直流线路问题立刻变脏电力电子器件引入了新的控制变量直流侧的电压、电流和功率与交流侧的电压、相角相互耦合而且换流站的运行模式定功率、定电压、定熄弧角等还会改变方程的结构。如果你想“先把交流算完再把直流结果代进去”迭代过程中大概率会来回震荡甚至直接发散。我最早接触这个课题时用的就是交替求解法交流系统算一轮把换流站交流母线功率提取出来折算成直流侧的注入量再去解直流网络方程算完再反馈回交流侧。理论上这个思路很直观但实际跑起来问题一大堆。最主要的矛盾是换流站两侧的方程本来就是联立的交直流接口处的功率、电压、控制量必须同时满足两侧关系你硬要拆开算就相当于把一个二次耦合问题降级成了两个循环嵌套的一阶逼近收敛性完全看初值脸色。统一迭代法的思路恰好相反——把换流站方程、直流网络方程、交流网络方程全部塞进同一个牛顿-拉夫逊迭代框架里联立求解。代价是雅可比矩阵规模变大但只要初始值给得靠谱收敛速度和稳定性远优于交替法。另一个让我当时忽略、后来吃了亏的细节是统一迭代法并不是简单地把直流节点“附加”到交流节点矩阵里就完事它要求你重新思考每个节点的功率平衡方程。换流站交流母线不再是普通PQ节点它的注入功率由直流侧传输功率和换流器损耗共同决定而这个关系本身就是一个隐式方程。你必须在牛顿迭代的每一步同时更新交流状态量和直流状态量否则“统一”就名不副实。说到底选统一迭代法不是因为它新潮而是因为从物理本质上讲交流系统和直流系统之间不存在时间尺度上的解耦它们的稳态解必须同时满足。下面我按实际建模顺序把这套方法从方程搭建到Matlab实现整个拆开讲。2. 模型搭建第一部分换流站稳态方程与直流网络方程2.1 换流站的基本方程交直流潮流里最核心的元件就是电压源换流器VSC或者电网换相换流器LCC。目前学术界和工程界讨论最多的其实是VSC-HVDC因为它能独立控制有功和无功而且不需要额外的换相电压支撑。VSC换流站的稳态模型通常用以下三个方程描述交流侧与直流侧的功率平衡关系 (P_{ac}P_{dc}P_{loss})换流器输出电压与直流电压的关系 (U_c \frac{\sqrt{3}}{2\sqrt{2}} M U_{dc}) 调制比为M换流变压器与电抗器上的电压降落关系 (U_c U_s - (R_c jX_c)I_c)这三个方程是基础但真正用的时候需要根据控制模式作出调整。比如定有功功率控制时(P_{dc})是已知量定直流电压控制时(U_{dc})是已知量定无功功率控制时交流侧无功注入是已知量。每种控制模式会改变雅可比矩阵中对应行的结构这也是统一迭代法实现时最容易出错的地方。2.2 直流网络方程直流网络相对简单没有频率、没有相角只有节点电压和注入电流的关系。用节点导纳矩阵 (Y_{dc}) 描述 (I_{dc}Y_{dc}U_{dc})每个直流节点的注入功率 (P_{dc}U_{dc}I_{dc})。这些方程在统一迭代法里会作为独立的功率失配方程进入牛顿迭代。注意直流侧没有无功功率的概念所以每个直流节点只提供一个功率方程。这意味着如果你把整个系统放在同一个牛顿法框架里每个直流节点贡献一个状态量 (U_{dc}) 和一个失配方程而换流站则根据控制模式贡献不同数量的状态量和方程。2.3 为什么要引入归一化和标幺值我在第一次写代码时直接用有名值结果交流侧电压是110kV量级直流侧电压是±200kV量级雅可比矩阵的条件数大得吓人迭代几步就出现数值溢出。后来老老实实全部转成标幺值问题立刻缓解。潮流计算里标幺值不只是一个“单位换算”问题它直接决定了牛顿法中海森矩阵的数值稳定性。建议交流侧以100MVA为基准直流侧同样以100MVA为基准电压基准取各自的额定电压。这样处理后所有变量都在0.8到1.2这个量级附近牛顿法的收敛行为会稳定得多。3. 模型搭建第二部分统一迭代法的失配方程与雅可比矩阵构造3.1 失配方程组的整体结构统一迭代法的本质就是构造一个整体的非线性方程组 (F(x)0)然后用牛顿-拉夫逊法迭代求解。这个方程组的未知量分为三块交流节点状态量交流母线电压幅值 (V) 和相角 (\theta)直流节点状态量直流母线电压 (U_{dc})换流站内部状态量换流器输出电压的相角 (\delta) 和幅值 (U_c)或者等效为调制比 (M) 和移相角对应地失配方程也分三块交流节点的有功、无功功率失配方程直流节点的有功功率失配方程换流站交直流接口处的功率平衡方程和控制方程这里有个容易被忽视的细节换流站内部状态量不是凭空多出来的。VSC换流站的交流侧电压 (U_c) 的幅值和相角是可调的它必须满足“换流变压器和电抗器上的电压降落”方程同时还要满足控制目标方程。这些方程和状态量必须在雅可比矩阵里一一对应否则矩阵就是奇异的。3.2 雅可比矩阵的分块结构统一迭代法的雅可比矩阵可以写成如下分块形式[ J \begin{bmatrix} J_{AC-AC} J_{AC-DC} \ J_{DC-AC} J_{DC-DC} \end{bmatrix} ]其中 (J_{AC-AC}) 是传统交流潮流的雅可比矩阵加上换流站功率对交流状态量的偏导数(J_{AC-DC}) 是交流失配方程对直流状态量的偏导数(J_{DC-AC}) 是直流失配方程对交流状态量的偏导数(J_{DC-DC}) 是直流失配方程对直流状态量的偏导数。很多人第一次写代码时会偷懒把交流雅可比矩阵当成常量只在迭代过程中更新右端项。这个做法在纯交流潮流里可能勉强能用在统一迭代法里绝对不行。因为换流站功率与直流电压直接相关而直流电压本身是迭代变量你必须把交流失配方程对直流电压的偏导数也计算出来。我踩过的坑就在这里一开始没算交叉偏导数矩阵接近奇异迭代始终在10的-2次方精度附近打转怎么都下不去。3.3 偏导数的推导示例以VSC换流站为例设交流母线电压为 (U_s \angle \theta_s)换流器输出电压为 (U_c \angle \theta_c)两者之间的等效阻抗为 (Z R jX)。从交流母线注入换流站的有功功率为[ P_s \frac{U_s U_c}{Z} \sin(\theta_s - \theta_c \alpha) - \frac{U_s^2}{Z} \sin \alpha ]这里 (\alpha \arctan(R/X))。这个功率既与 (U_s)、(\theta_s) 有关也与 (U_c)、(\theta_c) 有关。其中 (U_c) 又与直流电压 (U_{dc}) 满足 (U_c \frac{\sqrt{3}}{2\sqrt{2}} M U_{dc})。因此失配方程对 (U_{dc}) 的偏导数必须通过链式法则求得。这部分推导虽然繁琐但谁绕过去谁就会在矩阵奇异和收敛失败上栽跟头。3.4 控制模式的切换问题换流站是定有功功率还是定直流电压会直接影响雅可比矩阵里对应行的结构。定功率模式下换流站的有功功率失配方程是硬约束而直流电压是自由变量定电压模式下直流电压是已知量不再作为未知量而换流站的有功功率失配方程要替换成电压偏差方程。这个切换在编程时最容易产生索引错位。我的习惯是把换流站的控制模式定义成一个枚举变量比如1代表定有功功率2代表定直流电压3代表定无功功率在构造未知量索引和雅可比矩阵计算时都通过这个枚举变量来判断。这样虽然代码里分支判断多了一些但至少结构清晰调试时也容易定位问题。4. Matlab代码实现从数据输入到统一迭代求解4.1 数据结构设计写交直流潮流程序第一件事不是写牛顿迭代而是先设计好数据结构。我建议用Matlab的struct来组织整个系统的数据不要用零散的全局变量。下面是我常用的数据结构示例% 交流系统数据 system.ac.bus [ 1 1.06 0.0 0 0 0 0 1; 2 1.00 0.0 0 0 50 30 2; 3 1.00 0.0 0 0 60 40 2; ]; % bus数据格式: [节点编号, 电压幅值初值, 相角初值, 发电机有功, 发电机无功, 负荷有功, 负荷无功, 节点类型] % 节点类型: 1平衡节点, 2PQ节点, 3PV节点 % 换流站数据 system.vsc struct(... bus_ac, 3, ... % 交流侧接入节点 bus_dc, 4, ... % 直流侧节点编号 mode_p, 1, ... % 有功控制模式: 1定功率, 2定直流电压 p_ref, 20, ... % 有功功率参考值(MW) u_dc_ref, 1.0, ... % 直流电压参考值(标幺值) x_t, 0.15, ... % 换流变压器电抗(标幺值) r_c, 0.005, ... % 换流器等效电阻(标幺值) mode_q, 1); % 无功控制模式: 1定无功, 2定交流电压 % 直流网络数据 system.dc.bus [ 4 1.0 0; 5 1.0 0 ]; % 直流节点数据格式: [节点编号, 电压初值, 节点类型(0无控制)] system.dc.branch [ 4 5 0.01 0.1; ]; % 直流支路数据格式: [起始节点, 终止节点, 电阻, 电抗(直流无电抗, 这里填0)]4.2 初始化与初值设定统一迭代法对初值的要求比纯交流潮流要高。我的经验是交流母线电压初值取1.0∠0°直流电压初值取1.0换流器输出电压初值取交流母线电压的0.95倍左右相角初值取比交流母线相角滞后10°左右。这个初值策略在绝大多数VSC-HVDC算例里都能保证迭代收敛。如果初值给得太离谱比如直流电压初值给0.5雅可比矩阵可能会在早期迭代中出现奇异。% 初始化状态向量x % x的顺序: [交流节点相角(除平衡节点), 交流节点电压幅值(PQ节点), 直流节点电压, 换流站附加状态量] n_ac size(system.ac.bus, 1); n_pv_pq_ac ...; % 根据节点类型统计 n_dc size(system.dc.bus, 1); n_vsc length(system.vsc); x zeros(n_pv_pq_ac n_dc n_vsc * 2, 1); % 按顺序填充初值 % 交流相角初值 x(1:n_ac-1) 0.0; % 交流电压幅值初值(PQ节点) if n_pq 0 x(n_ac:n_acn_pq-1) 1.0; end % 直流电压初值 x(idx_dc_start:idx_dc_startn_dc-1) 1.0; % 换流站附加状态量初值 for k 1:n_vsc x(idx_vsc_start (k-1)*2) 0.95; % 换流器输出电压幅值 x(idx_vsc_start (k-1)*2 1) -0.1; % 换流器输出电压相角 end4.3 牛-拉夫逊迭代主循环统一迭代法的核心就是一个标准的牛顿法迭代循环计算失配量F计算雅可比矩阵J求解修正方程J·Δx -F更新x重复直到收敛。下面给出主循环框架max_iter 30; tol 1e-8; for iter 1:max_iter % 1. 计算失配方程 [F, ~] compute_mismatch(system, x); % 2. 计算雅可比矩阵 J compute_jacobian(system, x); % 3. 求解修正方程 dx -J \ F; % 4. 更新状态量 x x dx; % 5. 检查收敛 if norm(F, inf) tol fprintf(迭代收敛于第%d次\n, iter); break; end if iter max_iter error(迭代未收敛); end end这里我强烈建议用norm(F, inf)而不是norm(F, 2)来判断收敛。最大范数能直接反映最差节点的功率失配量二范数会把误差“平均”掉可能掩盖某个节点失配很大的问题。工程上通常要求功率失配量小于 (10^{-6}) 标幺值如果只是学习用途(10^{-6}) 到 (10^{-8}) 都可以接受。4.4 失配方程计算的详细实现失配方程是整个程序的核心也是容易出错的地方。下面给出compute_mismatch函数的简化实现思路function F compute_mismatch(system, x) % 从x中提取交流状态量 % 计算交流节点注入功率P_calc, Q_calc % 计算直流网络节点注入功率P_dc_calc % 计算换流站接口功率 % 交流节点功率失配 F_ac_p P_spec - P_calc; % 有功失配 F_ac_q Q_spec - Q_calc; % 无功失配 % 直流节点功率失配 F_dc_p P_dc_spec - P_dc_calc; % 换流站方程失配 % 这里需要根据控制模式组合方程 F_vsc compute_vsc_mismatch(system, x); % 组装 F [F_ac_p; F_ac_q; F_dc_p; F_vsc]; end这里有个关键点交流节点功率 (P_{calc})、(Q_{calc}) 的计算要考虑换流站的注入功率。换流站接入的交流母线其注入功率不再只是发电机和负荷的净值还要叠加换流站从交流侧吸收的功率。我在实现时是把换流站当成一个“可变功率注入源”每次迭代用当前状态量计算它的注入功率然后叠加到对应节点的功率失配里。这个处理方式比直接把换流站当成特殊节点更简洁而且不需要修改原始的交流潮流计算函数。4.5 换流站方程的详细实现VSC换流站的失配方程通常包括以下三类有功控制方程(P_{dc} - P_{ref} 0) 定功率模式 或 (U_{dc} - U_{dc,ref} 0) 定直流电压模式无功控制方程(Q_{ac} - Q_{ref} 0)定无功模式 或 (V_{ac} - V_{ac,ref} 0)定交流电压模式物理约束方程换流器交流侧电压、变压器阻抗压降、直流电压三者的关系第三类方程是统一迭代法特有的。它把换流器内部的状态量与交直流两侧的状态量联系起来是保证方程组完备性的关键。如果不写这个方程换流器输出电压幅值和相角就成了无法确定的自由变量雅可比矩阵必然奇异。function F_vsc compute_vsc_mismatch(system, x) F_vsc []; for k 1:length(system.vsc) vsc system.vsc(k); % 提取交流母线状态量 V_ac ...; theta_ac ...; % 提取直流节点电压 U_dc ...; % 提取换流器状态量 U_c ...; theta_c ...; % 计算交流注入功率 [P_ac, Q_ac] calc_vsc_ac_power(V_ac, theta_ac, U_c, theta_c, vsc); % 计算直流功率 P_dc U_dc * ...; % 由直流网络方程给出 % 有功平衡方程 F_vsc [F_vsc; P_ac - P_dc]; % 控制方程 if vsc.mode_p 1 F_vsc [F_vsc; P_dc - vsc.p_ref]; elseif vsc.mode_p 2 F_vsc [F_vsc; U_dc - vsc.u_dc_ref]; end % 无功控制方程 if vsc.mode_q 1 F_vsc [F_vsc; Q_ac - vsc.q_ref]; end end end4.6 雅可比矩阵的数值计算很多教材会花大量篇幅推导解析雅可比矩阵但实际用Matlab开发时我推荐先用数值差分验证解析结果或者直接采用数值雅可比矩阵来完成第一版功能。所谓数值雅可比就是对每个状态量加一个小扰动δ重新计算失配量然后用差分近似偏导数。function J compute_jacobian_numerical(system, x) n length(x); F0 compute_mismatch(system, x); J zeros(length(F0), n); delta 1e-7; for j 1:n x_pert x; x_pert(j) x_pert(j) delta; F_pert compute_mismatch(system, x_pert); J(:, j) (F_pert - F0) / delta; end end这种做法的优点是代码简洁、不易出错缺点是计算量大。对于小规模系统几十个节点以内完全够用但如果是大型系统建议在数值雅可比验证通过后再改写成解析雅可比。我个人的开发路径是先用数值雅可比跑通整体流程然后逐步把关键偏导数解析化每替换一块就用数值结果对比验证误差在1e-6以内就认为正确。这个习惯帮我省去了大量调试时间。5. 调试心得常见问题与排查技巧实录5.1 问题1迭代发散残差越来越大最常见的发散原因有两个初值不合适或者雅可比矩阵奇异。排查方法很简单在每次迭代后打印雅可比矩阵的条件数如果条件数在不断增大说明方程组本身有问题可能是指标缺失或者方程冗余。如果条件数一直正常但残差反复震荡大概率是初值距离真实解太远。我的处理顺序是先用扁平化系统把换流站替换成恒定功率注入跑一遍纯交流潮流验证交流网络本身没有建模错误。然后在纯交流潮流结果的基础上逐渐引入直流网络把换流站功率从恒定值改成迭代值。这样分步调试比一上来就上完整模型高效得多。5.2 问题2雅可比矩阵奇异雅可比矩阵奇异通常意味着未知量个数与方程个数不匹配或者存在冗余方程。检查方法计算秩或者行列式再看看哪些方程是线性相关的。最常见的情况是VSC换流站的无功控制方程与交流母线的电压幅值状态量没有形成正确的偏导关系。比如定交流电压模式下控制方程是 (V_{ac} - V_{ac,ref}0)但交流母线的电压幅值本身可能不是状态量如果该母线是平衡节点这时就必须引入换流站输出电压作为控制量否则方程无解。5.3 问题3收敛精度不够始终停在1e-4这种情况往往是失配方程中某个变量没参与迭代或者迭代过程中某处用了旧值没有更新。我调试时会在每次迭代后打印关键节点的电压幅值和功率观察它们是否在按照合理的趋势变化。如果发现某个变量一直不变检查它是否被错误地从状态向量中剔除或者雅可比矩阵中对应行全为零。5.4 问题4算例结果与文献对不上交直流潮流的结果验证我通常用两个方式一是把换流站的损耗设为零把直流网络电阻设为零这样直流侧总有功功率应该等于交流侧注入功率之差可以验证功率平衡二是把直流电压固定为1.0换流站定功率模式交流潮流结果应该与把换流站等效成恒功率负荷的结果一致。这两个“退化测试”能快速定位是模型错误还是数值错误。6. 一个完整算例三节点交流系统接两端直流网络为了让你能直接照猫画虎我给出一个最简单的算例。交流系统为三节点系统节点1是平衡节点节点2带有负荷节点3通过VSC换流站接入直流网络。直流网络是两端结构节点4通过直流线路接到节点5节点5连接另一个VSC换流站该换流站定直流电压控制支撑直流网络电压。系统参数如下交流基准容量100MVA交流节点1平衡节点电压1.06∠0°交流节点2PQ节点负荷50MWj30Mvar交流节点3PQ节点负荷60MWj40Mvar同时接入VSC换流站VSC1接入交流节点3定有功功率控制P_ref20MW定无功功率控制Q_ref0MvarVSC2接入交流节点5直流侧定直流电压控制U_dc_ref1.0直流线路节点4到节点5电阻R0.01标幺值运行统一迭代法正确结果应该是交流节点3的电压略低于1.0因为换流站从该节点吸取了20MW有功直流节点4的电压会略低于1.0因为有线路压降VSC1的换流器输出电压幅值在0.93到0.97之间相角比交流母线滞后几度。如果这些量级出现异常说明代码里有模型错误。这个算例虽然简单五脏俱全包含了交流网络、直流网络、两种控制模式定功率和定直流电压、换流站接口方程。跑通这个算例之后扩展到多端直流、多换流站系统就只是增加节点和分支的问题了。7. 从会写到会用几个提高效率的实用技巧7.1 把失配方程函数写成向量化形式Matlab的循环效率远不如向量化计算。如果你的系统规模不大循环无所谓但一旦节点数量上百建议把交流导纳矩阵计算、功率失配计算都写成矩阵运算形式。比如用稀疏矩阵存储导纳矩阵用矢量化公式一次性计算所有节点的注入功率可以提速一个数量级以上。7.2 使用匿名函数简化控制模式切换控制模式切换可以用匿名函数或者函数句柄数组来组织。比如定义两个函数句柄calc_p_error (P_dc, U_dc, ref) ...然后在主循环里根据模式调用不同的句柄。代码会更清爽调试时也更容易单步跟踪。7.3 用disp和fprintf记录迭代轨迹统一迭代法的调试比纯交流潮流更依赖迭代轨迹信息。我习惯在每次迭代后打印迭代次数、最大失配量、直流节点电压、换流站输出电压幅值。这样一旦发散从打印数据里能快速看出是交流侧先发散还是直流侧先发散。很多问题看一眼轨迹曲线就能定位不用逐行查代码。7.4 学会用退化测试验证代码退化测试是代码正确性的试金石。最简单的一个把直流线路电阻设为零两个换流站都设定有功功率控制其他什么都不变。这时直流网络相当于一个无损耗功率传输通道两侧换流站的有功功率应该完全一致交流系统的功率平衡也应该是精确的。如果算出来两侧功率有微小偏差说明直流方程的符号或者方向定义有问题。8. 一些基于经验的小忠告交直流潮流计算的难点不在数学而在建模。很多人一上来就盯着雅可比矩阵的公式推导其实真正花时间的是理解每个控制模式对应的方程结构。控制模式一变状态量的索引要变失配方程要变雅可比矩阵的稀疏结构也要变。如果一开始就设计成“控制模式可配置”的程序结构后面扩展多端直流、混合直流都会省力很多。另一个容易被低估的问题是单位制的统一。我见过有人在Matlab里用有名值又混入标幺值参数结果迭代出来一些诡异的数值查了一晚上才发现是单位混乱。建议整个程序从头到尾统一用标幺值只在输入输出边界做转换。直流电压基准和交流电压基准分别取各自的额定电压功率基准统一取系统基准容量这样最稳妥。最后再分享一个小技巧调试时不要同时验证多个新功能。我每次只改一个点比如先只验证“换流站定有功功率模式”跑通了再加“定直流电压模式”再加“无功控制模式”。每加一个模式都用退化测试验证结果没有回归。这个习惯看上去慢但实际总耗时最少。交直流潮流是一个经典但很容易写错的地方希望这整套流程能帮你少走弯路。
