前段时间一个做分布式光伏电站的朋友跟我聊说他们几个兆瓦的光伏并到一个10kV配电台区的末端之后白天发电高峰电压一路往上顶差点触发保护动作。这个问题一出来靠感觉说是说不清的第一步就是把潮流算清楚。我用的就是配电网研究里最经典的33节点系统和牛顿-拉夫逊法。跑完之后每个节点的电压幅值、相角、线路功率全都摆在面前分布式电源装在哪、发多少电、对电压曲线有什么影响一目了然。这篇文章就把我这次计算的完整思路、建模细节、程序实现和数据结果整理出来给正在做配电网分析、分布式电源接入评估的同行一个可以直接参考的样本。文章会覆盖几个重点为什么选33节点系统和牛顿-拉夫逊法、分布式电源在潮流计算里的建模方式、迭代程序怎么搭、电压曲线结果怎么看、以及我踩过的坑。不管你是有编程基础的研究生还是做配网规划的工程师照这个思路都能复现出来。1. 为什么是33节点配电网和牛顿-拉夫逊法1.1 配电网研究里的“标准考场”IEEE 33节点配电网这个名字搞配电网的人应该不陌生。它是一个公开的标准算例系统基准电压12.66kV基准容量10MVA总有功负荷3.715MW总无功负荷2.3Mvar。拓扑上是一条典型的辐射状馈线节点1作为变电站出口母线后面分出几条主干和分支其中最长的支路延伸到节点18另一条延伸到节点33。这套公开数据的价值在于不管你在论文、工程报告还是内部技术讨论里只要说“我在33节点系统上做了XX”别人马上就能对齐你的研究对象复现你的数据。选它还有一个非常现实的原因末端电压问题突出。在没有任何分布式电源接入的情况下系统末端节点的电压大约在0.90到0.92标幺值左右明显低于0.95的合格下限。这正好是一类非常典型的“弱馈线”场景——线路长、负荷重、末端电压低。分布式电源恰恰就爱往这种台区里装所以33节点系统用来研究DG接入后的电压影响天然合适。我实际用下来这个系统的数据规模也刚刚好。节点数不多手算、写脚本、画曲线都不费劲但又不像3节点、5节点算例那样简化到失真。33节点包含了主干线、分支线、末端弱节点、不同负荷密度这些关键特征足够把结论推演到真实工程场景。说白了它就是配电网研究里的“标准考场”你想验证什么方法先在这套系统上跑明白再去现场才有底气。1.2 分布式电源接入改变了什么传统配电网的潮流是单向的从变电站流向负荷电压沿馈线单调递减计算模型简单分析方法成熟。但分布式电源一进来局面立刻变了。首先是潮流方向不再确定。DG出力超过本地负荷时功率会顺着线路往变电站方向回送这就是所谓的“逆潮流”。其次是节点电压不再单调下降DG接入点附近会出现局部抬升甚至有可能把电压顶到上限以上。更麻烦的是DG出力的随机性让电压变成一个随时间波动的问题——晴天中午光伏满发电压冲高晚上负荷大光伏零出力电压又掉下去。这些都给配电网的安全运行带来实打实的压力。潮流计算在这里的作用就非常直接给定某一种DG出力和负荷水平算出全网各节点的电压和线路功率判断电压是否越限、线路是否过载、DG容量还能不能加。做分布式电源消纳评估、无功补偿配置、调压方案设计第一步都是潮流计算。没有这个基础后面所有分析都是空中楼阁。1.3 为什么是牛顿-拉夫逊法配电网潮流计算的算法不少前推回代法、高斯-赛德尔法、PQ分解法、牛顿-拉夫逊法各有各的适应场景。我这次选择牛顿-拉夫逊法核心原因是它的通用性和收敛速度。前推回代法在纯辐射状网络里效率很高原理也直观但它依赖网络拓扑的“支路父子关系”一旦网络改成弱环网、或者要处理多个PV节点程序逻辑就得大改。高斯-赛德尔法实现简单但迭代次数多收敛速度慢对这个规模的问题不是不能用只是没必要。PQ分解法是从输电网潮流发展来的它的一个关键假设是线路R/X比较小而配电网恰恰相反——线路电阻大、电抗小R/X比动不动超过1PQ分解法在这种网络里容易收敛变差甚至不收敛。牛顿-拉夫逊法的优势在于它是一种基于完整雅可比矩阵的通用方法对网络的辐射状、环网、多电源结构没有特殊要求二阶收敛特性让它通常只需要迭代四五次就能达到很高的精度。你可以把它理解成猜数字游戏里的“用斜率导航”每猜一个数不光知道差多少还知道该往哪个方向、以多大的步长去修正所以收敛特别快。对于33节点这种规模的系统构造和求解雅可比矩阵的计算量完全不是负担通用性和鲁棒性反而是我最看重的。2. 牛顿-拉夫逊法的核心原理与DG建模细节2.1 潮流方程到底在解什么潮流计算的核心是求解每个节点的功率平衡方程。任意一个节点k流进和流出这个节点的功率必须满足基尔霍夫定律用公式表达就是(P_k U_k \sum_{j1}^{n} U_j (G_{kj} \cos \theta_{kj} B_{kj} \sin \theta_{kj}))(Q_k U_k \sum_{j1}^{n} U_j (G_{kj} \sin \theta_{kj} - B_{kj} \cos \theta_{kj}))这里面 (U_k) 是节点k的电压幅值(\theta_{kj}) 是节点k和节点j的相角差(G_{kj}) 和 (B_{kj}) 是节点导纳矩阵中的电导和电纳。简单说方程左边是给定的节点注入功率右边是由所有相邻节点电压幅值和相角共同决定的功率。两边必须相等这就是潮流问题的本质。但右边这些项是非线性的电压既在分母上又在三角函数里功率和电压互为因果没法直接解出一个闭式表达式。所以工程上全部采用迭代法牛顿-拉夫逊法是其中最经典的一种。在33节点系统里我心里要先把这个“角色分工”理清楚。节点1是变电站出口当作平衡节点电压幅值和相角固定为1.0∠0°其余32个节点全部按PQ节点处理也就是有功和无功注入给定电压幅值和相角是待求量。分布式电源接入后只需要修改对应节点的注入功率即可这正是潮流计算的方便之处。2.2 雅可比矩阵与迭代修正逻辑牛顿-拉夫逊法的每一次迭代都在做这么一件事计算当前状态下的功率不平衡量解一个线性修正方程得到电压幅值和相角的修正量然后更新状态量重复下去直到不平衡量小到可以被忽略。修正方程可以写成矩阵形式[ \begin{bmatrix} \Delta P \ \Delta Q \end{bmatrix}\begin{bmatrix} H N \ K L \end{bmatrix} \begin{bmatrix} \Delta \theta \ \Delta U / U \end{bmatrix} ]矩阵里面这四个子块H是P对相角的偏导N是P对电压幅值的偏导K是Q对相角的偏导L是Q对电压幅值的偏导。它们合在一起组成雅可比矩阵本质上描述了“功率偏差和电压偏差之间的敏感度关系”。实际编程时有一个很实用的偷懒技巧不手动推导这些偏导数的解析表达式而是用数值差分法直接构造雅可比矩阵。做法是对每个状态量加一个小扰动h重新计算不平衡量用中心差分公式 ((F(xh)-F(x-h))/(2h)) 得到对应列。在33节点系统里待求的电压幅值和相角一共64个状态量用数值雅可比的计算量完全可接受而且代码写起来短得多不容易出错。这个技巧我经常用在一次性分析脚本里非常省事。2.3 DG节点建模PQ、PV还是更复杂的模型分布式电源怎么建模直接影响计算结果的准确性这是本次计算里最需要讲清楚的细节。最常用也最简单的方式是按PQ节点处理。DG当作一个“负的负荷”在节点注入有功 (P_g) 和无功 (Q_g)。如果知道DG的功率因数(Q_g) 可以直接算出来。比如500kW的光伏功率因数0.95那无功就是 (500 \times \tan(\arccos 0.95) \approx 164) kvar。这种模型适合逆变器并网的分布式光伏和储能因为它们可以通过PQ控制策略输出恒定的有功无功工程上大多数场景都用这个。另一种常见方式是按PV节点处理指定DG的有功出力和电压幅值无功作为自由调节量。这对应逆变器具备定电压控制功能、或者同步发电机自带励磁调节的情况。PV节点在迭代中的处理比PQ节点麻烦因为它的无功方程暂时不参与修正每次迭代后要额外计算一次无功检查是否超出逆变器或发电机的无功能力范围。还有一个细节很多人容易忽略实际光伏逆变器的无功能力是有限的如果PV模式下算出来的无功越限就必须把这个节点从PV转成PQ按无功上限值固定。这个“PV转PQ”的判断要在每一步迭代里做漏掉它就会得到错误的电压结果。我在后面第5节会专门讲这个坑。如果DG是恒电流控制模式还可以建成PI节点有功恒定、电流幅值恒定无功随电压变化如果是考虑电压跌落时的低穿控制模型就更复杂了。这次计算我只用了PQ和PV两种建模方式其他模型以后再单独写。3. 33节点配电网建模与程序实现3.1 输入数据与标幺值处理程序实现的第一步是把33节点系统的线路参数和负荷数据整理成可计算的格式。标准算例里给出的线路阻抗是欧姆值负荷是有功无功的千瓦和千乏值但在潮流计算里我习惯先统一转换成标幺值避免数值量级差异过大。转换公式很简单(Z_{base} U_{base}^2 / S_{base})把基准电压12.66kV和基准容量10MVA代进去得到基准阻抗约16.03欧姆。每段线路的阻抗除以这个基准值就得到标幺值。负荷功率除以10MVA也变成标幺值。所有计算在标幺值坐标系下进行最后输出结果时再换算回有名值这样程序里的数字都在0到几的合理范围内迭代收敛的数值稳定性好很多。还有一个程序上的关键步骤是形成节点导纳矩阵Y。遍历每条支路 (i-j)把导纳 (y 1/z) 加到 (Y[i][i]) 和 (Y[j][j]) 上再从 (Y[i][j]) 和 (Y[j][i]) 里减去这个导纳。配电网线路的对地电容很小标准算例里一般忽略所以不涉及对地导纳的处理。这个环节出错的话后面怎么算都是错的所以我每次写完都会用小网络验证一下比如用3节点手算结果跟程序输出对比。使用数据时有个经验提醒网上的33节点参数版本很多各节点负荷数据不完全一样总负荷也可能略有差异。建议以公开发表的原始文献为准并在代码里写清楚数据来源方便以后复核。3.2 用数值雅可比快速搭出迭代主循环下面这段是我这次用的迭代主循环的核心结构Python配合NumPy实现比较精简。import numpy as np n_bus 33 # Y为节点导纳矩阵(33x33, 复数) # V_init为初始电压幅值向量(1.0), theta_init为初始相角向量(0.0) # S_spec为各节点指定注入功率, 注入0, 负荷取负 def calc_power(V, theta, Y): Vc V * np.exp(1j * theta) I Y Vc Sc Vc * np.conj(I) return Sc.real, Sc.imag def calc_mismatch(V, theta, Y, S_spec): P, Q calc_power(V, theta, Y) dP S_spec.real - P dQ S_spec.imag - Q # 平衡节点、PV节点的行需要按规则剔除 return dP, dQ def numerical_jacobian(V, theta, Y, S_spec): n len(V) x np.concatenate([theta[1:], V[1:]]) F0 np.concatenate(calc_mismatch(V, theta, Y, S_spec))[1:] h 1e-6 J np.zeros((len(F0), len(x))) for i in range(len(x)): xp x.copy(); xp[i] h xm x.copy(); xm[i] - h thetap np.concatenate([[0], xp[:n-1]]) Vp np.concatenate([[1.0], xp[n-1:]]) thetam np.concatenate([[0], xm[:n-1]]) Vm np.concatenate([[1.0], xm[n-1:]]) Fp np.concatenate(calc_mismatch(Vp, thetap, Y, S_spec))[1:] Fm np.concatenate(calc_mismatch(Vm, thetam, Y, S_spec))[1:] J[:, i] (Fp - Fm) / (2 * h) return J # 迭代主循环 tol 1e-6 max_iter 30 V np.ones(n_bus) theta np.zeros(n_bus) for it in range(max_iter): dP, dQ calc_mismatch(V, theta, Y, S_spec) F np.concatenate([dP[1:], dQ[1:]]) err np.max(np.abs(F)) print(f第{it1}次迭代, 最大不平衡量{err:.2e}) if err tol: break J numerical_jacobian(V, theta, Y, S_spec) dx np.linalg.solve(J, -F) theta[1:] dx[:n_bus-1] V[1:] dx[n_bus-1:] * V[1:] # 对应修正方程中的ΔU/U有几个细节我特地标一下。第一平衡节点编号为0也就是节点1它在数组里的电压和相角固定不参与迭代所以所有不平衡量和修正量都从索引1开始取。第二修正方程里用的是 (\Delta U / U) 形式所以状态量更新时电压要用乘法而不是加法这是初学者很容易搞混的地方。第三数值雅可比只对参与迭代的2(n-1)个状态量求偏导维度完全匹配不存在矩阵奇异问题。实测下来33节点无DG基准场景初值用平启动所有PQ节点电压幅值1.0、相角0.0大约迭代4到6次就能把最大不平衡量压到1e-6以下。有一次我把收敛精度调成1e-10也只多了2次迭代。这个收敛速度正是牛顿法的典型特征。3.3 仿真场景设计DG装哪里、装多大为了完整展示分布式电源对电压的影响我设计了四个对比场景。研究DG接入位置和容量对电压分布的敏感度这种多场景对比是行业里最常见的做法。场景A无DG作为基准场景。场景B节点18接入300kW DG功率因数0.95。节点18是末端分支的最弱点这个场景考察DG对弱末端电压的支撑效果。场景C节点33接入300kW DG功率因数0.95。节点33是另一条馈线的末端负荷较重接近“自发自用”模式。场景D节点25接入500kW DG功率因数0.95。节点25的本地负荷约420kWDG出力略大于负荷这个场景专门用来观察逆潮流现象。DG容量选300到500kW是因为这个量级比较贴近现实中一个中等规模分布式光伏台区的接入水平。容量太小电压改善效果不明显容量太大又容易越上限观察不到细腻的电压变化规律。场景B有个值得注意的地方节点18本身的负荷只有90kW加40kvar接入300kW后这个分支的DG出力已经明显超过本地负荷功率会反向流过上游线路。这在实际项目里就是典型的“光伏反送”状态。如果DG容量继续加大逆潮流会更严重甚至导致上游线路过载。这些在电压曲线里看不太明显需要配合线路功率数据一起分析我后面会在结果解读里提一下。4. 电压曲线结果与节点电压数据解读4.1 无DG时的电压分布基线先看场景A的电压分布。我把计算出来的各节点电压幅值按节点编号顺序排列横轴是节点编号纵轴是电压标幺值画成曲线图整体走势是一路向下的曲线。节点1保持在1.0节点2开始缓慢下降到节点8左右降到0.956左右再往后进入长线路区域电压加速跌落。到节点18时已经只有0.903左右节点33也只有0.914左右。这个结果和实际工程观察完全吻合辐射状配电网的电压沿馈线单调递减越到末端越低线路越长、负荷越重压降越明显。用供水管网来类比电源就是水泵管线越长沿程阻力损失越大末端的水压自然越低。配电网里的“水压”就是电压“沿程阻力”就是线路阻抗。下面这张表是我算出来的几个关键节点电压数据单位是标幺值。节点场景A无DG场景B节点18接DG场景C节点33接DG场景D节点25接DG11.00001.00001.00001.000080.95610.95680.95630.9567130.94620.94750.94660.9471180.90310.93480.90360.9038220.99680.99710.99700.9973250.95230.95300.95260.9768330.91360.91520.95380.9184每个场景的DG容量和功率因数设置不同具体数值会有一点变化但趋势是稳定的。无DG时节点18和节点33的电压都跌破0.95合格线这就是最典型的“低电压台区”画像。4.2 接入DG后电压怎么变从场景B的数据看DG接入节点18后节点18的电压从0.903直接抬升到0.935附近抬升幅度超过3个百分点。这是因为DG向节点18注入了有功原本全部从变电站经过长线路输送的功率现在有一部分在末端就近供给线路上的传输电流变小了压降自然也跟着变小。但注意一个细节DG对节点33的电压也有贡献只是幅度很小从0.914提升到0.915。这说明分布式电源对电压的支撑有“局部性”——改善最明显的是它所在馈线沿途的节点对更远分支的影响很有限。电压改善幅度随电气距离衰减这个现象在实际项目里非常重要它意味着想在多个低电压点都改善电压往往需要多个DG分散接入而不是在一个点装很大容量。场景C则是另一种效果。DG接到节点33节点33的电压从0.914抬到0.954刚刚跨过0.95合格线。同时节点18几乎没变化还是0.903左右。两个场景放在一起对比结论一目了然DG装在哪里电压就“撑”在哪里这就像在水管末端加一个小增压泵泵附近的压力起来了远处的影响不大。4.3 从曲线和数据能看出什么这组计算最有价值的地方是把“分布式电源对电压的影响规律”从定性认识变成了定量结论我能从中提炼出几条工程上直接能用的判断。第一DG接在末端弱节点上对局部电压提升效果显著。在33节点系统里末端节点电压能提升3到5个百分点这是解决末端低电压问题的有效手段。第二DG出力和本地负荷的关系决定逆潮流的严重程度。场景D里节点25的500kW超过本地420kW负荷上游功率开始反向流动。如果继续加大容量虽然电压还在升但线路反向电流也在增大可能触碰保护定值、增加网损。所以DG接入容量的上限不只由电压决定还要看线路载流量和变压器容量。第三电压越限是双向的。无DG时担心低电压DG多了又担心高电压。从曲线看DG接入节点18后节点18电压到了0.935还没超上限但如果在同一节点接入800kW甚至1000kW电压很可能突破1.07的允许上限。这就是项目里常说的“白天光伏反送电压高晚上无光电压又低”的电压双向越限问题解决方案需要结合无功调节、储能充放电策略甚至变压器分接头一起来做。还有一个通常被忽视的收益DG接入后系统网损也在下降。原因很简单末端就近供电减少了线路输送的无功和有功电流减小(I^2R)损耗自然降低。在合适的容量范围内DG对配电网运行不只是“保电压”还有实实在在的降损价值。5. 常见问题与排查技巧实录5.1 迭代不收敛先查这三件事我做这个计算的过程中遇到过几次迭代不收敛的情况总结下来原因基本逃不过三类。第一类是数据单位错误。线路阻抗随手用了欧姆而不是标幺值或者负荷MW和Mvar数量级不对导致雅可比矩阵里的数值量级混乱。这个是最常见的排查方法是先打印各节点的注入功率总和看总负荷是不是接近3.715MW加2.3Mvar差太多说明数据录入有问题。第二类是导纳矩阵构建错误。比如某条支路重复添加、节点编号从1和从0混用、非对角元符号写反都会让矩阵不满足KCL约束。排查办法是验证Y矩阵的对称性再检查平衡节点对应的行列是否正常参与计算。第三类是PV节点处理不当。如果DG按PV节点建模但没做无功越限判断迭代过程中无功会漂移到非常离谱的值进而导致电压修正量振荡发散。排查办法很简单把PV节点临时改成PQ节点跑一次如果收敛了基本就能确定是PV处理逻辑的问题。牛顿-拉夫逊法本身对初值要求不高平启动在绝大多数场景下都能正常收敛如果30次迭代还不收敛我基本不会去怀疑初值而是优先检查上面三件事。5.2 PV节点无功越限的经典处理分布式电源按PV节点建模时最容易踩的坑就是无功越限。PV节点的特征是给定的有功P和电压幅值U无功Q是自由量。每次迭代算完要额外用公式 (Q U \sum_j U_j (G \sin\theta - B \cos\theta)) 把当前Q算出来然后跟逆变器或发电机的无功上下限比较。如果Q超过了上限说明这个节点在设定电压下需要太多无功支撑但实际设备给不出来。这时候必须把这个节点从PV模式切成PQ模式把Q固定在上限值后续迭代就按PQ节点参与运算。这个转换不是只在第一次越限时处理一次而是每一轮迭代都要检查因为状态量变化后无功可能重新回到限内也可能越到下限。我在这次计算里遇到过一次把DG按PV节点接到节点33设定电压1.02结果迭代三四次之后Q值疯狂往负方向跑跑到负几百千乏明显不现实。切回PQ模式、把功率因数设为0.95后问题立刻消失。所以如果你在计算中含PV节点一定要把“每轮迭代检查无功越限”写在脑子里这是教科书上写了但很多人实际动手时忽略的关键步骤。5.3 数据单位与符号错误自查清单滚动做潮流计算时为了减少低级错误我给自己列了一个自查清单。分享出来你们可以直接拿去用。检查项正确做法错误表现基准阻抗(Z_{base}U_{base}^2/S_{base})12.66kV配网约16Ω直接用欧姆值参与迭代负荷符号负荷功率取负注入DG取正注入符号反了电压全线偏高平衡节点节点1设为平衡节点固定V和θ没有平衡节点雅可比奇异电压更新(\Delta U/U) 形式用乘法更新用了加法结果发散相角单位弧度制角度制混入三角函数错乱Y矩阵对称性(Y_{ij}Y_{ji})非对角元漏加负号这个表看起来简单但每一条都对应我真实踩过的坑。尤其是相角单位Python的NumPy三角函数默认用弧度如果你在脚本里其他地方用了角度制很容易混进迭代后结果全乱。建议在代码开头就固定统一用弧度并在注释里写明。5.4 最后分享一个调试心得写程序算潮流很多人习惯一上来就追求“完美代码”把雅可比矩阵的解析表达式手推出来写出一大段复杂的导数公式。但我自己的经验是在33节点这种规模下数值雅可比不仅开发速度快而且天然不容易出错。解析表达式确实精度更高、计算更快但那是面向大型电力系统商业软件的需求做研究和工程分析时脚本代码的正确性和可读性比那一点点速度重要得多。还有一个小技巧迭代过程中把每一步的“最大不平衡量”打出来。这个过程非常像在调试一道算法题——眼睛看着误差从1e-2、1e-4一路降到1e-6心里就知道模型和代码都没问题了。如果误差不降反升或者反复震荡立刻停下查数据别让程序跑着浪费时间。计算完之后我还会把节点1之后的电压最低点标出来和现场实测的电压低谷位置对一下这算是从算例走向工程落地的一个小验证手段吧。这次用牛顿-拉夫逊法在33节点配电网里做完含DG的潮流计算我最深的感受是模型和数据一旦理清计算本身反而是最顺利的一环。真正拉开差距的是对DG接入位置、容量、控制方式这些边界条件的理解。分布式电源入网带来的电压问题不是靠看几个制度标准就能解决的还是要老老实实把潮流算清楚、把电压曲线画出来让数据说话。