第一次接触潮流计算时，我对着两节点电路图犯过嘀咕：这不过是一个电源、一条线路、一个负荷的简单电路，欧姆定律加基尔霍夫电压定律分分钟就能列方程，为什么还要专门搞一套 Gauss-Seidel 迭代？直到自己动手在 MATLAB 里把代码写出来，才意识到哪怕只有两个节点，母线 2 的电压也会顽固地出现在方程右边的分母里，而且是以共轭形式出现。这是个彻头彻尾的非线性问题，手算解析解不是不行，但换到几十上百个节点的真实电网就彻底没戏。这篇文章就从两节点系统入手，讲清楚 Gauss-Seidel 潮流计算的完整建模、MATLAB 实现、收敛过程，以及我实际调试中踩过的一堆坑。
1. 两节点也要迭代？先想清楚潮流问题到底在求什么
1.1 一个看起来“很简单”的两节点系统
先把这个两节点系统的算例参数摆出来：
- 母线 1：平衡节点（可以理解成大电网或者发电机），电压幅值固定为 V1 = 1.05 pu，相角为 0°
- 母线 2：PQ 节点，负荷吸收有功 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 + jQ2？
1.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^(k+1) = (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）表示。本算例的具体设定：
- 母线 1：V1 = 1.05 pu，θ1 = 0°，类型为平衡节点
- 母线 2：P2 = -0.8 pu，Q2 = -0.6 pu，类型为 PQ 节点（负号表示吸收功率）
- 线路阻抗：Z = 0.05 + j0.2 pu
如果用 100 MVA 作基准容量，那么 0.8 pu 的有功就是 80 MW，0.6 pu 的无功就是 60 Mvar，V = 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 已知，恒定） % 母线2：PQ节点（P2、Q2 已知，求电压幅值和相角） % 所有量均为标幺值（pu） clear; clc; close all; %% ---------- 1. 系统参数 ---------- % 线路阻抗（标幺值）：架空线电抗远大于电阻 Z = 0.05 + 1j*0.2; % 节点导纳矩阵 Ybus（2x2，忽略对地导纳） 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(maxIter+1,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(k+1) = 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.8+j0.6) pu\n'); %% ---------- 6. 收敛曲线 ---------- iterUsed = k; % 实际迭代次数 figure('Color','w'); subplot(1,2,1); plot(0:iterUsed, abs(V2_history(1:iterUsed+1)), '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:iterUsed+1))*180/pi, 'r-s', 'LineWidth',1.2, 'MarkerSize',5); xlabel('迭代次数 k'); ylabel('δ_2 (deg)'); title('PQ节点电压相角收敛过程'); grid on;3.3 单次迭代的逐步解读
用第一次迭代演示这个公式到底在算什么。初始值 V2 = 1 + j0：
- conj(S2) = -0.8 + j0.6
- conj(S2) / conj(V2) = -0.8 + j0.6（因为 conj(V2)=1）
- -Y21·V1 = -(-1/Z)·1.05 = 1.05/(0.05+j0.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 次迭代收敛。前几次迭代记录大致如下：
| 迭代次数 k | V2 实部 | V2 虚部 | |V2| (pu) | δ2 (°) | ΔV |
|---|---|---|---|---|---|
| 0 | 1.0000 | 0.0000 | 1.0000 | 0.0000 | - |
| 1 | 0.8900 | -0.1300 | 0.8994 | -8.308 | 1.70e-1 |
| 2 | 0.8531 | -0.1173 | 0.8612 | -7.831 | 3.90e-2 |
| 3 | 0.8453 | -0.1243 | 0.8544 | -8.365 | 1.05e-2 |
| 4 | 0.8426 | -0.1232 | 0.8516 | -8.314 | 2.90e-3 |
| 5 | 0.8420 | -0.1238 | 0.8512 | -8.370 | 8.50e-4 |
| 6 | 0.8418 | -0.1238 | 0.8508 | -8.370 | 2.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，很多情况下迭代次数能减少接近一半，但在重负荷场景下也会直接震荡给你看。这种“亲手让它发散一次再救回来”的经历，比任何理论分析都更能帮你建立对潮流算法的直觉。
