☰
三机九节点潮流计算:从导纳矩阵到牛顿-拉夫逊法MATLAB实现
2026/10/14 14:38:19 网站建设 项目流程

简介:这份资源面向电力系统专业学生、研究人员及工程师,提供一套基于MATLAB实现的IEEE 3机9节点潮流计算程序,适合用于理解潮流计算原理、开展教学演示或作为课程设计的基础代码。压缩包共3个文件,包含2个txt数据描述文件和1个m脚本文件,整体仅2KB,其中txt文件用于存放节点与支路参数,m文件则承载核心计算逻辑,结构精简、便于快速上手。目前已有4069人学习下载,说明该模型在电力系统分析教学中具有较高的参考价值。程序围绕牛顿-拉夫森迭代法展开,涵盖数据输入、网络建模、雅可比矩阵构建、迭代求解与结果输出等关键环节,读者可借此掌握节点电压幅值、相角及支路功率的求解流程,并理解功率平衡与电压稳定性等概念。代码还可作为扩展基础,适配更大规模系统或不同运行场景,是学习潮流计算算法与MATLAB编程实践的实用工具。

1. 三机九节点算例为什么成了电力仿真的“最小可用系统”

如果你正在做电力系统潮流计算,大概率绕不开 IEEE3机9节点这个算例。它只有3台发电机、9条母线、3个负荷,规模小到可以在纸上画出来,却是绝大多数调度、稳定、新能源并网研究的起点。原因很直接:它足够小,小到你能在几分钟内跑完一次潮流并肉眼核对结果;又足够典型,包含平衡节点、PV节点、PQ节点三种母线类型,包含变压器支路和并联补偿,能覆盖潮流计算里几乎所有关键机制。

我见过太多人卡在第一步:拿到一个三机九节点的数据表,却不知道怎么把它变成能跑的 MATLAB 程序。要么是导纳矩阵拼错,要么是雅可比矩阵维度对不上,要么是迭代不收敛却不知道从哪查。这篇笔记就按我自己的实操路径,把从数据到结果的全过程拆开讲清楚,包括参数怎么设、代码怎么写、不收敛时先看哪里。适合刚接触潮流计算的学生,也适合需要快速搭一个验证环境的工程师。

2. 从节点数据到导纳矩阵:三机九节点的建模底座

2.1 先搞清楚三机九节点里每个数字代表什么

三机九节点的标准数据一般包含三张表:母线数据、支路数据、发电机数据。母线表里每行给出节点编号、电压幅值初值、相角初值、有功负荷、无功负荷、发电机有功、发电机无功,以及节点类型。节点类型用 1 表示 PQ 节点,2 表示 PV 节点,3 表示平衡节点。三机九节点里,节点1通常是平衡节点,节点2和节点3是 PV 节点,其余节点是 PQ 节点。

支路表里每行给出首末端节点编号、电阻、电抗、对地电纳(总充电电纳)、变压器变比。如果变比等于1,说明是普通线路;如果变比不等于1,说明是变压器支路。三机九节点里有三条变压器支路,分别连接节点1-4、节点2-7、节点3-9,变比通常取1.0或者接近1.0的值,具体看算例版本。

发电机表里给出每台发电机所在节点、有功出力、电压幅值设定值、无功上下限。PV节点的电压幅值由发电机励磁维持,所以在潮流计算中作为已知量处理。

注意:不同教材和课程资料里的三机九节点数据在负荷和阻抗数值上可能有细微差别,但拓扑结构一致。建议先固定一套数据,跑通之后再换数据验证。

2.2 用节点支路关联矩阵拼出导纳矩阵

导纳矩阵 Ybus 是 n×n 的复数矩阵,n 是节点数。对于三机九节点,n=9。拼装规则很直接:对角元素是连接到该节点的所有支路导纳之和,加上该节点的对地并联导纳;非对角元素是节点 i 和节点 j 之间支路导纳的负值。

