前一阵帮人调试一套环形配电网的潮流程序对方从网上找的是前推回代法的代码在辐射网上跑得挺漂亮一换成环形网络就直接原地爆炸、迭代发散、报索引越界各种问题轮着来。这个事其实很典型很多初学者没太想清楚潮流算法的收敛性和网络拓扑有很强的关联辐射网用前推回代当然舒服但城市配网、输电网大量采用环形接线这时候就得老老实实上牛拉法也就是牛顿-拉夫逊法Newton-Raphson。这篇文章分享我用Matlab写的一个通用牛拉法潮流计算程序的核心思路输入只有两张表——节点数据表和支路数据表程序自动形成节点导纳矩阵不管辐射、单环、多环还是随便怎么连都能直接计算通用性不绑定具体算例。电力系统专业的同学、刚接触配网或输电网仿真的工程师以及想搞懂牛拉法代码细节的人都可以参考。1. 为什么环形网络绕不开牛拉法先看拓扑再选算法1.1 辐射网算法与环网的根本冲突很多教材在讲潮流时都会先提前推回代法因为它简单、占内存少、对配电网这种典型辐射状网络收敛速度飞快。前推回代的核心思想是顺着树走从根节点出发功率一层一层往下推再从末端节点往回代电压本质上是利用了树状网络每一个节点只有唯一供电路径的结构特性。但环形网络一出现这个结构就失效了。环网里两个节点之间存在两条甚至更多条供电路径功率流向不再唯一会出现环流和功率分叉。你用前推回代时根本不知道前是哪边、后是哪边硬要处理就得在环的某处把网络切开加入断口电压和注入功率的修正反复迭代逼近。这个思路理论上可行但实现复杂度成倍上升而且对多环网络切哪里、怎么切、怎么收敛都很头疼。我自己见过不少同学试图把前推回代强行扩展成能算环网的版本最后程序里塞满了对特定拓扑的断点处理逻辑换一个网络就得重新调。与其在树的逻辑里勉强补环不如直接换一套根本不依赖树枝结构的算法这就是牛拉法。1.2 牛拉法的数学内核把潮流问题变成多元求根问题牛拉法处理环网的思路非常直接先不管网络是树还是环把每个节点的功率平衡方程全部写出来然后求解这个非线性方程组。极坐标下节点注入功率的表达式是有功P_i V_i Σ V_k (G_ik cosθ_ik B_ik sinθ_ik)无功Q_i V_i Σ V_k (G_ik sinθ_ik - B_ik cosθ_ik)其中θ_ik θ_i - θ_kG和B来自节点导纳矩阵。每个PQ节点有P和Q两个方程每个PV节点只有一个P方程平衡节点的电压和相角给定不用参与迭代。把所有方程写到一起就是一个典型的F(X)0问题。牛拉法的解法是先给一组初值然后对F(X)做一阶泰勒展开得到线性修正方程组J·ΔX -F其中J就是雅可比矩阵。解出修正量ΔX更新状态再重复直到不平衡量小到可以接受。生活里有个差不多的例子猜一个人的体重先估一个数站上秤发现偏差根据多吃一顿大概涨多少这种敏感度信息来反向修正猜测多试几次就接近真实值。这里的敏感度就是雅可比矩阵。牛拉法之所以快是因为它在解附近是二阶收敛的通常几次迭代就能把不平衡量压到1e-8以下。1.3 为什么牛拉法天然适用于任意拓扑关键点在于牛拉法自始至终只需要两样东西节点导纳矩阵和功率平衡方程。导纳矩阵只是根据支路连接关系生成的环网和辐射网的区别在矩阵里不过是多几个非零元素没有任何树状结构的假设。所以只要你的程序里没有写死拓扑类型的判断逻辑牛拉法天然就能处理辐射网、单环网、双环网、多环网根本不需要为环单独做准备。这也是为什么我推荐核心判断标准看拓扑树状结构用前推回代能省事但只要网络带环直接上牛拉法省下的调试时间远远多于多写的那点代码。2. 数据接口怎么设计两张矩阵表描述任意环形网络2.1 节点与支路的输入约定要通用第一步就是输入格式不能绑死网络规模。我习惯用两个矩阵作为全部输入节点数量和支路数量完全取决于矩阵行数不用改任何代码。节点数据矩阵 bus每行八个字段字段含义bus(i,1)节点编号整数可以乱序bus(i,2)节点类型1为平衡节点2为PV节点3为PQ节点bus(i,3)注入有功功率P标幺值注入网络为正bus(i,4)注入无功功率Q标幺值PQ节点必须给定PV节点给初始值迭代中会重新计算bus(i,5)电压幅值初值V标幺值bus(i,6)电压相角初值θ弧度一般给0bus(i,7)PV节点无功上限Qmaxbus(i,8)PV节点无功下限Qmin支路数据矩阵 line每行七个字段字段含义line(k,1)首端节点编号line(k,2)末端节点编号line(k,3)电阻R标幺值line(k,4)电抗X标幺值line(k,5)线路对地充电电纳的一半B/2普通线路填实际值变压器支路填0line(k,6)变压器变比k普通线路填1变压器支路按高压侧/低压侧填写line(k,7)支路类型0为普通线路1为变压器这里说明一下变压器支路的变比定义k V_高压侧 / V_低压侧且变压器阻抗归算到低压侧即line矩阵中的R、X是低压侧的阻抗标幺值。这个约定必须统一否则导纳矩阵会错得很隐蔽。如果你暂时只算纯线路环网可以把变压器部分忽略所有支路的类型标志填0、变比填1即可。为什么用纯矩阵而不用GUI或者面向对象封装因为矩阵格式最简单可以直接在脚本里写死两个数组就开跑也方便后面用循环批量生成随机网络做测试。真要做成大工程矩阵格式也能作为底层数据源完全不冲突。2.2 自动生成节点导纳矩阵通用程序的地基就是自动形成节点导纳矩阵Y。给定bus和line程序循环每一条支路把导纳元素累加进矩阵对应位置。普通线路和变压器支路的累加规则不一样但都只需要局部信息跟网络整体拓扑无关。这段Matlab代码是完整可用的核心函数function Y formY(bus, line) % 根据节点表和支路表生成节点导纳矩阵 nb size(bus, 1); nl size(line, 1); Y zeros(nb, nb); % 先按节点数建满矩阵 for k 1:nl i line(k, 1); j line(k, 2); R line(k, 3); X line(k, 4); bHalf line(k, 5); tap line(k, 6); type line(k, 7); if type 0 % 普通线路等值π模型电纳取一半 y 1 / (R 1j * X); Y(i, i) Y(i, i) y 1j * bHalf; Y(j, j) Y(j, j) y 1j * bHalf; Y(i, j) Y(i, j) - y; Y(j, i) Y(j, i) - y; else % 变压器支路阻抗归算到j侧变比k高压/低压 yT 1 / (R 1j * X); Y(i, i) Y(i, i) yT / tap^2; Y(j, j) Y(j, j) yT; Y(i, j) Y(i, j) - yT / tap; Y(j, i) Y(j, i) - yT / tap; end end end注意普通线路里的bHalf是线路总充电电纳的一半对应等值π模型的半集中电容。如果是不计分布电容的短线路填0就行。2.3 节点类型统计与修正方程维度牛拉法的修正方程不是所有节点都参与的。平衡节点的电压幅值和相角都给定不迭代PV节点的电压幅值固定只有相角未知PQ节点两个未知量都参与。所以迭代前的统计工作很关键n size(bus, 1); slackIdx find(bus(:, 2) 1); pvIdx find(bus(:, 2) 2); pqIdx find(bus(:, 2) 3); nPv length(pvIdx); nPq length(pqIdx); % 修正方程维度 非平衡节点的相角数 PQ节点的电压幅值数 dim (n - 1) nPq;这个dim后面要用来初始化雅可比矩阵和右端项如果统计错了矩阵维度不匹配Matlab会直接报错。我建议在程序开头把n、nPv、nPq、dim都打印出来看一眼能提前暴露很多数据错误。3. 牛拉法求解核心不平衡量、雅可比矩阵与迭代主循环3.1 功率不平衡量的计算迭代的第一步是拿着当前电压幅值和相角用功率方程算一遍注入功率然后和给定值做差ΔP_i P_给定_i - P_计算_iΔQ_i Q_给定_i - Q_计算_i其中平衡节点不参与修正PV节点只保留ΔPPQ节点ΔP和ΔQ都保留。这个计算本身很简单两层循环遍历所有节点即可function [P, Q] calPower(V, theta, Y) G real(Y); B imag(Y); n length(V); P zeros(n, 1); Q zeros(n, 1); for i 1:n for k 1:n th theta(i) - theta(k); P(i) P(i) V(i) * V(k) * (G(i,k) * cos(th) B(i,k) * sin(th)); Q(i) Q(i) V(i) * V(k) * (G(i,k) * sin(th) - B(i,k) * cos(th)); end end end这里有个细节k等于i的对角项会自动包含进去因为夹角为0sin项为0cos项为1刚好对应G_ii V_i²和-B_ii V_i²不需要单独处理。3.2 雅可比矩阵四个子块的公式与代码雅可比矩阵是牛拉法程序里最核心、也最容易写错的部分。我在极坐标下采用ΔV/V作为电压修正变量这样做的好处是雅可比非对角元素恰好都是V_i V_j乘三角函数的形式公式很整齐。为了方便对照先把公式列出来这里的功率方程定义和上一节完全一致。雅可比矩阵分四个子块分别对应ΔP对Δθ、ΔP对ΔV/V、ΔQ对Δθ、ΔQ对ΔV/V子块i ≠ j有支路相连时i jH ∂P/∂θV_i V_j (-G_ij sinθ_ij B_ij cosθ_ij)-B_ii V_i² - Q_iN V·∂P/∂VV_i V_j (G_ij cosθ_ij B_ij sinθ_ij)G_ii V_i² P_iJ ∂Q/∂θV_i V_j (-G_ij cosθ_ij - B_ij sinθ_ij)P_i - G_ii V_i²L V·∂Q/∂VV_i V_j (G_ij sinθ_ij - B_ij cosθ_ij)Q_i - B_ii V_i²注意这里P_i和Q_i是当前迭代状态下用calPower算出的值不是给定的额定功率这一点漏了会直接导致结果错乱。组装雅可比矩阵的通用写法我推荐先建索引映射再按节点对循环填入。索引映射就是提前算好每个节点的theta行号和V行号function J buildJ(bus, Y, P, Q, V, theta, nPq) n size(bus, 1); dim (n - 1) nPq; J zeros(dim, dim); G real(Y); B imag(Y); % 建立节点到修正方程行号的映射 rowTheta zeros(n, 1); rowV zeros(n, 1); pos 1; for i 1:n if bus(i, 2) ~ 1 rowTheta(i) pos; pos pos 1; end end for i 1:n if bus(i, 2) 3 rowV(i) pos; pos pos 1; end end for i 1:n ti rowTheta(i); vi rowV(i); for k 1:n if i k % 对角元素 if ti 0 J(ti, ti) J(ti, ti) (-B(i,i) * V(i)^2 - Q(i)); % 对应H_ii end if ti 0 vi 0 J(ti, vi) J(ti, vi) (G(i,i) * V(i)^2 P(i)); % 对应N_ii end if vi 0 J(vi, vi) J(vi, vi) (Q(i) - B(i,i) * V(i)^2); % 对应L_ii end if vi 0 ti 0 J(vi, ti) J(vi, ti) (P(i) - G(i,i) * V(i)^2); % 对应J_ii end else if Y(i,k) ~ 0 tk rowTheta(k); vk rowV(k); th theta(i) - theta(k); Vij V(i) * V(k); % H_ik if ti 0 tk 0 J(ti, tk) J(ti, tk) Vij * (-G(i,k) * sin(th) B(i,k) * cos(th)); end % N_ik if ti 0 vk 0 J(ti, vk) J(ti, vk) Vij * (G(i,k) * cos(th) B(i,k) * sin(th)); end % J_ik if vi 0 tk 0 J(vi, tk) J(vi, tk) Vij * (-G(i,k) * cos(th) - B(i,k) * sin(th)); end % L_ik if vi 0 vk 0 J(vi, vk) J(vi, vk) Vij * (G(i,k) * sin(th) - B(i,k) * cos(th)); end end end end end end这段代码读起来比直接手写子矩阵长但它有一个明显好处完全由节点对循环决定网络怎么连都能算不需要在代码里写死任何拓扑结构。而且对角元素只需要叠加P_i和Q_i跟其他节点无关正好利用calPower的结果。3.3 迭代主循环组装右端项、解方程、更新状态主循环是标准的牛拉法流程。每次迭代做四件事算不平衡量、检查收敛、组装雅可比、解修正方程并更新tol 1e-8; maxIter 50; theta bus(:, 6); V bus(:, 5); for iter 1:maxIter [Pcal, Qcal] calPower(V, theta, Y); dP bus(:, 3) - Pcal; dQ bus(:, 4) - Qcal; % 平衡节点的dP不参与非PQ节点的dQ不参与 dP(slackIdx) 0; dQ(setdiff(1:n, pqIdx)) 0; % 组装右端项F F [dP(setdiff(1:n, slackIdx)); dQ(pqIdx)]; if norm(F, inf) tol fprintf(已收敛迭代次数: %d\n, iter); break; end % 组装雅可比 J buildJ(bus, Y, Pcal, Qcal, V, theta, nPq); % 用Matlab反斜杠求解 dX J \ F; % 更新状态量 nTheta n - 1; theta(setdiff(1:n, slackIdx)) theta(setdiff(1:n, slackIdx)) dX(1:nTheta); V(pqIdx) V(pqIdx) .* (1 dX(nTheta 1:end)); end这里用norm(F, inf)做收敛判据表示所有不平衡量的最大值小于容差就停止比只判断某一个节点更稳妥。迭代上限设50次正常环形网络平启动下5次左右就能收敛如果50次还不收敛通常意味着数据有问题或者网络本身接近电压崩溃点。4. 通用性打磨PV节点越限、编号重排与收敛控制4.1 PV节点无功越限自动转换前面处理的都是理想的PQ节点和PV节点。实际电网里PV节点通常代表发电机母线无功出力有上下限。当牛拉法迭代到某个PV节点时计算出的无功Qcal超出限值说明这个节点已经不能继续维持给定电压需要把它转成PQ节点把无功固定在限值上重新迭代。这个逻辑不做程序在真实算例里很容易不收敛。在每次迭代里加一段检查for i pvIdx if Qcal(i) bus(i, 7) || Qcal(i) bus(i, 8) fprintf(节点%d无功越限转为PQ节点Q%.4f\n, ... i, min(max(Qcal(i), bus(i,8)), bus(i,7))); bus(i, 2) 3; bus(i, 4) min(max(Qcal(i), bus(i,8)), bus(i,7)); end end % 越限后需要重新统计PQ、PV节点信息 pvIdx find(bus(:, 2) 2); pqIdx find(bus(:, 2) 3); nPv length(pvIdx); nPq length(pqIdx); dim (n - 1) nPq;越限转换后雅可比矩阵的维度会变化所以必须在循环体内更新统计量再进入下一次迭代。有些场景下后续迭代又可能让无功回到限值以内想恢复PV节点也是可以的但通常简化处理为先转PQ不再转回工程上够用。4.2 与节点编号无关内部只认矩阵行号通用性强还有一个很容易被忽视的指标节点编号能不能乱很多网上的程序隐式假设节点编号从1到n连续且有序一旦你的网络编号跳号或者不从1开始程序就崩了。这个问题在数据接口层就要处理。比较好的做法是程序入口就做一次重编号把bus矩阵中出现的节点编号统一映射到1到n的连续行号同时把line矩阵里的首末端编号也映射过去。简单做法是检查line中的节点编号是否都在bus表里存在然后直接用bus的行号代替节点编号参与所有计算。这样节点编号只是一个标签随便怎么编都不影响结果。我在formY函数里就是用bus行号直接索引的所以只要把line的首末端按实际编号换成bus行号就完成了脱敏。建议在main程序一开头做这个映射并且把映射后的矩阵回显出来方便检查数据有没有输错。4.3 收敛性控制与发散问题排查用平启动所有PV/PQ节点V1θ0在绝大多数环形网络上都能顺利收敛但遇到重负载、线路电阻电抗比很高、或者接近极限送电能力的算例牛拉法也会耍脾气。我常用的几个控制手段阻尼因子。把修正量乘以一个系数θ和V都变成走半步而不是走一步。代码里可以改成theta alpha * dTheta、V V .* (1 alpha * dV)alpha从0.5开始试收敛稳定后再慢慢调到1。代价是迭代次数变多但能救回很多难收敛的算例。迭代中监视不平衡量的变化趋势。如果norm(F, inf)在振荡而不是单调下降基本可以判定初值离解太远或者雅可比矩阵状态很差这时候与其硬迭代不如回头检查数据。检查是否出现孤立节点。如果某个节点没有任何支路连接Y阵该行全零雅可比矩阵奇异Matlab会警告矩阵接近奇异。这种情况收敛不了是正常的属于网络数据本身不合法。排查发散问题时我的经验是先缩到最小复现把网络节点数减到3个、只保留一条环路逐步加支路/加载荷看哪一步开始不收敛。比盯着大网络日志猜快得多。5. 环形算例实测迭代表现与三个容易踩的坑5.1 一个四节点环形网络的实际表现用一个简单的四节点环网做测试网络结构是1-2-3-4-1形成的闭环。数据如下节点类型PQ1平衡002PQ0.80.33PQ0.60.24PQ0.40.1支路分别是1-2R0.02X0.06、2-3R0.08X0.24、3-4R0.06X0.18、4-1R0.05X0.15四条线路都不计分布电容。平启动后迭代过程大致如下迭代次数最大不平衡量max|ΔP,ΔQ|收敛状态17.8e-1继续22.3e-2继续32.4e-4继续43.1e-8继续51e-9收敛从第二轮开始不平衡量几乎以二次速率下降这是牛拉法的典型特征。换到辐射网测试收敛行为也完全一样说明程序对拓扑确实不敏感。我还做过把节点编号打乱重排的测试只要映射逻辑正确收敛轨迹基本一致结果完全相同。5.2 三个我自己踩过的坑第一个坑是雅可比矩阵符号不统一。早期的程序里功率方程用的是从网络吸取功率为正的定义从网上抄来的雅可比公式却是注入为正混在一起迭代直接发散。这个问题的排查思路很简单找一个两节点手算例子把第一轮迭代的雅可比元素打印出来对照公式一项一项查。符号错误通常会在非对角元素上暴露得非常明显。第二个坑是变压器变比定义反了。变比到底是高压侧除以低压侧还是反过来不同的书定义不同。我做了一个含变压器的简单算例发现电压全都偏高反复核对后才发现是Y阵里的变比放反了。解决方法是把定义写死在注释里并且用一个非常简单的双节点变压器网络单独测试一次性确认Y阵的对角元素和互导纳符合手册值再往大网络上用。第三个坑是修改变量选择混乱。有的资料用ΔV做变量有的用ΔV/V做变量两者对应的雅可比公式不一样N和L的对角元素差一个V的倍数。只要初值偏离1比较远结果就会出问题。这个坑特别隐蔽因为形式上看程序能跑但收敛结果偏得莫名其妙。建议全程统一使用ΔV/V版本并保证电压更新的写法是V V .* (1 dV)别混用。还有一个容易被忽略的经验Matlab里求修正方程不要自己写高斯消元直接用反斜杠运算符J \ F。节点规模到几百上千的时候反斜杠会自动选择合适解法比手写消元更稳也省去一堆索引调试的麻烦。5.3 再补充一个技巧把雅可比矩阵的稀疏结构利用起来牛拉法在规模较小或中等网络时用full矩阵完全足够。但如果你打算拿这个程序去算上百节点的环形配网或输电网full矩阵的雅可比会越来越大内存和速度都不太好看。这时可以把Y阵改成sparse存储雅可比也按sparse分配其他逻辑完全不变Matlab反斜杠会自动切换到稀疏求解。改完以后大算例的速度提升非常明显代码改动量却很小只是把zeros改成spalloc而已。如果你以后打算把牛拉法扩展到更大的网络建议提前在数据量大的测试里验证一下Y阵的非零结构是否正确再上稀疏化。毕竟稀疏矩阵肉眼看起来不如full矩阵直观排错时要多留几个心眼。