第一次接触潮流计算时我对着两节点电路图犯过嘀咕这不过是一个电源、一条线路、一个负荷的简单电路欧姆定律加基尔霍夫电压定律分分钟就能列方程为什么还要专门搞一套 Gauss-Seidel 迭代直到自己动手在 MATLAB 里把代码写出来才意识到哪怕只有两个节点母线 2 的电压也会顽固地出现在方程右边的分母里而且是以共轭形式出现。这是个彻头彻尾的非线性问题手算解析解不是不行但换到几十上百个节点的真实电网就彻底没戏。这篇文章就从两节点系统入手讲清楚 Gauss-Seidel 潮流计算的完整建模、MATLAB 实现、收敛过程以及我实际调试中踩过的一堆坑。1. 两节点也要迭代先想清楚潮流问题到底在求什么1.1 一个看起来“很简单”的两节点系统先把这个两节点系统的算例参数摆出来母线 1平衡节点可以理解成大电网或者发电机电压幅值固定为 V1 1.05 pu相角为 0°母线 2PQ 节点负荷吸收有功 0.8 pu、无功 0.6 pu也就是注入功率 P2 -0.8 pu、Q2 -0.6 pu线路阻抗Z 0.05 j0.2 pu电抗远大于电阻符合高压架空线的典型特征如果按照中学物理的思路电源电动势 1.05线路阻抗已知负载功率已知似乎应该“一步就能算出来”。但这里有个关键区别电力系统里的负荷模型是恒功率不是恒阻抗。恒阻抗负载可以用欧姆定律直接算但恒功率负荷的意思是不管母线 2 电压降到多少它消耗的有功和无功都保持 0.8 j0.6 pu。这就把问题变成了一个耦合的非线性方程——我们要求解的是什么样的 V2 才能让节点 2 的注入功率恰好等于给定的 P2 jQ21.2 为什么解析法在这里卡壳把节点 2 的功率方程写出来就明白了。节点注入电流方向取“流入网络为正”那么从节点 1 流向节点 2 的线路电流是 (V1 - V2)/Z而节点 2 的注入电流是它的相反数I2 (V2 - V1)/Z又因为复功率的定义是 S V · conj(I)所以S2 V2 · conj(I2) V2 · conj((V2 - V1)/Z)把 V1 视为实数 1.05展开S2 V2 · [conj(V2) - V1] / conj(Z) (|V2|² - V1·conj(V2)) / conj(Z)这个式子同时包含 V2 和 conj(V2)也就是同时包含实部虚部两个独立变量拆开写就是两个耦合的二元二次方程。手算两节点系统还能硬解但这个方法不可扩展——真实电网动辄几百个节点每个节点的功率方程都和其他所有节点电压耦合解析法在规模面前毫无意义。1.3 Gauss-Seidel 的思路把方程变成迭代Gauss-Seidel 法做的事情很简单把节点电压方程改写成不动点迭代格式。对 PQ 节点 i核心迭代公式是V_i^(k1) (1/Y_ii) · [conj(S_i)/conj(V_i^(k)) - Σ(Y_ij · V_j)]其中已经更新过的节点用新值没更新的用旧值——这就是“Gauss-Seidel”和“Jacobi 法”的本质区别。每次用当前估计值算出新的电压估计反复迭代直到修正量足够小。两节点系统是这条思路的最小实现但已经包含全部核心要素导纳矩阵、共轭复功率、节点电压耦合。2. 节点三分法与两节点系统的建模细节2.1 PQ、PV、平衡节点各自要什么、给什么潮流计算里节点不是随便分类的每种节点对应一组已知量和待求量背后的工程含义完全不同节点类型已知量待求量工程角色PQ 节点有功 P、无功 Q电压幅值 V、相角 δ负荷节点、恒功率负荷绝大多数节点属于此类PV 节点有功 P、电压幅值 V无功 Q、相角 δ发电机节点有自动电压调节器AVR能维持机端电压平衡节点电压幅值 V、相角 δ有功 P、无功 Q参考节点承担全网功率不平衡量相位基准为什么必须有一个平衡节点因为全网功率损耗在潮流计算之前是未知的总要有一个节点吸收“发电总功率 - 负荷总功率 - 网损”这个差值。另外电压相角是相对量所有母线相角需要一个基准平衡节点就是这个基准通常取 V 1.0∠0° 或 1.05∠0°。2.2 本算例的参数设定与标幺值为了让大家直接复现我把参数统一用标幺值pu表示。本算例的具体设定母线 1V1 1.05 puθ1 0°类型为平衡节点母线 2P2 -0.8 puQ2 -0.6 pu类型为 PQ 节点负号表示吸收功率线路阻抗Z 0.05 j0.2 pu如果用 100 MVA 作基准容量那么 0.8 pu 的有功就是 80 MW0.6 pu 的无功就是 60 MvarV 1.05 pu 在 220 kV 系统里对应 231 kV。标幺值的意义在于把电压、功率、阻抗都归一化到 1 附近数值稳定性好而且三相系统用线电压标幺、三相功率标幺之后可以直接用单相等效电路建模非常方便。2.3 符号约定是最大的坑我在给研究生改代码时发现Gauss-Seidel 潮流代码跑不出来的原因十有八九是符号约定乱了。这里必须提前说清楚本代码中 P、Q 的正方向定义为“注入节点”方向所以负荷吸收功率时 P、Q 取负值复功率定义S V · conj(I)因此节点注入电流 I conj(S) / conj(V)注意是 conj(S) 除以 conj(V)不是 S 除以 V。这个错误会让迭代结果完全错误甚至发散有的教材把 P、Q 正方向定义为“流出节点”即负荷方向那公式里的符号会整体变号。写代码前先定好约定然后所有公式保持一致。我见过有人前五分钟用注入方向定义 S后五分钟又从负荷方向抄公式结果自然是灾难。2.4 节点导纳矩阵 Ybus 的构建两节点系统忽略对地电容后导纳矩阵非常简洁Ybus [1/Z, -1/Z; -1/Z, 1/Z]对角线元素是自导纳等于与该节点相连的所有支路导纳之和非对角线元素是互导纳等于连接两节点的支路导纳的负值。这个结论在没有变压器和移相器时成立。以后扩展到多节点系统Ybus 的构建规律是一样的自导纳累加、互导纳取负。3. Gauss-Seidel 迭代公式推导与 MATLAB 实现3.1 从节点电压方程到迭代格式对于节点 2节点电压方程是I2 Y21·V1 Y22·V2把 I2 conj(S2) / conj(V2) 代入conj(S2)/conj(V2) Y21·V1 Y22·V2移项得到V2 (1/Y22) · [conj(S2)/conj(V2) - Y21·V1]这就是 PQ 节点的 Gauss-Seidel 迭代式。注意 V1 是平衡节点电压恒定不参与迭代但如果系统里有多个 PQ 节点那么每个新算出的电压都要立刻用于后续节点计算这正是 Gauss-Seidel 和 Jacobi 的区别。3.2 完整可运行的 MATLAB 代码下面给出我调试过的完整代码。注释写得比较详细可以粘贴到 MATLAB 里直接运行%% 两节点电力系统 Gauss-Seidel 潮流计算 % 母线1平衡节点V1 已知恒定 % 母线2PQ节点P2、Q2 已知求电压幅值和相角 % 所有量均为标幺值pu clear; clc; close all; %% ---------- 1. 系统参数 ---------- % 线路阻抗标幺值架空线电抗远大于电阻 Z 0.05 1j*0.2; % 节点导纳矩阵 Ybus2x2忽略对地导纳 Y11 1/Z; Y12 -1/Z; Y21 -1/Z; Y22 1/Z; Ybus [Y11 Y12; Y21 Y22]; %% ---------- 2. 节点数据 ---------- % 平衡节点节点1 V1 1.05; % 幅值 theta1 0; % 相角rad V1_phasor V1 * exp(1j*theta1); % PQ节点节点2注入节点方向为正负荷为负 P2 -0.8; % 吸收有功 0.8 pu Q2 -0.6; % 吸收无功 0.6 pu S2 P2 1j*Q2; % 节点2的注入复功率 %% ---------- 3. 初始化 ---------- V2 1.0 1j*0.0; % 平启动flat start tol 1e-7; % 电压修正量收敛容忍度 maxIter 50; % 最大迭代次数 fprintf( 两节点系统 Gauss-Seidel 潮流计算 \n); fprintf(线路阻抗 Z %.3f j%.3f pu\n, real(Z), imag(Z)); fprintf(平衡节点 V1 %.2f∠%.1f°\n, abs(V1_phasor), angle(V1_phasor)*180/pi); fprintf(PQ节点 P2 %.2f pu, Q2 %.2f pu\n\n, P2, Q2); % 记录每次迭代的电压便于画收敛曲线 V2_history zeros(maxIter1,1); V2_history(1) V2; %% ---------- 4. Gauss-Seidel 主迭代 ---------- for k 1:maxIter % 核心迭代公式 % V2_new (1/Y22) * [ conj(S2)/conj(V2) - Y21*V1 ] % 其中 conj(S2) P2 - j*Q2 V2_new (1/Y22) * ( conj(S2)/conj(V2) - Y21*V1_phasor ); dV abs(V2_new - V2); V2 V2_new; V2_history(k1) V2; fprintf(k%2d | V2%.8f %.8fj | |V2|%.8f | δ2%9.6f° | ΔV%.2e\n, ... k, real(V2), imag(V2), abs(V2), angle(V2)*180/pi, dV); if dV tol fprintf(\n第 %d 次迭代收敛ΔV %.0e\n\n, k, tol); break; end if k maxIter error(达到最大迭代次数仍未收敛检查参数或减小负荷); end end %% ---------- 5. 结果输出与验证 ---------- fprintf( 计算结果 \n); fprintf(母线2电压幅值 |V2| %.6f pu\n, abs(V2)); fprintf(母线2电压相角 δ2 %.6f°\n, angle(V2)*180/pi); fprintf(母线2电压相量 V2 %.6f ∠ %.6f° pu\n\n, abs(V2), angle(V2)*180/pi); % 线路电流正方向 1-2及相关功率 I12 (V1_phasor - V2) / Z; S12 V1_phasor * conj(I12); % 节点1注入线路的功率 S21 V2 * conj(I12); % 节点2从线路接收的功率 fprintf(线路电流 I12 %.6f ∠ %.6f° pu\n, abs(I12), angle(I12)*180/pi); fprintf(节点1注入线路 S12 %.6f %.6fj pu\n, real(S12), imag(S12)); fprintf(节点2接收线路 S21 %.6f %.6fj pu\n, real(S21), imag(S21)); fprintf(验证S21 应约等于负荷功率 (0.8j0.6) pu\n); %% ---------- 6. 收敛曲线 ---------- iterUsed k; % 实际迭代次数 figure(Color,w); subplot(1,2,1); plot(0:iterUsed, abs(V2_history(1:iterUsed1)), b-o, LineWidth,1.2, MarkerSize,5); xlabel(迭代次数 k); ylabel(|V_2| (pu)); title(PQ节点电压幅值收敛过程); grid on; subplot(1,2,2); plot(0:iterUsed, angle(V2_history(1:iterUsed1))*180/pi, r-s, LineWidth,1.2, MarkerSize,5); xlabel(迭代次数 k); ylabel(δ_2 (deg)); title(PQ节点电压相角收敛过程); grid on;3.3 单次迭代的逐步解读用第一次迭代演示这个公式到底在算什么。初始值 V2 1 j0conj(S2) -0.8 j0.6conj(S2) / conj(V2) -0.8 j0.6因为 conj(V2)1-Y21·V1 -(-1/Z)·1.05 1.05/(0.05j0.2)这是由平衡节点通过互导纳注入耦合的部分括号内求和后乘以 1/Y22 Z最终得到 V2^(1) ≈ 0.8900 - j0.1300即第一次迭代就从 1.0∠0° 跳到了约 0.899∠-8.3°。大跨步一步到位之后每次修正量迅速衰减。这就是 Gauss-Seidel 法在简单系统中的典型表现——开始几步下降很快越往后越慢呈线性收敛。4. 收敛过程实测初始值、判据与松弛因子的那些事4.1 我的机器上跑出来的迭代记录在 MATLAB R2023b 下运行收敛容忍度设为 1e-7 时大约 9 到 10 次迭代收敛。前几次迭代记录大致如下迭代次数 kV2 实部V2 虚部|V2| (pu)δ2 (°)ΔV01.00000.00001.00000.0000-10.8900-0.13000.8994-8.3081.70e-120.8531-0.11730.8612-7.8313.90e-230.8453-0.12430.8544-8.3651.05e-240.8426-0.12320.8516-8.3142.90e-350.8420-0.12380.8512-8.3708.50e-460.8418-0.12380.8508-8.3702.00e-4可以看到修正量每个迭代大约缩小 3 到 4 倍这就是线性收敛的特征。最终收敛结果约在 |V2| 0.8508 pu、δ2 -8.37°也就是 V2 ≈ 0.8418 - j0.1238 pu。4.2 收敛判据的两个流派我代码里用的是电压修正量判据也就是 max|ΔV| tol。这个判据直观教学演示够用。但实际工程潮流程序更常用的是功率失配量判据定义如下对每个 PQ 节点ΔP P_spec - real(V · conj(Σ Y_ij·V_j)) ΔQ Q_spec - imag(V · conj(Σ Y_ij·V_j))当所有节点的 max|ΔP| 和 max|ΔQ| 都小于容忍度比如 0.01 MW / 0.01 Mvar时判定收敛。功率失配量直接对应物理量物理意义清晰。为什么还要提两种判据因为存在极端情况电压修正量很小但功率失配量仍然偏大或者反过来。一个严谨的代码最好两种判据同时检查或者至少理解两者的差异。教学代码用电压修正量没问题但心里要清楚这不是工程标准。4.3 松弛因子怎么用、怎么翻车Gauss-Seidel 法收敛速度不够快时可以引入松弛因子V_new_acc V_old ω · (V_new - V_old)ω 1 叫超松弛SORω 1 叫欠松弛。在高压输电网中ω 取 1.2 到 1.7 通常能显著减少迭代次数。但我的实际经验是系统重载或线路 R/X 比较大时超松弛很容易让迭代震荡甚至发散。调试期我强烈建议先用 ω 1确认收敛后再慢慢往上调。你可以把上面的代码改成omega 1.3; % 松弛因子 V2_new_raw (1/Y22) * ( conj(S2)/conj(V2) - Y21*V1_phasor ); V2_new V2 omega * (V2_new_raw - V2);如果发现迭代曲线在某个值附近来回振荡第一个该怀疑的就是 ω 太大了改成 0.9 试试。松弛因子不是越大越好它和系统条件数、线路参数都有关系。4.4 发散时的五步排查我自己调试潮流代码时踩过不少坑发散时按这个顺序查基本能定位问题符号约定P、Q 方向定义是不是前后不一致。负荷写成正值但公式里又按注入推导符号必错Ybus 构建互导纳漏了负号或者自导纳忘了加并联导纳。两节点系统一眼能看出来多节点就麻烦了负荷是否超过系统极限重负荷下潮流可能根本无解任何初始值都会发散。可以把负荷按比例降低比如乘 0.9验证初始值平启动通常是安全的但极端情况下要换一个有物理依据的初值松弛因子把 ω 降回 1 或 0.9排除超松弛导致的振荡5. 收敛结果怎么验证功率平衡才是硬道理5.1 用线路功率反推验证潮流算完不能直接信必须验证。最直接的办法是拿收敛的电压值回代计算线路功率看节点 2 接收到的功率是否等于给定的负荷功率。代码里已经算好了几个关键量I12 (V1 - V2) / Z这是从母线 1 流向母线 2 的线路电流S12 V1 · conj(I12)母线 1 注入线路的功率S21 V2 · conj(I12)母线 2 从线路接收的功率用我上面手推的最终电压值粗算S21 ≈ 0.8003 j0.5986 pu和给定的负荷功率 0.8 j0.6 pu 非常接近手算精度导致的微小差异是正常的计算机跑出来会完全一致。这说明潮流计算结果满足功率平衡迭代没有问题。5.2 功率不平衡量落在哪里问题就在哪里如果 S21 和给定 P2 jQ2 对不上差值就指向了问题源头实部对不上检查有功数据或电压实部计算虚部对不上检查无功数据或电压虚部计算两个都对不上大概率符号约定错了另外还可以检查全网功率平衡S12 应该等于 S21 加上线路损耗。线路损耗等于 I12²·Z这是另一条独立的验证路径。5.3 结果的工程解读算出来的 V2 ≈ 0.85 pu比 V1 1.05 pu 低了约 0.2 pu。这个压降哪里来的就是负荷电流流过线路阻抗产生的。高压线路电抗远大于电阻所以无功分量贡献了主要压降有功分量贡献了相角差。δ2 ≈ -8.4° 说明有功功率从母线 1 流向母线 2——电力系统里相角超前的节点向滞后节点送出有功这是基本规律。顺便说一个有趣的现象如果把 P2、Q2 从负值改成正值也就是母线 2 变成电源比如分布式光伏接入那 V2 非但不会跌落反而会升高。这就是“分布式电源接入导致电压抬升”的直观解释。一个两节点系统把负荷改成电源你就能亲手看到电压曲线从下降变成上升这对理解新能源并网电压问题特别有帮助。6. 从两节点到 N 节点这套代码的扩展方向6.1 通用化改造点两节点代码跑通之后往多节点扩展是很自然的下一步。核心改造点有四个Ybus 要用完整的构建逻辑遍历所有节点和支路自导纳累加、互导纳取负有变压器时要引入变比修正节点数据要结构化用矩阵或结构体数组存储节点编号、类型、P、Q、V 设定值等迭代逻辑要区分节点类型平衡节点跳过PQ 节点按基本公式更新PV 节点每次迭代后要把电压幅值拉回设定值所谓“钳位”并重新计算无功 Q收敛判据改为全节点失配量判据通用化后的迭代核心伪代码如下for iter 1:maxIter for i 2:nBus % 平衡节点不参与迭代 if 节点i是PQ节点 V(i) (1/Ybus(i,i)) * ( conj(S(i))/conj(V(i)) - sum(Ybus(i,j)*V(j)) ); elseif 节点i是PV节点 % 先求Q(i)再迭代V(i)最后把幅值钳位到设定值 Q(i) -imag( conj(V(i)) * sum(Ybus(i,j)*V(j)) ); V(i) (1/Ybus(i,i)) * ( conj(S(i))/conj(V(i)) - sum(Ybus(i,j)*V(j)) ); V(i) V_set(i) * V(i) / abs(V(i)); % 幅值钳位 end end if max(abs(deltaV)) tol; break; end end注意其中 sum 函数要实现“已更新的节点用新值”因为 MATLAB 的向量化 sum 没法区分新旧需要用循环或者循环内逐个累加。这是 Gauss-Seidel 和多节点扩展最容易出 bug 的地方。6.2 节点编号顺序会影响迭代次数这是 Gauss-Seidel 法一个很容易被忽略的特性节点编号顺序影响收敛速度。迭代时先更新的节点其新值会参与后续节点计算因此把强耦合节点排在前面能加快收敛。实际工程里常见的做法是按“PQ 节点在前、PV 节点在后”的顺序编号但这只是经验规则不是严格最优。6.3 Gauss-Seidel 法在今天的位置现代电力系统潮流计算的主流是牛顿-拉夫逊法和快速解耦法Gauss-Seidel 在大规模输电网里确实不够快。但它在几个场景里依然有价值教学场景代码短、逻辑直观是理解潮流计算的最佳入口小型系统节点数几十个以下时G-S 的迭代次数可以接受配电网计算配电网辐射状结构下常用的前推回代法本质上和 G-S 有血缘关系这个两节点系统是整个潮流计算领域的最小完整示例。我在实际教学和工程调试中的体会是把这段代码玩透亲手改变负荷、改阻抗、改松弛因子、观察收敛曲线的变化对潮流计算的理解会远超看十遍推导。后面再学牛顿法、快速解耦法、最优潮流都建立在同样的物理和数学直觉上。最后分享一个小技巧把代码里的 tol 改成 1e-4迭代次数会少两三次教学演示时实时感更强等你需要精细结果时再改回 1e-7。然后试试把松弛因子设成 1.4很多情况下迭代次数能减少接近一半但在重负荷场景下也会直接震荡给你看。这种“亲手让它发散一次再救回来”的经历比任何理论分析都更能帮你建立对潮流算法的直觉。
第一次接触潮流计算时我对着两节点电路图犯过嘀咕这不过是一个电源、一条线路、一个负荷的简单电路欧姆定律加基尔霍夫电压定律分分钟就能列方程为什么还要专门搞一套 Gauss-Seidel 迭代直到自己动手在 MATLAB 里把代码写出来才意识到哪怕只有两个节点母线 2 的电压也会顽固地出现在方程右边的分母里而且是以共轭形式出现。这是个彻头彻尾的非线性问题手算解析解不是不行但换到几十上百个节点的真实电网就彻底没戏。这篇文章就从两节点系统入手讲清楚 Gauss-Seidel 潮流计算的完整建模、MATLAB 实现、收敛过程以及我实际调试中踩过的一堆坑。1. 两节点也要迭代先想清楚潮流问题到底在求什么1.1 一个看起来“很简单”的两节点系统先把这个两节点系统的算例参数摆出来母线 1平衡节点可以理解成大电网或者发电机电压幅值固定为 V1 1.05 pu相角为 0°母线 2PQ 节点负荷吸收有功 0.8 pu、无功 0.6 pu也就是注入功率 P2 -0.8 pu、Q2 -0.6 pu线路阻抗Z 0.05 j0.2 pu电抗远大于电阻符合高压架空线的典型特征如果按照中学物理的思路电源电动势 1.05线路阻抗已知负载功率已知似乎应该“一步就能算出来”。但这里有个关键区别电力系统里的负荷模型是恒功率不是恒阻抗。恒阻抗负载可以用欧姆定律直接算但恒功率负荷的意思是不管母线 2 电压降到多少它消耗的有功和无功都保持 0.8 j0.6 pu。这就把问题变成了一个耦合的非线性方程——我们要求解的是什么样的 V2 才能让节点 2 的注入功率恰好等于给定的 P2 jQ21.2 为什么解析法在这里卡壳把节点 2 的功率方程写出来就明白了。节点注入电流方向取“流入网络为正”那么从节点 1 流向节点 2 的线路电流是 (V1 - V2)/Z而节点 2 的注入电流是它的相反数I2 (V2 - V1)/Z又因为复功率的定义是 S V · conj(I)所以S2 V2 · conj(I2) V2 · conj((V2 - V1)/Z)把 V1 视为实数 1.05展开S2 V2 · [conj(V2) - V1] / conj(Z) (|V2|² - V1·conj(V2)) / conj(Z)这个式子同时包含 V2 和 conj(V2)也就是同时包含实部虚部两个独立变量拆开写就是两个耦合的二元二次方程。手算两节点系统还能硬解但这个方法不可扩展——真实电网动辄几百个节点每个节点的功率方程都和其他所有节点电压耦合解析法在规模面前毫无意义。1.3 Gauss-Seidel 的思路把方程变成迭代Gauss-Seidel 法做的事情很简单把节点电压方程改写成不动点迭代格式。对 PQ 节点 i核心迭代公式是V_i^(k1) (1/Y_ii) · [conj(S_i)/conj(V_i^(k)) - Σ(Y_ij · V_j)]其中已经更新过的节点用新值没更新的用旧值——这就是“Gauss-Seidel”和“Jacobi 法”的本质区别。每次用当前估计值算出新的电压估计反复迭代直到修正量足够小。两节点系统是这条思路的最小实现但已经包含全部核心要素导纳矩阵、共轭复功率、节点电压耦合。2. 节点三分法与两节点系统的建模细节2.1 PQ、PV、平衡节点各自要什么、给什么潮流计算里节点不是随便分类的每种节点对应一组已知量和待求量背后的工程含义完全不同节点类型已知量待求量工程角色PQ 节点有功 P、无功 Q电压幅值 V、相角 δ负荷节点、恒功率负荷绝大多数节点属于此类PV 节点有功 P、电压幅值 V无功 Q、相角 δ发电机节点有自动电压调节器AVR能维持机端电压平衡节点电压幅值 V、相角 δ有功 P、无功 Q参考节点承担全网功率不平衡量相位基准为什么必须有一个平衡节点因为全网功率损耗在潮流计算之前是未知的总要有一个节点吸收“发电总功率 - 负荷总功率 - 网损”这个差值。另外电压相角是相对量所有母线相角需要一个基准平衡节点就是这个基准通常取 V 1.0∠0° 或 1.05∠0°。2.2 本算例的参数设定与标幺值为了让大家直接复现我把参数统一用标幺值pu表示。本算例的具体设定母线 1V1 1.05 puθ1 0°类型为平衡节点母线 2P2 -0.8 puQ2 -0.6 pu类型为 PQ 节点负号表示吸收功率线路阻抗Z 0.05 j0.2 pu如果用 100 MVA 作基准容量那么 0.8 pu 的有功就是 80 MW0.6 pu 的无功就是 60 MvarV 1.05 pu 在 220 kV 系统里对应 231 kV。标幺值的意义在于把电压、功率、阻抗都归一化到 1 附近数值稳定性好而且三相系统用线电压标幺、三相功率标幺之后可以直接用单相等效电路建模非常方便。2.3 符号约定是最大的坑我在给研究生改代码时发现Gauss-Seidel 潮流代码跑不出来的原因十有八九是符号约定乱了。这里必须提前说清楚本代码中 P、Q 的正方向定义为“注入节点”方向所以负荷吸收功率时 P、Q 取负值复功率定义S V · conj(I)因此节点注入电流 I conj(S) / conj(V)注意是 conj(S) 除以 conj(V)不是 S 除以 V。这个错误会让迭代结果完全错误甚至发散有的教材把 P、Q 正方向定义为“流出节点”即负荷方向那公式里的符号会整体变号。写代码前先定好约定然后所有公式保持一致。我见过有人前五分钟用注入方向定义 S后五分钟又从负荷方向抄公式结果自然是灾难。2.4 节点导纳矩阵 Ybus 的构建两节点系统忽略对地电容后导纳矩阵非常简洁Ybus [1/Z, -1/Z; -1/Z, 1/Z]对角线元素是自导纳等于与该节点相连的所有支路导纳之和非对角线元素是互导纳等于连接两节点的支路导纳的负值。这个结论在没有变压器和移相器时成立。以后扩展到多节点系统Ybus 的构建规律是一样的自导纳累加、互导纳取负。3. Gauss-Seidel 迭代公式推导与 MATLAB 实现3.1 从节点电压方程到迭代格式对于节点 2节点电压方程是I2 Y21·V1 Y22·V2把 I2 conj(S2) / conj(V2) 代入conj(S2)/conj(V2) Y21·V1 Y22·V2移项得到V2 (1/Y22) · [conj(S2)/conj(V2) - Y21·V1]这就是 PQ 节点的 Gauss-Seidel 迭代式。注意 V1 是平衡节点电压恒定不参与迭代但如果系统里有多个 PQ 节点那么每个新算出的电压都要立刻用于后续节点计算这正是 Gauss-Seidel 和 Jacobi 的区别。3.2 完整可运行的 MATLAB 代码下面给出我调试过的完整代码。注释写得比较详细可以粘贴到 MATLAB 里直接运行%% 两节点电力系统 Gauss-Seidel 潮流计算 % 母线1平衡节点V1 已知恒定 % 母线2PQ节点P2、Q2 已知求电压幅值和相角 % 所有量均为标幺值pu clear; clc; close all; %% ---------- 1. 系统参数 ---------- % 线路阻抗标幺值架空线电抗远大于电阻 Z 0.05 1j*0.2; % 节点导纳矩阵 Ybus2x2忽略对地导纳 Y11 1/Z; Y12 -1/Z; Y21 -1/Z; Y22 1/Z; Ybus [Y11 Y12; Y21 Y22]; %% ---------- 2. 节点数据 ---------- % 平衡节点节点1 V1 1.05; % 幅值 theta1 0; % 相角rad V1_phasor V1 * exp(1j*theta1); % PQ节点节点2注入节点方向为正负荷为负 P2 -0.8; % 吸收有功 0.8 pu Q2 -0.6; % 吸收无功 0.6 pu S2 P2 1j*Q2; % 节点2的注入复功率 %% ---------- 3. 初始化 ---------- V2 1.0 1j*0.0; % 平启动flat start tol 1e-7; % 电压修正量收敛容忍度 maxIter 50; % 最大迭代次数 fprintf( 两节点系统 Gauss-Seidel 潮流计算 \n); fprintf(线路阻抗 Z %.3f j%.3f pu\n, real(Z), imag(Z)); fprintf(平衡节点 V1 %.2f∠%.1f°\n, abs(V1_phasor), angle(V1_phasor)*180/pi); fprintf(PQ节点 P2 %.2f pu, Q2 %.2f pu\n\n, P2, Q2); % 记录每次迭代的电压便于画收敛曲线 V2_history zeros(maxIter1,1); V2_history(1) V2; %% ---------- 4. Gauss-Seidel 主迭代 ---------- for k 1:maxIter % 核心迭代公式 % V2_new (1/Y22) * [ conj(S2)/conj(V2) - Y21*V1 ] % 其中 conj(S2) P2 - j*Q2 V2_new (1/Y22) * ( conj(S2)/conj(V2) - Y21*V1_phasor ); dV abs(V2_new - V2); V2 V2_new; V2_history(k1) V2; fprintf(k%2d | V2%.8f %.8fj | |V2|%.8f | δ2%9.6f° | ΔV%.2e\n, ... k, real(V2), imag(V2), abs(V2), angle(V2)*180/pi, dV); if dV tol fprintf(\n第 %d 次迭代收敛ΔV %.0e\n\n, k, tol); break; end if k maxIter error(达到最大迭代次数仍未收敛检查参数或减小负荷); end end %% ---------- 5. 结果输出与验证 ---------- fprintf( 计算结果 \n); fprintf(母线2电压幅值 |V2| %.6f pu\n, abs(V2)); fprintf(母线2电压相角 δ2 %.6f°\n, angle(V2)*180/pi); fprintf(母线2电压相量 V2 %.6f ∠ %.6f° pu\n\n, abs(V2), angle(V2)*180/pi); % 线路电流正方向 1-2及相关功率 I12 (V1_phasor - V2) / Z; S12 V1_phasor * conj(I12); % 节点1注入线路的功率 S21 V2 * conj(I12); % 节点2从线路接收的功率 fprintf(线路电流 I12 %.6f ∠ %.6f° pu\n, abs(I12), angle(I12)*180/pi); fprintf(节点1注入线路 S12 %.6f %.6fj pu\n, real(S12), imag(S12)); fprintf(节点2接收线路 S21 %.6f %.6fj pu\n, real(S21), imag(S21)); fprintf(验证S21 应约等于负荷功率 (0.8j0.6) pu\n); %% ---------- 6. 收敛曲线 ---------- iterUsed k; % 实际迭代次数 figure(Color,w); subplot(1,2,1); plot(0:iterUsed, abs(V2_history(1:iterUsed1)), b-o, LineWidth,1.2, MarkerSize,5); xlabel(迭代次数 k); ylabel(|V_2| (pu)); title(PQ节点电压幅值收敛过程); grid on; subplot(1,2,2); plot(0:iterUsed, angle(V2_history(1:iterUsed1))*180/pi, r-s, LineWidth,1.2, MarkerSize,5); xlabel(迭代次数 k); ylabel(δ_2 (deg)); title(PQ节点电压相角收敛过程); grid on;3.3 单次迭代的逐步解读用第一次迭代演示这个公式到底在算什么。初始值 V2 1 j0conj(S2) -0.8 j0.6conj(S2) / conj(V2) -0.8 j0.6因为 conj(V2)1-Y21·V1 -(-1/Z)·1.05 1.05/(0.05j0.2)这是由平衡节点通过互导纳注入耦合的部分括号内求和后乘以 1/Y22 Z最终得到 V2^(1) ≈ 0.8900 - j0.1300即第一次迭代就从 1.0∠0° 跳到了约 0.899∠-8.3°。大跨步一步到位之后每次修正量迅速衰减。这就是 Gauss-Seidel 法在简单系统中的典型表现——开始几步下降很快越往后越慢呈线性收敛。4. 收敛过程实测初始值、判据与松弛因子的那些事4.1 我的机器上跑出来的迭代记录在 MATLAB R2023b 下运行收敛容忍度设为 1e-7 时大约 9 到 10 次迭代收敛。前几次迭代记录大致如下迭代次数 kV2 实部V2 虚部|V2| (pu)δ2 (°)ΔV01.00000.00001.00000.0000-10.8900-0.13000.8994-8.3081.70e-120.8531-0.11730.8612-7.8313.90e-230.8453-0.12430.8544-8.3651.05e-240.8426-0.12320.8516-8.3142.90e-350.8420-0.12380.8512-8.3708.50e-460.8418-0.12380.8508-8.3702.00e-4可以看到修正量每个迭代大约缩小 3 到 4 倍这就是线性收敛的特征。最终收敛结果约在 |V2| 0.8508 pu、δ2 -8.37°也就是 V2 ≈ 0.8418 - j0.1238 pu。4.2 收敛判据的两个流派我代码里用的是电压修正量判据也就是 max|ΔV| tol。这个判据直观教学演示够用。但实际工程潮流程序更常用的是功率失配量判据定义如下对每个 PQ 节点ΔP P_spec - real(V · conj(Σ Y_ij·V_j)) ΔQ Q_spec - imag(V · conj(Σ Y_ij·V_j))当所有节点的 max|ΔP| 和 max|ΔQ| 都小于容忍度比如 0.01 MW / 0.01 Mvar时判定收敛。功率失配量直接对应物理量物理意义清晰。为什么还要提两种判据因为存在极端情况电压修正量很小但功率失配量仍然偏大或者反过来。一个严谨的代码最好两种判据同时检查或者至少理解两者的差异。教学代码用电压修正量没问题但心里要清楚这不是工程标准。4.3 松弛因子怎么用、怎么翻车Gauss-Seidel 法收敛速度不够快时可以引入松弛因子V_new_acc V_old ω · (V_new - V_old)ω 1 叫超松弛SORω 1 叫欠松弛。在高压输电网中ω 取 1.2 到 1.7 通常能显著减少迭代次数。但我的实际经验是系统重载或线路 R/X 比较大时超松弛很容易让迭代震荡甚至发散。调试期我强烈建议先用 ω 1确认收敛后再慢慢往上调。你可以把上面的代码改成omega 1.3; % 松弛因子 V2_new_raw (1/Y22) * ( conj(S2)/conj(V2) - Y21*V1_phasor ); V2_new V2 omega * (V2_new_raw - V2);如果发现迭代曲线在某个值附近来回振荡第一个该怀疑的就是 ω 太大了改成 0.9 试试。松弛因子不是越大越好它和系统条件数、线路参数都有关系。4.4 发散时的五步排查我自己调试潮流代码时踩过不少坑发散时按这个顺序查基本能定位问题符号约定P、Q 方向定义是不是前后不一致。负荷写成正值但公式里又按注入推导符号必错Ybus 构建互导纳漏了负号或者自导纳忘了加并联导纳。两节点系统一眼能看出来多节点就麻烦了负荷是否超过系统极限重负荷下潮流可能根本无解任何初始值都会发散。可以把负荷按比例降低比如乘 0.9验证初始值平启动通常是安全的但极端情况下要换一个有物理依据的初值松弛因子把 ω 降回 1 或 0.9排除超松弛导致的振荡5. 收敛结果怎么验证功率平衡才是硬道理5.1 用线路功率反推验证潮流算完不能直接信必须验证。最直接的办法是拿收敛的电压值回代计算线路功率看节点 2 接收到的功率是否等于给定的负荷功率。代码里已经算好了几个关键量I12 (V1 - V2) / Z这是从母线 1 流向母线 2 的线路电流S12 V1 · conj(I12)母线 1 注入线路的功率S21 V2 · conj(I12)母线 2 从线路接收的功率用我上面手推的最终电压值粗算S21 ≈ 0.8003 j0.5986 pu和给定的负荷功率 0.8 j0.6 pu 非常接近手算精度导致的微小差异是正常的计算机跑出来会完全一致。这说明潮流计算结果满足功率平衡迭代没有问题。5.2 功率不平衡量落在哪里问题就在哪里如果 S21 和给定 P2 jQ2 对不上差值就指向了问题源头实部对不上检查有功数据或电压实部计算虚部对不上检查无功数据或电压虚部计算两个都对不上大概率符号约定错了另外还可以检查全网功率平衡S12 应该等于 S21 加上线路损耗。线路损耗等于 I12²·Z这是另一条独立的验证路径。5.3 结果的工程解读算出来的 V2 ≈ 0.85 pu比 V1 1.05 pu 低了约 0.2 pu。这个压降哪里来的就是负荷电流流过线路阻抗产生的。高压线路电抗远大于电阻所以无功分量贡献了主要压降有功分量贡献了相角差。δ2 ≈ -8.4° 说明有功功率从母线 1 流向母线 2——电力系统里相角超前的节点向滞后节点送出有功这是基本规律。顺便说一个有趣的现象如果把 P2、Q2 从负值改成正值也就是母线 2 变成电源比如分布式光伏接入那 V2 非但不会跌落反而会升高。这就是“分布式电源接入导致电压抬升”的直观解释。一个两节点系统把负荷改成电源你就能亲手看到电压曲线从下降变成上升这对理解新能源并网电压问题特别有帮助。6. 从两节点到 N 节点这套代码的扩展方向6.1 通用化改造点两节点代码跑通之后往多节点扩展是很自然的下一步。核心改造点有四个Ybus 要用完整的构建逻辑遍历所有节点和支路自导纳累加、互导纳取负有变压器时要引入变比修正节点数据要结构化用矩阵或结构体数组存储节点编号、类型、P、Q、V 设定值等迭代逻辑要区分节点类型平衡节点跳过PQ 节点按基本公式更新PV 节点每次迭代后要把电压幅值拉回设定值所谓“钳位”并重新计算无功 Q收敛判据改为全节点失配量判据通用化后的迭代核心伪代码如下for iter 1:maxIter for i 2:nBus % 平衡节点不参与迭代 if 节点i是PQ节点 V(i) (1/Ybus(i,i)) * ( conj(S(i))/conj(V(i)) - sum(Ybus(i,j)*V(j)) ); elseif 节点i是PV节点 % 先求Q(i)再迭代V(i)最后把幅值钳位到设定值 Q(i) -imag( conj(V(i)) * sum(Ybus(i,j)*V(j)) ); V(i) (1/Ybus(i,i)) * ( conj(S(i))/conj(V(i)) - sum(Ybus(i,j)*V(j)) ); V(i) V_set(i) * V(i) / abs(V(i)); % 幅值钳位 end end if max(abs(deltaV)) tol; break; end end注意其中 sum 函数要实现“已更新的节点用新值”因为 MATLAB 的向量化 sum 没法区分新旧需要用循环或者循环内逐个累加。这是 Gauss-Seidel 和多节点扩展最容易出 bug 的地方。6.2 节点编号顺序会影响迭代次数这是 Gauss-Seidel 法一个很容易被忽略的特性节点编号顺序影响收敛速度。迭代时先更新的节点其新值会参与后续节点计算因此把强耦合节点排在前面能加快收敛。实际工程里常见的做法是按“PQ 节点在前、PV 节点在后”的顺序编号但这只是经验规则不是严格最优。6.3 Gauss-Seidel 法在今天的位置现代电力系统潮流计算的主流是牛顿-拉夫逊法和快速解耦法Gauss-Seidel 在大规模输电网里确实不够快。但它在几个场景里依然有价值教学场景代码短、逻辑直观是理解潮流计算的最佳入口小型系统节点数几十个以下时G-S 的迭代次数可以接受配电网计算配电网辐射状结构下常用的前推回代法本质上和 G-S 有血缘关系这个两节点系统是整个潮流计算领域的最小完整示例。我在实际教学和工程调试中的体会是把这段代码玩透亲手改变负荷、改阻抗、改松弛因子、观察收敛曲线的变化对潮流计算的理解会远超看十遍推导。后面再学牛顿法、快速解耦法、最优潮流都建立在同样的物理和数学直觉上。最后分享一个小技巧把代码里的 tol 改成 1e-4迭代次数会少两三次教学演示时实时感更强等你需要精细结果时再改回 1e-7。然后试试把松弛因子设成 1.4很多情况下迭代次数能减少接近一半但在重负荷场景下也会直接震荡给你看。这种“亲手让它发散一次再救回来”的经历比任何理论分析都更能帮你建立对潮流算法的直觉。