对于普通线路,支路导纳 y = 1/(r + jx),对地电纳的一半分别加到两端节点的对角元素上。对于变压器支路,如果变比 k 不等于1,需要把变压器等效成 π 型或者直接用变比修正。常见做法是把变压器支路处理成:非标准变比侧的对角元素加上 y/k,标准变比侧加上 y,非对角元素为 -y/k。

下面是我常用的导纳矩阵拼装代码,直接可抄:

% 三机九节点导纳矩阵拼装 % 节点数 n = 9; % 初始化 Ybus Ybus = zeros(n, n); % 支路数据: [首端 末端 R X B 变比] % B 为总充电电纳,变比 k 非零 branch = [ 1 4 0.0 0.0576 0.0 1.0; 4 5 0.017 0.092 0.158 1.0; 5 6 0.039 0.170 0.358 1.0; 3 6 0.0 0.0586 0.0 1.0; 6 7 0.0119 0.1008 0.209 1.0; 7 8 0.0085 0.072 0.149 1.0; 8 2 0.0 0.0625 0.0 1.0; 8 9 0.032 0.161 0.306 1.0; 9 4 0.010 0.085 0.176 1.0; ]; for k = 1:size(branch, 1) i = branch(k, 1); j = branch(k, 2); R = branch(k, 3); X = branch(k, 4); B = branch(k, 5); tap = branch(k, 6); % 支路串联导纳 y = 1 / (R + 1j*X); if tap == 1.0 % 普通线路 Ybus(i, i) = Ybus(i, i) + y + 1j*B/2; Ybus(j, j) = Ybus(j, j) + y + 1j*B/2; Ybus(i, j) = Ybus(i, j) - y; Ybus(j, i) = Ybus(j, i) - y; else % 变压器支路,i 侧为非标准变比侧 Ybus(i, i) = Ybus(i, i) + y / (tap * conj(tap)); Ybus(j, j) = Ybus(j, j) + y; Ybus(i, j) = Ybus(i, j) - y / conj(tap); Ybus(j, i) = Ybus(j, i) - y / tap; end end

这段代码的逻辑说明:先初始化全零矩阵,然后逐条支路累加。普通线路的对角元素加串联导纳加半个充电电纳,非对角元素减串联导纳。变压器支路用变比的共轭修正,保证功率守恒。参数说明:R 和 X 是标幺值,B 是总充电电纳标幺值,tap 是变比。如果你的算例里变压器变比是1.0,那三条变压器支路会走普通线路分支,结果一样。

拼完 Ybus 后,建议立刻检查两个东西:一是 Ybus 是否对称,二是每行对角元素绝对值是否大于该行非对角元素绝对值之和。如果不对称,说明变比处理有问题;如果对角不占优,说明某条支路参数填错了。

2.3 节点分类和功率方程怎么写

潮流计算的核心方程是每个节点的功率平衡。对 PQ 节点,有功和无功都是已知量,电压幅值和相角是未知量。对 PV 节点,有功和电压幅值是已知量,无功和相角是未知量。对平衡节点,电压幅值和相角都已知,有功和无功是待求量。

三机九节点里,节点1是平衡节点,电压幅值设为1.04,相角设为0。节点2和节点3是 PV 节点,电压幅值分别设为1.025和1.025,有功出力根据发电机数据给定。其余节点是 PQ 节点,负荷有功无功按母线数据给定。

功率方程写成极坐标形式:

P_i = V_i * sum_j V_j * (G_ij * cos(theta_i - theta_j) + B_ij * sin(theta_i - theta_j))

Q_i = V_i * sum_j V_j * (G_ij * sin(theta_i - theta_j) - B_ij * cos(theta_i - theta_j))

其中 G_ij 和 B_ij 是 Ybus 的实部和虚部。这个方程对每个节点都成立,但已知量和未知量的划分不同。

3. 牛顿-拉夫逊法迭代:雅可比矩阵怎么拼、方程怎么解

3.1 为什么三机九节点适合用牛顿-拉夫逊法

潮流计算方法有好几种:高斯-赛德尔法、牛顿-拉夫逊法、PQ分解法。三机九节点规模小,高斯-赛德尔法也能收敛,但迭代次数多,而且对初值敏感。牛顿-拉夫逊法收敛速度快,通常4到6次迭代就能到1e-8精度,而且雅可比矩阵的结构清晰,适合用来理解潮流计算的本质。

