简介这份资源是电力系统稳态分析课程的完整设计报告面向电气工程及相关专业学生、课程设计选课者与需要复习潮流计算的读者聚焦“两机五节点”网络这一典型算例解决牛顿-拉夫逊法与PQ分解法的原理理解、程序实现与结果对比问题适合具备电路与电力系统基础、正在做课设或准备答辩的中高年级学生参考。压缩包共1个doc文件约271KB即整份报告书正文内含摘要、原理简介、MATLAB简介、程序与结果、总结及参考文献等章节。报告先比较高斯-赛德尔法、牛顿法与PQ分解法的收敛性与内存占用差异再给出设计资料与参数、牛顿法程序框图、MATLAB程序编写步骤与运行结果并附PQ法程序实现最后对比两种算法的迭代次数与精度表现。目前已有94人学习下载。读者可据此获得可直接借鉴的建模思路、雅可比矩阵构建与迭代代码框架、参数填写格式以及算法性能对比结论便于快速复现两机五节点潮流计算并完成自己的课设报告。1. 两机五节点系统与潮流计算任务拆解一份名为《两机五节点网络潮流计算方法牛拉法和 pq 法电力系统稳态分析课程设计报告书.doc》的资料通常同时装着三样东西直角坐标下的修正方程推导、一套能跑出k1,2,3的 MATLAB 主程序、以及一张只有五个节点却让不少人卡住一整天的等值电路图。卡住的地方多半不是公式本身而是把接线图翻译成节点导纳矩阵时节点顺序怎么排以及雅可比矩阵那 8 个分块往哪放。两机五节点指的是一台平衡机加一台电源、五条母线的小系统节点 1 是平衡节点节点 2 到节点 5 按 PQ 节点处理。规模不大但牛顿-拉夫逊法的全部环节一个不缺导纳阵生成、初值设定、雅可比矩阵组装、修正量求解、收敛判定、平衡节点功率回代。算得快、看得清正是它适合当课程设计的原因。适合读这份材料的人有三类做电力系统稳态分析课程设计的学生、需要复现小规模潮流程序验证算法的工程师、以及想对比牛拉法与 PQ 分解法收敛特性的从业者。下面的拆解从建模开始一路走到结果校验。2. 从系统接线图到节点导纳矩阵的建模2.1 节点编号为什么要把平衡节点挪到末尾原始接线图里节点 1 是平衡节点、节点 2 到 5 是 PQ 节点但程序第一行注释写着「为了使节点按照先 PQ、再 PV、最后平衡节点的次序编号以便与公式对照节点 1 与节点 5 对调」。这不是多此一举。牛顿-拉夫逊法在直角坐标下只对 PQ 节点列写有功、无功两个不平衡方程平衡节点的电压幅值和相角是已知量不参与迭代。如果平衡节点编号夹在中间求和不平衡量、组装雅可比矩阵时每次都要跳过某个下标循环里就得多写一层if出错概率成倍上升。把平衡节点放到最后一位PQ 节点就是连续的1:m求和可以写成for n1:5配合节点自身项的修正代码和教材公式一一对应。对调之后要同步改三处阻抗矩阵行列次序、节点功率给定向量S的元素次序、以及最终输出时把电压向量还原回原始编号。很多人程序跑出来的电压幅值看着「差不多」但角度符号全反了基本就是漏掉了最后这一步还原。2.2 支路阻抗填写与导纳矩阵生成两机五节点的支路数据是一组纯复数阻抗不含变压器变比和对地支路所以导纳矩阵可以完全由阻抗矩阵推导不需要额外处理非标准变比。按对调后的编号支路参数如下表。支路阻抗 (p.u.)串联导纳 y1/Z对应 Y 阵元素1-20.04 j0.122.5000 − j7.5000Y(1,2) −2.5 j7.51-40.08 j0.241.2500 − j3.7500Y(1,4) −1.25 j3.752-30.06 j0.181.6667 − j5.0000Y(2,3) −1.6667 j5.02-40.06 j0.181.6667 − j5.0000Y(2,4) −1.6667 j5.02-50.02 j0.065.0000 − j15.0000Y(2,5) −5.0 j15.03-40.01 j0.0310.000 − j30.000Y(3,4) −10.0 j30.03-50.08 j0.241.2500 − j3.7500Y(3,5) −1.25 j3.75生成代码分两步走先逐元素求倒数得到互导纳再按「对角为同行之和、非对角取负」组装节点导纳矩阵。% 节点阻抗矩阵零元素表示无直接支路连接 Z [0, 0.040.12i, 0, 0.080.24i, 0; 0.040.12i, 0, 0.060.18i, 0.060.18i, 0.020.06i; 0, 0.060.18i, 0, 0.010.03i, 0.080.24i; 0.080.24i, 0.060.18i, 0.010.03i, 0, 0; 0, 0.020.06i, 0.080.24i, 0, 0]; % 第一步逐元素取倒数得到互导纳 y零阻抗位置保持零 for m 1:5 for n 1:5 if Z(m,n) 0 y(m,n) 0; % 无支路互导纳为零 else y(m,n) 1 / Z(m,n); % 复数求逆注意用点除以外的标量除 end end end % 第二步对角元累加本行全部串联导纳非对角元取互导纳的相反数 for m 1:5 for n 1:5 if m ~ n Y(m,n) -y(m,n); else Y(m,n) sum(y(m,:)); % 自导纳 相连支路导纳之和 end end end G real(Y); % 电导矩阵 B imag(Y); % 电纳矩阵sum(y(m,:))这一句是关键。它默认把整行加了一遍而y(m,m)位置本来就是 0阻抗矩阵对角线上没有支路所以不需要再减去自阻抗项。如果哪天加了并联电容或者对地支路Z(m,m)就不为零这一行必须改成sum(y(m,:)) - y(m,m)再单独加并联支路导纳否则自导纳会算大一倍。跑出来的 Y 矩阵第一行是3.7500 - 11.2500i, -2.5000 7.5000i, 0, -1.2500 3.7500i, 0对角元 3.75 − j11.25 正好等于 y12 y14可以拿来当第一道自检。Y 矩阵是对称阵程序里可以顺手加一句norm(Y - Y.) 1e-10验证对称性非对称基本就是填表时漏了某条支路的反向。2.3 节点类型划分与初值设定节点 1 到 4 是 PQ 节点有功、无功注入给定节点 5 是平衡节点电压幅值 1.06、相角 0 给定。所有节点的有功无功注入按发电机为正、负荷为负填写功率向量写成% 节点注入功率PQ 节点给定额定值平衡节点填 0不参与不平衡量计算 S [-0.60-0.10i; % 节点1负荷 0.200.20i; % 节点2发电 本地负荷 -0.45-0.15i; % 节点3负荷 -0.40-0.05i; % 节点4负荷 0]; % 节点5平衡节点占位 % 平启动所有 PQ 节点电压实部取 1.0、虚部取 0平衡节点固定 U [1, 1, 1, 1, 1.06]; E real(U); F imag(U);平启动flat start在这个规模下通常能收敛因为各节点电压偏离额定值不远、相角差也不大。如果换成重负荷或者长距离输电线路的系统平启动可能反复振荡常见做法是先用高斯-赛德尔法迭代两三次把结果当作牛拉法的初值也就是所谓的「平启动 预热」组合。3. 直角坐标牛顿-拉夫逊法的雅可比矩阵与迭代实现3.1 功率不平衡量的节点级计算直角坐标下把电压写成U E jF节点注入功率的展开式是$$P_i E_i\sum_j (G_{ij}E_j - B_{ij}F_j) F_i\sum_j (G_{ij}F_j B_{ij}E_j)$$$$Q_i F_i\sum_j (G_{ij}E_j - B_{ij}F_j) - E_i\sum_j (G_{ij}F_j B_{ij}E_j)$$程序里把内层求和拆成两个临时数组Ai和Bi这两个数组在后面对角元素计算时还会被复用不要用完就丢。for m 1:4 Pt 0; Qt 0; for n 1:5 % 每个 n 单独算一项最后累加 Pt(n) E(m)*(G(m,n)*E(n) - B(m,n)*F(n)) ... F(m)*(G(m,n)*F(n) B(m,n)*E(n)); Qt(n) F(m)*(G(m,n)*E(n) - B(m,n)*F(n)) ... - E(m)*(G(m,n)*F(n) B(m,n)*E(n)); end dP(m) P(m) - sum(Pt); % 有功不平衡量 ΔP dQ(m) Q(m) - sum(Qt); % 无功不平衡量 ΔQ end不平衡量的符号约定要盯死dP P给定 - P计算。如果反过来写成计算值减给定值修正量方向就会反向迭代要么发散要么在两点之间来回跳。收敛判据用的是修正量模值的最大值C max(abs(dU))源码里取的阈值是1e-5。3.2 雅可比矩阵四个分块的对角与非对角元素雅可比矩阵按[H N; J L]排布对应[ΔP; ΔQ] [H N; J L]·[ΔF; ΔE]。非对角元素公式简单分块非对角元素表达式代码写法H∂P_i/∂F_j G_ij·F_i − B_ij·E_i-B(m,n)*E(m)G(m,n)*F(m)N∂P_i/∂E_j G_ij·E_i B_ij·F_iG(m,n)*E(m)B(m,n)*F(m)J∂Q_i/∂F_j −B_ij·F_i − G_ij·E_i-B(m,n)*F(m)-G(m,n)*E(m)L∂Q_i/∂E_j G_ij·F_i − B_ij·E_iG(m,n)*F(m)-B(m,n)*E(m)对角元素多出一项自身电压的贡献源码里的写法是借助前面已经算好的Ai、Bi做一次整体求和再修正for m 1:4 for n 1:5 Bi(n) G(m,n)*F(n) B(m,n)*E(n); Ai(n) G(m,n)*E(n) - B(m,n)*F(n); % 与上一节重复利用 end H(m,m) sum(Bi) - (B(m,m)*E(m) G(m,m)*F(m)) 2*G(m,m)*F(m); N(m,m) sum(Ai) - (G(m,m)*E(m) - B(m,m)*F(m)) 2*G(m,m)*E(m); J(m,m) -2*B(m,m)*F(m) sum(Ai) - (G(m,m)*E(m) - B(m,m)*F(m)); L(m,m) -2*B(m,m)*E(m) - (sum(Bi) - (B(m,m)*E(m) G(m,m)*F(m))); endsum(Bi) - (B(m,m)*E(m) G(m,m)*F(m))这一步等价于把 jm 的自项从整行和里剔出去得到的正是 Σ_{j≠m} 的贡献再加上 2·G_mm·F_m 就是完整的对角元。写成这样而不是直接展开公式好处是逐项对得上教材里的推导出题老师一眼能看出你确实理解了公式来源而不是从网上抄了一段。3.3 主循环、修正量求解与收敛判定雅可比矩阵按「对角块先填、再填非对角块」的顺序组装然后左除求解修正量。% 对角块填充索引映射 2m-1 → ΔP/ΔF2m → ΔQ/ΔE for m 1:4 JJ(2*m-1, 2*m-1) H(m,m); JJ(2*m-1, 2*m ) N(m,m); JJ(2*m, 2*m-1) J(m,m); JJ(2*m, 2*m ) L(m,m); end % 非对角块填充注意跳过 m n 的位置 for m 1:4 for n 1:4 if m ~ n H(m,n) -B(m,n)*E(m) G(m,n)*F(m); N(m,n) G(m,n)*E(m) B(m,n)*F(m); J(m,n) -B(m,n)*F(m) - G(m,n)*E(m); L(m,n) G(m,n)*F(m) - B(m,n)*E(m); JJ(2*m-1, 2*n-1) H(m,n); JJ(2*m-1, 2*n ) N(m,n); JJ(2*m, 2*n-1) J(m,n); JJ(2*m, 2*n ) L(m,n); end end end for m 1:4 PQ(2*m-1) dP(m); PQ(2*m) dQ(m); end dU JJ \ PQ.; % 推荐左除避免 inv 带来的数值误差 C max(abs(dU)); % 收敛判据修正量最大模值 for n 1:4 F(n) F(n) dU(2*n-1); % 虚部修正 E(n) E(n) dU(2*n); % 实部修正 end索引映射2m-1 → F方向、2m → E方向是最容易写反的地方。一旦把ΔF和ΔE的更新顺序对调程序照样能跑完但电压实部会越来越大、虚部趋近于零最后报出个看似收敛的荒谬结果。判断方法很简单正常情况下负荷节点的电压实部应该在 0.99~1.05 之间虚部是小负数。源码里用的是inv(JJ)*PQ从 8×8 的规模看差别不大但改成JJ \ PQ.有两个好处一是数值稳定性更好二是迭代次数一多显式求逆的累积误差会拖慢收敛。这一改动不影响任何结果数值属于可以无脑换的写法。按同一组数据跑迭代过程大致是第 1 次修正量最大模值约 0.108第 2 次降到 0.013第 3 次进入 1e-4 量级第 4 次即满足 1e-5 判据。平方收敛特性在这个小系统上体现得很直观——修正量每轮大约按上一轮结果的平方缩小。3.4 平衡节点注入功率的回代PQ 节点的电压算完平衡节点的功率还没着落。它由节点电压和导纳矩阵反推% 平衡节点电流I Σ Y(5,m)·U(m) for m 1:5 I(m) Y(5,m) * U(m); end S(5) U(5) * sum(conj(I)); % 复功率 S U·I*conj是共轭漏掉它算出来的就是视在功率而非有功无功分量平衡节点的功率会出现虚部符号翻转。回代得到的 S(5) 绝对值原则上等于全网损耗加上全部负荷可以用它做第一道能量守恒校验所有节点注入功率之和应为零忽略线路损耗时为近似成立含损耗时应满足sum(S) ≈ 0。4. PQ 分解法的解耦假设与程序改造4.1 P-θ 与 Q-V 强耦合的物理依据牛拉法每轮都要重算整个 2(n−1) 阶雅可比矩阵规模一上来存储和求逆的开销就顶不住。PQ 分解法抓住的是电力系统一条经验事实有功功率的平衡主要由节点电压相角决定无功功率的平衡主要由电压幅值决定两者的交叉敏感度弱。具体到公式上在极坐标雅可比矩阵中∂P/∂V和∂Q/∂θ这两块元素在高压电网中因为线路电抗远大于电阻X/R 比值大数值上比对角块小得多。把它们置零原来的 2n 阶方程就拆成了两个低阶子问题$$\Delta P/V B\Delta\theta, \qquad \Delta Q/V B\Delta V$$两个子问题各自阶数是 n_PQ 和 n_PQ总的求解规模降到原来的一半以下而且 B′、B″ 在迭代过程中是常数矩阵只需要在循环外做一次三角分解循环体内只做前代回代。这就是它成为目前计算速度最快的潮流算法之一的原因。4.2 B′ 与 B″ 矩阵的两种构成方式B′ 和 B″ 不是简单地把 Y 矩阵的虚部拆开具体构成有两种流派严格派B′ 取网络导纳矩阵虚部但构造时忽略对地支路和非标准变比的影响B″ 同样取虚部但保留对地支路。简化派XB 型B′ 直接用支路电抗的倒数-1/x_ij构成B″ 用导纳矩阵虚部构成。对两机五节点这种没有变压器、没有并联补偿的小系统两派结果几乎一样。程序里推荐用简化派因为支路电抗直接来自阻抗矩阵的虚部不用再处理 Y 矩阵里混在一起的接地项。% 用支路电抗倒数构成 B只对有支路的节点对赋值 Bp zeros(4,4); % PQ 节点阶数 for m 1:4 for n 1:4 if m ~ n imag(Z(m,n)) ~ 0 Bp(m,n) 1 / imag(Z(m,n)); % 非对角元素为支路电抗倒数 Bp(m,m) Bp(m,m) - 1 / imag(Z(m,n)); % 对角累加负值 end end endBp(m,m)的累加符号和非对角元相反这是由节点自导纳的定义决定的。如果对角元符号写反程序会在第二轮迭代就开始剧烈振荡dθ的模值不降反升。验证方法sum(Bp, 2)输出的每行对角线主导程度应该明显|B_ii| 大于同行其他元素绝对值之和不满足就说明符号写反了。4.3 共用同一套原始数据的程序改造牛拉法和 PQ 分解法共用完全相同的 Z 矩阵、Y 矩阵、S 向量、电压初值区别只在迭代内核。常见做法是写成两个独立的 m 文件或者一个文件里两个 function保证数据不分叉。改造时注意三点坐标形式变了牛拉法用直角坐标E jFPQ 法用极坐标V∠θ。初值要从直角坐标转换一次V abs(E jF)、theta angle(E jF)。修正量含义变了Δθ是弧度、ΔV是标量收敛判据要分别对待常用写法是max(abs(dtheta)) 1e-5与max(abs(dV./V)) 1e-5同时满足。B′、B″ 只在循环外算一次把它们放进内层循环会让 PQ 法丧失速度优势迭代次数看着差不多但单次迭代耗时接近牛拉法。4.4 两种方法的迭代行为对比以同一组数据实测牛拉法在 4 次内收敛PQ 分解法迭代次数略多但每次迭代的运算量只有牛拉法的一半左右总耗时反而更短。最终收敛结果必须一致这是交叉验证的核心。对比项牛顿-拉夫逊法直角坐标PQ 分解法极坐标迭代矩阵阶数8×82×PQ 节点数4×4 与 4×4 两个子问题矩阵是否随迭代变化每次迭代重新计算B′、B″ 为常数阵存储开销高雅可比元素约 64 个低两个子阵共 32 个典型收敛次数3~5 次6~10 次单次迭代耗时较大很小对初值敏感度较高病态系统需良好初值相对宽松收敛后的节点电压幅值标幺值与相角度如下这是对照两种算法的基准。节点电压实部电压虚部幅值相角 (°)10.9961−0.10441.0015−5.9821.0354−0.04771.0365−2.6431.0052−0.08451.0087−4.8141.0032−0.09011.0072−5.1351.060001.06000两种方法必须收敛到同一组数值误差在 1e-6 量级内。如果出现系统性偏差比如所有相角差零点几度先查 B′ 的对角符号再查极坐标与直角坐标的相位基准是否对齐——直角坐标的虚部是相对于实轴的分量极坐标的相角是相对于同步旋转参考轴的两者在初值转换时必须同源。5. 结果校验与收敛异常排查的实用手法5.1 功率平衡与支路潮流反查电压算出来只是第一步真正能说服人的是能量守恒。校验代码有两行% 全网复功率之和应接近零平衡节点回代后 S_total sum(S) ; % 近似为 0残差反映线路损耗计算是否正确 % 逐条支路潮流核对 for m 1:5 for n 1:5 if Y(m,n) ~ 0 m ~ n S_branch(m,n) U(m) * conj((U(m)-U(n)) * Y(m,n)); % m → n 方向潮流 end end endS_branch(m,n) S_branch(n,m)应该等于该支路的损耗|I|²·R是个小的正实数。如果算出负损耗说明阻抗矩阵的虚部符号或者电压虚部符号搞反了。5.2 不收敛时的排查顺序按这个顺序查能覆盖九成以上的问题先看修正量是震荡还是单调增。震荡通常是雅可比矩阵对角块符号错单调增通常是导纳矩阵自导纳算小了。检查功率给定值与节点类型的对应关系。上机时最常见的是把平衡节点的 S 值填成了非零程序不报错但结果全乱。打印第一轮的雅可比矩阵。对角占优程度可以肉眼判断如果非对角元素比对角还大说明节点编号顺序和公式里的索引对不上。把收敛判据从 1e-5 放到 1e-3 试跑。如果放宽带宽能收敛说明算法本身没问题只是初值或数据精度不足。检查 P 和 Q 单位是否统一为标幺值。混用有名值会让不平衡量差好几个数量级表现为「修正量第一轮就很大」。5.3 从 .doc 报告里提取可运行代码的技巧报告文档里的源码经常有换行符丢失、全角空格混入、if mn else ... end结构被压成一行的情况直接复制粘贴到 MATLAB 里会报语法错。稳妥的流程是先粘到支持正则替换的编辑器里把全角空格\u3000和中文标点批量换成半角再检查for/end配对是否为偶数最后把else单独拆行让if ... else ... end结构完整。数据部分注意一点报告里的 Y 矩阵经常只印到小数点后四位直接当输入用会引入截断误差。正确做法是从 Z 矩阵重新生成 Y把打印出的 Y 矩阵只当作校验用的参考值两者差异在 1e-6 以上就说明报告里的数据抄错了。最后给一个实用的自检习惯每次改完参数先把k 1那轮的dP、dQ打印出来和手算的功率不平衡量对一遍。手算时只需要代一个节点的两个方程两三分钟的事但能省掉后面反复跑程序的半小时。本文还有配套的精品资源点击获取