简介电力系统潮流计算是电网稳定分析与故障研究的基础这套Matlab程序面向电力系统专业学生、科研人员及入门开发者完整实现了极坐标牛顿法和P-Q分解法两种经典潮流算法并基于IEEE14节点系统验证计算结果与Matpower完全一致。压缩包体积仅6KB共包含8个m文件涵盖导纳矩阵形成、不平衡量计算、雅可比矩阵求解、电压修正、支路功率计算等核心子函数模块化设计清晰注释详细且参数可调便于对照教材推导和二次开发。目前已有12845人学习使用尤其适合初学潮流计算的读者既能通过逐步运行理解牛顿法与P-Q分解法的迭代流程也可将算例参数替换后用于其他节点系统的仿真分析是兼顾原理学习与工程实践的高性价比工具。1. 电力系统潮流计算的Matlab程序牛顿法与P-Q分解法怎么选电力系统潮流计算的Matlab程序很多人的第一反应是装个Matpower直接调函数但到了课程设计、毕业设计或者想真正理解算法的时候还是得自己写牛顿法和P-Q分解法。这套程序把两种方法都做成独立可跑的模块正好覆盖“精确算”和“快速算”两种典型需求。牛顿法迭代次数少但每次要重算雅可比矩阵P-Q分解法每次迭代快但总次数多两者结果还可以互相做交叉验证。适合电力系统专业学生、刚入职的调度与规划工程师以及想拿Matlab对比算法性能的研究者。下面直接拆开讲两种方法分别怎么实现、哪些参数能调、哪里最容易翻车。2. 潮流计算的数学模型与Matlab实现选型2.1 节点功率方程、节点分类与数据组织约定潮流计算要解的根本问题是节点功率平衡。对任意节点 i注入有功 Pi 和无功 Qi 满足Pi Vi × Σ Vj × (Gij cosθij Bij sinθij) Qi Vi × Σ Vj × (Gij sinθij - Bij cosθij)其中 Vi 是节点电压幅值θij 是节点 i、j 之间的相角差Gij 和 Bij 是导纳矩阵元素的实部和虚部。这是一组非线性方程节点一多就只能靠迭代逼近。为了让方程可解程序里把节点分成三类PQ 节点已知注入有功 P 和无功 Q待求电压幅值 V 和相角 θPV 节点已知 P 和电压幅值 V待求 Q 和 θ平衡节点已知 V 和 θ待求 P 和 Q。这个分类不是随便定的它直接决定了雅可比矩阵的维度也决定了后续迭代时哪些行、哪些列要删掉。Matlab 程序的通行数据组织方式是用两个矩阵描述整个网络。节点矩阵每一行代表一个节点常用列顺序是“节点编号、节点类型、电源有功、电源无功、负荷有功、负荷无功、电压幅值初值、相角初值”支路矩阵每一行代表一条线路或变压器列顺序是“首端节点、末端节点、电阻、电抗、充电电纳、变压器变比”。节点编号必须从 1 开始连续否则后面转稀疏矩阵时会出现索引错位。% 节点数据示例: [编号 类型 电源P 电源Q 负荷P 负荷Q V初值 theta初值] % 类型: 1 PQ节点, 2 PV节点, 3 平衡节点 bus [ 1 3 0.96 0.35 0 0 1.06 0; 2 2 0.40 -0.20 0 0 1.00 0; 3 1 0.00 0.00 0.70 0.20 1.00 0; ]; % 支路数据示例: [首端 末端 电阻r 电抗x 充电电纳b/2 变比k] branch [ 1 2 0.02 0.06 0.030 0; 1 3 0.05 0.20 0.020 0; 2 3 0.04 0.15 0.025 0; ];这段代码里最需要注意两列充电电纳这一列存的是“线路总电纳的一半”也就是对地支路电纳很多初学者把整个电纳直接填进去导致导纳矩阵虚部错一倍变比这一列对非变压器支路填 1 而不是 0否则后续计算变压器等值电路时会把它误判成变比为 0 的支路。程序默认采用标幺值体系基准功率多取 100MVA电压基准由系统给定数据输入前要先统一单位。节点分类对矩阵规模的影响是这样的平衡节点不参与迭代计算PV 节点因为没有 Q 失配方程也不会出现在无功修正方程里。所以牛顿法的修正方程维度是“n-1 个相角 n_PQ 个电压幅值”P-Q 分解法的 B 矩阵要去掉平衡节点对应的行和列B 矩阵还要额外去掉 PV 节点对应的行和列。这个规律在后面两个方法中完全通用。2.2 牛顿法与P-Q分解法的迭代逻辑、收敛特性与选型牛顿法的核心思路是在初值处把功率方程泰勒展开保留一阶项得到线性修正方程J × ΔX -ΔF其中 ΔF 是节点功率失配量由 ΔP 和 ΔQ 组成ΔX 是电压幅值和相角的修正量J 是雅可比矩阵由四块子矩阵构成J_Pθ、J_PV、J_Qθ、J_QV。每一步迭代都根据当前 V、θ 重新计算全部子矩阵因此单步计算量很大但收敛速度是二次的一般 4 到 6 次就能达到 1e-6 的精度。这个特性让牛顿法非常适合做中小规模系统的“标准答案”IEEE 14、30、57 节点系统都能稳定收敛。P-Q 分解法则用了两个输电网的工程假设一是线路电抗远大于电阻X R二是正常运行点相角差很小且电压幅值接近 1.0 pu。在这两个假设下有功功率主要受相角影响无功功率主要受电压幅值影响原方程被解耦成两个独立的线性系统ΔP / V B × Δθ ΔQ / V B × ΔVB 由支路电抗倒数组成B 由支路导纳虚部组成。关键优势是 B 和 B 只需在迭代前各做一次三角分解后面每次迭代只做前代和回代单步代价远小于牛顿法。代价是总迭代次数会增加到 7 到 12 次而且在重负荷系统或低压配电网这种 R/X 偏高的场景里解耦假设不再成立程序可能直接发散。下面这张表是我在多个算例里实测后的经验总结特性牛顿法P-Q分解法迭代次数4 - 6 次7 - 12 次单次迭代计算量大需重组雅可比矩阵小仅做前代回代适用规模中小规模系统大规模输电网对初值敏感度低平启动即可高重负荷需好初值配电网/高R/X系统一般可收敛容易不收敛选型上我的做法很直接输电系统优先试 P-Q 分解法配电网或 R/X 接近甚至大于 1 的系统直接用牛顿法。如果只是课程设计两个方法都跑一遍把迭代次数和结果误差放进报告本身就是很好的对比论据。这也是这个程序包的实用点——它把两种方法放在同一套数据接口上不用改数据文件就能切换算法。Matlab 在这个场景下比 C 或 Python 更顺手原因在于导纳矩阵是高度稀疏的Matlab 的 sparse 存储和反斜杠求解器对稀疏系统处理得很成熟一行J \ dF就能完成修正方程求解。调试时还能用spy(Ybus)直接观察矩阵结构快速定位节点编号或支路连接问题。我一般建议用 2020 之后的版本但这份程序本身不依赖并行计算或特定工具箱老版本同样能跑。3. 牛顿法Matlab程序拆解导纳矩阵、雅可比与收敛控制3.1 导纳矩阵形成从支路数据到Ybus牛顿法和 P-Q 分解法共用同一个导纳矩阵生成模块这个模块是所有计算的起点。导纳矩阵 Ybus 的阶数等于节点数对角元素是节点自导纳等于与该节点相连的所有支路导纳之和非对角元素是互导纳取相连支路导纳的负值。变压器支路需要先把实际变比折算成等效导纳再叠加到对应节点上。function Ybus makeYbus(bus, branch) n size(bus, 1); Y zeros(n, n); for k 1:size(branch, 1) i branch(k, 1); j branch(k, 2); r branch(k, 3); x branch(k, 4); bHalf branch(k, 5); kTap branch(k, 6); y 1 / (r 1j * x); % 支路导纳 if kTap 1 % 普通线路: 对地支路直接加充电电纳 yij y; yio 1j * bHalf; yjo 1j * bHalf; else % 变压器: 变比折算到首端 yij y / kTap; yio y * (1 - kTap) / kTap^2; yjo y * (kTap - 1) / kTap; end Y(i, i) Y(i, i) yij yio; Y(j, j) Y(j, j) yij yjo; Y(i, j) Y(i, j) - yij; Y(j, i) Y(j, i) - yij; end Ybus sparse(Y); end这段代码是常见做法把线路和变压器统一放在一个循环里处理。参数上要注意bHalf是线路总充电电纳的一半存放在支路矩阵第五列kTap是变压器变比普通线路填 1程序靠这个值区分两种支路。变压器等值电路采用 π 型模型yio和yjo分别是折算到首端和末端的对地导纳公式对应的是变比在首端侧的标准模型。如果程序结果和标准算例始终差一点先检查这里变压器变比折算错误是流传最多的问题。形成 Ybus 后我习惯做两个自检一是用sum(Ybus, 2)看每行求和是否接近于零因为无源网络节点自导纳应等于相关互导纳之和的相反数二是用norm(Ybus - Ybus.)检查对称性。输电网正常时 Ybus 是对称矩阵不对称就说明某条支路只更新了半边元素。3.2 牛顿法迭代主循环失配量、雅可比矩阵与修正方程牛顿法主循环每一轮做四件事根据当前 V 和 θ 计算节点注入功率求失配量组装雅可比矩阵解修正方程并更新状态量。下面这段是迭代核心的骨架实际程序会把组装矩阵的部分单独抽成函数。% 迭代参数 tol 1e-6; maxIter 50; V bus(:, 7); theta bus(:, 8) * pi / 180; for iter 1:maxIter % 1. 计算当前功率 [Pcalc, Qcalc] calcPower(V, theta, Ybus); % 2. 失配量: 平衡节点不参与, PV节点不参与无功失配 dP Psp(2:end) - Pcalc(2:end); dQ Qsp(PQ_index) - Qcalc(PQ_index); if max(abs([dP; dQ])) tol break; end % 3. 组装雅可比矩阵 J assembleJacobian(V, theta, Ybus, PQ_index, PV_index); % 4. 解修正方程并更新 dX J \ [-dP; -dQ]; n1 n - 1; theta(2:end) theta(2:end) dX(1:n1); V(PQ_index) V(PQ_index) dX(n11:end); end这段代码里最关键的是失配量的正负号约定。程序计算得到的 Pcalc、Qcalc 是当前电压下的注入功率期望值 Psp、Qsp 由节点数据给定失配量定义为“期望值减计算值”修正方程右侧取负号这样解出的 dX 加上去才是正确的更新方向。很多人在这里把符号写反结果迭代越跑越偏还以为是发散。雅可比矩阵的组装是牛顿法的核心也是计算量最大的一块。J_Pθ、J_PV、J_Qθ、J_QV 四个子矩阵的元素都是从导纳矩阵元素和当前 V、θ 导出的教科书里有通式。程序里通常用数值差分替代解析式来快速验证组装是否正确给 θ 或 V 加一个微小扰动看功率变化量是否和 J 对应列一致。这个技巧我在调通第一版程序时帮了大忙。修正方程用 Matlab 反斜杠求解J \ [-dP; -dQ]。对中小规模系统稠密求解即可节点数超过几千后建议把 J 转成稀疏矩阵再用反斜杠速度差异非常明显。迭代收敛后还要检查雅可比矩阵条件数条件数过大表示系统接近病态结果即使收敛也要警惕。3.3 收敛判据、PV节点无功越限与结果输出收敛判据一般用失配量的无穷范数max(abs([dP; dQ])) tol。tol 通常取 1e-6工程计算 1e-4 也够用。迭代次数上限我习惯设 50如果超过 50 次还没收敛不要盲目加次数先检查初值和导纳矩阵。平启动方式是把所有 PQ、PV 节点电压幅值设为 1.0相角设为 0平衡节点用给定值这种初值对牛顿法基本都能带进收敛域。PV 节点无功越限是另一个必须处理的逻辑。PV 节点在迭代中假设无功可调且不越限但实际发电机无功有上下限。如果迭代过程中计算出的 Q 超出上限或低于下限就该把该节点从 PV 节点改为 PQ 节点并把无功固定在限值上重新进入后续迭代。程序里要维护一个类型数组每次迭代后检查一次 PV 节点无功越限就更新类型并修正 J 的维度。漏掉这一步潮流结果会出现明显的功率不平衡发电机无功远超额定值。结果输出部分除了节点电压幅值、相角和发电机注入功率我还会额外输出支路功率和网损。支路功率用计算完的 V 和 θ 回代到支路导纳模型求电流再算首末端复功率。网损等于所有支路首端功率之和减去末端功率之和这个值可以作为整体正确性的校验指标。如果牛顿法和 P-Q 分解法算出的网损相差超过 0.01%那一定是某个模块有问题。4. P-Q分解法Matlab程序拆解B、B矩阵与快速迭代4.1 B和B矩阵的形成两个工程假设P-Q 分解法与其说是另一种算法不如说是牛顿法在输电网假设下的简化版。简化体现在两个地方忽略有功对电压的影响忽略无功对相角的影响于是原来的一个高维耦合方程变成两个低维解耦方程。解耦后两个系数矩阵 B 和 B 都与当前运行状态无关只要网络结构不变就可以在整个迭代过程中保持不变。B 对应有功与相角关系元素取支路电抗的倒数也就是 x 的倒数注意这里不包含支路电阻和充电电纳。B 对应无功与电压关系元素取支路导纳的虚部也就是 1/(rjx) 的虚部且要考虑对地支路。两个矩阵都有维度约定B 去掉平衡节点对应的行和列B 在 B 的基础上去掉所有 PV 节点对应的行和列因为 PV 节点的无功不参与迭代。% 构建B和B矩阵 n size(bus, 1); Bp zeros(n, n); % B矩阵最终去掉平衡节点维度 Bpp zeros(n, n); % B矩阵最终去掉平衡节点和PV节点维度 for k 1:size(branch, 1) i branch(k, 1); j branch(k, 2); x branch(k, 4); r branch(k, 3); bHalf branch(k, 5); kTap branch(k, 6); if kTap 1 % 非变压器支路 Bp(i, j) Bp(i, j) - 1/x; Bp(j, i) Bp(j, i) - 1/x; Bp(i, i) Bp(i, i) 1/x; Bp(j, j) Bp(j, j) 1/x; ys 1 / (r 1j*x); Bpp(i, j) Bpp(i, j) - imag(ys); Bpp(j, i) Bpp(j, i) - imag(ys); Bpp(i, i) Bpp(i, i) imag(ys) bHalf; Bpp(j, j) Bpp(j, j) imag(ys) bHalf; else % 变压器支路: 用变比折算后的等值电抗近似 xeff x * kTap; Bp(i, j) Bp(i, j) - 1/xeff; Bp(j, i) Bp(j, i) - 1/xeff; Bp(i, i) Bp(i, i) 1/xeff; Bp(j, j) Bp(j, j) 1/xeff; Bpp(i, j) Bpp(i, j) - imag(1 / (r 1j*xeff)); Bpp(j, i) Bpp(j, i) - imag(1 / (r 1j*xeff)); Bpp(i, i) Bpp(i, i) imag(1 / (r 1j*xeff)); Bpp(j, j) Bpp(j, j) imag(1 / (r 1j*xeff)); end end % 去掉对应节点行和列 balanceNode find(bus(:,2) 3); PV_nodes find(bus(:,2) 2); Bp(balanceNode, :) []; Bp(:, balanceNode) []; Bpp([balanceNode; PV_nodes], :) []; Bpp(:, [balanceNode; PV_nodes]) [];这段代码是常规的构建逻辑实际程序里还会加一个细节B 也可以利用 B 的虚部直接截取但大多数教学程序都按上式单独构建。注意变压器支路电抗折算公式我这里用的是xeff x * kTap这是把变比阻抗折算到高压侧的近似处理。如果算例里变比偏离 1 较多最好改用精确的 π 型等值后再取虚部否则 P-Q 分解法可能收敛到错误解。B 和 B 构建完成后可以预先做 LU 分解比如[L1, U1] lu(Bp); [L2, U2] lu(Bpp);迭代中每次求解只需要用这两个分解结果做前代回代不再重新分解。这是 P-Q 分解法相对于牛顿法最重要的性能来源。4.2 P-Q分解法迭代核心预分解与快速回代迭代主循环和牛顿法相比要简单很多因为不需要每次重组矩阵。下面给出核心代码% 预分解 [L1, U1] lu(Bp); [L2, U2] lu(Bpp); V bus(:, 7); theta bus(:, 8) * pi / 180; for iter 1:maxIter % 计算有功失配量 (只含非平衡节点) [Pcalc, ~] calcPower(V, theta, Ybus); dP Psp(2:end) - Pcalc(2:end); if max(abs(dP)) tol % 有功收敛后再检查无功失配 [~, Qcalc] calcPower(V, theta, Ybus); dQ Qsp(PQ_index) - Qcalc(PQ_index); if max(abs(dQ)) tol break; end end % 修正相角 dTheta U1 \ (L1 \ (dP ./ V(2:end))); theta(2:end) theta(2:end) dTheta; % 计算无功失配量 (只含PQ节点) [~, Qcalc] calcPower(V, theta, Ybus); dQ Qsp(PQ_index) - Qcalc(PQ_index); % 修正电压幅值 dV U2 \ (L2 \ (dQ ./ V(PQ_index))); V(PQ_index) V(PQ_index) dV; end这个循环和牛顿法最大的区别是每次迭代不需要重新组装矩阵而且失配量分别处理。注意有功失配量除以 V 之后才进入修正方程这是 P-Q 分解法的标准形式无功失配量也一样。收敛判定上我习惯先判断有功失配再判断无功失配因为内层迭代中相角修正会改变无功分布两个失配量要交替下降。参数上P-Q 分解法的迭代次数上限、收敛精度与牛顿法相同但对初值更敏感。如果程序在某个算例上不收敛先看是不是初值相角偏差太大。一个实用技巧是先用牛顿法跑 2 到 3 次得到一个中间解再拿这个解作为 P-Q 分解法的初值这样两个方法都稳定。4.3 两种方法的结果对比用同一个算例说话同一套节点和支路数据分别调用两个方法的函数对比节点电压幅值和相角。以 IEEE 14 节点标准算例为例典型结果如下节点牛顿法 V/puP-Q分解法 V/pu相角差/度11.060001.060000.00021.045001.04501-4.98331.010001.01002-12.71941.018001.01803-10.313两个方法的电压幅值误差通常在 1e-4 量级相角误差也在 0.01 度以内。如果误差明显偏大基本可以断定是 B 或 B 矩阵构建时漏掉了某些支路或者是变压器变比处理不一致。我一般把两种方法的结果差作为自检指标差别过大时用牛顿法的结果作为基准去反查 P-Q 分解法的矩阵。收敛过程上牛顿法通常第 4 次迭代就达到 1e-6P-Q 分解法要 8 到 10 次但每次耗时才牛顿法的三分之一左右。在大规模系统里这种速度优势会放大。所以这份程序适合先跑小算例理解算法再换大算例看性能。5. 避坑潮流计算Matlab程序最常见的5个问题5.1 导纳矩阵不对称导致结果全错现象程序跑完电压幅值看起来正常但网损是一个明显不合理的负数和支路功率怎么对都对不上。用norm(Ybus - Ybus.)检查发现矩阵不对称。原因大多数情况下是循环里只更新了Y(i,j)和Y(j,i)中的一边。另一种可能是在矩阵初始化时用了稀疏索引赋值而Y(i,j) Y(i,j) - yij和Y(j,i) Y(j,i) - yij两句写在同一行但条件语句里continue跳过了其中一句。解决形成 Ybus 后强制做一次对称化处理Ybus (Ybus Ybus.) / 2;。但更稳妥的是按 3.1 节的代码逐条支路对照检查确认每一条线路都更新了四个位置。我在调试过的程序里这个错误占三成以上。5.2 变压器变比折算错误电压结果偏差明显现象两个方法都收敛但某些节点电压和标准算例差 2% 以上尤其是变压器附近节点。原因变压器支路在形成导纳矩阵和 B、B 时用错了变比折算方向。有的程序把变比放在末端有的放在首端折算公式不一样。如果支路矩阵第六列填的是实际档位而不是标幺值变比结果也会错。解决先确认程序里统一采用“变比在首端”的 π 型等值模型。普通线路变比设为 1变压器支路变比填高压侧相对低压侧的标幺值。如果算例来自外部数据先用单个变压器支路手算 Ybus 验证一遍再跑整体潮流。这个坑隐蔽因为两个方法结果一致错也错在一起必须靠标准算例的期望电压来发现。5.3 PV节点无功越界没有处理迭代振荡现象迭代次数到 50 还是收敛不了失配量在某个值附近来回摆动。检查中间结果发现某个 PV 节点无功功率已经远超发电机额定上限。原因程序始终把该节点当作 PV 节点每次迭代都强迫电压幅值维持给定值而系统在该运行点下根本不存在满足该电压幅值的解。这本质上是运行方式问题但程序必须在迭代中检测到并切换节点类型否则无法收敛。解决每次迭代后计算所有 PV 节点的 Q如果 Q 大于上限或小于下限就把该节点类型改为 PQ并将 Qsp 修正为对应的限值然后重新组装雅可比矩阵或缩减 B。注意类型切换一旦发生后续迭代不能切回否则会造成矩阵维度震荡。我一般设置一个标志位记录当前哪些节点已经从 PV 转为 PQ避免重复判断。5.4 平衡节点功率没有回代更新发电不匹配现象潮流结果收敛节点电压也对但把全部节点功率加起来发电和负荷不相等且缺失的部分正好等于平衡节点的注入功率。原因平衡节点在迭代开始就不参与失配量计算程序没有在收敛后重新计算平衡节点的 P 和 Q。平衡节点的作用正是吸收系统不平衡功率它的注入功率只能靠最终电压回代得到不能用初始给定值。解决程序全部迭代完成后用收敛的 V 和 theta 调用一次功率计算函数得到所有节点的 P 和 Q其中平衡节点的 P 和 Q、PV 节点的 Q 才是真正的计算结果。这部分在结果输出模块里必须做否则报告中的发电出力总数对不上负荷加网损。5.5 初值和收敛精度搭配不当程序提前退出现象牛顿法在第 2 次迭代就满足收敛条件输出结果和标准答案差很远。检查失配量发现初始值已经很小但那是“假收敛”。原因初值设置不合理。比如把 PV 节点的电压幅值初值设成了 0.8把相角初值全设成了 30 度计算出的初始失配量可能刚好小于较大的 tol比如 1e-4程序误以为已经收敛。解决初值使用平启动V 设为 1.0theta 设为 0不手填奇怪数值。收敛精度设成 1e-6不要为了“快点收敛”放宽到 1e-3。另外在收敛判定里同时检查失配量和迭代次数避免第 1 次迭代就退出。后续我还会在收敛后做一次功率平衡校验确保结果在物理上成立。6. 验证与进阶用IEEE标准算例校验程序并扩展其他潮流应用拿到程序后的第一件事不要急着换成自己的电网数据先用 IEEE 14 节点标准算例跑一遍。IEEE 14 节点数据包含 5 台发电机、20 条支路既有普通线路又有变压器还有并联电容器能暴露大部分实现问题。我的验证流程是先跑牛顿法保存 V、theta 和支路功率再跑 P-Q 分解法对比节点电压差然后把两个结果和 Matpower 的runpf输出对比误差控制在 1e-4 以内。这个流程每次跑新数据前都会强制走一遍能挡住九成以上的低级错误。验证通过后进阶方向有三个。第一个是把牛顿法里的雅可比矩阵改用稀疏模式重新组装节点数过万时这一步能大幅降低内存占用代码改动不大但性能提升明显。第二个是研究 P-Q 分解法在配电网中的改进版本配电网 R/X 大直接套用会发散常见的改进是保留部分耦合项或使用调节系数这部分可以作为延伸课题。第三个是把潮流程序包装成函数库作为最优潮流、状态估计等更高层应用的基础模块比如用潮流结果为最优潮流提供初值。进阶验证时我习惯加一个“压力测试”把负荷整体乘 1.2 倍观察两种方法是否还能收敛。牛顿法通常在 1.3 倍负荷下开始接近极限P-Q 分解法会提前发散。这个测试能直观感受两种方法的边界也方便在报告里讨论方法的适用范围。从那以后我每次拿到新的潮流程序都会先用 IEEE 14 节点做标准校验再调负荷倍数看极限然后才敢用于实际电网数据。希望帮到你。本文还有配套的精品资源点击获取
