简介针对电力系统快速解耦法潮流计算与短路计算这份MATLAB算法讲解文档面向电力工程专业学生、科研人员及电网分析工程师系统拆解算法原理与程序实现。演示包共含1个doc文件大小仅384KB内容紧凑适合边看边练。文档从MATLAB矩阵运算与复数处理优势切入详细说明线路参数文件与节点状态文件的结构并通过循环自动识别最大节点编号动态生成节点导纳矩阵使程序可适配任意规模网络。为实现通用化程序内置seqencing()函数自动将PQ表排列为平衡节点、PV节点、PQ节点顺序并与支路参数对应同时给出潮流计算的完整迭代框架涵盖初值设置、收敛判断、修正方程求解、电压幅值及相角更新、结果输出等关键步骤。结合各函数作用与实际算例读者能快速理解快速解耦法的化简假设、计算流程及MATLAB实现细节进而提升电力系统建模与仿真能力。目前该资源已有214人学习下载。1. 为什么电力系统计算里还留着快速解耦法写电力系统潮流计算程序时多数人第一反应是牛顿-拉夫逊法迭代几次就收敛精度也不错。但实际节点规模一过千每次迭代都要重新计算雅可比矩阵内存和耗时的代价相当明显。快速解耦法是牛顿法在高压输电网条件下的工程近似它把有功-功角和无功-电压两组问题强行拆开用两个固定不变的常数矩阵反复迭代单次迭代开销小一个量级收敛速度只比牛顿法慢一两轮。这就是它在课程设计、毕业设计甚至在线计算场景里一直没被淘汰的原因。本文基于 matlab 实现电力系统快速解耦法潮流计算及短路计算程序先讲清楚 B 和 B 两个矩阵的来源与构造方法再给出完整的潮流迭代主循环代码最后在潮流初值上叠加故障分量完成三相短路和不对称短路计算。整个过程只依赖 matlab 基础函数不需要额外工具箱按步骤跑通后直接替换数据即可用于实际算例。2. 快速解耦法的B和B从牛顿法到定雅可比矩阵的简化路径2.1 牛顿法修正方程是怎么变成两组独立方程的完整的牛顿-拉夫逊法潮流修正方程为[ \begin{bmatrix} \Delta P \ \Delta Q \end{bmatrix}\begin{bmatrix} H N \ M L \end{bmatrix} \begin{bmatrix} \Delta \delta \ \Delta V / V \end{bmatrix} ]其中 H、N、M、L 是雅可比矩阵的四个分块N 反映无功对功角的耦合M 反映有功对电压的耦合。在高压输电网中线路电抗远大于电阻支路两端功角差通常在 10° 以内这时 cos(θij)≈1、sin(θij)≈0、Gij 远小于 Bij带入雅可比矩阵各项后N 和 M 两个耦合块趋近于零H 块近似为常数矩阵L 块也近似为常数矩阵。修正方程由此拆成两组独立方程[ \Delta P / V -B \Delta \delta ][ \Delta Q / V -B \Delta V ]这就是快速解耦法也称 P-Q 分解法的核心。注意这里的 ΔP 和 ΔQ 都除以了电压 V这是为了配合 B 和 B 的常数化处理。实际编程时直接计算 ΔP/V 和 ΔQ/V 作为残差向量即可不要漏掉这个 V。2.2 B和B矩阵的物理含义和形成规则B 对应有功-功角子问题它的元素定义为非对角元 B(i,j) -1/xij其中 xij 为支路电抗对角元 B(i,i) 与节点 i 相连的所有支路电抗倒数之和不含接地支路B 对应无功-电压子问题它的元素定义为非对角元 B(i,j) -1/xij有变压器时考虑变比折算对角元 B(i,i) 与节点 i 相连的支路电抗倒数之和加上该节点对地导纳含线路充电电容 B/2构造时的维数规则是两者最容易搞混的地方B 保留除平衡节点之外的所有节点即 PQ 节点加 PV 节点B 只保留 PQ 节点。原因是 PV 节点的电压幅值由励磁调节器维持给定值不需要参与无功迭代平衡节点的功角和电压都给定完全退出迭代。2.3 matlab中构造B和B的函数骨架直接基于支路数据逐条累加构造比先求完整节点导纳矩阵再从中抠子矩阵更直观也方便处理 B 不含接地支路这一要求。function [Bp, Bpp] build_B_matrices(node, branch) % node: 节点编号, 类型码(1平衡/2PV/3PQ), 有功, 无功, 电压初值 % branch: 首端, 末端, R, X, B/2, 变比 n size(node, 1); type node(:, 2); bp_idx find(type ~ 1); % 除平衡节点外全部参与B bpp_idx find(type 3); % 只有PQ节点参与B Bp zeros(length(bp_idx)); Bpp zeros(length(bpp_idx)); % 支路非对角元和对角元累加 for k 1:size(branch, 1) i branch(k, 1); j branch(k, 2); x branch(k, 4); tap branch(k, 6); % tap取1表示无变比 bi find(bp_idx i); bj find(bp_idx j); qi find(bpp_idx i); qj find(bpp_idx j); if ~isempty(bi) ~isempty(bj) Bp(bi, bj) Bp(bi, bj) - 1 / x; Bp(bj, bi) Bp(bj, bi) - 1 / x; Bp(bi, bi) Bp(bi, bi) 1 / x; Bp(bj, bj) Bp(bj, bj) 1 / x; end if ~isempty(qi) ~isempty(qj) Bpp(qi, qj) Bpp(qi, qj) - 1 / x / tap; Bpp(qj, qi) Bpp(qj, qi) - 1 / x / tap; Bpp(qi, qi) Bpp(qi, qi) 1 / x / tap; Bpp(qj, qj) Bpp(qj, qj) 1 / x / tap; end end % B需要补充对地支路线路充电电容B/2 for k 1:size(branch, 1) i branch(k, 1); j branch(k, 2); b_half branch(k, 5); qi find(bpp_idx i); qj find(bpp_idx j); if ~isempty(qi), Bpp(qi, qi) Bpp(qi, qi) b_half; end if ~isempty(qj), Bpp(qj, qj) Bpp(qj, qj) b_half; end end % 转稀疏矩阵大系统内存差距很大 Bp sparse(Bp); Bpp sparse(Bpp); end参数说明node 矩阵第二列的类型码决定了哪些节点进入哪个矩阵branch 第四列是电抗 X 而非阻抗模值填入时不要混入 R。变比 tap 在 B 中按 1/tap 折算在 B 中经典做法是不计非标准变比直接按线路电抗处理即可。这个函数把 R 完全丢弃这正是快速解耦法的前提条件。3. matlab实现快速解耦法潮流计算迭代主循环和收敛控制3.1 数据组织和节点编号方案先把数据格式定死后续迭代逻辑才不用反复改。常见做法是用两个矩阵node 存节点信息branch 存支路信息。节点编号建议重新排序平衡节点排第 1PV 节点其次PQ 节点放最后。这样 B 正好对应 node 矩阵后面的一段索引写代码时能少很多查找操作。矩阵列含义node1节点编号重新编号后node2类型1 平衡、2 PV、3 PQnode3注入有功 P发电机为正负荷为负node4注入无功 Q仅 PQ 节点有意义node5电压幅值初值branch1,2首端、末端节点编号branch3,4支路电阻 R、电抗 X标幺值branch5单端充电电容 B/2branch6变压器变比非变压器支路为 1对中小型系统初值一律取平启动所有 PQ 节点电压幅值 1.0所有节点功角 0PV 节点电压用给定值。快速解耦法对初值不敏感不需要像连续潮流那样做斜坡启动。迭代前要把潮流计算里会用到的完整导纳矩阵 Ybus 构造出来用于每次迭代计算节点注入功率。3.2 迭代主循环残差计算、常数矩阵分解、修正量更新快速解耦法省时间的核心在于 B 和 B 在迭代前只分解一次迭代中只做前代和回代。matlab 中推荐对两个矩阵做预先分解迭代循环内直接用分解对象求解避免每次迭代重新执行耗时的inv。% 构造Ybus用于功率残差计算 Y full_ybus(node, branch); % 自写函数按支路数据累加导纳 G real(Y); B imag(Y); n size(node, 1); V node(:, 5); Theta zeros(n, 1); Pspec node(:, 3); Qspec node(:, 4); pv find(node(:, 2) 2); pq find(node(:, 2) 3); slack find(node(:, 2) 1); bp_order [pv; pq]; % B的节点排列顺序 [Bp, Bpp] build_B_matrices(node, branch); dBP decomposition(Bp); dBPP decomposition(Bpp); max_iter 30; tol 1e-4; for k 1:max_iter % ---------- 有功迭代 ---------- dP zeros(length(bp_order), 1); for r 1:length(bp_order) i bp_order(r); Pcal V(i) * sum(V .* (G(i, :) * cos(Theta(i) - Theta) ... B(i, :) * sin(Theta(i) - Theta))); dP(r) (Pspec(i) - Pcal) / V(i); end dTheta dBP \ dP; for r 1:length(bp_order) Theta(bp_order(r)) Theta(bp_order(r)) dTheta(r); end % ---------- 无功迭代 ---------- dQ zeros(length(pq), 1); for r 1:length(pq) i pq(r); Qcal V(i) * sum(V .* (G(i, :) * sin(Theta(i) - Theta) ... - B(i, :) * cos(Theta(i) - Theta))); dQ(r) (Qspec(i) - Qcal) / V(i); end dV dBPP \ dQ; for r 1:length(pq) V(pq(r)) V(pq(r)) dV(r); end % ---------- 收敛检查 ---------- if max(abs(dP)) tol max(abs(dQ)) tol break; end end代码里有三个关键点。第一dTheta和dV解出后的放回位置必须与bp_order和pq的排列一致初学者最常见的错误就是把解向量顺序和节点编号混用导致修正量加错节点残差永远不收敛。第二Pcal和Qcal的符号约定与 node 矩阵第三、四列的注入功率符号必须一致——发电机为正、负荷为负Pcal 是节点所有连接支路流出的功率之和正值表示该节点实际向网络注入功率。第三decomposition返回的对象直接用反斜杠求解即可matlab 会自动选择合适的前代回代策略不要手动拆 L 和 U。3.3 无功越限处理PV节点转PQ节点实际系统中 PV 节点发电机的无功出力有上下限。迭代过程中若某 PV 节点的无功越限说明该节点电压幅值已无法维持给定值需要把它改为 PQ 节点无功取越限值电压由迭代结果决定。这个转换对 B 的影响是必须重新分解一次 B因为矩阵的维数变了。% 在每次无功迭代结束后检查PV节点无功是否越限 for i pv Qcal V(i) * sum(V .* (G(i, :) * sin(Theta(i) - Theta) ... - B(i, :) * cos(Theta(i) - Theta))); if Qcal node(i, 6) % node第6列存无功上限 node(i, 2) 3; % 改为PQ节点 node(i, 4) node(i, 6); % 无功锁定在上限 converted true; elseif Qcal node(i, 7) % node第7列存无功下限 node(i, 2) 3; node(i, 4) node(i, 7); converted true; end end if converted pq find(node(:, 2) 3); [~, Bpp] build_B_matrices(node, branch); dBPP decomposition(Bpp); end越限判断必须放在无功迭代结束、电压更新完之后不能在有功迭代里做因为此时 PV 节点的无功还没有更新。另外 PV 转 PQ 后B 不需要重新分解因为 B 本来就同时包含 PV 和 PQ 两批节点矩阵维数和内容都不变只有 B 需要重建。反过来 PQ 转 PV 的情况在常规迭代中基本不会出现程序里可以不写。3.4 收敛判据怎么选快速解耦法常见的收敛判据是检查不平衡功率相对值 ΔP/V 和 ΔQ/V 的最大绝对值。阈值从 0.01 到 1e-6 都有人用实际建议取 1e-4 p.u.。理由有三条0.01 p.u. 在 100MW 基准下约等于 1MW对潮流计算明显偏粗1e-6 p.u. 会把迭代次数推到 10 次以上快速解耦法省下的单次迭代时间会被多余的迭代次数吃掉而 1e-4 p.u. 对应的节点电压误差通常在 0.0001 p.u. 以内对后续短路计算和静态安全分析完全够用。判断收敛时用最大绝对值max(abs(dP))比用二范数更合理因为范数对个别残差大的节点不敏感一个节点出问题可能被其他正常节点平均掉导致误判为收敛。4. 短路计算在潮流初值上叠加故障分量4.1 短路计算的总体流程短路计算分两步先跑上一章的潮流程序得到故障前稳态各节点电压的幅值和相角都已知再在这个初值上叠加故障分量。故障分量法的核心是节点阻抗矩阵 Zbus短路电流只与故障点的自阻抗 Zkk 有关故障后各节点电压变化量只与互阻抗 Zik 有关计算量远小于重新求解全网方程。用快速解耦法跑潮流给短路计算提供初值比用牛顿法更合适的一点是短路计算中往往需要对多种预想故障重复计算快速解耦法每次都能在较少的迭代内给出足够精确的初值整体计算时间明显更短。4.2 用节点阻抗矩阵算三相短路三相短路是完全对称的故障类型不需要拆序网络。核心算法为故障电流 If Vk0 / (Zkk Zf)故障后节点电压 Vi Vi0 - Zik × If其中 Vk0 是故障点故障前电压相量来自潮流结果Zf 是过渡电阻金属性短路取 0。matlab 中小系统直接对 Ybus 求逆得到 Zbus 即可几百节点以内计算量完全可接受。% 假设潮流已跑完V0为节点电压幅值向量Theta0为相角向量 V0_phase V0 .* exp(1j * Theta0); % 转化为复数相量 Y full_ybus(node, branch); Zbus inv(Y); k 1; % 故障点节点编号 Vk0 V0_phase(k); Zkk Zbus(k, k); If Vk0 / Zkk; % 金属性三相短路电流 V_after V0_phase - Zbus(:, k) * If; % 故障后各节点电压相量这里的If是复数代表短路电流的幅值和相位。V_after是故障后全部节点的电压相量要查看某节点电压幅值只需取abs(V_after(i))相角取angle(V_after(i))。如果故障点在线路中间而非节点上需要把故障点作为新增节点并入网络再重新形成 Zbus课程设计中常用近似做法是把故障点放到线路两端节点之一或者把线路按故障比例拆成两条支路加入支路表。4.3 不对称短路的序网络连接不对称短路要拆成正序、负序和零序三个序网络。正序网络就是正常网络的参数负序网络在静态分析中近似取与正序相同零序网络需要单独输入零序阻抗参数。三个序网在故障点的连接方式取决于故障类型下面这张表是程序实现时直接查的故障类型序网连接方式正序故障电流表达式三相短路仅正序网络If1 Vk0 / Zkk1单相接地三序网络串联If1 Vk0 / (Zkk1 Zkk2 Zkk0)两相短路正序与负序并联If1 Vk0 / (Zkk1 Zkk2)两相接地正序与负序-零序并联支路串联If1 Vk0 / (Zkk1 (Zkk2 ∥ Zkk0))matlab 实现时需要分别构造正序、负序和零序的 Zbus。单相接地时短路点 A 相电流等于 3 倍正序电流B 相和 C 相电流为零需要通过对称分量变换矩阵把序量转换回 abc 三相量。这里最容易出错的是变换矩阵元素的正负号建议对照电力系统分析教科目录中的对称分量变换矩阵逐项核对或者用单机无穷大系统手算结果做验证。故障后的支路电流如果也要输出把短路点电流折算到该节点注入电流中重新做一次全网求解即可。5. 让程序在大算例上站稳验证手段和提速技巧5.1 用IEEE标准节点数据和MATPOWER做对照手算小系统只能验证接线逻辑真正检验快速解耦法实现对不对要用 IEEE 标准算例。常见做法是找 IEEE 9 节点或 14 节点系统数据把支路参数转成标幺值按第 3.1 节的数据格式填入然后用 MATLAB 自带基础函数跑一遍与已知结果的电压幅值逐节点对比。校核项检查方法可接受误差电压幅值逐节点对比 0.001 p.u.电压相角逐节点对比 0.01°平衡节点功率与全网损耗对比 0.1%PV节点无功与基准潮流对比 0.5 Mvar标幺值 0.005动辄对比几十个节点时直接打印全表肉眼检查不现实建议用max(abs(V_calc - V_ref))一条语句输出最大误差单独看这一个数就够了。误差超过 0.001 p.u. 再逐节点定位重点检查是哪个区域的节点偏差最大通常能直接指向某条支路的参数填错或 B 矩阵缺失了对地支路。5.2 不收敛的三个排查方向快速解耦法不收敛绝大多数不是迭代逻辑问题而是数据或参数问题。第一r/x 比值过大。线路电阻远大于电抗时有功只跟功角有关的假设不成立B 丢弃电阻的简化直接导致修正方向错误。这类问题典型出现在 10kV 及以下电压等级的配电网中输电网 220kV 以上基本不会遇到。解决方法有两条改用保留电阻项的 XB 法或 BX 法改进型快速解耦法或者直接换牛顿法。第二B 矩阵构造或节点顺序出错。正常情况下 B 的对角元为正、非对角元为负。如果出现正的非对角元说明某条支路的电抗符号填反了如果某个 PQ 节点的对角元为 0则说明该节点没有任何支路连接数据里存在孤岛节点。调试时把 B 打印出来检查每一行非零元的位置是否与该节点的实际连接关系一致能排查掉大多数矩阵构造错误。第三PV 转 PQ 后没有重新分解 B。PV 转 PQ 后 PQ 集合变大B 的维数增大若不重新生成并分解矩阵程序会在下一次无功迭代时因矩阵阶数不匹配直接报错或者使用旧矩阵导致计算混乱。把越限判断放在无功迭代结束后、下一个迭代周期开始前用标志位控制只在必要时重新分解。5.3 提速与验证的小技巧大系统下形成 B 和 B 时尽量用稀疏矩阵填充避免用全零矩阵逐元素赋值再转稀疏。构造过程中至少要先预分配非零元素个数或者直接用sparse(i, j, x, m, n)三参数形式从坐标向量生成矩阵内存开销差别很大。迭代部分用decomposition预分解后反复回代对几百节点的系统提速非常明显。短路计算模块的验证推荐构造一个单机无穷大系统一台发电机经单回线接无穷大母线手算出线路末端三相短路电流再和程序结果对比。这个系统没有 PV 与 PQ 的复杂交互潮流结果也足够简单能直接检验 Zbus 形成和叠加定理的相位处理是否正确。最后检查时把每次迭代的dP和dQ最大残差打印出来观察变化趋势。正常收敛时两个残差单调下降如果曲线先降后升说明某次节点类型变换后矩阵没有同步更新直接去查越限处理分支是否有converted标志位遗漏的情况。本文还有配套的精品资源点击获取 MATLAB实现电力系统快速解耦法潮流计算与短路分析
简介针对电力系统快速解耦法潮流计算与短路计算这份MATLAB算法讲解文档面向电力工程专业学生、科研人员及电网分析工程师系统拆解算法原理与程序实现。演示包共含1个doc文件大小仅384KB内容紧凑适合边看边练。文档从MATLAB矩阵运算与复数处理优势切入详细说明线路参数文件与节点状态文件的结构并通过循环自动识别最大节点编号动态生成节点导纳矩阵使程序可适配任意规模网络。为实现通用化程序内置seqencing()函数自动将PQ表排列为平衡节点、PV节点、PQ节点顺序并与支路参数对应同时给出潮流计算的完整迭代框架涵盖初值设置、收敛判断、修正方程求解、电压幅值及相角更新、结果输出等关键步骤。结合各函数作用与实际算例读者能快速理解快速解耦法的化简假设、计算流程及MATLAB实现细节进而提升电力系统建模与仿真能力。目前该资源已有214人学习下载。1. 为什么电力系统计算里还留着快速解耦法写电力系统潮流计算程序时多数人第一反应是牛顿-拉夫逊法迭代几次就收敛精度也不错。但实际节点规模一过千每次迭代都要重新计算雅可比矩阵内存和耗时的代价相当明显。快速解耦法是牛顿法在高压输电网条件下的工程近似它把有功-功角和无功-电压两组问题强行拆开用两个固定不变的常数矩阵反复迭代单次迭代开销小一个量级收敛速度只比牛顿法慢一两轮。这就是它在课程设计、毕业设计甚至在线计算场景里一直没被淘汰的原因。本文基于 matlab 实现电力系统快速解耦法潮流计算及短路计算程序先讲清楚 B 和 B 两个矩阵的来源与构造方法再给出完整的潮流迭代主循环代码最后在潮流初值上叠加故障分量完成三相短路和不对称短路计算。整个过程只依赖 matlab 基础函数不需要额外工具箱按步骤跑通后直接替换数据即可用于实际算例。2. 快速解耦法的B和B从牛顿法到定雅可比矩阵的简化路径2.1 牛顿法修正方程是怎么变成两组独立方程的完整的牛顿-拉夫逊法潮流修正方程为[ \begin{bmatrix} \Delta P \ \Delta Q \end{bmatrix}\begin{bmatrix} H N \ M L \end{bmatrix} \begin{bmatrix} \Delta \delta \ \Delta V / V \end{bmatrix} ]其中 H、N、M、L 是雅可比矩阵的四个分块N 反映无功对功角的耦合M 反映有功对电压的耦合。在高压输电网中线路电抗远大于电阻支路两端功角差通常在 10° 以内这时 cos(θij)≈1、sin(θij)≈0、Gij 远小于 Bij带入雅可比矩阵各项后N 和 M 两个耦合块趋近于零H 块近似为常数矩阵L 块也近似为常数矩阵。修正方程由此拆成两组独立方程[ \Delta P / V -B \Delta \delta ][ \Delta Q / V -B \Delta V ]这就是快速解耦法也称 P-Q 分解法的核心。注意这里的 ΔP 和 ΔQ 都除以了电压 V这是为了配合 B 和 B 的常数化处理。实际编程时直接计算 ΔP/V 和 ΔQ/V 作为残差向量即可不要漏掉这个 V。2.2 B和B矩阵的物理含义和形成规则B 对应有功-功角子问题它的元素定义为非对角元 B(i,j) -1/xij其中 xij 为支路电抗对角元 B(i,i) 与节点 i 相连的所有支路电抗倒数之和不含接地支路B 对应无功-电压子问题它的元素定义为非对角元 B(i,j) -1/xij有变压器时考虑变比折算对角元 B(i,i) 与节点 i 相连的支路电抗倒数之和加上该节点对地导纳含线路充电电容 B/2构造时的维数规则是两者最容易搞混的地方B 保留除平衡节点之外的所有节点即 PQ 节点加 PV 节点B 只保留 PQ 节点。原因是 PV 节点的电压幅值由励磁调节器维持给定值不需要参与无功迭代平衡节点的功角和电压都给定完全退出迭代。2.3 matlab中构造B和B的函数骨架直接基于支路数据逐条累加构造比先求完整节点导纳矩阵再从中抠子矩阵更直观也方便处理 B 不含接地支路这一要求。function [Bp, Bpp] build_B_matrices(node, branch) % node: 节点编号, 类型码(1平衡/2PV/3PQ), 有功, 无功, 电压初值 % branch: 首端, 末端, R, X, B/2, 变比 n size(node, 1); type node(:, 2); bp_idx find(type ~ 1); % 除平衡节点外全部参与B bpp_idx find(type 3); % 只有PQ节点参与B Bp zeros(length(bp_idx)); Bpp zeros(length(bpp_idx)); % 支路非对角元和对角元累加 for k 1:size(branch, 1) i branch(k, 1); j branch(k, 2); x branch(k, 4); tap branch(k, 6); % tap取1表示无变比 bi find(bp_idx i); bj find(bp_idx j); qi find(bpp_idx i); qj find(bpp_idx j); if ~isempty(bi) ~isempty(bj) Bp(bi, bj) Bp(bi, bj) - 1 / x; Bp(bj, bi) Bp(bj, bi) - 1 / x; Bp(bi, bi) Bp(bi, bi) 1 / x; Bp(bj, bj) Bp(bj, bj) 1 / x; end if ~isempty(qi) ~isempty(qj) Bpp(qi, qj) Bpp(qi, qj) - 1 / x / tap; Bpp(qj, qi) Bpp(qj, qi) - 1 / x / tap; Bpp(qi, qi) Bpp(qi, qi) 1 / x / tap; Bpp(qj, qj) Bpp(qj, qj) 1 / x / tap; end end % B需要补充对地支路线路充电电容B/2 for k 1:size(branch, 1) i branch(k, 1); j branch(k, 2); b_half branch(k, 5); qi find(bpp_idx i); qj find(bpp_idx j); if ~isempty(qi), Bpp(qi, qi) Bpp(qi, qi) b_half; end if ~isempty(qj), Bpp(qj, qj) Bpp(qj, qj) b_half; end end % 转稀疏矩阵大系统内存差距很大 Bp sparse(Bp); Bpp sparse(Bpp); end参数说明node 矩阵第二列的类型码决定了哪些节点进入哪个矩阵branch 第四列是电抗 X 而非阻抗模值填入时不要混入 R。变比 tap 在 B 中按 1/tap 折算在 B 中经典做法是不计非标准变比直接按线路电抗处理即可。这个函数把 R 完全丢弃这正是快速解耦法的前提条件。3. matlab实现快速解耦法潮流计算迭代主循环和收敛控制3.1 数据组织和节点编号方案先把数据格式定死后续迭代逻辑才不用反复改。常见做法是用两个矩阵node 存节点信息branch 存支路信息。节点编号建议重新排序平衡节点排第 1PV 节点其次PQ 节点放最后。这样 B 正好对应 node 矩阵后面的一段索引写代码时能少很多查找操作。矩阵列含义node1节点编号重新编号后node2类型1 平衡、2 PV、3 PQnode3注入有功 P发电机为正负荷为负node4注入无功 Q仅 PQ 节点有意义node5电压幅值初值branch1,2首端、末端节点编号branch3,4支路电阻 R、电抗 X标幺值branch5单端充电电容 B/2branch6变压器变比非变压器支路为 1对中小型系统初值一律取平启动所有 PQ 节点电压幅值 1.0所有节点功角 0PV 节点电压用给定值。快速解耦法对初值不敏感不需要像连续潮流那样做斜坡启动。迭代前要把潮流计算里会用到的完整导纳矩阵 Ybus 构造出来用于每次迭代计算节点注入功率。3.2 迭代主循环残差计算、常数矩阵分解、修正量更新快速解耦法省时间的核心在于 B 和 B 在迭代前只分解一次迭代中只做前代和回代。matlab 中推荐对两个矩阵做预先分解迭代循环内直接用分解对象求解避免每次迭代重新执行耗时的inv。% 构造Ybus用于功率残差计算 Y full_ybus(node, branch); % 自写函数按支路数据累加导纳 G real(Y); B imag(Y); n size(node, 1); V node(:, 5); Theta zeros(n, 1); Pspec node(:, 3); Qspec node(:, 4); pv find(node(:, 2) 2); pq find(node(:, 2) 3); slack find(node(:, 2) 1); bp_order [pv; pq]; % B的节点排列顺序 [Bp, Bpp] build_B_matrices(node, branch); dBP decomposition(Bp); dBPP decomposition(Bpp); max_iter 30; tol 1e-4; for k 1:max_iter % ---------- 有功迭代 ---------- dP zeros(length(bp_order), 1); for r 1:length(bp_order) i bp_order(r); Pcal V(i) * sum(V .* (G(i, :) * cos(Theta(i) - Theta) ... B(i, :) * sin(Theta(i) - Theta))); dP(r) (Pspec(i) - Pcal) / V(i); end dTheta dBP \ dP; for r 1:length(bp_order) Theta(bp_order(r)) Theta(bp_order(r)) dTheta(r); end % ---------- 无功迭代 ---------- dQ zeros(length(pq), 1); for r 1:length(pq) i pq(r); Qcal V(i) * sum(V .* (G(i, :) * sin(Theta(i) - Theta) ... - B(i, :) * cos(Theta(i) - Theta))); dQ(r) (Qspec(i) - Qcal) / V(i); end dV dBPP \ dQ; for r 1:length(pq) V(pq(r)) V(pq(r)) dV(r); end % ---------- 收敛检查 ---------- if max(abs(dP)) tol max(abs(dQ)) tol break; end end代码里有三个关键点。第一dTheta和dV解出后的放回位置必须与bp_order和pq的排列一致初学者最常见的错误就是把解向量顺序和节点编号混用导致修正量加错节点残差永远不收敛。第二Pcal和Qcal的符号约定与 node 矩阵第三、四列的注入功率符号必须一致——发电机为正、负荷为负Pcal 是节点所有连接支路流出的功率之和正值表示该节点实际向网络注入功率。第三decomposition返回的对象直接用反斜杠求解即可matlab 会自动选择合适的前代回代策略不要手动拆 L 和 U。3.3 无功越限处理PV节点转PQ节点实际系统中 PV 节点发电机的无功出力有上下限。迭代过程中若某 PV 节点的无功越限说明该节点电压幅值已无法维持给定值需要把它改为 PQ 节点无功取越限值电压由迭代结果决定。这个转换对 B 的影响是必须重新分解一次 B因为矩阵的维数变了。% 在每次无功迭代结束后检查PV节点无功是否越限 for i pv Qcal V(i) * sum(V .* (G(i, :) * sin(Theta(i) - Theta) ... - B(i, :) * cos(Theta(i) - Theta))); if Qcal node(i, 6) % node第6列存无功上限 node(i, 2) 3; % 改为PQ节点 node(i, 4) node(i, 6); % 无功锁定在上限 converted true; elseif Qcal node(i, 7) % node第7列存无功下限 node(i, 2) 3; node(i, 4) node(i, 7); converted true; end end if converted pq find(node(:, 2) 3); [~, Bpp] build_B_matrices(node, branch); dBPP decomposition(Bpp); end越限判断必须放在无功迭代结束、电压更新完之后不能在有功迭代里做因为此时 PV 节点的无功还没有更新。另外 PV 转 PQ 后B 不需要重新分解因为 B 本来就同时包含 PV 和 PQ 两批节点矩阵维数和内容都不变只有 B 需要重建。反过来 PQ 转 PV 的情况在常规迭代中基本不会出现程序里可以不写。3.4 收敛判据怎么选快速解耦法常见的收敛判据是检查不平衡功率相对值 ΔP/V 和 ΔQ/V 的最大绝对值。阈值从 0.01 到 1e-6 都有人用实际建议取 1e-4 p.u.。理由有三条0.01 p.u. 在 100MW 基准下约等于 1MW对潮流计算明显偏粗1e-6 p.u. 会把迭代次数推到 10 次以上快速解耦法省下的单次迭代时间会被多余的迭代次数吃掉而 1e-4 p.u. 对应的节点电压误差通常在 0.0001 p.u. 以内对后续短路计算和静态安全分析完全够用。判断收敛时用最大绝对值max(abs(dP))比用二范数更合理因为范数对个别残差大的节点不敏感一个节点出问题可能被其他正常节点平均掉导致误判为收敛。4. 短路计算在潮流初值上叠加故障分量4.1 短路计算的总体流程短路计算分两步先跑上一章的潮流程序得到故障前稳态各节点电压的幅值和相角都已知再在这个初值上叠加故障分量。故障分量法的核心是节点阻抗矩阵 Zbus短路电流只与故障点的自阻抗 Zkk 有关故障后各节点电压变化量只与互阻抗 Zik 有关计算量远小于重新求解全网方程。用快速解耦法跑潮流给短路计算提供初值比用牛顿法更合适的一点是短路计算中往往需要对多种预想故障重复计算快速解耦法每次都能在较少的迭代内给出足够精确的初值整体计算时间明显更短。4.2 用节点阻抗矩阵算三相短路三相短路是完全对称的故障类型不需要拆序网络。核心算法为故障电流 If Vk0 / (Zkk Zf)故障后节点电压 Vi Vi0 - Zik × If其中 Vk0 是故障点故障前电压相量来自潮流结果Zf 是过渡电阻金属性短路取 0。matlab 中小系统直接对 Ybus 求逆得到 Zbus 即可几百节点以内计算量完全可接受。% 假设潮流已跑完V0为节点电压幅值向量Theta0为相角向量 V0_phase V0 .* exp(1j * Theta0); % 转化为复数相量 Y full_ybus(node, branch); Zbus inv(Y); k 1; % 故障点节点编号 Vk0 V0_phase(k); Zkk Zbus(k, k); If Vk0 / Zkk; % 金属性三相短路电流 V_after V0_phase - Zbus(:, k) * If; % 故障后各节点电压相量这里的If是复数代表短路电流的幅值和相位。V_after是故障后全部节点的电压相量要查看某节点电压幅值只需取abs(V_after(i))相角取angle(V_after(i))。如果故障点在线路中间而非节点上需要把故障点作为新增节点并入网络再重新形成 Zbus课程设计中常用近似做法是把故障点放到线路两端节点之一或者把线路按故障比例拆成两条支路加入支路表。4.3 不对称短路的序网络连接不对称短路要拆成正序、负序和零序三个序网络。正序网络就是正常网络的参数负序网络在静态分析中近似取与正序相同零序网络需要单独输入零序阻抗参数。三个序网在故障点的连接方式取决于故障类型下面这张表是程序实现时直接查的故障类型序网连接方式正序故障电流表达式三相短路仅正序网络If1 Vk0 / Zkk1单相接地三序网络串联If1 Vk0 / (Zkk1 Zkk2 Zkk0)两相短路正序与负序并联If1 Vk0 / (Zkk1 Zkk2)两相接地正序与负序-零序并联支路串联If1 Vk0 / (Zkk1 (Zkk2 ∥ Zkk0))matlab 实现时需要分别构造正序、负序和零序的 Zbus。单相接地时短路点 A 相电流等于 3 倍正序电流B 相和 C 相电流为零需要通过对称分量变换矩阵把序量转换回 abc 三相量。这里最容易出错的是变换矩阵元素的正负号建议对照电力系统分析教科目录中的对称分量变换矩阵逐项核对或者用单机无穷大系统手算结果做验证。故障后的支路电流如果也要输出把短路点电流折算到该节点注入电流中重新做一次全网求解即可。5. 让程序在大算例上站稳验证手段和提速技巧5.1 用IEEE标准节点数据和MATPOWER做对照手算小系统只能验证接线逻辑真正检验快速解耦法实现对不对要用 IEEE 标准算例。常见做法是找 IEEE 9 节点或 14 节点系统数据把支路参数转成标幺值按第 3.1 节的数据格式填入然后用 MATLAB 自带基础函数跑一遍与已知结果的电压幅值逐节点对比。校核项检查方法可接受误差电压幅值逐节点对比 0.001 p.u.电压相角逐节点对比 0.01°平衡节点功率与全网损耗对比 0.1%PV节点无功与基准潮流对比 0.5 Mvar标幺值 0.005动辄对比几十个节点时直接打印全表肉眼检查不现实建议用max(abs(V_calc - V_ref))一条语句输出最大误差单独看这一个数就够了。误差超过 0.001 p.u. 再逐节点定位重点检查是哪个区域的节点偏差最大通常能直接指向某条支路的参数填错或 B 矩阵缺失了对地支路。5.2 不收敛的三个排查方向快速解耦法不收敛绝大多数不是迭代逻辑问题而是数据或参数问题。第一r/x 比值过大。线路电阻远大于电抗时有功只跟功角有关的假设不成立B 丢弃电阻的简化直接导致修正方向错误。这类问题典型出现在 10kV 及以下电压等级的配电网中输电网 220kV 以上基本不会遇到。解决方法有两条改用保留电阻项的 XB 法或 BX 法改进型快速解耦法或者直接换牛顿法。第二B 矩阵构造或节点顺序出错。正常情况下 B 的对角元为正、非对角元为负。如果出现正的非对角元说明某条支路的电抗符号填反了如果某个 PQ 节点的对角元为 0则说明该节点没有任何支路连接数据里存在孤岛节点。调试时把 B 打印出来检查每一行非零元的位置是否与该节点的实际连接关系一致能排查掉大多数矩阵构造错误。第三PV 转 PQ 后没有重新分解 B。PV 转 PQ 后 PQ 集合变大B 的维数增大若不重新生成并分解矩阵程序会在下一次无功迭代时因矩阵阶数不匹配直接报错或者使用旧矩阵导致计算混乱。把越限判断放在无功迭代结束后、下一个迭代周期开始前用标志位控制只在必要时重新分解。5.3 提速与验证的小技巧大系统下形成 B 和 B 时尽量用稀疏矩阵填充避免用全零矩阵逐元素赋值再转稀疏。构造过程中至少要先预分配非零元素个数或者直接用sparse(i, j, x, m, n)三参数形式从坐标向量生成矩阵内存开销差别很大。迭代部分用decomposition预分解后反复回代对几百节点的系统提速非常明显。短路计算模块的验证推荐构造一个单机无穷大系统一台发电机经单回线接无穷大母线手算出线路末端三相短路电流再和程序结果对比。这个系统没有 PV 与 PQ 的复杂交互潮流结果也足够简单能直接检验 Zbus 形成和叠加定理的相位处理是否正确。最后检查时把每次迭代的dP和dQ最大残差打印出来观察变化趋势。正常收敛时两个残差单调下降如果曲线先降后升说明某次节点类型变换后矩阵没有同步更新直接去查越限处理分支是否有converted标志位遗漏的情况。本文还有配套的精品资源点击获取