简介面向电力系统初学者的MATLAB潮流计算入门示例基于三机九节点IEEE标准节点系统演示从网络拓扑定义、节点功率平衡方程建立到牛顿-拉弗森迭代求解的完整流程。这份压缩包仅含1个m脚本文件大小仅2KB属于精简的纯代码实现无额外数据文件便于直接阅读与修改调试。目前已有854人学习使用适合电力专业学生、科研人员及工程入门者对照理论教材逐行理解潮流计算。脚本重点覆盖节点导纳矩阵构建、电压初值设定、迭代收敛判据及结果输出等关键环节能帮助读者快速掌握利用MATLAB求解电力系统稳态潮流的方法验证节点电压、线路功率等物理量并为后续扩展到IEEE14节点、39节点等更复杂系统提供可直接套用的思路。1. 三机九节点IEEE节点潮流计算这套基准算例在解决什么三机九节点这个名字在电力系统仿真里出现得比多数人预想得更频繁它对应的标准称呼是IEEE 9节点测试系统3台发电机、9条母线、3台变压器和6条输电线路负荷集中在母线5、6、8。很多入门者以为它只是教学模型实际它把多机协调、多电压等级和分布式无功补偿都压进了一个能看清矩阵结构的尺度里。潮流计算在这个网络上要做的事情很明确通过迭代求解每个节点的电压幅值和相角进而得到全网功率分布。对想验证一套“电力系统潮流计算matlab”代码的人来说这是能同时检验节点编号、导纳矩阵组装和牛顿迭代的最小系统。下面按实际排错的顺序从节点分类一直写到收敛性调整。2. 三机九节点的母线分类与潮流计算基准数据2.1 母线类型决定潮流计算方程数潮流计算必须先区分三类节点因为每个节点上已知量和未知量的组合完全不同。松弛母线提供全网相位基准电压幅值和相角固定有功和无功要等迭代收敛后才能算出来PV母线给定有功出力和电压幅值相角待求无功由迭代结果回代PQ母线给定有功和无功需求电压幅值和相角都待求。三机九节点正好把三类母线凑齐了母线1接平衡机母线2和3接发电机母线4到9属于输电网络节点其中5、6、8带负荷。列方程时PQ节点同时写有功偏差和无功偏差两个方程PV节点因为电压幅值已知只写有功偏差平衡节点不参与迭代计算。因此本算例的未知量数量为 6 个 PQ 节点乘 2 加上 2 个 PV 节点乘 1一共 14 个方程数量同样是 14 个矩阵不会出现欠定或超定。节点类型已知量待求量参与方程三机九节点中的母线松弛/平衡节点V、相角P、Q不参与迭代1PV/发电机节点P、VQ、相角ΔP2、3PQ/负荷节点P、QV、相角ΔP、ΔQ4 9这三种节点的“已知/待求”对调是新手最常写错的地方。例如把PV节点当成PQ节点去给无功初值雅可比矩阵的列编号就会错位迭代初期数值波动剧烈很难在默认迭代次数内收敛。2.2 基准容量与标幺值参数归算的关键参数三机九节点标准数据默认基于 S_base 100 MVA。发电机侧母线1、2、3的有名电压等级分别是16.5 kV、18 kV、13.8 kV母线4到9是230 kV。从表面看电压标幺值都在1.0附近但实际有名值完全不同归算时用错电压等级会导致变压器阻抗偏差十倍以上。有名值转标幺值使用 Z_pu Z(Ω) × S_base / (V_base²)其中 V_base 用该母线所在电压等级的基准值。变压器支路的变比在三机九节点里按标幺值是1:1这并不表示没有变比而是两侧基准电压已经吸收了16.5/230、18/230、13.8/230的有名变比。如果去复制别人代码时把变比写成13.9或12.78潮流结果中的无功分布会明显偏移。2.3 三机九节点常用母线数据与MATLAB初始赋值表2-2是后续代码使用的母线数据负荷功率全部按吸收方向取正值发电机有功为注入方向。这里只列潮流计算真正需要的初值Qg在PV和平衡节点上迭代后会变化初值先填0。母线号类型电压设定/初值(pu)PG(pu)PD(pu)QD(pu)1松弛1.0400.716*002PV1.0251.630003PV1.0250.850004PQ1.0000005PQ1.00001.2500.5006PQ1.00000.9000.3007PQ1.0000008PQ1.00001.0000.3509PQ1.000000*标记表示松弛节点PG是收敛后回代的量输入时给一个参考值即可。PQ母线第4列的1.000是平启动初值不是给定的电压约束。MATLAB里初始化节点数据和支路数据常见做法是用矩阵直接存储每行含义写进注释% 节点数据行格式[编号, 类型, 电压幅值, 相角, Pg, Qg, Pd, Qd] % 类型标记1PQ2PV3平衡节点 bus [ 1 3 1.040 0 0.716 0 0 0; 2 2 1.025 0 1.630 0 0 0; 3 2 1.025 0 0.850 0 0 0; 4 1 1.000 0 0 0 0 0; 5 1 1.000 0 0 0 1.250 0.500; 6 1 1.000 0 0 0 0.900 0.300; 7 1 1.000 0 0 0 0 0; 8 1 1.000 0 0 0 1.000 0.350; 9 1 1.000 0 0 0 0 0; ];这段代码把相角初值全部置0PQ节点电压幅值置1.0也就是所谓的平启动。平启动适合三机九节点这种参数不极端的小系统通常10次迭代以内就能落入1e-6精度如果从其它潮流结果文件导入初值反而可能带来人工写入的电压越限。3. 节点导纳矩阵与功率偏差方程的MATLAB构建3.1 支路数据到Y矩阵的组装规则节点导纳矩阵是潮流计算的第一个硬门槛。对角线元素 Yii 是连接在节点i上所有支路导纳之和再加上该节点对地电纳非对角线元素 Yij 是连接节点i和j的支路导纳取负。变压器支路还要按变比折算否则1、2、3号母线与高压网络之间会多出几十Mvar的误差。三机九节点的支路数据一般按9条支路给出其中3条是变压器支路。线路充电电纳容易踩坑不少资料把表头写成“B/2”实际值是线路全电纳的一半组装时如果直接当成全电纳加进对角线相当于少了一半并联无功支撑最后的电压会整体偏高0.01~0.02 pu。下面表3-1直接给全电纳B代码里除以2加到两端避免歧义。支路首端末端R(pu)X(pu)全电纳B(pu)变比kT1140.00000.057601.0L1450.01700.09200.31601.0L2560.03900.17000.71601.0T3360.00000.058601.0L3670.01190.10080.41801.0L4780.00850.07200.29801.0T2270.00000.062501.0L5890.03200.16100.61201.0L6940.01000.08500.35201.0变比k在标准数据里按标幺值都是1.0不代表变压器没有变比只是计算时已经归算。如果哪天拿到非1变比参数需要引入变压器π型等值电路这时三机九节点的模板就不够用了。3.2 用MATLAB把三机九节点支路表组装成Y矩阵按支路表逐条填充Y矩阵代码写成独立函数方便后面牛顿迭代复用function Yb build_Y(branch, n) Yb zeros(n, n); for k 1:size(branch, 1) n1 branch(k, 1); n2 branch(k, 2); r branch(k, 3); x branch(k, 4); b branch(k, 5); tap branch(k, 6); if tap 0 tap 1; end y 1 / (r 1j * x); % 首端导纳按变比折算末端不折算 Yb(n1, n1) Yb(n1, n1) y / tap^2; Yb(n2, n2) Yb(n2, n2) y; % 互导纳是 -y/tap Yb(n1, n2) Yb(n1, n2) - y / tap; Yb(n2, n1) Yb(n2, n1) - y / tap; % 充电电纳分两侧各加一半 if b 0 Yb(n1, n1) Yb(n1, n1) 1j * b / 2; Yb(n2, n2) Yb(n2, n2) 1j * b / 2; end end Yb sparse(Yb); end代码中变量tap对应表3-1的变比列。这里没有写完整的变压器补偿支路因为三机九节点标准数据标幺变比就是1足够覆盖算例。第5列b用的是全电纳如果读者手里的数据表写的是“B/2”记得把代码里的加法改成1j * b或者把数据乘以2。组装完成后建议先打印对角线元素观察量级。母线4的对角导纳大约是若干个0.1~0.3级别虚数的叠加数值量级正常再往下做如果出现NaN或无穷大多半是某个支路阻抗x被填成了0这在变压器支路上最容易发生。3.3 功率偏差方程和迭代变量的对应关系潮流计算的核心不是直接解电压而是构造“给定功率”与“计算功率”的偏差。节点注入有功和无功计算式通常写成P_i V_i × Σ V_j × (G_ij × cos(δ_ij) B_ij × sin(δ_ij))Q_i V_i × Σ V_j × (G_ij × sin(δ_ij) - B_ij × cos(δ_ij))其中 G_ij、B_ij 来自Y矩阵实部和虚部δ_ij δ_i - δ_j。对PQ节点给定值来自负荷和发电机功率对PV节点只写有功偏差ΔP_i (Pg_i - Pd_i) - P_iΔQ_i (Qg_i - Qd_i) - Q_i只有PQ节点需要写出ΔQ方程PV节点的Q在迭代中是被计算出来再回代给定值的。到这里牛顿-拉夫森法的右端项和未知量映射已经明确未知量是PQ节点的V和相角、PV节点的相角右端项是这些偏差。4. 牛顿-拉夫森法求解IEEE 9节点潮流计算迭代细节4.1 极坐标雅可比矩阵四个子块的解析表达式牛顿法在每一步迭代中要解一个线性方程组系数矩阵就是雅可比矩阵。极坐标形式下未知量分成相角修正量Δδ和电压修正量ΔV/V两段所以雅可比矩阵天然分成四块H、N、M、L。子块偏导对象表达式(i≠j)在迭代中的作用HΔP对Δδ-V_i·V_j·(G_ij·sinδ_ij - B_ij·cosδ_ij)反映有功和相角的强耦合NΔP对ΔV/V-V_i·V_j·(G_ij·cosδ_ij B_ij·sinδ_ij)有功对电压的灵敏度MΔQ对ΔδV_i·V_j·(G_ij·cosδ_ij B_ij·sinδ_ij)无功对相角的灵敏度LΔQ对ΔV/V-V_i·V_j·(G_ij·sinδ_ij - B_ij·cosδ_ij)无功和电压的强耦合对角元素稍有不同计算时用 H_ii -Q_i - B_ii·V_i²、N_ii P_i G_ii·V_i²、M_ii P_i - G_ii·V_i²、L_ii Q_i - B_ii·V_i²。注意这里Q_i、P_i是当前迭代点的计算注入功率不是给定功率组装时要区分清楚。4.2 牛顿-拉夫森主循环的MATLAB代码和参数说明主循环采用常见的极坐标牛顿-拉夫森写法先算功率残差再组装雅可比解方程后更新相角和电压% 平启动 V ones(9, 1); V([1 2 3]) [1.040; 1.025; 1.025]; theta zeros(9, 1); id_PQ find(bus(:, 2) 1); id_PV find(bus(:, 2) 2); id_SL find(bus(:, 2) 3); nPQ length(id_PQ); nPV length(id_PV); nTheta nPQ nPV; P_spec bus(:, 5) - bus(:, 7); Q_spec bus(:, 6) - bus(:, 8); tol 1e-6; max_iter 20; for iter 1:max_iter % 计算当前注入功率 Pc zeros(9, 1); Qc zeros(9, 1); for i 1:9 for k 1:9 Gik real(Y(i, k)); Bik imag(Y(i, k)); th theta(i) - theta(k); Pc(i) Pc(i) V(i) * V(k) * (Gik * cos(th) Bik * sin(th)); Qc(i) Qc(i) V(i) * V(k) * (Gik * sin(th) - Bik * cos(th)); end end dP P_spec - Pc; dQ Q_spec - Qc; dP(id_SL) []; % 去掉平衡节点 dQ(id_PV) []; % 去掉PV节点 dQ(id_SL) []; mismatch [dP; dQ]; if max(abs(mismatch)) tol fprintf(converged at iteration %d\n, iter); break; end % 组装雅可比 J zeros(nTheta nPQ, nTheta nPQ); mapT zeros(9, 1); mapV zeros(9, 1); idAP [id_PQ; id_PV]; mapT(idAP) (1:nTheta); mapV(id_PQ) (1:nPQ); for i idAP rP mapT(i); for k 1:9 if k i || mapT(k) 0 continue; end Gik real(Y(i, k)); Bik imag(Y(i, k)); th theta(i) - theta(k); c V(i) * V(k); J(rP, mapT(k)) J(rP, mapT(k)) - c * (Gik * sin(th) - Bik * cos(th)); if mapV(k) 0 J(rP, nTheta mapV(k)) J(rP, nTheta mapV(k)) ... - c * (Gik * cos(th) Bik * sin(th)); end if mapV(i) 0 J(nTheta mapV(i), mapT(k)) J(nTheta mapV(i), mapT(k)) ... c * (Gik * cos(th) Bik * sin(th)); if mapV(k) 0 J(nTheta mapV(i), nTheta mapV(k)) J(nTheta mapV(i), nTheta mapV(k)) ... - c * (Gik * sin(th) - Bik * cos(th)); end end end % 对角项 Gi real(Y(i, i)); Bi imag(Y(i, i)); J(rP, mapT(i)) J(rP, mapT(i)) (-Qc(i) - Bi * V(i)^2); if mapV(i) 0 J(rP, nTheta mapV(i)) J(rP, nTheta mapV(i)) (Pc(i) Gi * V(i)^2); J(nTheta mapV(i), mapT(i)) J(nTheta mapV(i), mapT(i)) (Pc(i) - Gi * V(i)^2); J(nTheta mapV(i), nTheta mapV(i)) J(nTheta mapV(i), nTheta mapV(i)) (Qc(i) - Bi * V(i)^2); end end dx J \ mismatch; theta(idAP) theta(idAP) dx(1:nTheta); V(id_PQ) V(id_PQ) .* (1 dx(nTheta 1:end)); end这段代码的要点在于变量编号mapT把PQ和PV节点统一编成相角修正列mapV只给PQ节点编电压修正列两段列在雅可比中首尾相接。第14行到第17行删除残差时顺序必须和列编号一致否则求解出的修正量被错位加到错误的母线上典型表现是迭代两三步后残差反而增大。收敛后松弛母线的电压和相角固定它的 Pc(1)、Qc(1) 就是整网不平衡功率的提供值。4.3 收敛判据与初值设置的直接影响三机九节点算例中1e-6的收敛精度对应约0.1 kW/Mvar级别的残差对通用潮流分析已经足够。如果只是看电压分布1e-4也能用但如果后面要做静态安全分析、N-1校核建议保持1e-6或更严因为线路潮流的微小误差在切除支路后会被放大。初值方面平启动基本不需要改除非负荷数据被改成重载工况。若是将总负荷提高50%平启动仍能收敛只是迭代次数从7次左右增加到12次左右当负荷超过系统输送极限就会出现相角拉不开、残差震荡的问题这时先降低负荷确认代码正确性再逐步加载荷找临界点。这种做法比直接调大迭代次数更有诊断价值。5. 潮流计算结果的验证技巧与不收敛修正5.1 用松弛节点功率和线路负载率反向验证计算结果迭代收敛不等于数据正确。第一个要验证的是全网有功平衡松弛节点有功应等于总负荷加网络损耗再减去其它发电机有功。三机九节点总负荷是3.15 pu第二、三台发电机合计2.48 pu所以松弛节点的P应在0.68 pu附近再叠加线路损耗通常落在0.70~0.75 pu区间。如果算出来是负数或超过1.0先查节点数据里的PG、PD填反了没有。接着可以按支路首端潮流公式核算线路负载率% 支路i-j的首端潮流 Pij V(i)^2 * gij - V(i) * V(j) * ... (gij * cos(delta(i) - delta(j)) bij * sin(delta(i) - delta(j))); load_rate sqrt(Pij^2 Qij^2) / S_max;S_max 取该线路热稳定极限三机九节点中没有统一给定工程上常用100 MVA基准下按1.0~1.5 pu假设。负载率超过0.8的支路优先复核。实际的IEEE 9节点基准状态下线路负载率多在0.3到0.6之间如果某条支路达到0.95以上多半是母线编号或变比数据错位。5.2 三种常见不收敛场景和调整顺序现象优先检查位置第一次迭代就出现NaNY矩阵中有0阻抗支路或导纳稀疏矩阵被错误稀疏化残差前几步下降之后停在1e-3量级PV节点的无功越限或残差行与列编号映射不一致迭代次数超过20仍不收敛负荷太重或线路电纳符号加反如果现象属于第一类用full(diag(Y))打印对角线能直接看到NaN来自哪条母线。第二类情况在入门级代码里最隐蔽dQ删除PV行时没有同步删除雅可比对应行导致矩阵不对称。第三类则先用轻载数据验证程序比如把全部负荷乘0.5能收敛就说明算法骨架没问题。提示调试牛顿法时每次都打印max(abs(mismatch))。正常收敛会看到残差大致按二次速率下降例如从1e-2掉到1e-4再掉到1e-8如果收敛曲线变平先怀疑无功越限处理再怀疑数据符号这是定位潮流计算问题最快的手段。本文还有配套的精品资源点击获取