简介这份 Word 文档是《基于Matlab的IEEE14节点潮流计算程序》的配套说明书面向电力系统稳态分析课程的学习者、课程设计及大作业撰写者帮助其理解并复现 IEEE14 节点潮流计算。内容按论文格式组织含目录、摘要、正文与参考文献系统梳理牛顿—拉夫逊法与 P-Q 分解法的原理、实现步骤和流程并给出节点数据与网络拓扑说明程序部分还介绍了将普通支路等效为变比为一的变压器支路、自动判断节点类型等简化思路便于读者对照代码理解导纳矩阵形成与迭代求解过程。资源包共 1 个 docx 文件约 1.04MB结构完整、便于直接查阅和改写。已有 3442 人学习适合作为课程大作业模板、算法对比参考及潮流计算入门材料。1. 潮流计算14 个节点的电压为什么不能一次解出只能一步步逼近打开一份 IEEE14 节点数据很多人的第一反应是列 14 个方程解 14 个电压一次就能解出来。可真跑起来就会发现这组方程是电压的二次型还带三角函数属于非线性代数方程组根本没有求根公式。Matlab 能做的只有逼近给一组初值算出偏差修正一次再算偏差。牛顿-拉夫逊法和 P-Q 分解法的差别本质上就是两件事怎么修正以及修正时愿意丢掉多少物理上不重要的东西。想明白这一点再看那套潮流计算程序的每一行才知道它为什么那样写。还有一件事更容易被忽略算完电压只是中间结果真正要看的是支路功率、网损和平衡节点的出力而这些量全是从同一组电压反推出来的。2. 牛顿-拉夫逊法雅可比矩阵里那四块偏导到底在修正什么2.1 未知量计数为什么方程数一定等于未知量数潮流计算里最容易被跳过的一步是把未知量和方程数数清楚。一个 n 节点的系统按约束条件分三类平衡节点slack给定电压幅值和相角功率完全由计算结果决定PV 节点给定有功和电压幅值无功待求PQ 节点给定有功和无功电压幅值和相角都待求。设 PV 节点 m 个、平衡节点 1 个那么 PQ 节点就是 n-1-m 个。未知量分两组相角 θ 除了平衡节点外全部待求共 n-1 个电压幅值 V 只有 PQ 节点待求共 n-1-m 个。合计 2n-2-m 个。方程也分两组有功不平衡方程对除平衡节点外的所有节点列写共 n-1 个无功不平衡方程只对 PQ 节点列写共 n-1-m 个。合计也是 2n-2-m 个。这个数字对上不是巧合而是节点分类的直接后果。PV 节点为什么不列无功方程因为它的无功是自由量列出来反而多一个未知数。平衡节点为什么两组方程都不列因为它的 V 和 θ 都已知注入功率是被动算出来的。我见过不少初学时的错误都是在这里搞混的要么把 PV 节点的无功方程也塞进雅可比矩阵矩阵维度对不上Matlab 抛出Matrix dimensions must agree要么忘了把平衡节点的行列剔除解出来相角全偏。提示写程序时不要硬编码节点编号和类型用索引向量idx_p参与有功方程的节点和idx_q参与无功方程的节点来切片后面切换算例时不用改代码。2.2 极坐标功率方程电压的二次型从哪来节点注入功率用极坐标写出来是两组求和式。记 θ_ij θ_i - θ_jY_ij G_ij jB_ij则P_i V_i · Σ_j V_j (G_ij·cosθ_ij B_ij·sinθ_ij) Q_i V_i · Σ_j V_j (G_ij·sinθ_ij - B_ij·cosθ_ij)式子里 V_i 乘在外面里面的求和又带一个 V_j所以对电压而言是二次的θ 出现在 cos 和 sin 里所以对相角是非线性的。这两条性质决定了后面所有的事二次型让方程可以求偏导得到解析的雅可比矩阵三角函数的周期性则意味着初值选得离谱时可能收敛到别的解分支上去。求和范围 j 包含 i 本身把 j i 的项单独拎出来cosθ_ii 1、sinθ_ii 0就得到 P_i V_i²G_ii Σ_{j≠i}(…)这个拆分在后面推对角元素时要用。功率不平衡量就是给定值减去计算值% Sbus 为节点给定注入复功率V/th 为当前迭代的电压 S V .* conj(Ybus * V); % 由电压反算各节点注入功率 dP real(Sbus) - real(S); % 有功不平衡量 dQ imag(Sbus) - imag(S); % 无功不平衡量这里有个细节值得注意Ybus * V得到的是节点注入电流再乘conj(V)才是复功率。写成V .* (Ybus*V)少了共轭算出来的无功符号会反电压越迭代越离谱属于典型的翻车现场。2.3 四块矩阵 H、N、M、L偏导的物理含义把不平衡量按泰勒级数展开、只保留一次项得到修正方程[ΔP] [H N] [Δθ ] [ΔQ] [M L] [ΔV/V]雅可比矩阵的四块分别是 ∂P/∂θ、V·∂P/∂V、∂Q/∂θ、V·∂Q/∂V。非对角元素i ≠ j由微分直接得到分块表达式物理含义H_ijV_iV_j(G_ij·sinθ_ij − B_ij·cosθ_ij)节点 j 相角变化对节点 i 有功的影响N_ijV_iV_j(G_ij·cosθ_ij B_ij·sinθ_ij)节点 j 电压变化对节点 i 有功的影响M_ij−V_iV_j(G_ij·cosθ_ij B_ij·sinθ_ij)节点 j 相角变化对节点 i 无功的影响L_ijV_iV_j(G_ij·sinθ_ij − B_ij·cosθ_ij)节点 j 电压变化对节点 i 无功的影响对照一下就能看出 H_ij L_ij、N_ij −M_ij。这不是凑出来的巧合而是因为 P 和 Q 的表达式只差一个 90° 的相移。这个对称性在后面 P-Q 分解法里会被用到。对角元素i j如果老老实实对求和式逐项求导计算量会翻好几倍。实际程序里都用另一套等价写法把求和项替换成 P_i 和 Q_i 本身H_ii ∂P_i/∂θ_i −Q_i − B_ii·V_i² N_ii V_i·∂P_i/∂V_i P_i G_ii·V_i² M_ii ∂Q_i/∂θ_i P_i − G_ii·V_i² L_ii V_i·∂Q_i/∂V_i Q_i − B_ii·V_i²推导过程不复杂关键一步是把 Σ_{j≠i}V_iV_j(G_ij·sinθ_ij − B_ij·cosθ_ij) 认出来等于 Q_i B_ii·V_i²。这么一换每个对角元素从一次 O(n) 求和降成一次乘加。别小看这一处n 上千的时候雅可比矩阵的组装时间能差出好几倍这是实打实的优化不是写法偏好。2.4 为什么修正量写成 ΔV/V 而不是 ΔV修正方程里电压用的是 ΔV/V也就是相对修正量。这个写法有两种解释一是让雅可比矩阵各元素的量纲和数量级更接近条件数更好高斯消元时不容易被舍入误差放大二是让 N、L 两块的结构和 H、M 对齐代码里写起来一致。代价是更新电压时得多乘一次% dx 的前 np 项是 Δθ后 nq 项是 ΔV/V th(idx_p) th(idx_p) dx(1:np); V(idx_q) V(idx_q) .* (1 dx(np1:end)); % 相对修正乘回去如果这里写成V(idx_q) V(idx_q) dx(...)量纲就错了第一次迭代可能还看不出来迭代两三次之后电压直接冲到 0 或者 5 以上然后振荡发散。这类错误在调试时特别难定位因为矩阵求解本身没有任何报错。收敛判据一般取 max(|ΔP|, |ΔQ|) 1e-5 到 1e-6标幺值判据放在算完不平衡量、还没组装雅可比矩阵之前能省掉最后一次没必要的矩阵分解。牛顿法的收敛是平方级的误差大致按 e → e² 缩小所以从 1e-2 到 1e-6 通常只要三四次迭代。这也是它在教科书里被当成标杆的原因——单次迭代贵但迭代次数少得离谱。3. P-Q 分解法把 2n 阶雅可比拆成两个常数矩阵省下的到底是什么3.1 物理前提高压网里有功跟相角走无功跟幅值走牛顿法每次迭代都要重建并分解一个 2n-2-m 阶的雅可比矩阵复杂度是 O(n³)。电网规模一大这一步就成了瓶颈。P-Q 分解法的出发点不是数学上的化简技巧而是一条电网的物理事实在高压输电网里线路电抗远大于电阻节点电压相角差很小有功功率的流动主要取决于电压相角差无功功率的分布主要取决于电压幅值差。把这条事实翻译成雅可比矩阵的语言就是 N 块和 M 块的元素数值远小于 H 块和 L 块。既然小那就直接扔掉认为 N M 0。原来的 2n-2-m 阶方程立刻分裂成一个 n-1 阶的有功方程和一个 n-1-m 阶的无功方程两组轮流迭代。这一步叫有功-无功分解是第一次近似。扔掉 N、M 之后H 和 L 的元素仍然是电压幅值和相角的函数每步还得重算。所以还有第二次近似这才是省时间的关键。3.2 第二次近似H 和 L 怎么变成常数矩阵利用三个条件线路两端相角差很小cosθ_ij ≈ 1、sinθ_ij ≈ 0节点电压幅值接近 1标幺值下与节点无功相对应的并联导纳远小于该节点自导纳的虚部即 Q_i ≪ B_ii·V_i²。代进 H_ij 的表达式V_iV_j(G_ij·sinθ_ij − B_ij·cosθ_ij) 在 sinθ_ij → 0、cosθ_ij → 1、V → 1 之后退化成 −B_ij。于是 H 矩阵就变成了 −B一个只由节点导纳矩阵虚部构成的常数矩阵。L 块同理。这一步做完系数矩阵在迭代过程中彻底不变只需要在最开始做一次三角分解LU之后每次迭代只做前代和回代复杂度从 O(n³) 降到 O(n²)。这是 P-Q 分解法真正的收益所在不是迭代次数少而是每次迭代便宜得多。3.3 B 和 B两个矩阵的构造规则为什么不一样实际程序里有两个系数矩阵习惯叫 B 和 B阶数不同构造来源也不同ΔP_i / V_i Σ_j B_ij · Δθ_j ΔQ_i / V_i Σ_j B_ij · ΔV_j B_ij −1/x_ij (i ≠ j)B_ii Σ_{j≠i} 1/x_ij 只取串联电抗忽略并联支路和变比 B_ij −B_ij (i ≠ j)B_ii −B_ii 取自 Ybus 虚部含并联补偿B 只保留串联电抗的倒数把并联支路、线路电阻、变压器变比全部忽略。原因是有功潮流主要由串联电抗决定并联支路消耗的有功极少留着反而破坏常数矩阵的性质。B 不一样它必须保留并联补偿因为并联电容器和电抗器直接注入无功去掉之后无功平衡就错了迭代会明显变慢甚至不收敛。B 的阶数是 n-1不含平衡节点的行列B 的阶数是 n-1-m平衡节点和 PV 节点的行列都要剔除。两个矩阵的阶数不同这也是程序中容易出错的地方切片用的索引向量不是同一个。注意不同教材对 B、B 的符号约定不完全一致有的把负号吸进方程右边有的留在矩阵里。混用两套约定表现是第一次迭代电压就跳到 2.0 以上。判断方法很简单看 B 的对角元是不是正的在 −1/x 约定下非对角为负、对角的和对角线为正。3.4 迭代顺序为什么先解相角再解幅值P-Q 分解法每轮迭代做两件事先用当前电压算有功不平衡量解出相角修正量并更新 θ再用更新后的相角算无功不平衡量解出电压修正量并更新 V。顺序反过来也能收敛但慢一些因为无功不平衡量对相角更敏感先修相角能让下一步的无功不平衡量更接近真实值。[Lp, Up] lu(Bp); % 常数矩阵只分解一次 [Lq, Uq] lu(Bq); while iter 50 S V .* conj(Ybus*V); dP (real(Sbus) - real(S)) ./ V; % 注意先除以 V dth Up \ (Lp \ dP(idx_p)); % 前代回代O(n²) th(idx_p) th(idx_p) dth; S V .* conj(Ybus*V); % 相角变了功率要重算 dQ (imag(Sbus) - imag(S)) ./ V; dV Uq \ (Lq \ dQ(idx_q)); V(idx_q) V(idx_q) dV; if max(abs([dP(idx_p); dQ(idx_q)])) tol, break; end iter iter 1; end代码里有个容易漏的点更新完相角之后必须重新算一次 S不能沿用上一次的。因为 ΔP 和 ΔQ 都依赖同一组 (V, θ)相角一改无功不平衡量就变了直接复用会让收敛路径变得很奇怪迭代次数平白多出好几次。收敛判据仍然用功率不平衡量不是用 Δθ 和 ΔV 的模。这一点很重要——P-Q 分解法有两次近似修正量本身不再精确对应真实的牛顿步长拿修正量当判据会提前停机结果精度不够。4. MATLAB 实现节点导纳矩阵怎么组、节点类型怎么自动判4.1 数据表设计节点表、支路表、发电机表三张表潮流程序的数据组织基本是固定的三张表节点表存编号、类型、负荷 P/Q、给定电压幅值支路表存首末节点、电阻、电抗、充电电纳、变比发电机表存接入节点、有功出力、电压设定值。三张表的行顺序不要求一致靠节点编号关联。节点类型自动判断的核心思路是从设备表反推约束挂了发电机的节点至少是 PV因为电压可控平衡节点单独指定其余默认是 PQ。写成代码就是几行nb size(bus, 1); gen_bus unique(gen(:, 1)); % 有发电机接入的节点集合 type 3 * ones(nb, 1); % 先全部按 PQ 处理3 PQ type(gen_bus) 2; % 挂机组 - PV2 PV type(1) 1; % 平衡节点1 slack一般取 1 号节点 idx_p find(type ~ 1); % 参与有功方程的节点 idx_q find(type 3); % 参与无功方程的节点 np numel(idx_p); nq numel(idx_q);这么写的好处是换算例时只改数据文件不动代码。前提是数据表里的发电机接入节点必须准确漏填一台机组的接入节点那个节点会被误判成 PQ程序照样能跑但结果里那个节点的电压幅值会偏离给定值——这种错误不报异常只能靠对比参考结果发现。提示PV 节点在迭代过程中如果无功越限严格来说要转成 PQ 节点再继续迭代。教学用的算例一般不做这层检查但如果拿这份程序去算重载工况不做限值处理会得到一个物理上不可实现的解。4.2 导纳矩阵组装把普通支路当作变比 1 的变压器项目里有一句话很关键为了简化导纳矩阵的程序把普通支路当作变比为一的变压器支路计算。这不是偷懒而是把两种元件合并到同一段代码里支路表多一列变比普通线路填 1代码只写一遍。变压器支路的等效是变比 k 折算在 i 侧串联导纳 y 1/(r jx)则Y_ii y/k² j·b/2 Y_jj y j·b/2 Y_ij - y/k Y_ji - y/k当 k 1 时y/k² y、y/k y整套式子退化成普通线路的 Π 型等效。写进 Matlab 就是一个循环nb size(bus, 1); Ybus zeros(nb, nb); for k 1:size(branch, 1) i branch(k, 1); j branch(k, 2); z branch(k, 3) 1i * branch(k, 4); % 串联阻抗 r jx y 1 / z; % 串联导纳 b branch(k, 5); % 线路总充电电纳 kt branch(k, 6); % 变比普通线路填 1 if kt 0, kt 1; end % 容错漏填变比按 1 处理 Ybus(i, i) Ybus(i, i) y / kt^2 1i * b / 2; Ybus(j, j) Ybus(j, j) y 1i * b / 2; Ybus(i, j) Ybus(i, j) - y / kt; Ybus(j, i) Ybus(j, i) - y / kt; end Ybus Ybus diag(shunt); % 叠加节点并联补偿电容器为正电抗器为负几个参数的含义要说清楚。kt是非标准变比标幺值下通常写成 1 ± 0.0x 的形式它折算在哪一侧决定了y/kt²落在 i 还是 j。如果搞反了导纳矩阵不再对称Ybus(i,j) ≠ Ybus(j,i)迭代照样能收敛但网损和电压分布会偏属于最难查的一类错误。b/2是把线路总充电电纳对半分到两端这是 Π 型等效的标准做法。如果数据表里给的是单侧电纳这里就不能再除以 2。shunt是节点级的并联补偿跟支路无关所以放在循环外面单独叠加。这个顺序不能反因为支路的Y_ii是累加节点补偿也是累加两者互不影响但分开写更清楚。4.3 牛顿法主循环矩阵组装和求解的分界主循环的结构是固定的算不平衡量 → 判收敛 → 组装雅可比 → 解修正方程 → 更新状态量。雅可比矩阵的组装建议单独写成一个函数输入是当前的 V、θ 和 Ybus输出是 (2n-2-m) 阶方阵。这样主循环只剩十几行逻辑一目了然。V ones(nb, 1); th zeros(nb, 1); % 平启动 V(type 2) bus(type 2, 5); % PV 节点电压幅值固定 for iter 1:30 S V .* conj(Ybus * V); dP real(Sbus) - real(S); dQ imag(Sbus) - imag(S); mis [dP(idx_p); dQ(idx_q)]; % 只取待解方程 if max(abs(mis)) 1e-6, break; end % 先判收敛再组装 J form_jacobian(V, th, Ybus, idx_p, idx_q); dx J \ mis; % 解修正方程 th(idx_p) th(idx_p) dx(1:np); V(idx_q) V(idx_q) .* (1 dx(np1:end)); endV(type 2) bus(type 2, 5)这一行是 PV 节点的电压约束必须在迭代开始前设好而且迭代过程中不能再被修正。如果忘了这一行PV 节点的电压会从 1.0 慢慢漂到给定值附近——看起来也能收敛但收敛到的不是给定电压的解结果对不上参考数据。J \ mis用的是 Matlab 的左除内部走 LU 分解。如果矩阵接近奇异Matlab 会给警告但不报错解出来的修正量可能大得离谱。遇到这种情况先检查是不是有孤立节点导纳矩阵某一行全零或者两个节点之间的支路参数填成了零阻抗。4.4 结果后处理平衡节点功率和网损从哪来迭代收敛后拿到的是全部节点的 V 和 θ剩下的事都是反推S V .* conj(Ybus * V); % 各节点注入复功率 loss sum(S); % 全网损耗 所有节点注入之和 Slack S(1); % 平衡节点出力sum(S)等于全网损耗这个等式成立是因为所有节点的注入功率加总内部支路上的流动两两抵消只剩下电阻上的有功损耗和电抗上的无功消耗。这个关系可以拿来验算如果算出来的网损是负的说明导纳矩阵或者电压相位有问题基本不可能是数据本身如此。5. 排查潮流不收敛、结果对不上时先看哪几处5.1 不收敛残差在 1e-2 附近反复横跳现象迭代次数打满上限不平衡量在同一个量级上来回振荡既不下降也不发散。原因通常是初值或者约束条件不自洽。最常见的是平衡节点选错——如果指定的平衡节点其实是个负荷节点它给定的注入功率和电压固定这个约束直接冲突方程无解。第二常见的是 PV 节点给定的有功超过了网络传输能力物理上就不存在可行解。解决先检查平衡节点的类型标记和注入功率是否同时被使用平衡节点的 Sbus 应该被忽略再把 PV 节点的有功出力调回合理范围试一次。如果还是不收敛把迭代上限调大、打印每次迭代的最大不平衡量看它是单调下降还是震荡——单调下降只是慢震荡才是真的有问题。5.2 收敛了但电压偏低PV 节点被误判成 PQ现象程序正常收敛但某些节点的电压幅值明显低于数据表里给的给定值偏差在 0.05 到 0.1 之间。原因发电机表里漏了某个节点的接入记录或者节点编号写成了字符串导致unique匹配失败那个节点被当成 PQ 处理电压自然不受约束。解决在组装索引之前打印一次type向量和数据表逐个核对。更稳妥的做法是加一句断言PV 节点数量应该等于发电机接入节点数减一去掉平衡节点数量对不上直接报错退出别让它悄悄算下去。5.3 结果和参考值差一点点变压器变比折算侧反了现象所有节点电压都能对上三位有效数字但网损偏大或偏小几个百分点个别支路的功率方向甚至反了。原因变比 k 折算在 i 侧还是 j 侧搞混了。导纳矩阵不再对称但依然是个合法的复数矩阵迭代照样收敛只是收敛到了另一组解。解决加一条自检——不含变压器的算例Ybus必须严格对称max(max(abs(Ybus - Ybus.)))应该是 0 或者机器精度量级。如果明显不为零说明变比处理有问题。5.4 P-Q 法迭代次数比牛顿法还多B 用错了矩阵现象P-Q 分解法跑了二十几次才收敛比牛顿法的四五次多出一个数量级。原因B 应该只由串联电抗构成如果图省事直接拿-imag(Ybus)当 B并联支路和电阻都被带进来常数矩阵的近似前提被破坏收敛速度会退化成接近线性甚至更差。解决单独写一个构造 B 的循环只读支路表的电抗列忽略充电电纳和变比。B 反过来必须用完整的Ybus虚部两者不能共用同一段代码。5.5 相角单位混用sin 吃进去的是度现象第一次迭代修正量巨大电压直接跳到 3 以上或者 Matlab 报出NaN。原因数据表里相角按度存储或者某处用了sind/cosd另一处又用了sin/cos。Matlab 的三角函数只认弧度混用不会有任何警告。解决整个程序内部统一用弧度只在输入输出的边界做一次转换。最省事的做法是在数据读入之后就乘pi/180之后再也不碰度。6. 算法取舍迭代次数少不等于算得快收敛快不等于算得准两种算法在 IEEE14 这个规模上跑出来的差别其实不大真正拉开差距是在几百个节点以上。把关键维度摆在一起看维度牛顿-拉夫逊法P-Q 分解法系数矩阵每步重建的雅可比矩阵2n-2-m 阶两个常数矩阵 B、Bn-1 阶与 n-1-m 阶单次迭代开销组装 LU 分解O(n³)只有前代回代O(n²)迭代次数少平方收敛IEEE14 常见 4~5 次多线性收敛通常是牛顿法的 1.5~2 倍收敛可靠性初值差时仍可能收敛鲁棒性好依赖 R ≪ X重载或电缆占比高时变差内存占用大小这里有个容易被误读的地方P-Q 分解法的迭代次数更多但总时间往往更短。原因是单次迭代的成本差了一个数量级多跑几次完全划得来。反过来如果只盯着迭代次数这一列就会得出牛顿法更快的错误结论——这是我在看别人写的对比报告时见过最多的误判。第二个容易误读的是精度。两种算法收敛到的是同一组解因为收敛判据用的都是功率不平衡量只要容差设成一样最终电压的差别在 1e-6 量级可以认为完全一致。所谓P-Q 法精度略低实际指的是它在重载工况下可能收敛得慢、甚至不收敛而不是收敛之后解不准。第三个取舍是通用性。牛顿法的雅可比矩阵是严格推导出来的不依赖电网的特殊性质配电网络、含大量电缆的系统、R/X 比接近 1 的场合都能用。P-Q 分解法的两次近似建立在高压输电网的物理特征上一旦这些前提不成立收敛性会明显退化。所以工业级的潮流程序里通常两种都实现让用户按场景选或者干脆做成自动切换。这份说明书把公式推导、节点数据拓扑、两套源码和结果对照放在一起价值不在于给出一个能跑的脚本而在于让每一条近似都能在结果里找到对应的痕迹——改了哪个参数、收敛路径怎么变、网损差在哪里都能顺着章节倒推回去。配套的程序和这份文档一起看比自己从零推一遍雅可比矩阵要省很多时间。从那以后我每次拿到一份 IEEE14 节点数据都会先把导纳矩阵对着对称性和参考值核一遍再开始迭代——这一步花不了两分钟但能省掉后面几个小时的排查。希望帮到你。本文还有配套的精品资源点击获取