简介面向电气工程及相关专业学生的电力系统稳态分析课程设计资料聚焦两机五节点网络的潮流计算适合用于课程设计报告撰写或潮流算法编程入门。资料以PDF文件形式呈现共1个文件压缩包大小约310KB虽体量精简但内容覆盖完整目前已有203人学习下载。文档围绕牛顿拉夫逊法和PQ分解法展开先介绍潮流分析原理与常用方法对比再结合MATLAB环境讲解节点导纳矩阵形成、雅可比矩阵构建、迭代修正等关键步骤并给出两机五节点系统的具体计算程序与运行结果。通过对照高斯-赛德尔法可帮助理解不同潮流算法的收敛速度、适用场景及优缺点进而掌握电力系统稳态分析的核心技能。整份资料结构清晰包含目录、原理简介、程序框图、源码实现和结果分析适合用于课程报告撰写或作为MATLAB潮流计算编程的参考模板。1. 两机五节点潮流计算为什么是牛拉法和PQ法的合适试验场两机五节点这个规模看起来不大却正好卡在“算得明白”和“算得有内容”的交界处。它有两个电源节点、三个负荷节点既能看出平衡节点如何自动吸收全系统不平衡功率也能让牛拉法和PQ分解法在迭代次数、雅可比矩阵稠密度、初值敏感度上显出真实差异。比三节点算例信息量足又比IEEE 14节点系统容易把每一步修正过程打印出来逐行核对。电力系统稳态分析课程设计选这个题目核心工作就是两件事把节点功率方程对电压相角和幅值的偏导关系写对再把两种非线性迭代方法在MATLAB里跑通并解释清楚收敛行为。下文按“方程推导 → 矩阵形成 → 代码实现 → 收敛排错 → 结果验收”的顺序展开中途给出一套可以直接抄进课程设计报告的参数和MATLAB函数骨架。2. 牛拉法极坐标修正方程与两机五节点雅可比矩阵构造2.1 节点类型划分与两机五节点的基准值选择潮流计算的第一步不是写迭代而是先把节点分类和基准值定死。两机五节点系统里节点1设平衡节点节点2设PV节点其余三个设为PQ节点。基准功率取100 MVA基准电压取220 kV全部数据转成标幺值后进入同一个计算框架。节点类型给定条件说明1平衡节点V1.06∠0°有功无功不预先给定承担系统剩余不平衡功率2PV节点P_G0.50 puV1.00无功在-0.200.80 pu内可调3PQ节点P_L0.55 puQ_L0.25 pu恒功率负荷4PQ节点P_L0.50 puQ_L0.20 pu恒功率负荷5PQ节点P_L0.30 puQ_L0.10 pu恒功率负荷注意这里给的是负荷功率程序里需要换算成节点注入功率即注入为负值。支路数据按线路的串联电阻、电抗和半充电电容给出支路r (pu)x (pu)b/2 (pu)1-20.050.200.0201-30.030.150.0152-30.040.180.0102-40.040.200.0122-50.050.220.0123-40.060.250.0204-50.060.200.015b/2表示线路对地电容等值导纳的一半放在支路两端各一次所以一条支路对整个系统贡献的总对地导纳是2×(b/2)b。这个细节在构造节点导纳矩阵时最容易错课程设计报告里建议单独写一段说明。2.2 极坐标功率方程与注入不平衡量的计算极坐标潮流方程把每个节点的注入有功和无功写成电压幅值与相角的函数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来自节点导纳矩阵的实部与虚部。迭代时需要计算“给定注入”与“当前电压下计算注入”之差这就是功率不平衡量ΔP和ΔQ。用MATLAB可以非常简洁地写这段计算function [P, Q] calc_power(V, theta, Y) % 根据当前电压幅值和相角计算各节点注入功率 n length(V); E V .* exp(1i * theta); % 节点电压相量 I Y * E; % 节点注入电流相量 S E .* conj(I); % 复功率 P real(S); Q imag(S); end这里直接用了相量形式的复功率公式SV·conj(I)结果和逐项累加功率方程完全一致但代码更短适合在报告中作为“功率计算模块”单独贴出。注意输入V是幅值向量theta是弧度制的相角向量Y是n×n节点导纳矩阵。2.3 牛拉法雅可比矩阵四个子块的表达式牛拉法极坐标形式的修正方程写成块矩阵[ΔP; ΔQ] [H N; J L] × [Δθ; ΔV/V]这是“ΔV/V”形式代码里最后再把ΔV/V乘回V。四个子块的元素必须和功率方程一一对应求偏导下面是教材级结论直接用于MATLAB组装子块非对角元 (i≠k)对角元 (ik)H∂P/∂θ−V_i V_k (G_ik sinθ_ik − B_ik cosθ_ik)−Q_i − B_ii V_i²N∂P/∂V×VV_i V_k (G_ik cosθ_ik B_ik sinθ_ik)P_i G_ii V_i²J∂Q/∂θV_i V_k (G_ik cosθ_ik B_ik sinθ_ik)P_i − G_ii V_i²L∂Q/∂V×VV_i V_k (G_ik sinθ_ik − B_ik cosθ_ik)Q_i − B_ii V_i²注意代码实现时符号约定要统一。这里采用ΔPHΔθNΔV/V、ΔQJΔθLΔV/V的形式那么不平衡量必须按“给定注入功率 − 计算注入功率”来定义二者符号绑定不能只取一个负号。2.4 牛拉法迭代流程与收敛判据的组织顺序牛拉法单次迭代的流程固定为计算ΔP、ΔQ → 判断是否小于容差 → 组装四个子块 → 解修正方程 → 更新θ和V。程序里用idx记录除平衡节点外的所有节点用npq记录PQ节点这样修正方程的维度和映射关系才清楚。ns find(bus(:,1) ~ 1); % 非平衡节点 npq find(bus(:,1) 3); % PQ节点 npv find(bus(:,1) 2); % PV节点 Jac [H(ns,ns) N(ns,npq); J(npq,ns) L(npq,npq)]; dF [dP(ns); dQ(npq)]; dx Jac \ dF; dTheta zeros(n,1); dVnorm zeros(n,1); dTheta(ns) dx(1:length(ns)); dVnorm(npq) dx(length(ns)1:end); theta theta dTheta; V V V .* dVnorm;这段代码的要点在于矩阵索引。H的列必须对应所有参与θ迭代的节点而N的列只对应PQ节点因为PV节点的电压幅值不迭代。求解之后更新V时用的是逐元素乘法即VVV.*dVnorm反过来如果在组装前给N乘了一次V这里就不能再乘。收敛判据我一般用max(abs(dP(ns)))和max(abs(dQ(npq)))同时小于1e-6单纯看ΔP会漏掉无功不收敛的情况。3. 从牛拉法到PQ分解B与B矩阵的形成与迭代格式3.1 解耦假设成立的两条线路条件PQ分解法不是凭空出现的简化它基于电力系统高压输电网的两个物理事实正常运行时节点电压相角差很小所以sinθ_ik≈θ_ik、cosθ_ik≈1线路电抗远大于电阻无功功率对相角不敏感有功功率对电压幅值也不敏感。把这两个条件代入牛拉法修正方程N和J两个子块就可以忽略修正方程解耦成有功迭代和无功迭代两部分。需要强调的是这种解耦在配电网络里往往不成立因为配电网r/x比值高强行用PQ法会造成迭代发散或收敛到错误解。这一点课程设计答辩时经常被问到要提前准备好说法。3.2 B与B矩阵的元素构成与符号约定PQ分解法核心是两个常数矩阵B和B。这里直接给出与MATLAB代码匹配的形成规则矩阵参与元素忽略元素用途B串联支路电抗形成的导纳虚部电阻、对地电容用于ΔP/BB的基础上叠加对地电容导纳电阻用于ΔQ/符号约定上B和B都取节点导纳矩阵虚部的相反数所以对角元为正、非对角元为负。下面这段代码同时构造两个矩阵注意B2在B1基础上把每条线路两端节点的对地电容2b2加回对角function [B1, B2] build_B_matrix(branch, n) % 构造PQ分解法需要的B矩阵和B矩阵 B1 zeros(n); for k 1:size(branch,1) i branch(k,1); j branch(k,2); r branch(k,3); x branch(k,4); b2 branch(k,5); y_series 1/(r 1i*x); b -imag(y_series); % 串联导纳虚部取反 B1(i,i) B1(i,i) b; B1(j,j) B1(j,j) b; B1(i,j) B1(i,j) - b; B1(j,i) B1(j,i) - b; end B2 B1; for k 1:size(branch,1) i branch(k,1); j branch(k,2); b2 branch(k,5); B2(i,i) B2(i,i) - 2*b2; % 对地支路叠加 B2(j,j) B2(j,j) - 2*b2; end end这里为什么B2要减2b2而不是加因为B1定义的是“串联导纳虚部的相反数”而节点导纳矩阵对角元里对地电容的虚部是正的取相反数之后就变成负贡献所以要减掉。很多资料直接写B-Im(Y)从数学上等价但当你说不清符号时不如按代码逐条支路累加来得直观。3.3 PQ分解法的两步交替迭代格式PQ分解法的迭代流程是先固定V解有功方程得到θ再用新的θ计算无功不平衡量并解无功方程得到V然后回到有功迭代。两个方程交替推进直到ΔP和ΔQ都满足容差。步骤操作方程1解有功修正Δθ B⁻¹ (ΔP./V)2更新相角θ θ Δθ3计算Q和ΔQ复用calc_power4解无功修正ΔV B⁻¹ (ΔQ./V)5更新幅值V(PQ) V(PQ) ΔV注意无功迭代只更新PQ节点的电压幅值PV节点电压幅值始终固定在给定值。收敛判断上PQ法通常把有功容差和无功容差分开展示因为两步的计算量本来就不同。这道流程里需要预先从B1、B2中删去平衡节点对应的行和列B2还需要删去PV节点对应的行和列。3.4 在两机五节点规模下PQ法省在哪里牛拉法每轮要组装并求解一个(n-1nPQ)阶的稠密雅可比矩阵而PQ分解法只需在迭代前对B(n-1阶)和B(nPQ阶)做一次三角分解后续反复用。两机五节点系统里这个优势不明显但换到IEEE 118节点单次迭代计算量差距可以到一个数量级。课程设计报告里可以把两种方法的迭代次数和单次迭代浮点运算量分开统计结论会更扎实。另外PQ法对初值的要求略宽松这来自常数矩阵的近似处理但代价是坏初值下不收敛时的报错信息更隐蔽不容易定位是哪一步解耦假设出了问题。4. 用MATLAB把牛拉法和PQ法写成可复现的完整函数4.1 节点与支路数据表组织方式为了让代码既能跑两机五节点又能迁移到其他算例数据组织采用数组而非硬编码。节点数据每行固定为[类型, P_spec, Q_spec, V_spec, theta_spec]类型1为平衡节点、2为PV节点、3为PQ节点支路数据每行固定为[起点, 终点, r, x, b/2]。下面数据块直接复制就能跑bus [ 1 0.00 0.00 1.06 0.00; 2 0.50 0.00 1.00 0.00; 3 -0.55 -0.25 1.00 0.00; 4 -0.50 -0.20 1.00 0.00; 5 -0.30 -0.10 1.00 0.00; ]; branch [ 1 2 0.05 0.20 0.020; 1 3 0.03 0.15 0.015; 2 3 0.04 0.18 0.010; 2 4 0.04 0.20 0.012; 2 5 0.05 0.22 0.012; 3 4 0.06 0.25 0.020; 4 5 0.06 0.20 0.015; ];PV节点的Q_spec初始值可以先给0迭代中计算出的实际Q在越限检查后再处理。V_spec和theta_spec是潮流解的初始猜测两机五节点直接用平启动即可。4.2 导纳矩阵构建函数 buildY导纳矩阵的构建是整个潮流程序的地基。常见做法是逐条支路循环把串联导纳累加到自导纳和互导纳再把半充电电容加到两端节点的自导纳function Y buildY(branch, n) Y zeros(n); for k 1:size(branch,1) i branch(k,1); j branch(k,2); r branch(k,3); x branch(k,4); b2 branch(k,5); y 1/(r 1i*x); Y(i,i) Y(i,i) y 1i*b2; Y(j,j) Y(j,j) y 1i*b2; Y(i,j) Y(i,j) - y; Y(j,i) Y(j,i) - y; end end参数说明Y(i,i)累加的是支路串联导纳和本端半电容Y(i,j)减去的是串联导纳。因为半电容b2在两端各出现一次所以实际对地总导纳是2b2与第2章节点表里的说明一致。4.3 牛拉法主循环牛拉法函数要接收bus、branch和收敛容差tol返回最终V、theta和迭代次数。核心循环加上雅可比组装大约四十行function [V, theta, iter] nr_powerflow(bus, branch, tol) n size(bus,1); Y buildY(branch, n); V bus(:,4); theta bus(:,5); npq find(bus(:,1) 3); npv find(bus(:,1) 2); ns find(bus(:,1) ~ 1); iter 0; while iter 50 [P, Q] calc_power(V, theta, Y); dP bus(:,2) - P; dQ bus(:,3) - Q; dP(1) 0; dQ(1) 0; dQ(npv) 0; if max(abs(dP(ns))) tol max(abs(dQ(npq))) tol break; end H zeros(n); N zeros(n); J zeros(n); L zeros(n); G real(Y); Bm imag(Y); for i 1:n for k 1:n if i ~ k c V(i)*V(k); st sin(theta(i)-theta(k)); ct cos(theta(i)-theta(k)); H(i,k) -c*(G(i,k)*st - Bm(i,k)*ct); N(i,k) c*(G(i,k)*ct Bm(i,k)*st); J(i,k) c*(G(i,k)*ct Bm(i,k)*st); L(i,k) c*(G(i,k)*st - Bm(i,k)*ct); end end H(i,i) -Q(i) - Bm(i,i)*V(i)^2; N(i,i) P(i) G(i,i)*V(i)^2; J(i,i) P(i) - G(i,i)*V(i)^2; L(i,i) Q(i) - Bm(i,i)*V(i)^2; end Jac [H(ns,ns) N(ns,npq); J(npq,ns) L(npq,npq)]; dF [dP(ns); dQ(npq)]; dx Jac \ dF; dTheta zeros(n,1); dVnorm zeros(n,1); dTheta(ns) dx(1:length(ns)); dVnorm(npq) dx(length(ns)1:end); theta theta dTheta; V V V .* dVnorm; iter iter 1; end end逻辑说明dQ(npv)0保证了PV节点的无功不平衡量不进入修正方程但PV节点的无功功率仍会在每次迭代中计算出来用于后续越限检查。组装雅可比时用了全n阶矩阵再索引两机五节点规模下清晰优先大规模系统建议直接按ns和npq的维度预先分配内存。4.4 PQ分解法主循环PQ分解法主循环比牛拉法短但结构上更讲究“常数矩阵只分解一次”。小规模算例可以直接用反斜杠运算符工程上更稳妥的是先做LU分解function [V, theta, iter] pq_powerflow(bus, branch, tol) n size(bus,1); Y buildY(branch, n); [B1, B2] build_B_matrix(branch, n); V bus(:,4); theta bus(:,5); npq find(bus(:,1) 3); npv find(bus(:,1) 2); ns find(bus(:,1) ~ 1); B1 B1(ns,ns); B2 B2(npq,npq); [L1,U1] lu(B1); [L2,U2] lu(B2); iter 0; while iter 100 [P, Q] calc_power(V, theta, Y); dP bus(:,2) - P; dP(1) 0; if max(abs(dP(ns))) tol max(abs(dQ(npq))) tol break; end dTheta zeros(n,1); dTheta(ns) U1 \ (L1 \ (dP(ns)./V(ns))); theta theta dTheta; [P, Q] calc_power(V, theta, Y); dQ bus(:,3) - Q; dQ(1) 0; dQ(npv) 0; dV zeros(n,1); dV(npq) U2 \ (L2 \ (dQ(npq)./V(npq))); V V dV; iter iter 1; end end这里用dP./V而不是直接dP对应修正方程两边的对称化处理。PQ法迭代上限设100次而不是牛拉法的50次因为其收敛速度通常更慢尤其在接近重负荷极限时。LU分解只做一次循环内只做前代回代这就是PQ法计算量优势的代码形态。5. 牛拉法与PQ法的收敛对比、初值与PV节点越限处理5.1 两种方法的迭代行为对比在两机五节点轻负荷算例下牛拉法一般48次收敛PQ法需要612次。这个数字会随负荷水平和收敛容差变化课程设计里不要写死某个值而是用循环统计并给出对比表。单次迭代的计算量差异更值得分析牛拉法要组装4个子块并求解一次稠密线性方程组PQ法则只有常数矩阵的前代回代。指标牛拉法PQ分解法每次迭代成本组装分解n-1nPQ阶雅可比两次低阶前代回代典型迭代次数48612对r/x比敏感度低高r/x过高时发散对初值敏感度较敏感略宽松5.2 初值选择与收敛失败的排查顺序潮流不收敛时先看初值再调参数。平启动是V1.0∠0°、θ0°对两机五节点够用。如果电压初始给成1.1或0.9牛拉法仍能收敛但PV节点无功越限的路径会变化。排查顺序一般是检查导纳矩阵符号 → 单步打印ΔP、ΔQ看是否振荡 → 检查输出功率符号负荷是负注入→ 最后再怀疑解耦假设。一条很实用的小技巧是把容差临时放宽到1e-4看残余量是否单调下降。5.3 PV节点无功越限的处理逻辑PV节点无功功率超出给定上下限时该节点必须转为PQ节点把无功固定在上限或下限电压幅值变为待求量。这个处理要在每次迭代后执行且转换后B2和参与无功迭代的节点集合都会变化if bus(i,1) 2 if Q(i) Qmax || Q(i) Qmin bus(i,1) 3; bus(i,3) min(max(Q(i), Qmin), Qmax); % 重新提取npq、重新分解B2 end end转换后需要重新计算npq索引并重新对B2做LU分解。因为PV节点转PQ会改变无功方程的行列数如果不重新分解后续的无功修正方程维数就会对不上。在课程设计报告中建议把越限处理拆成独立函数主循环只负责检测和调用。6. 结果自洽性检查与迁移到IEEE算例的验收技巧6.1 用功率平衡和节点电压校验收敛结果潮流程序跑通不等于结果可信。收敛后第一件事是核算全系统功率平衡所有节点注入有功之和应等于系统总网损所有节点注入无功之和应等于总无功损耗减去线路充电功率。MATLAB里用sum(real(V.conj(YV)))直接能得到复功率总和实部应为正值且等于网损。电压幅值若有节点低于0.9或高于1.1多半是数据符号有问题或初值不合理而不是算法问题。更底层的验证是对雅可比矩阵做数值差分。对每个待求变量施加微小扰动重新调用calc_power比较差分结果和解析雅可比矩阵元素。这个方法能快速定位符号错误或下标错位尤其在把代码从5节点迁移到其他系统时价值极高。6.2 从五节点迁移到IEEE 14节点的改动清单换算例时只需要替换bus和branch数据但有三处必须重查节点导纳矩阵的稀疏性会变buildY函数中全零矩阵应改为sparse光伏节点数量增多雅可比矩阵的行列索引不再能靠ndgrid遍历建议统一用稀疏索引变压器支路的变比和移相器会影响非对角元需要在buildY里增加一个变比参数。把这些点检查完牛拉法和PQ法的主循环可以不做任何修改直接跑IEEE 14节点系统这也是这套代码设计成函数式而非脚本式的原因。本文还有配套的精品资源点击获取 基于MATLAB的两机五节点潮流计算：牛拉法与PQ分解法实现与对比
简介面向电气工程及相关专业学生的电力系统稳态分析课程设计资料聚焦两机五节点网络的潮流计算适合用于课程设计报告撰写或潮流算法编程入门。资料以PDF文件形式呈现共1个文件压缩包大小约310KB虽体量精简但内容覆盖完整目前已有203人学习下载。文档围绕牛顿拉夫逊法和PQ分解法展开先介绍潮流分析原理与常用方法对比再结合MATLAB环境讲解节点导纳矩阵形成、雅可比矩阵构建、迭代修正等关键步骤并给出两机五节点系统的具体计算程序与运行结果。通过对照高斯-赛德尔法可帮助理解不同潮流算法的收敛速度、适用场景及优缺点进而掌握电力系统稳态分析的核心技能。整份资料结构清晰包含目录、原理简介、程序框图、源码实现和结果分析适合用于课程报告撰写或作为MATLAB潮流计算编程的参考模板。1. 两机五节点潮流计算为什么是牛拉法和PQ法的合适试验场两机五节点这个规模看起来不大却正好卡在“算得明白”和“算得有内容”的交界处。它有两个电源节点、三个负荷节点既能看出平衡节点如何自动吸收全系统不平衡功率也能让牛拉法和PQ分解法在迭代次数、雅可比矩阵稠密度、初值敏感度上显出真实差异。比三节点算例信息量足又比IEEE 14节点系统容易把每一步修正过程打印出来逐行核对。电力系统稳态分析课程设计选这个题目核心工作就是两件事把节点功率方程对电压相角和幅值的偏导关系写对再把两种非线性迭代方法在MATLAB里跑通并解释清楚收敛行为。下文按“方程推导 → 矩阵形成 → 代码实现 → 收敛排错 → 结果验收”的顺序展开中途给出一套可以直接抄进课程设计报告的参数和MATLAB函数骨架。2. 牛拉法极坐标修正方程与两机五节点雅可比矩阵构造2.1 节点类型划分与两机五节点的基准值选择潮流计算的第一步不是写迭代而是先把节点分类和基准值定死。两机五节点系统里节点1设平衡节点节点2设PV节点其余三个设为PQ节点。基准功率取100 MVA基准电压取220 kV全部数据转成标幺值后进入同一个计算框架。节点类型给定条件说明1平衡节点V1.06∠0°有功无功不预先给定承担系统剩余不平衡功率2PV节点P_G0.50 puV1.00无功在-0.200.80 pu内可调3PQ节点P_L0.55 puQ_L0.25 pu恒功率负荷4PQ节点P_L0.50 puQ_L0.20 pu恒功率负荷5PQ节点P_L0.30 puQ_L0.10 pu恒功率负荷注意这里给的是负荷功率程序里需要换算成节点注入功率即注入为负值。支路数据按线路的串联电阻、电抗和半充电电容给出支路r (pu)x (pu)b/2 (pu)1-20.050.200.0201-30.030.150.0152-30.040.180.0102-40.040.200.0122-50.050.220.0123-40.060.250.0204-50.060.200.015b/2表示线路对地电容等值导纳的一半放在支路两端各一次所以一条支路对整个系统贡献的总对地导纳是2×(b/2)b。这个细节在构造节点导纳矩阵时最容易错课程设计报告里建议单独写一段说明。2.2 极坐标功率方程与注入不平衡量的计算极坐标潮流方程把每个节点的注入有功和无功写成电压幅值与相角的函数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来自节点导纳矩阵的实部与虚部。迭代时需要计算“给定注入”与“当前电压下计算注入”之差这就是功率不平衡量ΔP和ΔQ。用MATLAB可以非常简洁地写这段计算function [P, Q] calc_power(V, theta, Y) % 根据当前电压幅值和相角计算各节点注入功率 n length(V); E V .* exp(1i * theta); % 节点电压相量 I Y * E; % 节点注入电流相量 S E .* conj(I); % 复功率 P real(S); Q imag(S); end这里直接用了相量形式的复功率公式SV·conj(I)结果和逐项累加功率方程完全一致但代码更短适合在报告中作为“功率计算模块”单独贴出。注意输入V是幅值向量theta是弧度制的相角向量Y是n×n节点导纳矩阵。2.3 牛拉法雅可比矩阵四个子块的表达式牛拉法极坐标形式的修正方程写成块矩阵[ΔP; ΔQ] [H N; J L] × [Δθ; ΔV/V]这是“ΔV/V”形式代码里最后再把ΔV/V乘回V。四个子块的元素必须和功率方程一一对应求偏导下面是教材级结论直接用于MATLAB组装子块非对角元 (i≠k)对角元 (ik)H∂P/∂θ−V_i V_k (G_ik sinθ_ik − B_ik cosθ_ik)−Q_i − B_ii V_i²N∂P/∂V×VV_i V_k (G_ik cosθ_ik B_ik sinθ_ik)P_i G_ii V_i²J∂Q/∂θV_i V_k (G_ik cosθ_ik B_ik sinθ_ik)P_i − G_ii V_i²L∂Q/∂V×VV_i V_k (G_ik sinθ_ik − B_ik cosθ_ik)Q_i − B_ii V_i²注意代码实现时符号约定要统一。这里采用ΔPHΔθNΔV/V、ΔQJΔθLΔV/V的形式那么不平衡量必须按“给定注入功率 − 计算注入功率”来定义二者符号绑定不能只取一个负号。2.4 牛拉法迭代流程与收敛判据的组织顺序牛拉法单次迭代的流程固定为计算ΔP、ΔQ → 判断是否小于容差 → 组装四个子块 → 解修正方程 → 更新θ和V。程序里用idx记录除平衡节点外的所有节点用npq记录PQ节点这样修正方程的维度和映射关系才清楚。ns find(bus(:,1) ~ 1); % 非平衡节点 npq find(bus(:,1) 3); % PQ节点 npv find(bus(:,1) 2); % PV节点 Jac [H(ns,ns) N(ns,npq); J(npq,ns) L(npq,npq)]; dF [dP(ns); dQ(npq)]; dx Jac \ dF; dTheta zeros(n,1); dVnorm zeros(n,1); dTheta(ns) dx(1:length(ns)); dVnorm(npq) dx(length(ns)1:end); theta theta dTheta; V V V .* dVnorm;这段代码的要点在于矩阵索引。H的列必须对应所有参与θ迭代的节点而N的列只对应PQ节点因为PV节点的电压幅值不迭代。求解之后更新V时用的是逐元素乘法即VVV.*dVnorm反过来如果在组装前给N乘了一次V这里就不能再乘。收敛判据我一般用max(abs(dP(ns)))和max(abs(dQ(npq)))同时小于1e-6单纯看ΔP会漏掉无功不收敛的情况。3. 从牛拉法到PQ分解B与B矩阵的形成与迭代格式3.1 解耦假设成立的两条线路条件PQ分解法不是凭空出现的简化它基于电力系统高压输电网的两个物理事实正常运行时节点电压相角差很小所以sinθ_ik≈θ_ik、cosθ_ik≈1线路电抗远大于电阻无功功率对相角不敏感有功功率对电压幅值也不敏感。把这两个条件代入牛拉法修正方程N和J两个子块就可以忽略修正方程解耦成有功迭代和无功迭代两部分。需要强调的是这种解耦在配电网络里往往不成立因为配电网r/x比值高强行用PQ法会造成迭代发散或收敛到错误解。这一点课程设计答辩时经常被问到要提前准备好说法。3.2 B与B矩阵的元素构成与符号约定PQ分解法核心是两个常数矩阵B和B。这里直接给出与MATLAB代码匹配的形成规则矩阵参与元素忽略元素用途B串联支路电抗形成的导纳虚部电阻、对地电容用于ΔP/BB的基础上叠加对地电容导纳电阻用于ΔQ/符号约定上B和B都取节点导纳矩阵虚部的相反数所以对角元为正、非对角元为负。下面这段代码同时构造两个矩阵注意B2在B1基础上把每条线路两端节点的对地电容2b2加回对角function [B1, B2] build_B_matrix(branch, n) % 构造PQ分解法需要的B矩阵和B矩阵 B1 zeros(n); for k 1:size(branch,1) i branch(k,1); j branch(k,2); r branch(k,3); x branch(k,4); b2 branch(k,5); y_series 1/(r 1i*x); b -imag(y_series); % 串联导纳虚部取反 B1(i,i) B1(i,i) b; B1(j,j) B1(j,j) b; B1(i,j) B1(i,j) - b; B1(j,i) B1(j,i) - b; end B2 B1; for k 1:size(branch,1) i branch(k,1); j branch(k,2); b2 branch(k,5); B2(i,i) B2(i,i) - 2*b2; % 对地支路叠加 B2(j,j) B2(j,j) - 2*b2; end end这里为什么B2要减2b2而不是加因为B1定义的是“串联导纳虚部的相反数”而节点导纳矩阵对角元里对地电容的虚部是正的取相反数之后就变成负贡献所以要减掉。很多资料直接写B-Im(Y)从数学上等价但当你说不清符号时不如按代码逐条支路累加来得直观。3.3 PQ分解法的两步交替迭代格式PQ分解法的迭代流程是先固定V解有功方程得到θ再用新的θ计算无功不平衡量并解无功方程得到V然后回到有功迭代。两个方程交替推进直到ΔP和ΔQ都满足容差。步骤操作方程1解有功修正Δθ B⁻¹ (ΔP./V)2更新相角θ θ Δθ3计算Q和ΔQ复用calc_power4解无功修正ΔV B⁻¹ (ΔQ./V)5更新幅值V(PQ) V(PQ) ΔV注意无功迭代只更新PQ节点的电压幅值PV节点电压幅值始终固定在给定值。收敛判断上PQ法通常把有功容差和无功容差分开展示因为两步的计算量本来就不同。这道流程里需要预先从B1、B2中删去平衡节点对应的行和列B2还需要删去PV节点对应的行和列。3.4 在两机五节点规模下PQ法省在哪里牛拉法每轮要组装并求解一个(n-1nPQ)阶的稠密雅可比矩阵而PQ分解法只需在迭代前对B(n-1阶)和B(nPQ阶)做一次三角分解后续反复用。两机五节点系统里这个优势不明显但换到IEEE 118节点单次迭代计算量差距可以到一个数量级。课程设计报告里可以把两种方法的迭代次数和单次迭代浮点运算量分开统计结论会更扎实。另外PQ法对初值的要求略宽松这来自常数矩阵的近似处理但代价是坏初值下不收敛时的报错信息更隐蔽不容易定位是哪一步解耦假设出了问题。4. 用MATLAB把牛拉法和PQ法写成可复现的完整函数4.1 节点与支路数据表组织方式为了让代码既能跑两机五节点又能迁移到其他算例数据组织采用数组而非硬编码。节点数据每行固定为[类型, P_spec, Q_spec, V_spec, theta_spec]类型1为平衡节点、2为PV节点、3为PQ节点支路数据每行固定为[起点, 终点, r, x, b/2]。下面数据块直接复制就能跑bus [ 1 0.00 0.00 1.06 0.00; 2 0.50 0.00 1.00 0.00; 3 -0.55 -0.25 1.00 0.00; 4 -0.50 -0.20 1.00 0.00; 5 -0.30 -0.10 1.00 0.00; ]; branch [ 1 2 0.05 0.20 0.020; 1 3 0.03 0.15 0.015; 2 3 0.04 0.18 0.010; 2 4 0.04 0.20 0.012; 2 5 0.05 0.22 0.012; 3 4 0.06 0.25 0.020; 4 5 0.06 0.20 0.015; ];PV节点的Q_spec初始值可以先给0迭代中计算出的实际Q在越限检查后再处理。V_spec和theta_spec是潮流解的初始猜测两机五节点直接用平启动即可。4.2 导纳矩阵构建函数 buildY导纳矩阵的构建是整个潮流程序的地基。常见做法是逐条支路循环把串联导纳累加到自导纳和互导纳再把半充电电容加到两端节点的自导纳function Y buildY(branch, n) Y zeros(n); for k 1:size(branch,1) i branch(k,1); j branch(k,2); r branch(k,3); x branch(k,4); b2 branch(k,5); y 1/(r 1i*x); Y(i,i) Y(i,i) y 1i*b2; Y(j,j) Y(j,j) y 1i*b2; Y(i,j) Y(i,j) - y; Y(j,i) Y(j,i) - y; end end参数说明Y(i,i)累加的是支路串联导纳和本端半电容Y(i,j)减去的是串联导纳。因为半电容b2在两端各出现一次所以实际对地总导纳是2b2与第2章节点表里的说明一致。4.3 牛拉法主循环牛拉法函数要接收bus、branch和收敛容差tol返回最终V、theta和迭代次数。核心循环加上雅可比组装大约四十行function [V, theta, iter] nr_powerflow(bus, branch, tol) n size(bus,1); Y buildY(branch, n); V bus(:,4); theta bus(:,5); npq find(bus(:,1) 3); npv find(bus(:,1) 2); ns find(bus(:,1) ~ 1); iter 0; while iter 50 [P, Q] calc_power(V, theta, Y); dP bus(:,2) - P; dQ bus(:,3) - Q; dP(1) 0; dQ(1) 0; dQ(npv) 0; if max(abs(dP(ns))) tol max(abs(dQ(npq))) tol break; end H zeros(n); N zeros(n); J zeros(n); L zeros(n); G real(Y); Bm imag(Y); for i 1:n for k 1:n if i ~ k c V(i)*V(k); st sin(theta(i)-theta(k)); ct cos(theta(i)-theta(k)); H(i,k) -c*(G(i,k)*st - Bm(i,k)*ct); N(i,k) c*(G(i,k)*ct Bm(i,k)*st); J(i,k) c*(G(i,k)*ct Bm(i,k)*st); L(i,k) c*(G(i,k)*st - Bm(i,k)*ct); end end H(i,i) -Q(i) - Bm(i,i)*V(i)^2; N(i,i) P(i) G(i,i)*V(i)^2; J(i,i) P(i) - G(i,i)*V(i)^2; L(i,i) Q(i) - Bm(i,i)*V(i)^2; end Jac [H(ns,ns) N(ns,npq); J(npq,ns) L(npq,npq)]; dF [dP(ns); dQ(npq)]; dx Jac \ dF; dTheta zeros(n,1); dVnorm zeros(n,1); dTheta(ns) dx(1:length(ns)); dVnorm(npq) dx(length(ns)1:end); theta theta dTheta; V V V .* dVnorm; iter iter 1; end end逻辑说明dQ(npv)0保证了PV节点的无功不平衡量不进入修正方程但PV节点的无功功率仍会在每次迭代中计算出来用于后续越限检查。组装雅可比时用了全n阶矩阵再索引两机五节点规模下清晰优先大规模系统建议直接按ns和npq的维度预先分配内存。4.4 PQ分解法主循环PQ分解法主循环比牛拉法短但结构上更讲究“常数矩阵只分解一次”。小规模算例可以直接用反斜杠运算符工程上更稳妥的是先做LU分解function [V, theta, iter] pq_powerflow(bus, branch, tol) n size(bus,1); Y buildY(branch, n); [B1, B2] build_B_matrix(branch, n); V bus(:,4); theta bus(:,5); npq find(bus(:,1) 3); npv find(bus(:,1) 2); ns find(bus(:,1) ~ 1); B1 B1(ns,ns); B2 B2(npq,npq); [L1,U1] lu(B1); [L2,U2] lu(B2); iter 0; while iter 100 [P, Q] calc_power(V, theta, Y); dP bus(:,2) - P; dP(1) 0; if max(abs(dP(ns))) tol max(abs(dQ(npq))) tol break; end dTheta zeros(n,1); dTheta(ns) U1 \ (L1 \ (dP(ns)./V(ns))); theta theta dTheta; [P, Q] calc_power(V, theta, Y); dQ bus(:,3) - Q; dQ(1) 0; dQ(npv) 0; dV zeros(n,1); dV(npq) U2 \ (L2 \ (dQ(npq)./V(npq))); V V dV; iter iter 1; end end这里用dP./V而不是直接dP对应修正方程两边的对称化处理。PQ法迭代上限设100次而不是牛拉法的50次因为其收敛速度通常更慢尤其在接近重负荷极限时。LU分解只做一次循环内只做前代回代这就是PQ法计算量优势的代码形态。5. 牛拉法与PQ法的收敛对比、初值与PV节点越限处理5.1 两种方法的迭代行为对比在两机五节点轻负荷算例下牛拉法一般48次收敛PQ法需要612次。这个数字会随负荷水平和收敛容差变化课程设计里不要写死某个值而是用循环统计并给出对比表。单次迭代的计算量差异更值得分析牛拉法要组装4个子块并求解一次稠密线性方程组PQ法则只有常数矩阵的前代回代。指标牛拉法PQ分解法每次迭代成本组装分解n-1nPQ阶雅可比两次低阶前代回代典型迭代次数48612对r/x比敏感度低高r/x过高时发散对初值敏感度较敏感略宽松5.2 初值选择与收敛失败的排查顺序潮流不收敛时先看初值再调参数。平启动是V1.0∠0°、θ0°对两机五节点够用。如果电压初始给成1.1或0.9牛拉法仍能收敛但PV节点无功越限的路径会变化。排查顺序一般是检查导纳矩阵符号 → 单步打印ΔP、ΔQ看是否振荡 → 检查输出功率符号负荷是负注入→ 最后再怀疑解耦假设。一条很实用的小技巧是把容差临时放宽到1e-4看残余量是否单调下降。5.3 PV节点无功越限的处理逻辑PV节点无功功率超出给定上下限时该节点必须转为PQ节点把无功固定在上限或下限电压幅值变为待求量。这个处理要在每次迭代后执行且转换后B2和参与无功迭代的节点集合都会变化if bus(i,1) 2 if Q(i) Qmax || Q(i) Qmin bus(i,1) 3; bus(i,3) min(max(Q(i), Qmin), Qmax); % 重新提取npq、重新分解B2 end end转换后需要重新计算npq索引并重新对B2做LU分解。因为PV节点转PQ会改变无功方程的行列数如果不重新分解后续的无功修正方程维数就会对不上。在课程设计报告中建议把越限处理拆成独立函数主循环只负责检测和调用。6. 结果自洽性检查与迁移到IEEE算例的验收技巧6.1 用功率平衡和节点电压校验收敛结果潮流程序跑通不等于结果可信。收敛后第一件事是核算全系统功率平衡所有节点注入有功之和应等于系统总网损所有节点注入无功之和应等于总无功损耗减去线路充电功率。MATLAB里用sum(real(V.conj(YV)))直接能得到复功率总和实部应为正值且等于网损。电压幅值若有节点低于0.9或高于1.1多半是数据符号有问题或初值不合理而不是算法问题。更底层的验证是对雅可比矩阵做数值差分。对每个待求变量施加微小扰动重新调用calc_power比较差分结果和解析雅可比矩阵元素。这个方法能快速定位符号错误或下标错位尤其在把代码从5节点迁移到其他系统时价值极高。6.2 从五节点迁移到IEEE 14节点的改动清单换算例时只需要替换bus和branch数据但有三处必须重查节点导纳矩阵的稀疏性会变buildY函数中全零矩阵应改为sparse光伏节点数量增多雅可比矩阵的行列索引不再能靠ndgrid遍历建议统一用稀疏索引变压器支路的变比和移相器会影响非对角元需要在buildY里增加一个变比参数。把这些点检查完牛拉法和PQ法的主循环可以不做任何修改直接跑IEEE 14节点系统这也是这套代码设计成函数式而非脚本式的原因。本文还有配套的精品资源点击获取