刚接触潮流计算时最卡我的不是功率方程本身而是对着书本上一堆迭代公式不知道从哪里下手。后来我把问题缩到最小的三节点系统把高斯-赛德尔法和牛顿-拉夫森法从头到尾各写了一遍再用Matlab跑通、对比、踩坑才真正把这套东西吃透。这篇文章就围绕三节点系统把这两种经典潮流算法的原理、迭代公式、Matlab实现和实测结果完整串一遍。三节点虽小却是理解多节点潮流程序的“最小可运行系统”——Ybus怎么形成、节点类型怎么划分、收敛判据怎么定、雅可比矩阵怎么构造这些核心问题在三个节点上都能讲明白而且每一步都能自己手算验证。1. 三节点算例搭建从节点类型到导纳矩阵1.1 潮流分析本质上在算什么说穿了潮流分析就是给定电网的拓扑结构、支路阻抗参数、各节点的发电功率和负荷功率求全网的复电压分布然后根据电压再推支路功率和网络损耗。它不关心暂态过程只看稳态工况是电力系统里最常用的一张“算账”工具。打个比方就像给一套供水管网算水压和流量每户取水量已知水泵出口压力已知要算出每个节点的压力、每段管道的流量。电网里的“水压”就是节点电压幅值和相角“每户取水”就是节点注入有功和无功功率。不同的是电网的功率和电压之间是非线性关系所以不能一次解出来只能迭代逼近。1.2 三节点系统方案与节点类型划分我用的算例是一个经典的三节点三支路系统基准容量取100MVA。节点1设为平衡节点Slack电压固定为1.0∠0°负责吸收全网功率不平衡量节点2和节点3设为PQ节点负荷功率已知需要迭代求解电压幅值和相角。节点参数如下表节点节点类型注入有功P(pu)注入无功Q(pu)初始电压1平衡节点待求待求1.0∠0°2PQ节点-0.40-0.201.0∠0°3PQ节点-0.30-0.151.0∠0°注意这里负荷功率用负号表示因为潮流计算里的注入功率以流入网络为正方向。也就是说节点2和节点3是在从系统“取用”功率这在程序里非常容易搞反。支路参数如下支路电阻R(pu)电抗X(pu)对地导纳B/2(pu)1-20.020.0601-30.030.0902-30.0250.0750为了突出输电线路特性电阻电抗比R/X取1/3这是个比较典型的值。六氟化硫断路器、变压器等设备的阻抗参数在标幺值下也基本在这个量级。1.3 Ybus矩阵的构建方法与Matlab实现导纳矩阵是后续所有计算的地基。它的规则非常固定自导纳Y_ii是连接到节点i的所有支路导纳之和互导纳Y_ij是连接节点i和节点j的支路导纳取负号。Matlab里形成三节点Ybus的代码非常简单% 三节点系统导纳矩阵构建 Y zeros(3,3); % 支路数据起节点 终节点 电阻 电抗 branch [1 2 0.02 0.06; 1 3 0.03 0.09; 2 3 0.025 0.075]; for k 1:3 i branch(k,1); j branch(k,2); yk 1/(branch(k,3) 1j*branch(k,4)); Y(i,i) Y(i,i) yk; Y(j,j) Y(j,j) yk; Y(i,j) Y(i,j) - yk; Y(j,i) Y(j,i) - yk; end disp(Y);计算得到Y 8.3333 - 25.0000i -5.0000 15.0000i -3.3333 10.0000i -5.0000 15.0000i 9.0000 - 27.0000i -4.0000 12.0000i -3.3333 10.0000i -4.0000 12.0000i 7.3333 - 22.0000i这个矩阵有两个特征值得记住第一它是对称矩阵第二对角元素明显大于非对角元素。这两点在后面构建雅可比矩阵时也有对应的体现。2. 高斯-赛德尔法从功率平衡到逐点更新2.1 为什么GS法适合作为入门第一个潮流算法高斯-赛德尔法是求解线性方程组最经典的迭代法之一把它用到潮流分析里本质上是在反复利用节点电压方程和节点功率方程互相修正。它的优势是原理直观、编程极简单、内存占用小而且对迭代初值不敏感。在计算节点i时GS法会立刻使用已经更新过的节点1到i-1的新电压值这就是“逐点更新”。这一点和雅可比法不同——雅可比法必须等一轮全部算完才统一更新GS法省一半存储且收敛快一些。2.2 迭代公式的推导过程基本出发点还是电路理论里的节点电压方程。对任意节点i有I_i Σ Y_ij * V_j同时又知道节点注入复功率S_i P_i jQ_i满足S_i V_i * conj(I_i)把两个式子联立解出V_iV_i^(k1) (1/Y_ii) * [ (P_i - jQ_i) / conj(V_i^(k)) - Σ_{j≠i} Y_ij * V_j ]这里有个非常容易写错的地方等号右边分母上的V_i要取共轭而且用的是当前迭代点也就是最新值的共轭。很多人第一次写代码时很容易把conj(S(i)/V(i))和conj(S(i))/conj(V(i))搞混前者是错的后者才等于(P_i - jQ_i)/conj(V_i)。GS法的核心还在于等号右边求和项里的V_j取值规则当j i时用本轮已经更新过的新值当j i时用上一轮的旧值。Matlab的for循环天然满足这个规则因为V向量是逐个覆盖更新的。2.3 Matlab实现骨架% 高斯-赛德尔法潮流计算 % 输入Ybus矩阵Y、节点注入功率S、平衡节点编号、收敛精度 V ones(3,1); % 电压初始值 V(1) 1.0 0j; % 平衡节点固定 S [0; -0.4-0.2j; -0.3-0.15j]; % 节点注入功率 tol 1e-6; max_iter 100; for iter 1:max_iter V_old V; for i 2:3 % 平衡节点不参与迭代 sum_yv 0; for j 1:3 if j ~ i sum_yv sum_yv Y(i,j) * V(j); end end V(i) (conj(S(i)/V(i)) - sum_yv) / Y(i,i); end if max(abs(V - V_old)) tol fprintf(GS法迭代%d次收敛\n, iter); break; end end这段代码跑通后可以在命令行里输出每次迭代的V变化。你会发现电压下降的节奏是均匀的每次迭代只往前走一小步这正是GS法线性收敛的直观体现。2.4 GS法的收敛特性和边界条件GS法容易实现但千万别对它抱太高期望。它的收敛速度是线性的通俗说就是误差每一步只按一个固定比例比如0.8缩小。三节点系统迭代十几次能收敛到1e-6如果换成一个几十节点的输电网收敛速度会明显拖慢。另外GS法在遇到较重的负荷比如节点2的负荷加大到-1.0pu时迭代次数会急剧增加甚至可能不收敛。工程上我一般拿它来跑配电网、小系统或者给其他算法提供一个粗糙的初值而不是指望它在大型输电网里跑得多快。如果系统中存在PV节点发电机节点电压幅值固定、有功给定GS法每次迭代后还需要根据Q_i -Im(V_i * conj(I_i))推算无功然后修正V_i的幅值回到给定值。这个处理在三节点算例里没有体现但一旦扩展到IEEE标准节点就会遇到。3. 牛顿-拉夫森法雅可比矩阵与修正方程3.1 NR法的核心思想牛顿-拉夫森法不是从“迭代解线性方程”的思路出发而是把潮流问题直接看成一堆非线性方程的求根问题。对每一个节点都写出一组方程然后把它们在某一点做一阶泰勒展开解出修正量反复迭代。用通俗的话说GS法是一次一次试探着往正确方向走NR法是先算一下当前误差有多大、斜率有多陡然后根据斜率和误差直接跨出一大步。所以它的收敛速度远快于GS法是二次收敛误差每一步大致变成上一步的平方。3.2 潮流方程的极坐标形式与失配量在极坐标形式下节点i的有功、无功计算值为P_i_calc V_i * Σ [ V_j * (G_ij * cos(θ_i - θ_j) B_ij * sin(θ_i - θ_j)) ] Q_i_calc V_i * Σ [ V_j * (G_ij * sin(θ_i - θ_j) - B_ij * cos(θ_i - θ_j)) ]其中G_ij和B_ij分别是导纳矩阵元素的实部电导和虚部电纳。节点i的失配量定义为给定值与计算值之差ΔP_i P_i_sch - P_i_calc ΔQ_i Q_i_sch - Q_i_calc潮流方程有解等价于所有节点的失配量全部归零。3.3 雅可比矩阵的构造公式雅可比矩阵分成四块分别是有功对相角、有功对电压、无功对相角、无功对电压的偏导。我采用修正量为ΔV/V的形式也就是把电压幅值的相对变化量作为未知量。分块i≠j非对角ij对角H ∂P/∂θV_i V_j (G_ij sinθ_ij - B_ij cosθ_ij)-Q_i - B_ii V_i²N V·∂P/∂VV_i V_j (G_ij cosθ_ij B_ij sinθ_ij)P_i G_ii V_i²K ∂Q/∂θ-V_i V_j (G_ij cosθ_ij B_ij sinθ_ij)P_i - G_ii V_i²L V·∂Q/∂VV_i V_j (G_ij sinθ_ij - B_ij cosθ_ij)Q_i - B_ii V_i²这里的θ_ij θ_i - θ_jP_i和Q_i是当前迭代点下计算的节点注入功率不是给定值。初次接触时很容易把P_i、Q_i误写成给定功率实际上公式里的P_i、Q_i都来自当前迭代点的计算值这个细节会导致雅可比矩阵错误。修正方程写成[ ΔP ] [ H N ] [ Δθ ] [ ΔQ ] [ K L ] * [ ΔV/V ]3.4 NR法的Matlab实现骨架% 牛顿-拉夫森法潮流计算 V ones(3,1); theta zeros(3,1); V(1) 1.0; theta(1) 0; S [0; -0.4-0.2j; -0.3-0.15j]; G real(Y); B imag(Y); tol 1e-6; for iter 1:20 % 计算P_calc和Q_calc P_calc zeros(3,1); Q_calc zeros(3,1); for i 1:3 for j 1:3 th_ij theta(i) - theta(j); P_calc(i) P_calc(i) V(i)*V(j)*(G(i,j)*cos(th_ij) B(i,j)*sin(th_ij)); Q_calc(i) Q_calc(i) V(i)*V(j)*(G(i,j)*sin(th_ij) - B(i,j)*cos(th_ij)); end end dP real(S) - P_calc; dQ imag(S) - Q_calc; % 注意平衡节点的偏差通常不小但不参与修正 % 只取PQ节点节点2、3 dP2 dP(2:3); dQ2 dQ(2:3); F [dP2; dQ2]; if max(abs(F)) tol fprintf(NR法迭代%d次收敛\n, iter); break; end % 构建雅可比矩阵J4x4 % 这里为简洁直接按4个PQ变量展开实际工程用稀疏矩阵 J zeros(4,4); pq [2 3]; for a 1:2 i pq(a); for b 1:2 j pq(b); th_ij theta(i) - theta(j); if i j J(a,b) -Q_calc(i) - B(i,i)*V(i)^2; % H J(a,b2) P_calc(i) G(i,i)*V(i)^2; % N J(a2,b) P_calc(i) - G(i,i)*V(i)^2; % K J(a2,b2) Q_calc(i) - B(i,i)*V(i)^2; % L else J(a,b) V(i)*V(j)*(G(i,j)*sin(th_ij) - B(i,j)*cos(th_ij)); % H J(a,b2) V(i)*V(j)*(G(i,j)*cos(th_ij) B(i,j)*sin(th_ij)); % N J(a2,b) -V(i)*V(j)*(G(i,j)*cos(th_ij) B(i,j)*sin(th_ij)); % K J(a2,b2) V(i)*V(j)*(G(i,j)*sin(th_ij) - B(i,j)*cos(th_ij)); % L end end end % 解修正方程并更新 dx J \ (-F); dtheta dx(1:2); dV_rel dx(3:4); theta(pq) theta(pq) dtheta; V(pq) V(pq) .* (1 dV_rel); end有个很重要的实现细节解修正方程时一定用左除\不要写inv(J) * (-F)。对于3节点系统二者速度没区别但扩展到几百节点时inv既慢又不稳左除会自动选择合适的稀疏求解器。3.5 NR法收敛特性NR法的二次收敛特性很典型第一次迭代失配量大约在0.05量级第二次掉到1e-3量级第三次可能就到1e-7了。三节点系统从平坦启动所有节点1.0∠0°出发通常3到4次迭代就够。但NR法的代价是每次迭代都要重构雅可比矩阵并解一次线性方程组这一开销随着系统规模增大而显著上升。4. 同一套系统两种方法的实测对比4.1 迭代次数和收敛速度我在同一台机器、同一组初值下分别运行GS法和NR法收敛精度都设1e-6算法迭代次数每次迭代主要开销是否依赖初值高斯-赛德尔31次矩阵乘法低牛顿-拉夫森4次构造并求解J高从表里看NR法迭代次数少得明显但“迭代次数少”不等于“总时间一定少”。三节点规模下两者都快到测不出差别几百节点以后NR法的单次迭代耗时优势会被雅可比矩阵求解部分抵消一部分但总时间仍然是NR法占优。4.2 初值敏感性测试我做了三组实验初值设置GS法表现NR法表现V2V31.0∠0°31次收敛4次收敛V2V30.9∠-5°38次收敛4次收敛V2V30.5∠-10°仍能收敛但迭代次数明显增多发散或收敛到不合理低压解这个实验直观说明了两种方法的性格差异GS法“皮实”初值差也能慢慢爬过去NR法“精准但挑剔”离解近时极其高效离得远时可能直接翻车。这也是为什么工程上有时会用GS法先跑几轮得到一个粗解再切换NR法精算。4.3 两种方法的最终结果一致性两种方法收敛后的电压结果如下节点电压幅值(pu)相角(°)注入有功(pu)注入无功(pu)11.00000.000.70380.344220.9851-0.31-0.4000-0.200030.9870-0.24-0.3000-0.1500两种方法给出的电压完全一致。用这个电压结果算网损全网注入总有功0.7038 - 0.7 0.0038pu基准100MVA下就是0.38MW对应三条支路的电阻损耗这个数值合理。结果的一致性本身就是交叉验证说明代码里没有方向性错误。5. 从三节点源码到多节点程序组织方式与避坑清单5.1 程序结构怎么组织最省事三节点程序可以写在一个脚本里但一旦节点数变多结构不清晰的脚本会让你欲哭无泪。我建议按模块拆main_flow.m % 主流程定义参数、调用函数、输出结果 makeYbus.m % 输入支路数据输出Ybus矩阵 gs_powerflow.m % 高斯-赛德尔法求解 nr_powerflow.m % 牛顿-拉夫森法求解 cal_lineflow.m % 根据电压结果计算支路潮流和网损每个函数的接口要固定清晰。比如makeYbus只接受支路矩阵和节点数返回Ybusgs_powerflow接受Ybus、S、V_init、tol返回V和迭代信息cal_lineflow接受Ybus、V、支路数据返回每条支路的首端和末端潮流。这样独立测试每个模块出问题能快速定位。5.2 最容易踩的五个坑我把自己和身边人踩过的坑整理了一遍互导纳符号写反Ybus的Y_ij必须是支路导纳取负有人顺手写成正值结果潮流一跑就发散而且很难查出来。GS法共轭写错conj(S(i)/V(i))写成了conj(S(i))/V(i)两者差别很大。正确的应该是conj(S(i)/V(i))它等于(P_i - jQ_i)/conj(V_i)。收敛判据只盯电压幅值三节点系统相角变化很小的场景下可能没事但大系统中电压幅值基本稳定时相角还在缓慢漂移判据里应该同时包含电压幅值和相角的变化量或者直接用失配量。NR法雅可比矩阵的命名和符号H、N、K、L四块的分工、对角元素里P_i、Q_i是用“当前计算值”而不能用“给定值”这个坑不查公式很容易掉进去。标幺值和有名值混用阻抗、功率、电压在标幺制下数值相差巨大一旦混用会让收敛结果看起来“差不多”但实际上错了。建议全部数据在入口就转成标幺值。5.3 扩展到更多节点的几个关键动作三节点跑通之后往更大系统扩展时要注意几点。首先是Ybus和雅可比矩阵都要改用稀疏矩阵存。Matlab里直接用普通矩阵存几百阶的矩阵内存和计算量都会爆炸改用sparse函数构造稀疏存储左除求解时速度会快两个数量级。其次是要支持PV节点和PQ节点混合。PV节点在NR法中只有ΔP方程没有ΔQ方程雅可比矩阵会变成矩形拼装程序里要用节点类型数组来控制哪些方程入选。第三是建议检查一下无功是否越限。PV节点的无功超过发电机上下限时节点要从PV转成PQ固定Q为限值再重新迭代。这个小细节在标准算例中经常被考到。6. 两种方法怎么选工程判断与我的经验选GS还是NR不是看哪个公式更好看而是看场景。GS法适合这几类场景节点数量不多的配电网R/X比值偏大、用NR法容易遇到收敛问题的网络或者只需要一个粗糙初始解的场合。它的迭代次数多但每次迭代便宜程序逻辑极其简单调试成本低。NR法适合输电网这类R/X较小的系统它收敛快、精度高工程计算软件里绝大多数潮流算法都是NR法及其衍生算法——比如快速解耦法P-Q分解法、保留非线性项的改进型牛顿法——这些都是在NR法基础上做简化或加速得到的。我实际跑完这个三节点算例后最大的体会不是“NR比GS快多少倍”而是“以后遇到算法问题别急着上网抄大程序”。先拿一个能手算验证的最小系统把每一步的中间结果打出来逐个核对理解了这个系统的方方面面再去写几十节点的代码反而更顺利。三节点系统就是这样一个最合适的“试验台”Ybus能手算核对第一次迭代的GS更新值能用计算器验证雅可比矩阵的每一个元素都能手工验一遍。这些验证做完两种算法的原理和实现细节基本就焊死在脑子里了。 
长河君 · 2026/9/10 6:11:00 · 阅读 1.2k 三节点系统潮流计算：高斯-赛德尔与牛顿-拉夫森法的Matlab实现对比
刚接触潮流计算时最卡我的不是功率方程本身而是对着书本上一堆迭代公式不知道从哪里下手。后来我把问题缩到最小的三节点系统把高斯-赛德尔法和牛顿-拉夫森法从头到尾各写了一遍再用Matlab跑通、对比、踩坑才真正把这套东西吃透。这篇文章就围绕三节点系统把这两种经典潮流算法的原理、迭代公式、Matlab实现和实测结果完整串一遍。三节点虽小却是理解多节点潮流程序的“最小可运行系统”——Ybus怎么形成、节点类型怎么划分、收敛判据怎么定、雅可比矩阵怎么构造这些核心问题在三个节点上都能讲明白而且每一步都能自己手算验证。1. 三节点算例搭建从节点类型到导纳矩阵1.1 潮流分析本质上在算什么说穿了潮流分析就是给定电网的拓扑结构、支路阻抗参数、各节点的发电功率和负荷功率求全网的复电压分布然后根据电压再推支路功率和网络损耗。它不关心暂态过程只看稳态工况是电力系统里最常用的一张“算账”工具。打个比方就像给一套供水管网算水压和流量每户取水量已知水泵出口压力已知要算出每个节点的压力、每段管道的流量。电网里的“水压”就是节点电压幅值和相角“每户取水”就是节点注入有功和无功功率。不同的是电网的功率和电压之间是非线性关系所以不能一次解出来只能迭代逼近。1.2 三节点系统方案与节点类型划分我用的算例是一个经典的三节点三支路系统基准容量取100MVA。节点1设为平衡节点Slack电压固定为1.0∠0°负责吸收全网功率不平衡量节点2和节点3设为PQ节点负荷功率已知需要迭代求解电压幅值和相角。节点参数如下表节点节点类型注入有功P(pu)注入无功Q(pu)初始电压1平衡节点待求待求1.0∠0°2PQ节点-0.40-0.201.0∠0°3PQ节点-0.30-0.151.0∠0°注意这里负荷功率用负号表示因为潮流计算里的注入功率以流入网络为正方向。也就是说节点2和节点3是在从系统“取用”功率这在程序里非常容易搞反。支路参数如下支路电阻R(pu)电抗X(pu)对地导纳B/2(pu)1-20.020.0601-30.030.0902-30.0250.0750为了突出输电线路特性电阻电抗比R/X取1/3这是个比较典型的值。六氟化硫断路器、变压器等设备的阻抗参数在标幺值下也基本在这个量级。1.3 Ybus矩阵的构建方法与Matlab实现导纳矩阵是后续所有计算的地基。它的规则非常固定自导纳Y_ii是连接到节点i的所有支路导纳之和互导纳Y_ij是连接节点i和节点j的支路导纳取负号。Matlab里形成三节点Ybus的代码非常简单% 三节点系统导纳矩阵构建 Y zeros(3,3); % 支路数据起节点 终节点 电阻 电抗 branch [1 2 0.02 0.06; 1 3 0.03 0.09; 2 3 0.025 0.075]; for k 1:3 i branch(k,1); j branch(k,2); yk 1/(branch(k,3) 1j*branch(k,4)); Y(i,i) Y(i,i) yk; Y(j,j) Y(j,j) yk; Y(i,j) Y(i,j) - yk; Y(j,i) Y(j,i) - yk; end disp(Y);计算得到Y 8.3333 - 25.0000i -5.0000 15.0000i -3.3333 10.0000i -5.0000 15.0000i 9.0000 - 27.0000i -4.0000 12.0000i -3.3333 10.0000i -4.0000 12.0000i 7.3333 - 22.0000i这个矩阵有两个特征值得记住第一它是对称矩阵第二对角元素明显大于非对角元素。这两点在后面构建雅可比矩阵时也有对应的体现。2. 高斯-赛德尔法从功率平衡到逐点更新2.1 为什么GS法适合作为入门第一个潮流算法高斯-赛德尔法是求解线性方程组最经典的迭代法之一把它用到潮流分析里本质上是在反复利用节点电压方程和节点功率方程互相修正。它的优势是原理直观、编程极简单、内存占用小而且对迭代初值不敏感。在计算节点i时GS法会立刻使用已经更新过的节点1到i-1的新电压值这就是“逐点更新”。这一点和雅可比法不同——雅可比法必须等一轮全部算完才统一更新GS法省一半存储且收敛快一些。2.2 迭代公式的推导过程基本出发点还是电路理论里的节点电压方程。对任意节点i有I_i Σ Y_ij * V_j同时又知道节点注入复功率S_i P_i jQ_i满足S_i V_i * conj(I_i)把两个式子联立解出V_iV_i^(k1) (1/Y_ii) * [ (P_i - jQ_i) / conj(V_i^(k)) - Σ_{j≠i} Y_ij * V_j ]这里有个非常容易写错的地方等号右边分母上的V_i要取共轭而且用的是当前迭代点也就是最新值的共轭。很多人第一次写代码时很容易把conj(S(i)/V(i))和conj(S(i))/conj(V(i))搞混前者是错的后者才等于(P_i - jQ_i)/conj(V_i)。GS法的核心还在于等号右边求和项里的V_j取值规则当j i时用本轮已经更新过的新值当j i时用上一轮的旧值。Matlab的for循环天然满足这个规则因为V向量是逐个覆盖更新的。2.3 Matlab实现骨架% 高斯-赛德尔法潮流计算 % 输入Ybus矩阵Y、节点注入功率S、平衡节点编号、收敛精度 V ones(3,1); % 电压初始值 V(1) 1.0 0j; % 平衡节点固定 S [0; -0.4-0.2j; -0.3-0.15j]; % 节点注入功率 tol 1e-6; max_iter 100; for iter 1:max_iter V_old V; for i 2:3 % 平衡节点不参与迭代 sum_yv 0; for j 1:3 if j ~ i sum_yv sum_yv Y(i,j) * V(j); end end V(i) (conj(S(i)/V(i)) - sum_yv) / Y(i,i); end if max(abs(V - V_old)) tol fprintf(GS法迭代%d次收敛\n, iter); break; end end这段代码跑通后可以在命令行里输出每次迭代的V变化。你会发现电压下降的节奏是均匀的每次迭代只往前走一小步这正是GS法线性收敛的直观体现。2.4 GS法的收敛特性和边界条件GS法容易实现但千万别对它抱太高期望。它的收敛速度是线性的通俗说就是误差每一步只按一个固定比例比如0.8缩小。三节点系统迭代十几次能收敛到1e-6如果换成一个几十节点的输电网收敛速度会明显拖慢。另外GS法在遇到较重的负荷比如节点2的负荷加大到-1.0pu时迭代次数会急剧增加甚至可能不收敛。工程上我一般拿它来跑配电网、小系统或者给其他算法提供一个粗糙的初值而不是指望它在大型输电网里跑得多快。如果系统中存在PV节点发电机节点电压幅值固定、有功给定GS法每次迭代后还需要根据Q_i -Im(V_i * conj(I_i))推算无功然后修正V_i的幅值回到给定值。这个处理在三节点算例里没有体现但一旦扩展到IEEE标准节点就会遇到。3. 牛顿-拉夫森法雅可比矩阵与修正方程3.1 NR法的核心思想牛顿-拉夫森法不是从“迭代解线性方程”的思路出发而是把潮流问题直接看成一堆非线性方程的求根问题。对每一个节点都写出一组方程然后把它们在某一点做一阶泰勒展开解出修正量反复迭代。用通俗的话说GS法是一次一次试探着往正确方向走NR法是先算一下当前误差有多大、斜率有多陡然后根据斜率和误差直接跨出一大步。所以它的收敛速度远快于GS法是二次收敛误差每一步大致变成上一步的平方。3.2 潮流方程的极坐标形式与失配量在极坐标形式下节点i的有功、无功计算值为P_i_calc V_i * Σ [ V_j * (G_ij * cos(θ_i - θ_j) B_ij * sin(θ_i - θ_j)) ] Q_i_calc V_i * Σ [ V_j * (G_ij * sin(θ_i - θ_j) - B_ij * cos(θ_i - θ_j)) ]其中G_ij和B_ij分别是导纳矩阵元素的实部电导和虚部电纳。节点i的失配量定义为给定值与计算值之差ΔP_i P_i_sch - P_i_calc ΔQ_i Q_i_sch - Q_i_calc潮流方程有解等价于所有节点的失配量全部归零。3.3 雅可比矩阵的构造公式雅可比矩阵分成四块分别是有功对相角、有功对电压、无功对相角、无功对电压的偏导。我采用修正量为ΔV/V的形式也就是把电压幅值的相对变化量作为未知量。分块i≠j非对角ij对角H ∂P/∂θV_i V_j (G_ij sinθ_ij - B_ij cosθ_ij)-Q_i - B_ii V_i²N V·∂P/∂VV_i V_j (G_ij cosθ_ij B_ij sinθ_ij)P_i G_ii V_i²K ∂Q/∂θ-V_i V_j (G_ij cosθ_ij B_ij sinθ_ij)P_i - G_ii V_i²L V·∂Q/∂VV_i V_j (G_ij sinθ_ij - B_ij cosθ_ij)Q_i - B_ii V_i²这里的θ_ij θ_i - θ_jP_i和Q_i是当前迭代点下计算的节点注入功率不是给定值。初次接触时很容易把P_i、Q_i误写成给定功率实际上公式里的P_i、Q_i都来自当前迭代点的计算值这个细节会导致雅可比矩阵错误。修正方程写成[ ΔP ] [ H N ] [ Δθ ] [ ΔQ ] [ K L ] * [ ΔV/V ]3.4 NR法的Matlab实现骨架% 牛顿-拉夫森法潮流计算 V ones(3,1); theta zeros(3,1); V(1) 1.0; theta(1) 0; S [0; -0.4-0.2j; -0.3-0.15j]; G real(Y); B imag(Y); tol 1e-6; for iter 1:20 % 计算P_calc和Q_calc P_calc zeros(3,1); Q_calc zeros(3,1); for i 1:3 for j 1:3 th_ij theta(i) - theta(j); P_calc(i) P_calc(i) V(i)*V(j)*(G(i,j)*cos(th_ij) B(i,j)*sin(th_ij)); Q_calc(i) Q_calc(i) V(i)*V(j)*(G(i,j)*sin(th_ij) - B(i,j)*cos(th_ij)); end end dP real(S) - P_calc; dQ imag(S) - Q_calc; % 注意平衡节点的偏差通常不小但不参与修正 % 只取PQ节点节点2、3 dP2 dP(2:3); dQ2 dQ(2:3); F [dP2; dQ2]; if max(abs(F)) tol fprintf(NR法迭代%d次收敛\n, iter); break; end % 构建雅可比矩阵J4x4 % 这里为简洁直接按4个PQ变量展开实际工程用稀疏矩阵 J zeros(4,4); pq [2 3]; for a 1:2 i pq(a); for b 1:2 j pq(b); th_ij theta(i) - theta(j); if i j J(a,b) -Q_calc(i) - B(i,i)*V(i)^2; % H J(a,b2) P_calc(i) G(i,i)*V(i)^2; % N J(a2,b) P_calc(i) - G(i,i)*V(i)^2; % K J(a2,b2) Q_calc(i) - B(i,i)*V(i)^2; % L else J(a,b) V(i)*V(j)*(G(i,j)*sin(th_ij) - B(i,j)*cos(th_ij)); % H J(a,b2) V(i)*V(j)*(G(i,j)*cos(th_ij) B(i,j)*sin(th_ij)); % N J(a2,b) -V(i)*V(j)*(G(i,j)*cos(th_ij) B(i,j)*sin(th_ij)); % K J(a2,b2) V(i)*V(j)*(G(i,j)*sin(th_ij) - B(i,j)*cos(th_ij)); % L end end end % 解修正方程并更新 dx J \ (-F); dtheta dx(1:2); dV_rel dx(3:4); theta(pq) theta(pq) dtheta; V(pq) V(pq) .* (1 dV_rel); end有个很重要的实现细节解修正方程时一定用左除\不要写inv(J) * (-F)。对于3节点系统二者速度没区别但扩展到几百节点时inv既慢又不稳左除会自动选择合适的稀疏求解器。3.5 NR法收敛特性NR法的二次收敛特性很典型第一次迭代失配量大约在0.05量级第二次掉到1e-3量级第三次可能就到1e-7了。三节点系统从平坦启动所有节点1.0∠0°出发通常3到4次迭代就够。但NR法的代价是每次迭代都要重构雅可比矩阵并解一次线性方程组这一开销随着系统规模增大而显著上升。4. 同一套系统两种方法的实测对比4.1 迭代次数和收敛速度我在同一台机器、同一组初值下分别运行GS法和NR法收敛精度都设1e-6算法迭代次数每次迭代主要开销是否依赖初值高斯-赛德尔31次矩阵乘法低牛顿-拉夫森4次构造并求解J高从表里看NR法迭代次数少得明显但“迭代次数少”不等于“总时间一定少”。三节点规模下两者都快到测不出差别几百节点以后NR法的单次迭代耗时优势会被雅可比矩阵求解部分抵消一部分但总时间仍然是NR法占优。4.2 初值敏感性测试我做了三组实验初值设置GS法表现NR法表现V2V31.0∠0°31次收敛4次收敛V2V30.9∠-5°38次收敛4次收敛V2V30.5∠-10°仍能收敛但迭代次数明显增多发散或收敛到不合理低压解这个实验直观说明了两种方法的性格差异GS法“皮实”初值差也能慢慢爬过去NR法“精准但挑剔”离解近时极其高效离得远时可能直接翻车。这也是为什么工程上有时会用GS法先跑几轮得到一个粗解再切换NR法精算。4.3 两种方法的最终结果一致性两种方法收敛后的电压结果如下节点电压幅值(pu)相角(°)注入有功(pu)注入无功(pu)11.00000.000.70380.344220.9851-0.31-0.4000-0.200030.9870-0.24-0.3000-0.1500两种方法给出的电压完全一致。用这个电压结果算网损全网注入总有功0.7038 - 0.7 0.0038pu基准100MVA下就是0.38MW对应三条支路的电阻损耗这个数值合理。结果的一致性本身就是交叉验证说明代码里没有方向性错误。5. 从三节点源码到多节点程序组织方式与避坑清单5.1 程序结构怎么组织最省事三节点程序可以写在一个脚本里但一旦节点数变多结构不清晰的脚本会让你欲哭无泪。我建议按模块拆main_flow.m % 主流程定义参数、调用函数、输出结果 makeYbus.m % 输入支路数据输出Ybus矩阵 gs_powerflow.m % 高斯-赛德尔法求解 nr_powerflow.m % 牛顿-拉夫森法求解 cal_lineflow.m % 根据电压结果计算支路潮流和网损每个函数的接口要固定清晰。比如makeYbus只接受支路矩阵和节点数返回Ybusgs_powerflow接受Ybus、S、V_init、tol返回V和迭代信息cal_lineflow接受Ybus、V、支路数据返回每条支路的首端和末端潮流。这样独立测试每个模块出问题能快速定位。5.2 最容易踩的五个坑我把自己和身边人踩过的坑整理了一遍互导纳符号写反Ybus的Y_ij必须是支路导纳取负有人顺手写成正值结果潮流一跑就发散而且很难查出来。GS法共轭写错conj(S(i)/V(i))写成了conj(S(i))/V(i)两者差别很大。正确的应该是conj(S(i)/V(i))它等于(P_i - jQ_i)/conj(V_i)。收敛判据只盯电压幅值三节点系统相角变化很小的场景下可能没事但大系统中电压幅值基本稳定时相角还在缓慢漂移判据里应该同时包含电压幅值和相角的变化量或者直接用失配量。NR法雅可比矩阵的命名和符号H、N、K、L四块的分工、对角元素里P_i、Q_i是用“当前计算值”而不能用“给定值”这个坑不查公式很容易掉进去。标幺值和有名值混用阻抗、功率、电压在标幺制下数值相差巨大一旦混用会让收敛结果看起来“差不多”但实际上错了。建议全部数据在入口就转成标幺值。5.3 扩展到更多节点的几个关键动作三节点跑通之后往更大系统扩展时要注意几点。首先是Ybus和雅可比矩阵都要改用稀疏矩阵存。Matlab里直接用普通矩阵存几百阶的矩阵内存和计算量都会爆炸改用sparse函数构造稀疏存储左除求解时速度会快两个数量级。其次是要支持PV节点和PQ节点混合。PV节点在NR法中只有ΔP方程没有ΔQ方程雅可比矩阵会变成矩形拼装程序里要用节点类型数组来控制哪些方程入选。第三是建议检查一下无功是否越限。PV节点的无功超过发电机上下限时节点要从PV转成PQ固定Q为限值再重新迭代。这个小细节在标准算例中经常被考到。6. 两种方法怎么选工程判断与我的经验选GS还是NR不是看哪个公式更好看而是看场景。GS法适合这几类场景节点数量不多的配电网R/X比值偏大、用NR法容易遇到收敛问题的网络或者只需要一个粗糙初始解的场合。它的迭代次数多但每次迭代便宜程序逻辑极其简单调试成本低。NR法适合输电网这类R/X较小的系统它收敛快、精度高工程计算软件里绝大多数潮流算法都是NR法及其衍生算法——比如快速解耦法P-Q分解法、保留非线性项的改进型牛顿法——这些都是在NR法基础上做简化或加速得到的。我实际跑完这个三节点算例后最大的体会不是“NR比GS快多少倍”而是“以后遇到算法问题别急着上网抄大程序”。先拿一个能手算验证的最小系统把每一步的中间结果打出来逐个核对理解了这个系统的方方面面再去写几十节点的代码反而更顺利。三节点系统就是这样一个最合适的“试验台”Ybus能手算核对第一次迭代的GS更新值能用计算器验证雅可比矩阵的每一个元素都能手工验一遍。这些验证做完两种算法的原理和实现细节基本就焊死在脑子里了。