PQ分解法虽然更快,但它利用了高压电网的物理特性,对三机九节点这种小系统优势不明显,而且实现起来多一层近似,反而容易在参数设置上翻车。所以我一般先用牛顿-拉夫逊法跑通,确认结果正确后再考虑换算法。

3.2 雅可比矩阵的四个子块和拼装顺序

牛顿-拉夫逊法的核心是求解修正方程:

[ΔP; ΔQ] = J * [Δθ; ΔV/V]

其中 J 是雅可比矩阵,包含四个子块:H = ∂P/∂θ,N = ∂P/∂(V/V),M = ∂Q/∂θ,L = ∂Q/∂(V/V)。注意这里用的是 ΔV/V 而不是 ΔV,这样雅可比矩阵的元素形式更整齐。

对于三机九节点,未知量包括所有 PQ 节点的电压幅值和相角,以及所有 PV 节点的相角。平衡节点的幅值和相角都不参与迭代。假设有 nPQ 个 PQ 节点,nPV 个 PV 节点,那么未知量总数是 2nPQ + nPV。三机九节点里,PQ 节点有6个(节点4到9),PV 节点有2个(节点2和3),所以未知量是 26 + 2 = 14 个。雅可比矩阵是 14×14。

雅可比矩阵的元素公式在教科书里都有,但拼装顺序容易搞错。我习惯按节点编号顺序先排所有 PQ 节点的相角,再排所有 PV 节点的相角,最后排所有 PQ 节点的电压幅值。这样雅可比矩阵的结构是:

J = [H N; M L]

其中 H 是 (nPQ+nPV) × (nPQ+nPV),N 是 (nPQ+nPV) × nPQ,M 是 nPQ × (nPQ+nPV),L 是 nPQ × nPQ。

下面是对应的雅可比矩阵拼装代码:

% 雅可比矩阵拼装 % 输入: Ybus, V, theta, PQ节点索引, PV节点索引 % 输出: J 雅可比矩阵, mismatch 不平衡量 % 节点分类 pq = [4 5 6 7 8 9]; % PQ节点编号 pv = [2 3]; % PV节点编号 slack = 1; % 平衡节点 nbus = 9; npq = length(pq); npv = length(pv); % 计算注入功率 P = zeros(nbus, 1); Q = zeros(nbus, 1); for i = 1:nbus for j = 1:nbus P(i) = P(i) + V(i)*V(j)*(real(Ybus(i,j))*cos(theta(i)-theta(j)) + imag(Ybus(i,j))*sin(theta(i)-theta(j))); Q(i) = Q(i) + V(i)*V(j)*(real(Ybus(i,j))*sin(theta(i)-theta(j)) - imag(Ybus(i,j))*cos(theta(i)-theta(j))); end end % 给定功率 Psp = [0; 1.63; 0.85; 0; 0; 0; 0; 0; 0]; % 节点有功给定 Qsp = [0; 0; 0; 0; 0; 0; 0; 0; 0]; % 节点无功给定,PV和平衡节点不参与 % 不平衡量 dP = Psp - P; dQ = Qsp - Q; % 雅可比矩阵初始化 J = zeros(2*npq + npv, 2*npq + npv); % 排列顺序: [dtheta_pq; dtheta_pv; dV_pq] % 先填 H 块: dP/dtheta for i = 1:npq for j = 1:npq if i == j J(i, j) = -Q(pq(i)) - imag(Ybus(pq(i), pq(i)))*V(pq(i))^2; else J(i, j) = V(pq(i))*V(pq(j))*(real(Ybus(pq(i),pq(j)))*sin(theta(pq(i))-theta(pq(j))) - imag(Ybus(pq(i),pq(j)))*cos(theta(pq(i))-theta(pq(j)))); end end end % 后续块按同样方式填充,此处省略重复结构

