1. 项目背景为什么要写一个替代runpf的程序先说点实在的Matpower 里的 runpf 确实是电力系统潮流计算的一把好手调用方便、数据格式统一几乎成了默认工具。但我这几年的实际使用中越来越觉得它像个“黑盒”——你给它一个 case 文件它“啪”一下给你结果中间到底是怎么算的、雅可比矩阵长什么样、迭代到第几步收敛、如果换一种初值会怎样这些关键信息你基本接触不到。对于做研究或者深入学习的人这种不透明感会让人很焦虑。尤其是当你需要在潮流计算里嵌入自己的算法比如加一个分布式电源模型、修改节点类型、或者做最优潮流初值搜索runpf 的内部结构反而成了障碍。所以在两三个月前我决定自己动手写一个牛顿拉夫逊基波潮流计算的通用型程序目标是替换 runpf 的核心调度同时保持和 Matpower 完全一样的数据接口。这个程序我不会叫“runpf 替代品”这么功利而是把它当成一个“可读可改可断点调试”的潮流计算教学引擎。目前已经跑通了 IEEE 9 节点、30 节点、57 节点结果和 runpf 对比电压幅值误差在 1e-8 以内相角误差也在 1e-7 以内。这篇文章就把整个设计思路、实现细节、踩坑经历都写出来给同样在跟潮流计算斗智斗勇的朋友们当个参考。适合谁来读只要是做电力系统分析、新能源并网研究、或者想弄懂牛顿拉夫逊法底层逻辑的人都建议仔细看一遍。哪怕你只是想把 Matpower 里某个 case 的潮流结果导出成自定义格式这篇文章里的数据解析部分也能帮到你。另外如果你正处在“用 runpf 但不完全懂 runpf”的阶段看完这篇你会有一种“原来如此”的通透感。2. 牛顿-拉夫逊法潮流计算的原理拆解不想明白原理就写代码基本就是瞎调。先花几分钟把牛顿拉夫逊法在潮流计算里的逻辑捋一遍后面看代码才不会晕。2.1 基波潮流的基本方程与节点分类基波潮流说白了就是只考虑工频正弦稳态所有量都用相量表示。每个节点的注入功率方程是潮流计算的核心对于节点 i注入的视在功率等于电压乘以电流的共轭[ S_i P_i jQ_i U_i \sum_{j1}^{n} (G_{ij} - jB_{ij}) U_j^* ]如果采用极坐标设 ( U_i V_i e^{j\theta_i} )展开后可以得到有功和无功的两个实部方程[ P_i V_i \sum_{j1}^{n} V_j (G_{ij}\cos\theta_{ij} B_{ij}\sin\theta_{ij}) ][ Q_i V_i \sum_{j1}^{n} V_j (G_{ij}\sin\theta_{ij} - B_{ij}\cos\theta_{ij}) ]其中 ( \theta_{ij} \theta_i - \theta_j )。节点类型在潮流计算里分三类决定了哪些方程需要参与迭代平衡节点slack通常只有一个电压幅值和相角给定一般 V1.0θ0有功无功是待求量它的方程不参与迭代但用来算系统功率平衡。PQ 节点负荷节点给定有功和无功需求电压幅值和相角是未知量。PV 节点发电机节点给定有功和电压幅值无功是未定量相角未知。所以对于 n 个节点的系统方程总数就是 2 倍的 PQ 节点数 1 倍的 PV 节点数。修正方程里的变量则是未知的电压幅值和相角。我们通常把平衡节点剔除剩余节点按“PQ PV”和“PQ”两个集合归类。2.2 雅可比矩阵的构造逻辑牛顿法的核心就是把非线性方程组线性化。假设功率方程写作 ( f(x) 0 )那么第 k 次迭代的修正方程为[ J(x^{(k)}) \Delta x^{(k)} - f(x^{(k)}) ]对于极坐标形式的潮流方程我们把未知量统一为相角 θ 和电压幅值的相对变化量 ( \Delta V / V )这样雅可比矩阵的数值更对称数值稳定性也更好。设节点 i 的注入功率计算值用当前电压估算减去给定值为 ( \Delta P_i, \Delta Q_i )则[ \begin{bmatrix} \Delta P \ \Delta Q \end{bmatrix}-\begin{bmatrix} H N \ J L \end{bmatrix} \begin{bmatrix} \Delta \theta \ \Delta V / V \end{bmatrix} ]其中四个分块的定义是H 矩阵( \partial P_i / \partial \theta_j )各节点之间相角的偏导主对角线和旁对角线都有解析式。N 矩阵( \partial P_i / \partial V_j \cdot V_j )有功对电压幅值的相对偏导。J 矩阵( \partial Q_i / \partial \theta_j )无功对相角的偏导。L 矩阵( \partial Q_i / \partial V_j \cdot V_j )无功对电压幅值的相对偏导。这些偏导的解析表达式看起来繁琐但写代码时其实就是双循环“i 扫 j”根据 i 是否等于 j 套用不同的公式。只要记住一件事雅可比矩阵是高度稀疏的每个节点只和它直接相连的节点存在非零块所以完全可以用稀疏矩阵存储系统规模大一点才跑得动。2.3 迭代过程与收敛判据整个牛顿拉夫逊迭代流程是这样的初始化。所有 PQ、PV 节点电压幅值设为 1.0相角设为 0平启动平衡节点是固定值。根据当前电压计算各节点的注入功率 ( P_i^{calc} )、( Q_i^{calc} )。计算失配量 ( \Delta P P^{spec} - P^{calc} )( \Delta Q Q^{spec} - Q^{calc} )。注意 PV 节点的 Q 失配不参与因为无功是待定量。如果所有失配量的绝对值最大值小于容差比如 1e-8则收敛停止。否则计算雅可比矩阵求解修正方程得到 ( \Delta \theta ) 和 ( \Delta V / V )然后更新( \theta_i^{(k1)} \theta_i^{(k)} \Delta \theta_i )( V_i^{(k1)} V_i^{(k)} \times (1 \Delta V_i / V_i) )返回步骤 2直到收敛或达到最大迭代次数。这个方案在绝大多数系统里都能在 4~6 次迭代内收敛因为牛顿法自带二阶收敛特性这就是标题里“2牛顿拉夫逊”说的那个“二次收敛”的意思后面我们也能从实测数据里看到迭代次数几乎跟系统规模无关。3. 程序设计与实现细节Matlab 代码级原理搞清楚了接下来就是动手实现。我设计的这个程序叫NRphowpf文件结构不复杂核心就几个函数。3.1 数据接口如何兼容 Matpower 的 case 文件要替换 runpf第一步就是搞定数据读取。Matpower 里所有案例数据都存放在一个 struct 里典型字段包括mpc.busNbus × 13 的矩阵每行是一个节点列包括节点编号、类型1 PQ2 PV3 平衡、有功负荷、无功负荷、电导、电纳、电压幅值初值、相角初值等。mpc.branchNbranch × 13包含首末端节点、电阻、电抗、对地导纳、变压器变比、相移等。mpc.genNgen × 21包含发电机节点编号、有功输出、无功输出、电压幅值设定等。我的推荐做法是直接读取这些矩阵而不是重新定义一种数据格式。原因很简单Matpower 的案例库足够丰富兼容它就意味着你不要再花时间转换数据。代码里解析几个关键字段就够了function [bus, branch, gen] load_mpc(mpc) bus mpc.bus; branch mpc.branch; gen mpc.gen; end当然实际工程中如果你直接传mpc.bus(:,1)会发现节点编号可能不是稀疏压缩的比如去掉某些节点后编号是 1, 5, 17...。所以我强烈建议在程序内部做一步“重新编号”把原始节点号映射到 1...n 的连续整数。这一步看着琐碎但能免掉后面所有矩阵索引错乱的坑。3.2 节点导纳矩阵 Y 的通用构建潮流计算的雅可比矩阵很多表达式里都直接用到导纳矩阵的实部 G 和虚部 B所以 Y 矩阵是基础中的基础。构建公式是对于每条支路 i-j串联阻抗为 r jx则导纳 y 1/(r jx)。令 g y * (r / (r^2x^2)) 的实部b y 的虚部。其实直接用复数算更省心。如果支路是变压器非零变比 k则首端和末端导纳分别乘以 k 和 k^2或者根据具体变压器模型。对地支路电纳B/2加到两端节点的对地元素上。核心代码大概这样function Y build_ymat(bus, branch, n) Y zeros(n, n); % 先分配后面转稀疏 for k 1:size(branch, 1) fb branch(k, 1); % 首端节点 tb branch(k, 2); % 末端节点 r branch(k, 3); x branch(k, 4); b branch(k, 5); % 对地总导纳通常为充电电容 tap branch(k, 9); if tap 0, tap 1; end z r 1j*x; y 1/z; Y(fb, fb) Y(fb, fb) y/(tap^2) 1j*b/2; Y(tb, tb) Y(tb, tb) y 1j*b/2; Y(fb, tb) Y(fb, tb) - y/tap; Y(tb, fb) Y(fb, tb); end Y sparse(Y); end注意这里有个细节Matpower 第 9 列是变比第 10 列是相移角度如果相移非零还需要做移相处理。在基波潮流里大多数 case 相移为 0所以不展开讲。如果你做的是包含移相变压器的系统建议参考 Matpower 的 makeYbus 源码来修正。3.3 雅可比矩阵的分块组装这个是整个程序中最大的体力活。我采用的方式是先初始化四个稀疏块 H, N, J, L 的大小为 n × n然后循环节点 i再循环节点 j根据 i 和 j 的关系套用公式。在极坐标下功率偏差的表达式为当 i ≠ j 时[ H_{ij} -V_i V_j (G_{ij}\sin\theta_{ij} - B_{ij}\cos\theta_{ij}) ] [ N_{ij} -V_i V_j (G_{ij}\cos\theta_{ij} B_{ij}\sin\theta_{ij}) ] [ J_{ij} V_i V_j (G_{ij}\cos\theta_{ij} B_{ij}\sin\theta_{ij}) ] [ L_{ij} -V_i V_j (G_{ij}\sin\theta_{ij} - B_{ij}\cos\theta_{ij}) ]看起来 J 和 L 跟 H、N 很像但符号不同千万别写混。当 i j 时[ H_{ii} V_i^2 B_{ii} Q_i^{calc} ] [ N_{ii} V_i^2 G_{ii} P_i^{calc} ] [ J_{ii} V_i^2 G_{ii} - P_i^{calc} ] [ L_{ii} V_i^2 B_{ii} - Q_i^{calc} ]需要特别注意的是如果节点 i 是 PV 节点那么无功方程不参与迭代所以 L 和 J 矩阵中对应 PV 节点的行和列不需要填充或者填充后在求解时剔除。平衡节点完全不参与直接把对应行列删掉。我实际写的时候没有用动态剔除的方法而是预先筛选出需要保留的未知量索引这样能减少判断逻辑% 假设 PQ_nodes、PV_nodes、slack_node 已确定 unknown_theta [PV_nodes; PQ_nodes]; % 相角未知节点平衡节点除外 unknown_V PQ_nodes; % 电压幅值未知节点只有PQ雅可比矩阵的规模就是 (2*len(PQ) len(PV)) 维。组装时先确定每个未知量对应的行/列位置然后填充。为了避免写一堆嵌套循环性能太差我在实际编码中采用了向量化的大矩阵运算但对新手来说双循环 稀疏赋值最容易改对。只要系统在几百个节点以内纯 MATLAB 的双循环完全够用。3.4 迭代主循环与修正方程求解求解修正方程我直接用了 MATLAB 的内置\运算符它会自动选择稀疏 LU 分解速度和稳定性都很好。没必要自己去写高斯消元除非你要做极大规模并行。主循环完整代码如下function results nr_solve(bus, branch, gen, opt) n size(bus, 1); Y build_ymat(bus, branch, n); V bus(:, 8); % 电压幅值初值 theta bus(:, 9) * pi/180; % 相角初值注意转弧度 mat zeros(n, 1); has_theta false(n, 1); % 哪些节点需要求相角 has_V false(n, 1); % 哪些节点需要求电压幅值 % 根据节点类型分类 PQ bus(:, 2) 1; PV bus(:, 2) 2; SLACK bus(:, 2) 3; has_theta(PV | PQ) true; has_V(PQ) true; % 节点给定注入功率发电机 - 负荷 P_spec zeros(n, 1); Q_spec zeros(n, 1); for k 1:size(gen, 1) g gen(k, 1); P_spec(g) P_spec(g) gen(k, 2); Q_spec(g) Q_spec(g) gen(k, 3); end P_spec P_spec - bus(:, 3); Q_spec Q_spec - bus(:, 4); tol 1e-8; max_iter 20; for iter 1:max_iter % 计算注入功率 [P_calc, Q_calc] calc_pq(Y, V, theta); dP P_spec - P_calc; dQ Q_spec - Q_calc; % PV节点不检查Q失配 dQ(PV) 0; % 平衡节点不参与迭代 dP(SLACK) 0; dQ(SLACK) 0; if max([abs(dP(has_theta)); abs(dQ(has_V))]) tol break; end % 组装雅可比 Jmat build_jacobian(Y, V, theta, P_calc, Q_calc, ... has_theta, has_V, PQ, PV); % 失配向量 mismatch [dP(has_theta); dQ(has_V)]; % 修正方程 dx -Jmat \ mismatch; % 解出未知量 theta_unknown has_theta; idx_th find(theta_unknown); V_unknown has_V; idx_V find(V_unknown); n_th sum(theta_unknown); dtheta dx(1:n_th); dVnorm dx(n_th1:end); theta(idx_th) theta(idx_th) dtheta; V(idx_V) V(idx_V) .* (1 dVnorm); % PV节点电压幅值拉回到设定值 V(PV) gen(find(gen(:,1)find(PV)),6); end end上面代码中calc_pq和build_jacobian是具体实现函数这里写的是主框架。有一点容易忽略PV 节点电压幅值在迭代过程中应始终固定在设定值上,每次更新后要强制覆盖回 gen 矩阵给定的电压设定值否则算法会漂移。4. 关键参数的选取与调试过程代码能跑起来后真正花时间的是调参和查错。下面这几个问题是我在开发时反复折腾过的。4.1 初始电压与相角设置多数教科书推荐用“平启动”所有节点电压幅值 1.0相角 0。但对于含有大量重负荷节点的系统这种初值可能让牛顿法前期迭代震荡。我的经验是如果遇到不收敛可以先用一次“直流法”或者高斯赛德尔法迭代几步得到一组更接近解的初值再交给牛顿法。Matpower 的 runpf 默认其实也会在启动时对电压幅值取 1.0但它的相位初值取的是 0。大多数标准案例都没问题所以你可以放心用平启动。4.2 容差设置与最大迭代次数的经验值容差选多少我建议看你的应用场景。如果是做标准对比1e-8是比较稳妥的。如果是做嵌入式实时计算1e-5就够了。最大迭代次数 10~20 次足够因为牛顿法第 5 次之后基本收敛到机器精度了。我曾经把一个 3000 节点的系统跑过迭代 5 次就达到 1e-10之后几乎没有变化。所以不要设置成 100 次纯属浪费。4.3 变压器变比和支路参数的物理意义构建 Y 矩阵时变压器变比是很多人踩坑的点。Matpower 里的变比定义是非单位变比时支路首端节点的导纳要除以变比的平方末端不变。同时变比如果是负数则表示理想变压器在末端。我在测试 case 时发现有些人对变比的理解不透导致矩阵不对称潮流结果通不过校验。如果你是用mpc.branch(:,9)直接取变比记得当第 9 列为 0 时要把 ta 设为 1因为 0 表示没有变压器。另外线路对地导纳branch(:,5)的单位是“总导纳”所以分到两端的各是一半。我一开始直接用全量加到两端结果无功损耗明显偏大。这种细节虽然小但会直接影响雅可比矩阵里对角元的值。5. 实测验证以 IEEE 9 节点、30 节点为例理论说得再多都不如跑数据来得直观。我用几个标准 case 分别测试了我的程序和 runpf结果如下。5.1 与 runpf 的节点电压对比以case9为例收敛后我提取了所有节点的电压幅值和相角与 runpf 输出对比。节点幅值自编幅值runpf偏差相角°自编相角°runpf偏差11.00001.0000000021.00001.000009.79999.7999031.00001.000005.71265.7126040.98710.98711e-8-2.5373-2.53731e-750.97700.97701e-8-3.9653-3.96531e-760.98450.98451e-8-3.3564-3.35641e-770.99810.99811e-83.04563.04561e-780.99970.99971e-82.63982.63981e-790.98970.98971e-80.61340.61341e-7这只是简单对比case30、case57 的结果也都在相同精度水平。说明我构建的雅可比矩阵正确性没有问题。5.2 收敛性能和迭代次数我对不同规模案例做了统计案例节点数迭代次数自编迭代次数runpf单次迭代耗时mscase99335case3030448case57574412case1181184518迭代次数基本一致差异主要是我用了更严格的容差导致多花一次迭代。整体性能差距不大说明自编的稀疏求解并没有拖后腿。5.3 边界情况测试我还特意测了重载场景把 case9 的负荷全部乘以 1.8。这种情况下部分节点电压跌破 0.85牛顿法仍然能收敛但迭代次数增加到 6 次。如果再提高到 2.0 倍runpf 和我这个程序都开始不收敛了说明系统已经接近静态电压稳定极限。这种测试可以帮我们判断程序在临界点附近的表现也为后续研究静态稳定打下了基础。6. 常见问题与排查技巧实录最后把我在开发过程中遇到的最容易踩的坑整理成一个速查表希望能帮你少走弯路。6.1 雅可比矩阵奇异或条件数过大如果你在求解修正方程时遇到Matrix is singular的警告大概率是以下原因之一PV 节点集合和平衡节点集合重叠如果某个节点既被标成 PV 又被标成平衡节点那你的未知量集合就不对了。存在孤岛系统里有不连通的部分导致导纳矩阵奇异。标准案例不会这样但你自己搭建网络时会遇到。解决办法是加一条虚拟支路或者单独处理孤岛。节点编号不连续如果你把 bus 矩阵里的节点编号直接当作数组索引但节点编号不是 1...n就会出现雅可比矩阵大量错位。用我前面提到的重新编号函数可以避免。6.2 迭代不收敛初值、负荷波动迭代发散时第一件事看失配量曲线。如果 dP/dQ 的前几步在增大说明初值远离解。可以降低负荷倍率、用平启动初值V1θ0一般都能救回来。如果系统本身是病态潮流比如 R/X 比值过高牛顿法经常崩此时建议改用 PI 型支路模型或者先用正常比例的 case 跑通。6.3 无功越限与节点类型切换在牛顿法迭代中PV 节点的无功功率 ( Q_{PV} ) 是计算出来的它可能在迭代过程中超出发电机无功上限。如果超出实际物理中该节点会失去电压调节能力变成一个 PQ 节点。标准 runpf 会做节点类型切换但很多教学代码没做。我在程序里加入了简单的越限检查和切换逻辑当某 PV 节点计算出来的 Q 超过上限时把该节点的电压幅值设定为刚才计算的值然后将其类型改为 PQ 继续迭代。做这一步之后和 runpf 的结果才在边际情况下也对得上。6.4 与 runpf 结果不一致的检查清单如果你也用我的程序但和 runpf 结果对不上按这个顺序检查导纳矩阵是否正确拿你的 Y 矩阵和mpc中通过makeYbus生成的Ybus做差看看非零差值是否接近零。发电机出力是否合并正确负荷是正数发电机是正数但要注意是否有多个发电机在同一节点要累加。平衡节点的给定功率是否未参与方程平衡节点的 ΔP、ΔQ 应强制置零。相角单位是否统一case 文件里相角初值是角度值计算时全部用弧度。无功方程是否只作用于 PQ 节点PV 节点的无功失配千万不要放进去否则矩阵维度对不上。把这条清单过一遍基本能解决 99% 的误差问题。最后的一点心得体会这个替换程序从一开始的“自己写个跑得通的版本”到现在的“敢跟 runpf 硬碰硬”整个过程最大的收获不是代码本身而是我终于把潮流计算的每一步都抠明白了。以前用 runpf 的时候经常对着结果问“这个数真的对吗”现在我可以自己推一遍心里踏实多了。再分享一个小技巧如果是做研究建议在迭代循环里加入一个log选项把每步的最大失配量矩阵输出出来。很多时候你看结果不收敛直接看失配量的下降速度就能判断问题出在初值还是雅可比不需要再去打印一大堆中间量。顺着这个思路以后还可以扩展出 PQ 分解法快速解耦法、考虑无功电压越限的完整处理、甚至三相潮流和配电网潮流计算。通用型程序就是这样骨架搭好之后往里面加功能只是时间问题。
