MATLAB潮流计算课程设计：节点导纳矩阵与牛顿-拉夫逊法

简介：一份围绕电力系统潮流计算的课程设计文档，重点讲解基于MATLAB的牛顿—拉夫逊法潮流计算实现，适合电气工程专业学生完成算法类课程设计或初步接触潮流计算时参考。资源为单个doc文档，共1个文件，压缩包约346KB，内容涵盖设计目的与要求、题目分析、节点导纳矩阵构建、雅可比矩阵形成、迭代求解流程、流程图与源程序，以及手工迭代计算过程等完整章节，便于按步骤理解算法从模型到代码的落地过程。文档还专门介绍了变压器的∏型等值电路、节点电压方程、MATLAB矩阵运算特点，并给出潮流计算流程图、源程序及运行结果，帮助读者直接对照实现。同时，针对传统潮流计算程序依赖大量手动输入、界面不直观的痛点，说明了MATLAB在矩阵运算与可视化分析上的优势。课程设计说明书还包含摘要、关键词、总结与参考文献等模块，结构完整，可兼作报告写作模板。对需要撰写课程设计报告、完成上机调试或准备答辩的读者具有明确参考价值。已有161人学习浏览，说明该资料在同类课程设计中具备一定实用性和关注度。

1. 潮流计算为什么值得自己写一遍：MATLAB 选型与课程设计的真实门槛

很多教科书把潮流计算讲成一组漂亮的非线性方程组，但真正动手做课程设计时才会发现，难的不是牛顿-拉夫逊法本身，而是从一张线路表到可收敛程序的完整链路。这份基于 MATLAB 的电力系统潮流计算课程设计，本质上是一个六节点、七支路的经典算例，要求同时完成手算迭代和程序实现，覆盖了节点导纳矩阵构建、变压器 Π 型等值电路、雅可比矩阵迭代求解、支路功率与平衡节点功率计算等完整环节。传统 C 语言方案在矩阵运算和复数处理上需要写大量底层代码，而 MATLAB 的矩阵天然支持复数运算，调试迭代过程时可以直接观察变量变化，这让它成为绝大多数电气工程专业学生完成潮流计算课程设计的首选环境。适合正在做课程设计需要完整思路参考的本科生，也适合刚接触电力系统分析、想把算法细节落实到代码里的研究生，甚至对牛顿法收敛边界感兴趣的从业者也能从中看到一些工程实现层面的细节。

2. 节点导纳矩阵：从 Π 型等值电路到可执行的 MATLAB 函数

2.1 为什么必须先把变压器折算成 Π 型等值电路

实际电力系统中变压器的变比往往不等于 1，如果直接按原始参数建立节点电压方程，变比会出现在理想变压器的约束方程里，让节点导纳矩阵的对称性被破坏，程序实现时也会多出不少分支判断。课程设计里给出的线路数据中，支路 1-2 的变比为 1.025，支路 4-3 的变比为 1.100，这就是典型的需要折算的场景。

双绕组变压器在不计励磁支路时，可以用阻抗串联理想变压器来等效。设理想变压器变比为 k，变压器阻抗为 Z_T，推导后得到 Π 型等值电路的三个支路导纳为：

y_ij = y_T / k y_i0 = y_T * (1 - k) / k^2 y_j0 = y_T * (k - 1) / k

其中 y_T = 1 / Z_T。注意这三条式子中，k 的位置不能记错，否则程序跑出来的结果会和手工计算结果完全对不上。实际编程时，最稳妥的做法是写一个独立函数处理变压器支路，把折算后的三条支路导纳返回给主程序，再并入整体导纳矩阵。

2.2 用 MATLAB 函数实现 Y 矩阵构建

以这份课程设计的六节点系统为例，线路参数包括电阻 R、电抗 X 和变比 Tap Ratio，基准容量为 100 MVA。构建节点导纳矩阵的常见做法是逐支路遍历，先求串联导纳，再按支路类型决定是否做变压器折算。下面这个函数可以直接用于该算例：