这段代码的关键点:不平衡量的符号是给定值减计算值,雅可比矩阵的元素符号要和这个定义一致。如果符号搞反,迭代会发散。参数说明:Psp 和 Qsp 是节点给定功率,PV 节点的 Qsp 不参与不平衡量计算,因为无功是待求量。雅可比矩阵的维度是 14×14,填完后用左除求解修正量。

3.3 迭代收敛判据和最大迭代次数怎么设

收敛判据一般用不平衡量的无穷范数,也就是 max(abs(dP), abs(dQ)) 小于某个阈值。阈值通常取 1e-6 到 1e-8。三机九节点用 1e-8 就足够,迭代4到6次就能到。最大迭代次数设 20 到 50 都行,我一般设 30,超过就说明有问题,需要查数据或初值。

迭代过程就是:算不平衡量,判断是否收敛,不收敛就拼雅可比矩阵,解修正方程,更新电压幅值和相角,再回到算不平衡量。循环直到收敛或达到最大次数。

注意:更新电压幅值时用的是 ΔV/V,所以实际更新是 V = V * (1 + ΔV/V)。相角直接加 Δθ。平衡节点的幅值和相角不更新。

4. 不收敛、结果不对、雅可比矩阵报错:排查清单

4.1 迭代发散或者振荡

现象:不平衡量在几次迭代后越来越大,或者来回振荡不下降。

原因:最常见的是雅可比矩阵符号错误,或者不平衡量定义反了。另一个常见原因是导纳矩阵拼错,特别是变压器支路的变比处理。还有一种可能是初值给得太离谱,比如电压幅值全设成1.0但相角全设成0,对三机九节点来说其实够用,但如果负荷特别重,可能需要平启动之外的处理。

解决:先检查不平衡量符号。给定功率减计算功率,如果写成计算减给定,迭代方向就反了。再检查雅可比矩阵的 H 块对角元素,应该是 -Q_i - B_ii * V_i^2。如果写成正号,也会发散。最后检查 Ybus 是否对称,不对称就查变压器支路。

4.2 收敛后结果和标准答案对不上

现象:迭代收敛了,但节点电压幅值和相角跟教材或标准结果差很多。

原因:可能是节点类型搞错了,比如把 PV 节点当成 PQ 节点算。也可能是负荷数据或发电机出力填错了。还有一种情况是变比方向搞反了,变压器支路的非标准变比侧和标准变比侧弄反,导致功率分布不对。

解决:先核对节点类型。三机九节点里节点1是平衡节点,节点2和3是 PV 节点,其余是 PQ 节点。再核对发电机有功出力和负荷有功无功。最后检查变压器支路的首末端,非标准变比侧通常是发电机侧。

4.3 雅可比矩阵维度对不上

现象:拼雅可比矩阵时报错,说维度不一致。

原因:未知量个数算错了。三机九节点里 PQ 节点6个,PV 节点2个,未知量是 2*6+2=14。如果误把平衡节点也算进去,或者把 PV 节点的电压幅值也当成未知量,维度就会多出来。

解决:明确未知量列表:所有 PQ 节点的相角和幅值,所有 PV 节点的相角。平衡节点的幅值和相角都不参与。PV 节点的幅值是已知量,不参与迭代。

4.4 程序跑得慢或者内存不够

现象:三机九节点这么小的系统,跑起来却卡顿。

原因:可能是用了全矩阵存储但没预分配,或者循环里反复调用函数。MATLAB 里循环不预分配数组会动态扩容,小系统看不出来,但养成习惯不好。

解决:雅可比矩阵和 Ybus 都提前用 zeros 初始化。不平衡量也用 zeros 预分配。三机九节点规模小,预分配后跑一次不到0.1秒。

4.5 平衡节点功率算出来不对

现象:节点电压和相角都收敛了,但平衡节点的有功无功跟预期差很多。

原因:平衡节点的功率是在迭代收敛后回代算出来的,如果 Ybus 或电压相角有误差,平衡节点功率会放大误差。另外,如果平衡节点电压幅值设得不对,比如设成1.0而不是1.04,平衡节点功率也会偏。

