电流注入型牛拉法：配电网潮流计算的统一建模与实现

这段时间我在做一套配电网潮流计算程序。一开始按教材老路走的是功率注入型牛拉法程序写了大半突然发现要面对一堆恒电流负荷、电流源型分布式电源的时候代码越改越别扭一会儿要把恒阻抗负荷折算进导纳矩阵一会儿又要单独算恒电流负荷的等效功率残差逻辑散得到处都是。后来我索性把整个框架换成了电流注入型牛拉法这才算把负荷模型和分布式电源统一进了一套方程里。这篇博文就把我从数学模型推导到程序落地的完整过程理一遍中间包含了我踩过的几个坑以及最终整理出来的一套可运行实现。这个内容适合三类人看一是课程设计或者研究生阶段需要自己手写潮流程序的同学二是做配电网、微电网平台开发需要把负荷模型和电源模型灵活接进程序的工程师三是已经会用功率注入型牛拉法、想搞清楚两种写法差在哪的进阶用户。我会把“为什么这么设计”和“实际代码怎么落”放在同等位置不搞纸上谈兵。1. 为什么把程序从功率注入型换成电流注入型1.1 牛拉法在潮流计算里解决的本质问题潮流计算的目标其实很朴素给定网络拓扑、线路参数、各节点的注入功率或者电压约束求全网各节点的电压幅值和相角。有了电压之后线路潮流、变压器潮流、网损、平衡节点出力这些量都可以顺着一股脑算出来。但问题在于节点功率平衡方程是非线性的。一个节点的注入功率取决于本节点电压和所有相邻节点电压的乘积电压之间还带着三角函数关系。这个非线性方程组没有解析解只能迭代。牛顿-拉夫逊法的思路是先猜一组初值把非线性方程在当前点做一阶泰勒展开形成一个线性方程组解出修正量更新状态再重复。这个“拉”的过程本质上是在“拉”电压状态量让不平衡量残差逐步逼近零。残差怎么定义直接决定了雅可比矩阵长什么样也决定了程序的数据流怎么组织。1.2 教材主流的功率注入型牛拉法到底哪里不够用绝大部分教材讲牛拉法潮流用的都是功率注入型。在极坐标下状态量是电压幅值V和相角θ残差是节点有功不平衡量ΔP和无功不平衡量ΔQ。对PQ节点两个方程对PV节点一个有功方程加一个电压幅值约束。算法成熟MATPOWER这类工具也是这么实现的收敛性经受了大量工程检验。但我在实际项目里遇到了两个比较头疼的问题。第一个是负荷模型。配电网的负荷通常不是纯恒功率很多地方用ZIP模型描述也就是一部分是恒阻抗一部分是恒电流一部分是恒功率。功率注入型程序里恒阻抗负荷需要折算成等效导纳并进导纳矩阵恒电流负荷需要按当前电压折算成等效功率放进P、Q给定值。问题来了每迭代一步电压都在变恒电流负荷折算出来的等效功率也在变你得在残差计算里动态更新这部分代码写起来很绕。第二个是分布式电源。很多并网逆变器在稳态分析里的控制特性是电流源型的比如PQ控制模式下直接给电流指令或者采用恒电流保护限流特性。用功率注入型去模拟电流源型电源等于把一个天然用电流表达的东西硬掰成功率表达每次算完还得验证等效功率对不对。这两个问题有一个共同解放弃“功率不平衡量”作为残差改用“电流不平衡量”。这就是电流注入型牛拉法。1.3 电流注入型的适用边界电流注入型牛拉法并不是要取代功率注入型它更适合特定场景。适合的场景包括配电网三相潮流计算因为配电网负荷模型复杂ZIP模型很常见含大量电流控制型逆变器、储能变流器的微电网和主动配电网分析需要同时把恒功率负荷、恒电流负荷、恒阻抗负荷、电流源型电源统一建模的程序框架。不太适合的场景也很明确大电网在线潮流计算功率注入型结合稀疏技术已经做到极致没必要换还有就是刚接触潮流计算的初学者功率注入型的物理直觉更直接先学那个再学电流注入型会更顺。我的经验是如果模型的复杂性上来了程序的架构成本反而比数学推导更值得重视。电流注入型虽然雅可比推导麻烦一点但数据流是统一的所有设备模型最后都汇到“电流规格值”和“导纳矩阵”这两个入口扩展起来特别舒服。2. 电流注入型牛拉法的数学骨架2.1 直角坐标下的电压、电流、功率转化电流注入型牛拉法我建议用直角坐标实现原因后面会说。先约定基本变量。节点电压写为复数和直角坐标两种形式[ V_i e_i jf_i ]节点导纳矩阵元素[ Y_{ij} G_{ij} jB_{ij} ]节点注入功率 (S_i P_i jQ_i)电源注入为正负荷注入为负。根据复数功率公式 (S VI^*)给定功率对应的注入电流规格值为[ I_i^{spec} \left( \frac{S_i}{V_i} \right)^* \frac{P_i - jQ_i}{V_i^*} \frac{P_i - jQ_i}{e_i - jf_i} ]把分母实数化之后实部和虚部分别是[ I_i^{spec, re} \frac{P_i e_i Q_i f_i}{e_i^2 f_i^2} ][ I_i^{spec, im} \frac{P_i f_i - Q_i e_i}{e_i^2 f_i^2} ]这个公式是电流注入型的核心。注意这里 (P_i)、(Q_i) 是给定注入功率不是负荷功率建模的时候如果输入的是负荷吸收功率要取相反数。另一方面根据节点电压方程网络中流动的电流是[ I_i^{net} \sum_j Y_{ij} V_j \sum_j (G_{ij} jB_{ij})(e_j jf_j) ]展开后[ I_i^{net, re} \sum_j (G_{ij} e_j - B_{ij} f_j) ][ I_i^{net, im} \sum_j (G_{ij} f_j B_{ij} e_j) ]电流注入型的残差定义就是规格电流减去网络电流[ \Delta I_i^{re} I_i^{spec, re} - I_i^{net, re} ][ \Delta I_i^{im} I_i^{spec, im} - I_i^{net, im} ]2.2 三类节点的残差构造PQ、PV与平衡节点现在进入关键点三类节点分别怎么列方程。PQ节点最直接。状态量是 (e_i) 和 (f_i)给定的是 (P_i) 和 (Q_i)所以两个电流残差方程都用上。PV节点稍微复杂。给定的是有功 (P_i) 和电压幅值 (V_i^{spec})无功 (Q_i) 是未知的所以 (I_i^{spec, im}) 里的 (Q_i) 没法直接给。如果硬要用电流残差就得把 (Q_i) 设成状态量多一个未知数多一个方程程序复杂度上升。工程上更常用的做法是PV节点保留两个功率类残差[ \Delta P_i P_i^{spec} - P_i^{calc} ][ \Delta V_i^2 (V_i^{spec})^2 - (e_i^2 f_i^2) ]其中 (P_i^{calc}) 是直角坐标下的有功表达式[ P_i^{calc} e_i \sum_j (G_{ij} e_j - B_{ij} f_j) f_i \sum_j (G_{ij} f_j B_{ij} e_j) ]这样PV节点的两个方程分别来自有功平衡和电压幅值约束跟功率注入型牛拉法在形式上接上了。平衡节点不参与迭代。它的 (e)、(f) 是给定值潮流算完之后再用电压和导纳矩阵算出它的注入功率。算完之后一定要检查平衡节点功率是否在合理范围内这是验证全局功率平衡最直接的信号。还有一点需要提醒如果有节点的电压幅值约束和PQ给定量都出现但系统里可能出现平衡节点那么总未知量个数等于 (2 \times (N-1))方程数也是 (2 \times (N-1))一除以二刚好对上。2.3 雅可比矩阵的实用组装方式电流注入型雅可比矩阵的解析推导比功率注入型繁琐主要是电流规格值表达式的分母里有 (e_i^2f_i^2)求偏导要套除法法则。我在第一版程序里用了一个比较省事的办法有限差分。所谓的有限差分就是对每一个状态量 (x_k)加一个小扰动 (h)重新算一遍残差然后用[ J_{:,k} \approx \frac{F(x_k h) - F(x_k)}{h} ]来计算雅可比矩阵第 (k) 列。扰动量的经验值是 (10^{-6}) 到 (10^{-7}) 标幺值具体大小要保证电压幅值在0.95~1.1之间时残差变化量不会因为浮点误差被吃掉。用有限差分的最大好处是残差函数写完之后雅可比矩阵自动就有了不用手动推导一堆偏导公式不容易错。代价是每次迭代要多算 (2(N-1)) 次残差对几万节点的系统来说不划算但几百节点以内的教学和验证系统完全够用。后面需要提速再针对残差表达式的每一项做解析求导机器学习里叫“先跑通再优化”电力系统里这叫“先验证模型再优化算法”。当然如果你的系统对性能要求高建议直接用解析雅可比。在这里给两条关键结论第一对PQ节点的电流残差来说网络电流部分对 (e_k)、(f_k) 的偏导就是导纳矩阵元素的实部和虚部这正是电流注入型适合稀疏化的原因。第二规格电流部分的偏导最终只出现在对角线块上非对角线元素全部来自网络电流项。这个稀疏结构比功率注入型更规整。2.4 一次牛顿迭代的完整计算链条把上面的方程合在一起一次牛顿迭代的步骤如下。第一步给定当前所有节点的 (e)、(f)更新每个节点的电压平方和电流规格值。第二步计算残差向量 (F)此时PQ节点往里放电流残差实部和虚部PV节点放有功残差和电压幅值平方残差平衡节点不放。第三步用有限差分或者解析求导组装雅可比矩阵 (J)。第四步解线性方程组[ J \Delta x -F ]第五步更新状态量[ x x \Delta x ]第六步检查收敛判据。常用的判据是残差向量最大绝对值小于 (10^{-6}) 标幺或者计算迭代步前后状态量变化量小于某个阈值。我习惯两者都看因为残差很小但状态量还在缓慢漂移的情况在病态网络里出现过。这就是电流注入型牛拉法的完整骨架。公式不多但写代码的时候要把节点类型、变量索引、残差顺序完全对齐不然雅可比矩阵的行列对应关系很容易对不上。3. 程序实现里最容易出错的部分3.1 数据模型节点和支路怎么设计才不坑写潮流程序第一个容易翻车的点是节点和支路的数据结构。很多初学者喜欢把节点属性揉在一个大字典里节点类型、给定功率、电压都给串起来结果程序一复杂查错查到崩溃。我最终采用的方式是节点用一个简单的类或者命名元组字段包括节点编号、节点类型、给定有功、给定无功、给定电压幅值、电压初值、是否是平衡节点。支路数据用表格存每条支路包含起点、终点、电阻、电抗、对地电纳、变压器变比。对齐方式用“节点编号从0开始”的数组下标这样导纳矩阵的索引和列表下标一致省去一堆减一操作。实际项目里还会遇到双绕组变压器、三绕组变压器、移相器这些特殊设备。第一版不建议全支持先把普通线路和双绕组变压器跑通后面的设备类型在支路上打一个type字段慢慢扩展。3.2 导纳矩阵从支路数据到Y矩阵的累加逻辑导纳矩阵形成的核心是累加。每条支路遍历一次把贡献累加到自导纳和互导纳上。普通线路用 (\pi) 型等值模型自导纳累加(Y_{ii} y_{series} j\frac{b}{2})(Y_{jj} y_{series} j\frac{b}{2})互导纳累加(Y_{ij} -y_{series})(Y_{ji} -y_{series})其中 (y_{series} 1 / (r jx))。变压器支路要用变比折算。如果变比 (k) 定义为起点侧对终点侧的电压比那么[ Y_{ii} \frac{y_{series}}{k^2}, \quad Y_{jj} y_{series}, \quad Y_{ij} -\frac{y_{series}}{k}, \quad Y_{ji} -\frac{y_{series}}{k} ]这个公式在不同教材里变比定义可能不一样我一开始因为变比定义确认不到位算出来的结果怎么都对不上MATPOWER。后来统一在程序开头的注释里写明“变比定义为PSS/E风格的起点侧变比”再没出过这类问题。恒阻抗负荷可以直接并进导纳矩阵对角元。如果负荷吸收功率是 (P_{load} jQ_{load})那么等效并联导纳[ y_{load} \frac{P_{load} - jQ_{load}}{V_{spec}^2} ]加在对应节点的对角元上。这一步做好之后恒阻抗负荷就彻底从迭代变量里消失了它消耗的功率会自动体现在网络方程里。这是电流注入型框架最舒服的一点。3.3 变量编号与稀疏矩阵求解的坑变量编号是最容易被忽视的性能瓶颈。最简单的做法是把所有节点的 (e)、(f) 放进一个大数组迭代时只更新自由节点的位置。这个逻辑正确但如果做稀疏矩阵必须注意编号顺序对稀疏分解性能影响很大。实际工程里如果节点数超过几百个一定不要用稠密矩阵做 LU 分解内存和计算量都会爆炸。MATLAB里直接构造sparse矩阵Python里用scipy.sparse求解用sparse.linalg.spsolve。配电网节点编号如果按照馈线支路顺序排稀疏分解的填充会很友好如果编号特别乱可以用AMD或最小度排序算法预处理一下。还有一个细节残差向量的顺序必须和雅可比矩阵行顺序一致而且和状态量数组的顺序一致。我自己的习惯是状态量数组先排所有PQ节点的 (e)、(f)再排所有PV节点的 (e)、(f)残差顺序和它一一对应。这样虽然PV节点的方程类型不一样但矩阵组装逻辑可以统一走同一个索引表。3.4 阻尼步长与收敛控制牛拉法在初值不好或者系统接近电压崩溃点时很容易震荡甚至发散。我这个程序里加了一个最简单的阻尼策略每次求出的修正量 (\Delta x)不一定全量加上而是先乘一个阻尼因子 (\lambda)。先试 (\lambda 1)如果残差范数变小就接受这一步如果变大了就把 (\lambda) 减半再试最多试几次。这就是回溯线搜索几行代码的事能解决很多发散问题。阻尼因子对配电网特别重要因为配电网R/X比高导纳矩阵对角占优程度不如输电网牛顿法容易走冤枉路。加阻尼之后急剧震荡的情况明显减少。收敛判据我建议同时监控两个量残差最大绝对值和状态量修正量最大绝对值。工业级程序常用 (10^{-6}) 这个量级如果只是验证算法 (10^{-5}) 也够用。4. 一个可运行的算例实现4.1 测试系统数据为了让代码可以直接跑我设计了一个三节点系统覆盖三种节点类型。系统基准容量取100 MVA。节点1是平衡节点电压给定 (1.06\angle 0^\circ)节点2是PV节点出力 (50) MW电压给定 (1.0) p.u.节点3是PQ节点负荷 (100) MW、(50) Mvar这个负荷吸收功率所以在程序里注入功率用负数。三条支路参数如下起点终点电阻(p.u.)电抗(p.u.)对地电纳(p.u.)120.030.090.03130.020.060.03230.040.120.02这个系统参数比较温和用平启动应该三四次迭代就能收敛。4.2 导纳矩阵生成函数先写一个导纳矩阵函数输入节点数和支路列表输出复导纳矩阵。重点是把变压器和非变压器支路分开处理。import numpy as np def build_ybus(n_bus, lines): Y np.zeros((n_bus, n_bus), dtypecomplex) for line in lines: i, j, r, x, b, k line if r 0 and x 0: y_series 0.0 else: z complex(r, x) y_series 1.0 / z if k is None or k 0: # 普通线路 y_shunt complex(0, b / 2.0) Y[i, i] y_series y_shunt Y[j, j] y_series y_shunt Y[i, j] -y_series Y[j, i] -y_series else: # 变压器线路 Y[i, i] y_series / (k * k) Y[j, j] y_series Y[i, j] -y_series / k Y[j, i] -y_series / k return Y这个函数里没处理恒阻抗负荷的并入因为负荷信息在节点表里放在残差函数之外一起管理会更清晰。如果需要并入恒阻抗负荷在组装导纳矩阵之后再循环一次节点数组把负荷导纳累加到对角元即可。4.3 残差计算函数残差函数是整个程序的灵魂。输入是当前所有节点的 (e)、(f) 数组以及节点定义输出是自由节点对应的残差向量。def calc_residual(e, f, buses, Y): n len(buses) V_sq e ** 2 f ** 2 Ie_net np.zeros(n) If_net np.zeros(n) for i in range(n): Ie_net[i] np.sum(Y[i].real * e - Y[i].imag * f) If_net[i] np.sum(Y[i].real * f Y[i].imag * e) F_list [] for i in range(n): bus_type buses[i][type] if bus_type PQ: P buses[i][P] Q buses[i][Q] Ie_spec (P * e[i] Q * f[i]) / V_sq[i] If_spec (P * f[i] - Q * e[i]) / V_sq[i] F_list.append(Ie_spec - Ie_net[i]) F_list.append(If_spec - If_net[i]) elif bus_type PV: P_spec buses[i][P] V_spec buses[i][V_spec] P_calc (e[i] * Ie_net[i] f[i] * If_net[i]) F_list.append(P_spec - P_calc) F_list.append(V_spec ** 2 - V_sq[i]) elif bus_type SLACK: # 平衡节点不加入残差 pass return np.array(F_list)注意这里PQ节点的 (P) 和 (Q) 是节点注入功率。如果定义节点数据时给的是负荷功率千万记得取负数。程序跑不出来先别急着查雅可比先把符号搞对这是最容易犯的低级错误。4.4 雅可比矩阵与主循环用有限差分组装雅可比矩阵然后写主迭代循环。def build_jacobian(e, f, buses, Y): n len(buses) free_idx [i for i in range(n) if buses[i][type] ! SLACK] m len(free_idx) x_all np.zeros(2 * n) x_all[0:n] e x_all[n:2 * n] f x0 np.concatenate([e[free_idx], f[free_idx]]) F0 calc_residual(e, f, buses, Y) J np.zeros((m, m)) h 1e-7 for k in range(m): x_pert x0.copy() x_pert[k] h e_pert e.copy() f_pert f.copy() e_pert[free_idx] x_pert[0:m] f_pert[free_idx] x_pert[m:2 * m] F_pert calc_residual(e_pert, f_pert, buses, Y) J[:, k] (F_pert - F0) / h return J, free_idx主循环里每次迭代组装雅可比解方程更新状态检查残差。代码省略了收敛判据和阻尼步长但骨架很清楚先算残差、再算雅可比、然后解线性方程组。求解用np.linalg.solve节点多就换成scipy.sparse.linalg.spsolve。4.5 跑出来的结果怎么验证跑通三节点之后第一件事不是直接换大系统而是做几个必要的自检。第一看收敛曲线。正常的牛拉法迭代次数在3到5次残差应该是单调下降的。如果残差来回震荡多半是变量索引对不上、雅可比行顺序错位、或者阻尼策略没加。第二看平衡节点功率。算完之后用下面公式计算平衡节点注入功率S_bal complex(e[0], f[0]) * np.sum(Y[0, :] * (e 1j * f)).conj()对于三节点算例平衡节点发出的功率应该等于总负荷功率加上全网网损。如果平衡节点功率异常巨大或者为负先查节点功率符号有没有搞反。第三把PQ节点电压和MATPOWER的结果对比。MATPOWER里有现成的runpf同一套数据填进去跑一遍电压幅值误差在 (10^{-5}) 以内就说明核心逻辑没问题。我拿三节点这个简单系统做过实验平启动条件下(e)、(f) 初值分别取1.0和0.0迭代4次左右残差就降到 (10^{-7}) 以下结果和MATPOWER完全对得上。5. 调试中踩过的坑与排查思路5.1 初值、发散与“电流分母爆掉”电流注入型有一个独特的坑规格电流表达式的分母里有 (e_i^2 f_i^2)也就是电压幅值的平方。如果初值电压设成0第一轮迭代直接除零如果某些节点电压在迭代过程中掉到0附近量级也会爆炸。我一开始用全0作为初值结果一跑就崩。后来改成平启动所有节点 (e1.0)(f0.0)平衡节点用给定值这才稳定。这里多说一句某些病态系统里平启动也会出问题可以试试冷启动加更大的阻尼或者用带修边的两步启动法。排查发散问题的时候把每一步残差范数打出来是最直接的。如果残差前半段下降后半段突然上扬基本可以断定是雅可比矩阵数值出了问题如果一开始就乱跳重点查初值和符号。5.2 PV节点无功越限的自动转换这是一个特别容易被忽略的工程细节。PV节点的无功出力 (Q_i^{calc}) 在迭代过程中是隐式算出来的但发电机无功是有上下限的。迭代收敛之后一定要检查每个PV节点的无功功率是否在给定范围内。如果越限比如上机发电机的无功已经顶到极限那这个节点实际上不能再维持给定的电压幅值应该自动转成PQ节点并把无功固定在限值上然后重新迭代。我在程序里实现了一个两层循环外层循环负责节点类型转换内层循环做牛拉法迭代。每次牛拉法收敛后检查所有PV节点的无功如果越限就改类型、置一个type_changed标记重新做一轮牛拉法。这个流程虽然多算几次但结果才是物理上可信的。5.3 高R/X比系统的收敛应对配电网的R/X比通常比输电网高很多线路电阻和电抗差不多甚至电阻更大牛拉法的收敛域会明显变窄。这时候单纯靠平启动加阻尼有时候还是不够。我试过几个办法。第一个是提高导纳矩阵的数值稳定性把线路电纳并入对角元时注意标幺基准一致第二个是缩小伪随机初值范围比如让 (e) 在0.95到1.05之间、(f) 在-0.05到0.05之间而不是从很偏的值出发第三个是用连续潮流或者逐步加载的思路先跑一个轻负荷到重负荷的延续过程用上一档结果做下一档初值。配电网还有一种常见做法是把恒功率负荷的一部分按恒阻抗折算进导纳矩阵降低迭代过程中的非线性度。反正电流注入型框架里恒阻抗并入和恒功率进规格电流都是现成的搭配起来非常灵活。5.4 正确性验证的几条铁律程序写到最后验证手段不能只靠一个三节点算例。我在团队里的习惯是搞一个回归测试集把IEEE 4节点、5节点、14节点、30节点系统都放进去每次改完代码全量跑一遍。任何一个系统电压结果和MATPOWER对不上肯定是改坏了。另外两个物理校验必须做。一个是全网功率平衡。全网总发电减去总负荷必须等于总网损误差在 (10^{-5}) 标幺以内。这里总发电要把平衡节点功率也算进去。另一个是线路潮流双向一致。从节点i流向节点j的功率和从节点j流向节点i的功率两者之差必须等于线路损耗。如果不一致说明线路潮流公式和导纳矩阵定义不匹配。个人经验是潮流计算程序真正的门槛往往不在牛拉法公式本身而在数据约定和索引管理。只要把“功率符号”“变比定义”“节点编号”这三件事在代码注释里写死后面所有调试都会顺很多。电流注入型牛拉法因为残差物理意义直观、负荷模型接入自然一旦跨过初期的数学适应期实际工程里比功率注入型更耐造。