function Y = formYbus(nb, branch) % branch 每行: [from, to, R, X, TapRatio] % TapRatio = 1 表示普通线路；否则按变压器处理 Y = zeros(nb, nb); for k = 1:size(branch, 1) from = branch(k, 1); to = branch(k, 2); R = branch(k, 3); X = branch(k, 4); Tap = branch(k, 5); z = R + 1j * X; % 串联阻抗，标幺值 y = 1 / z; % 串联导纳 if Tap ~= 1.0 % 变压器 Π 型等值电路折算 y_from_to = y / Tap; y_shunt_from = y * (1 - Tap) / Tap^2; y_shunt_to = y * (Tap - 1) / Tap; Y(from, from) = Y(from, from) + y_from_to + y_shunt_from; Y(to, to) = Y(to, to) + y_from_to + y_shunt_to; Y(from, to) = Y(from, to) - y_from_to; Y(to, from) = Y(to, from) - y_from_to; else % 普通线路：两端直接加串联导纳 Y(from, from) = Y(from, from) + y; Y(to, to) = Y(to, to) + y; Y(from, to) = Y(from, to) - y; Y(to, from) = Y(to, from) - y; end end end

关键参数说明：from和to是支路两端节点编号，R、X必须是标幺值，如果原始数据给的是有名值，需要先按基准容量和基准电压归算；Tap只有在变压器支路上才不等于 1，普通线路置 1 即可。自导纳的累加逻辑是核心，对角元 Y(from, from) 要叠加本支路的所有关联导纳，包括串联部分和变压器折算后的对地部分；互导纳则取负值。这段代码可以处理任意节点数的系统，只需要改nb和branch数据矩阵。

2.3 六节点系统的 Y 矩阵数据准备

课程设计给出的线路表需要整理成 MATLAB 可直接读取的矩阵。支路数据共有 7 条，其中有两条带变比的变压器支路，其余为普通线路。我把线路表整理成如下形式：

注意支路 1 的 R 为 0，这意味着该支路是纯电抗支路，导纳是纯虚数，在雅可比矩阵中对应的偏导数项不会出现实部耦合，这会让某些子块的计算简化。支路 6 同样如此，且变比为 1.100，是最容易出错的一条支路，因为它的 k 值参与 Π 型等值电路折算，任何一处符号写反都会导致矩阵不对称。运行formYbus后，可以用full(Y)查看完整矩阵，检查对角元是否为该节点所有关联支路导纳之和，非对角元是否为负的互导纳。如果矩阵不对称，优先检查变压器折算公式里的(1-Tap)/Tap^2和(Tap-1)这两项的符号。

3. 牛顿-拉夫逊法极坐标迭代：雅可比矩阵的构造与收敛控制

3.1 极坐标形式的节点功率方程

节点电压用极坐标表示，即 U_i = U_i∠δ_i，功率方程分为有功和无功两部分。对于 PQ 节点，已知 P_i 和 Q_i，待求量为电压幅值 U_i 和相角 δ_i；对于 PV 节点，已知 P_i 和 U_i，待求量为 δ_i 和注入无功 Q_i；平衡节点电压幅值和相角都给定，用于功率平衡。课程设计选择了极坐标牛顿-拉夫逊法，因为相比直角坐标，极坐标下待求方程数量更少，PV 节点的处理也更直接。

有功和无功的不平衡量表达式为：

ΔP_i = P_ispec - P_i_cal = P_ispec - U_i * Σ(U_j * (G_ij*cosδ_ij + B_ij*sinδ_ij)) ΔQ_i = Q_ispec - Q_i_cal = Q_ispec - U_i * Σ(U_j * (G_ij*sinδ_ij - B_ij*cosδ_ij))

其中 δ_ij = δ_i - δ_j，G_ij 和 B_ij 分别是导纳矩阵元素的实部和虚部。迭代过程中，每次更新 δ 和 U，重新计算所有节点的注入功率，然后求不平衡量，判断是否满足精度要求。

3.2 雅可比矩阵各子块的物理含义

极坐标牛顿法的雅可比矩阵分为四个子块：H（有功对相角偏导）、N（有功对电压幅值偏导）、J（无功对相角偏导）、L（无功对电压幅值偏导）。编程时最容易混淆的是 N 和 L 子块中电压幅值项的处理，因为偏导结果中会出现 U_j 或 U_i，不同写法差一个系数。我采用的约定是：H 和 N 对应 ΔP，J 和 L 对应 ΔQ，其中 N 和 L 子块中的偏导项乘以 U_j，这样修正方程中的变量就是 Δδ 和 ΔU/U 或者 ΔU，具体取决于编程约定。

下面的代码展示了雅可比矩阵的核心计算片段，采用了直接按公式计算每个元素的写法，便于和教材对照：