解决:收敛后重新计算所有节点的注入功率,包括平衡节点。检查平衡节点电压幅值是否和算例给定一致。三机九节点里平衡节点电压通常设1.04,相角设0。

5. 从能跑到好用:三机九节点程序的验证和扩展技巧

跑通一次潮流不算完,能验证结果、能改数据、能扩展到其他算例,才算真正把这个程序用起来。我一般会做三件事:核对平衡节点功率、检查线路损耗、换一组负荷数据再跑一次。

核对平衡节点功率是最直接的验证。三机九节点收敛后,平衡节点的有功应该等于所有负荷有功加上线路有功损耗,再减去其他发电机有功。如果差得远,说明 Ybus 或电压结果有问题。线路损耗可以用支路电流算,也可以用两端节点电压和支路导纳算。我习惯用支路功率求和的方式,把所有支路的有功损耗加起来,和平衡节点有功减发电机有功加负荷有功对比,误差应该在1e-6以内。

检查线路损耗还能发现另一个问题:如果某条支路损耗特别大,可能是这条支路的电阻填错了,或者变比方向反了。三机九节点里线路损耗不大,但变压器支路如果变比处理不对,损耗会明显异常。

换一组负荷数据再跑一次,是验证程序鲁棒性的好办法。把负荷有功无功乘以1.2,看程序还能不能收敛。如果收敛,说明程序没问题;如果不收敛,可能是负荷过重导致电压崩溃,这时候需要调整初值或者用更稳健的算法。三机九节点在负荷乘以1.2时通常还能收敛,乘以1.5可能就接近极限了。

下面是一个验证平衡节点功率的代码片段:

% 收敛后验证平衡节点功率 % 计算所有节点注入功率 P_calc = zeros(nbus, 1); Q_calc = zeros(nbus, 1); for i = 1:nbus for j = 1:nbus P_calc(i) = P_calc(i) + V(i)*V(j)*(real(Ybus(i,j))*cos(theta(i)-theta(j)) + imag(Ybus(i,j))*sin(theta(i)-theta(j))); Q_calc(i) = Q_calc(i) + V(i)*V(j)*(real(Ybus(i,j))*sin(theta(i)-theta(j)) - imag(Ybus(i,j))*cos(theta(i)-theta(j))); end end % 平衡节点功率 P_slack = P_calc(1); Q_slack = Q_calc(1); % 总负荷 P_load_total = sum(Psp(pq)); Q_load_total = sum(Qsp(pq)); % 发电机总有功(不含平衡节点) P_gen_total = sum(Psp(pv)); % 线路总有功损耗 P_loss = P_slack + P_gen_total - P_load_total; fprintf('平衡节点有功: %.6f\n', P_slack); fprintf('线路总有功损耗: %.6f\n', P_loss);

这段代码的逻辑:收敛后重新计算所有节点注入功率,平衡节点的注入功率就是平衡节点功率。线路总有功损耗等于平衡节点有功加其他发电机有功减总负荷有功。如果 P_loss 是负数或者特别大,说明结果有问题。参数说明:Psp 和 Qsp 是给定功率向量,pq 和 pv 是节点索引。三机九节点里 P_loss 通常在0.05到0.1之间,具体看算例数据。

扩展方面,三机九节点程序稍作修改就能跑其他算例。把节点数、支路数据、节点类型改成对应算例的,雅可比矩阵拼装部分不用大改,因为牛顿-拉夫逊法的结构是通用的。我一般会把节点数据和支路数据写成独立的输入文件,程序里只负责读数据和算潮流,这样换算例时只改数据文件,不动代码。

最后一个技巧:用稀疏矩阵。三机九节点规模小,稠密矩阵够用,但如果扩展到几十个节点,雅可比矩阵会变得很大,这时候用 sparse 存储能省内存。MATLAB 里把 zeros 换成 sparse 就行,左除会自动用稀疏求解器。不过三机九节点没必要,稠密矩阵反而更快。

我自己的习惯是:每次改完数据或代码,先跑一次三机九节点,看平衡节点功率和线路损耗是否合理,再跑目标算例。三机九节点就像一把尺子,能量出程序有没有基本问题。希望帮到你。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询