% Y = G + jB，U 为幅值向量，theta 为相角向量 % 计算 H 子块对角元与非对角元 for i = 1:n for j = 1:n if i == j H(i,i) = -U(i)^2 * B(i,i) - sum_J; % sum_J 为当前节点注入无功功率相关的累加项 else H(i,j) = U(i)*U(j)*(G(i,j)*sin(theta(i)-theta(j)) ... - B(i,j)*cos(theta(i)-theta(j))); end end end

这段代码没有把完整的 sum_J 计算展开，实际编写时需要先在循环外求各节点注入有功和无功，再回代到 H、N、J、L 的每个元素。这种做法代码量大一些，但每一步都能和教材公式对应，排错时方便用 disp 打印中间矩阵。另一种做法是采用稀疏矩阵和向量化计算，适合节点数很多的系统，但课程设计场景下可读性优先。

3.3 迭代主循环与收敛判据

主循环的整体流程是：初始化电压幅值和相角，PQ 节点电压幅值取 1.0，PV 节点取给定值，相角全部取 0；然后循环计算不平衡量 ΔP、ΔQ，检查是否小于收敛精度；不满足则构造雅可比矩阵，求解修正方程得到 Δδ 和 ΔU，更新状态变量，进入下一轮。典型收敛精度取 1e-6，对应课程设计中"迭代到电压变化小于精度"的要求。

求解修正方程时，直接对雅可比矩阵做左除即可，即dx = J \ dpq，MATLAB 会选择合适的稀疏求解器。需要注意，雅可比矩阵的维度不是固定的，它取决于 PQ 节点和 PV 节点的数量。如果系统有 n 个节点，其中有 m 个 PQ 节点，平衡节点 1 个，那么待求的相角数量为 n-1，待求的电压幅值数量为 m，雅可比矩阵维度为 (n-1+m) × (n-1+m)。

4. 六节点系统实战：手工迭代与 MATLAB 程序的结果对照

4.1 系统给定参数与节点分类

课程设计给定的系统中共有 6 个节点，其中节点 5 的注入有功 P5 = 50.16 MW，同时存在负荷数据 55.0+j13.0、50.0+j5.0、30.0+j18.0。按照潮流计算的一般设置，节点 1 作为平衡节点，电压幅值和相角固定；其余节点中，给定有功和无功的为 PQ 节点，给定有功和电压幅值的为 PV 节点。实际编程时需要先确定哪些节点是 PV 节点，因为 PV 节点的无功是待求量，不能作为已知条件。

从题目数据推测，节点 5 带有发电机，可能被设定为 PV 节点，而带负荷的节点为 PQ 节点。这个分类直接决定雅可比矩阵的维度，如果分类错误，修正方程的维度对不上，程序会直接报错或者迭代发散。工程上常用做法是先把所有节点默认为 PQ 节点，再根据给定的数据类型逐个修改节点类型标志位。

4.2 手工迭代至少两次的计算步骤

课程设计要求至少手工迭代 2 次，这不仅是计算过程，更是对算法理解的检验。手工迭代时一般使用简化的系统模型，将变压器折算到 Π 型等值电路后，建立 Y 矩阵，然后按牛顿法的流程逐步计算。

第一次迭代的步骤如下：设所有节点电压初值为 1.0∠0°，代入功率方程，计算各节点注入功率与给定值的偏差 ΔP、ΔQ；然后构造雅可比矩阵，解方程求得 Δδ 和 ΔU；更新电压值后进入第二次迭代。手工计算时通常只保留 3 到 4 位小数，所以结果和程序计算的差异在半次迭代后就可能出现，这是正常现象，程序用的双精度显然更精确。需要注意的是，手工迭代时雅可比矩阵里的元素是随电压变化而变化的，不能复用第一次算出的矩阵，这是最常见的计算错误。

4.3 程序计算结果与手算结果的偏差分析

程序计算时，收敛后的典型输出包括各节点电压幅值、相角、节点注入功率以及各支路功率分布。以这个六节点系统为例，由于存在变压器变比，某些节点的电压幅值会偏离 1.0 较多，相角差则主要集中在线路阻抗较大的支路上。

对比手工迭代与程序结果时，一个实用方法是把手算第 2 次迭代后的电压值打印出来，与程序第 2 次迭代的中间结果做对比。如果偏差在 0.01 以内，说明手算过程和程序逻辑一致；如果偏差很大，优先检查 Y 矩阵中变压器折算部分，再用YY = full(Y)打印矩阵，逐元素核对自导纳和互导纳。此外要注意程序迭代过程中每轮的 ΔP 和 ΔQ 是否单调递减，如果某轮突然增大，通常说明雅可比矩阵构造有问题或者初值选择不当。

4.4 支路功率与网损的计算

收敛后还需要计算各支路功率以及全网功率损耗，这部分是课程设计评分时容易被忽略的要点。支路功率的计算公式为：从节点 i 流向节点 j 的功率 S_ij = U_i * conj(I_ij)，其中 I_ij 是支路电流。对于变压器支路，电流计算要考虑 Π 型等值电路中的对地支路，不能直接用串联导纳乘电压差。

MATLAB 里可以用如下片段计算支路功率：

for k = 1:size(branch,1) from = branch(k,1); to = branch(k,2); Tap = branch(k,5); y = 1/(branch(k,3) + 1j*branch(k,4)); if Tap == 1 Iij = y * (V(from) - V(to)); Iji = y * (V(to) - V(from)); else % 折算后的 Π 型等值电路 y12 = y / Tap; y10 = y*(1-Tap)/Tap^2; y20 = y*(Tap-1)/Tap; Iij = y12*(V(from) - V(to)) + y10*V(from); Iji = y12*(V(to) - V(from)) + y20*V(to); end end

这部分代码中 V 是复数电压向量，可以直接用 U.exp(1jtheta) 构造。支路损耗等于 S_ij + S_ji 的实部，全网网损则是所有支路损耗之和，平衡节点功率可以通过节点注入功率之和来校验，理论上全网注入功率之和应等于总负荷加网损。

5. 收敛失败定位：初值、雅可比矩阵病态与 MATLAB 调试技巧

5.1 初值选择对收敛性的实际影响

牛顿-拉夫逊法的收敛性对初值比较敏感，这个在理论上很明确，但课程设计里通常不会有人告诉你什么是"合理的初值"。实际经验是：对于 110kV 及以上电压等级的系统，电压幅值初值取 1.0，相角取 0，绝大多数情况可以收敛；但如果系统中有变压器变比偏离 1 较多的支路，或者重负荷节点，平直初值可能让迭代发散。一个实用做法是先用高斯-塞德尔法迭代几次得到一个粗略解，再用这个解作为牛顿法的初值，这也是很多教科书隐含推荐的组合。

遇到迭代发散时，不要急着改雅可比矩阵，先检查前两轮迭代的电压修正量方向。如果修正量符号交替变化且幅度增大，通常是初值离解太远；如果修正量单调增大，大概率是雅可比矩阵元素符号错误。另外一个快速排查方法是用 MATLAB 的cond函数检查雅可比矩阵的条件数，条件数超过 1e10 时矩阵接近奇异，即使迭代收敛也会出现振荡。

5.2 常见 MATLAB 实现错误与排查顺序

综合多份课程设计代码，我发现最容易出错的位置集中在三处：变压器 Π 型等值电路的符号处理、雅可比矩阵中 N 和 L 子块的 U_j 因子遗漏、平衡节点在雅可比矩阵中的行与列删除。前两类错误会导致迭代不收敛或收敛到错误解，第三类错误则直接导致矩阵维度不匹配。排查时建议按如下顺序进行：

如果是 PV 节点无功越限，需要在迭代过程中将其转换为 PQ 节点，重新给定无功值，这也是完整的潮流计算程序必须具备的功能，很多课程设计代码没有处理这一点，导致结果虽然收敛但物理上不合理。

5.3 快速验证程序正确性的小技巧

一个实用的验证技巧是：把变压器所有变比都设为 1，此时系统退化为纯线路网络，程序结果应该与直接用普通线路公式计算的结果完全一致。另外，可以把收敛后的电压代入功率方程，计算各节点注入功率，与给定值对比，偏差应在 1e-6 量级。最后再用平衡节点功率校核全网功率平衡，即所有节点的注入功率之和等于零（在标幺值下，总注入等于总负荷加网损），这一步能同时验证支路功率计算和网损计算的正确性。如果网损出现负值，几乎可以断定变压器支路的功率方向计算错误，值得回头检查 Π 型等值电路中电流方向的约定是否一致。

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