简介电力系统潮流计算是电力系统分析中的核心环节用于确定给定运行条件下各节点电压、支路电流及功率分布。这套PPT面向电气工程专业学生、电力系统分析初学者及需进行上机计算的工程师以计算机计算方法为主线依次讲解开式网、简单闭式网及复杂电力系统的潮流计算数学模型其中开式网可采用代数方法直接求解简单闭式网需考虑环路并运用回路电流法处理。资源压缩包内为1个PPT演示文稿大小约1.66MB内容围绕第11章“潮流计算计算机计算”展开涵盖节点电压方程、平衡节点/PV节点/PQ节点分类、功率平衡方程及电压、功率上下限等约束条件同时对高斯迭代法、高斯-赛得尔法和牛顿-拉夫逊法等求解非线性代数方程组的算法进行逐步推演并涉及PV节点无功越界处理、平衡节点不参与迭代等实现细节便于读者对照公式理解编程思路与收敛过程。已有184人浏览学习可作为课程复习、自学备考或编写潮流计算程序的入门参考。1. 电力系统潮流计算要解决的不只是解方程看到电力系统分析教材里那份“潮流计算.ppt”很多人的第一反应是又要面对一堆矩阵和迭代。实际做过电网断面计算后瓶颈通常不在解方程而在把节点类型、负荷方向和导纳数据组织成一套不自洽的输入再让迭代算法恰好收敛到物理上可行的解。潮流计算本身可以说得很简单给定发电机出力和负荷分布求出全网每个节点的电压幅值、相角进而算出每条线路和变压器的功率流动。它的结果服务于电压越限评估、网损统计、保护定值配合、新能源消纳分析这些具体业务。写代码、调参数、判断结果是否合理比背数学公式更能决定你能不能把一次潮流算明白。2. 潮流计算的数学模型从节点功率方程到雅可比矩阵2.1 三类节点决定方程个数和可解性潮流计算第一步不是列方程而是给每个节点定类型。电力系统分析里最常见的分类只有三种平衡节点、PV 节点、PQ 节点。平衡节点通常选在骨干电厂或者主网枢纽它的电压幅值和相角是已知的有功和无功是待求的PV 节点是发电机节点有功出力和电压幅值给定无功是待求的PQ 节点是负荷节点或无功可调但电压不控制的节点有功和无功给定电压幅值和相角待求。节点类型已知量待求量典型对象平衡节点VθV、θP、Q主网骨干电源PV 节点P、VQ、θ发电机组PQ 节点P、QV、θ负荷、新能源场站为什么必须有一个平衡节点核心原因是电力系统的网损在求解之前是未知的。所有发电机出力之和必须覆盖所有负荷加上网损而网损只有算出各支路电流之后才知道所以必须留一个节点不指定有功让它吸收整个系统的功率差额。如果把所有发电机有功都固定方程个数会超过未知量个数潮流无解。写成数学形式对每个非平衡节点列有功方程对每个 PQ 节点再列无功方程。系统有 n 个节点、m 个 PV 节点则未知量个数为 (n-1) 个相角加上 (n-1-m) 个 PQ 节点电压幅值方程式个数也正好是 (n-1)(n-1-m)方程闭合可以求解。2.2 极坐标下的功率失配方程与雅可比四个分块在极坐标下节点导纳矩阵写成 Y G jB节点 i 的有功和无功注入为P_i V_i Σ V_j (G_ij cosθ_ij B_ij sinθ_ij)Q_i V_i Σ V_j (G_ij sinθ_ij - B_ij cosθ_ij)其中 θ_ij θ_i - θ_j。潮流求解的实质是让残差趋近于零定义功率失配方程ΔP_i P_i_spec - P_i(V,θ) 0ΔQ_i Q_i_spec - Q_i(V,θ) 0对 PV 节点只列 ΔP对平衡节点两者都不列对 PQ 节点两者都列。牛顿-拉夫逊法的思路是把上述非线性方程在当前点做一阶泰勒展开用线性方程组解出修正量再更新电压和相角反复迭代直到残差满足精度。这里的雅可比矩阵按极坐标变量可以分为四个分块H 对应有功对相角的偏导N 对应有功对电压幅值的偏导K 对应无功对相角的偏导L 对应无功对电压幅值的偏导。非对角元的形式比较规整对角元需要额外叠加本节点的注入项。实际写程序时最容易错的就是对角元符号因为不同教材对 θ_ij 的正负定义不同照抄公式前务必和自己在潮流方程里用的夹角定义对齐。2.3 牛拉法和 PQ 分解法的选用边界牛顿-拉夫逊法具备二阶收敛速度正常情况下 3 到 5 次迭代就能把残差压到 1e-6 以下代价是每次迭代都要重新计算并分解雅可比矩阵。对于几百节点的小系统这个成本可以忽略对于上万节点的大电网稀疏分解的开销会明显影响计算速度。PQ 分解法是在牛拉法基础上做的近似。输电网支路的电阻远小于电抗节点电压幅值变化对有功潮流影响很小相角差变化对无功潮流影响也很小于是可以把雅可比矩阵固定成两个常数矩阵分别在 P-θ 迭代和 Q-V 迭代里使用。这样一来迭代次数会从几次增加到十几次甚至几十次但每轮迭代不再需要重新组装矩阵和重新分解总耗时往往比牛拉法低一个量级。算法单轮成本典型迭代次数适用场景牛顿-拉夫逊法高每轮重算并分解 J3~5 次任意 r/x 比值、中小规模PQ 分解法低只做三角回代10~30 次高压输电网、大规模在线分析选择哪种算法不是一个性能判断题而是收敛性和计算资源的权衡。下面两章分别给出两种算法的最小可运行实现。3. 用 Python 从零实现牛拉法潮流计算3.1 先建导纳矩阵支路参数表与 build_y潮流程序的第一份输入不是迭代算法而是支路参数。我用一个四节点系统作为教学算例包含 1 个平衡节点、1 个 PV 节点和 2 个 PQ 节点。所有数值都用标幺值100 MVA 基准下可以直接复现。支路首端末端rp.u.xp.u.对地电纳 b/2p.u.L1120.0050.100.02L2130.0100.150.03L3240.0150.200.04L4340.0080.100.02节点数据为节点 1 是平衡节点V1.06θ0节点 2 是 PV 节点P0.50V1.04节点 3 是 PQ 节点P-0.70Q-0.30节点 4 是 PQ 节点P-0.55Q-0.20。负荷的 P 和 Q 都带负号这是约定俗成的方向注入网络为正值消耗为负值。import numpy as np def build_y(n, branches): y np.zeros((n, n), dtypecomplex) for i, j, r, x, b_half in branches: z r 1j * x # 支路阻抗 ys 1.0 / z # 支路导纳 y[i, i] ys 1j * b_half # 首端自导纳加对地电纳 y[j, j] ys 1j * b_half # 末端自导纳加对地电纳 y[i, j] - ys # 互导纳为负 y[j, i] - ys return y这段代码把一个支路表变成节点导纳矩阵。自导纳由与该节点相连的所有支路导纳之和再加上该支路两端各一半的对地电纳互导纳等于支路导纳的负数。直接使用 dense 矩阵只适合教学和小系统验证工程实现应该用 scipy.sparse 的 CSC 格式保存后续稀疏 LU 分解速度会快几个数量级。3.2 牛拉法迭代核心数值雅可比版本我一般先把收敛流程跑通再用解析雅可比替换。下面的版本用中心差分组装雅可比矩阵逻辑最直观不容易被符号约定带偏。搞清楚流程后再按 2.2 节的四个分块改成解析实现。def calc_power(y, v, th): vc v * np.exp(1j * th) sc vc * np.conj(y vc) return sc.real, sc.imag def newton_lf(y, ref, pv, pq, p_spec, q_spec, v_spec, eps1e-6, max_iter15): n y.shape[0] v np.ones(n) th np.zeros(n) v[list(ref) list(pv)] [v_spec[i] for i in list(ref) list(pv)] theta_idx [i for i in range(n) if i not in ref] # 需要求相角的节点 v_idx list(pq) # 需要求电压幅值的节点 def mismatch(theta, volt): p, q calc_power(y, volt, theta) d [] for i in theta_idx: d.append(p_spec[i] - p[i]) for i in v_idx: d.append(q_spec[i] - q[i]) return np.array(d) x0 np.concatenate([th[theta_idx], v[v_idx]]) for it in range(max_iter): th_full np.zeros(n); th_full[theta_idx] x0[:len(theta_idx)] v_full np.ones(n); pv_and_ref set(pv) | set(ref); v_full[list(ref) list(pv)] [v_spec[i] for i in list(ref) list(pv)] v_full[v_idx] x0[len(theta_idx):] r mismatch(th_full, v_full) if np.max(np.abs(r)) eps: break J np.zeros((len(x0), len(x0))) h 1e-6 for k in range(len(x0)): xp x0.copy(); xp[k] h thp np.zeros(n); thp[theta_idx] xp[:len(theta_idx)] vp np.ones(n); vp[list(ref) list(pv)] [v_spec[i] for i in list(ref) list(pv)] vp[v_idx] xp[len(theta_idx):] xm x0.copy(); xm[k] - h thm np.zeros(n); thm[theta_idx] xm[:len(theta_idx)] vm np.ones(n); vm[list(ref) list(pv)] [v_spec[i] for i in list(ref) list(pv)] vm[v_idx] xm[len(theta_idx):] J[:, k] (mismatch(thp, vp) - mismatch(thm, vm)) / (2 * h) dx np.linalg.solve(J, -r) x0 x0 dx return v_full, th_full这里把整个状态向量拆成两块所有非平衡节点的相角放在前面所有 PQ 节点的电压幅值放在后面。每次计算残差时再从这两块拼回完整的电压向量。之所以不直接对导纳矩阵求导是为了先避开雅可比符号约定等程序能够稳定收敛后再把中心差分替换成解析式。中心差分每轮迭代要多算很多次功率只适合教学验证不能直接用在大规模算例上。np.linalg.solve(J, -r)这一步是解线性方程组-r 是牛顿法的修正方向。要注意残差量级和 eps 的单位都是标幺值100 MVA 基准下 1e-6 对应 0.1 kW已经远高于工程精度需求。3.3 收敛判据和四个关键参数参数常见取值说明收敛容差 eps1e-6 p.u.指最大功率失配不指电压变化量最大迭代次数10~15牛拉法正常 3~5 次收敛初始电压平启动 1.0∠0°适用于绝大多数输电网断面初值电压幅值PV 和平衡节点按给定值其他节点取 1.0提示判断收敛时不要同时看残差和电压变化量。潮流方程是功率残差方程只看 max(|ΔP|, |ΔQ|) 是否小于 eps不要在迭代早期混入有功功率不平衡判断否则会导致迭代提前终止在非收敛点。对小系统而言平启动初值几乎不会出问题。但系统重载、存在长距离重负荷线路时牛拉法也可能发散。这时候要检查的不只是初值还包括导纳矩阵里是否漏掉了对地电纳、并联电容器等无功补偿设备。4. PQ 分解法输电网潮流计算的快速通道4.1 两个近似和两个矩阵B′、B″PQ 分解法的基础是高压输电网中支路电抗远大于电阻电压幅值和相角两个问题近似解耦。有功功率主要和相角差有关无功功率主要和电压幅值有关。于是牛拉法里耦合的 2x2 分块雅可比被拆成两个独立的常系数矩阵。第一个矩阵 B′ 用于 P-θ 迭代由支路电抗倒数组成忽略电阻和对地电纳。第二个矩阵 B″ 用于 Q-V 迭代保留了对地电纳的影响通常直接取导纳矩阵虚部的负值。二者只需要在迭代开始前组装并分解一次后续每轮迭代只做两次线性回代计算量大减。这种近似不是对所有系统都成立。配电网 r/x 比值较大P-V 和 Q-θ 的耦合不可忽略PQ 分解法很容易发散或收敛到错误解。所以该方法只建议在高压输电网、主网分析这类场景使用。对比项牛拉法PQ 分解法雅可比矩阵每轮重算、重分解常数矩阵分解一次迭代次数3~5 次10~30 次有功与无功迭代联立求解交替求解对 r/x 比值的依赖不敏感较敏感4.2 PQ 分解法迭代代码与牛拉法的行为差异def pq_lf(y, ref, pv, pq, p_spec, q_spec, v_spec, eps1e-6, max_iter50): n y.shape[0] bpp -np.imag(y) # 取导纳矩阵虚部并取负作为 Q-V 矩阵 bp np.zeros((n, n)) for i in range(n): for j in range(n): if i ! j and abs(y[i, j]) 1e-12: bp[i, j] -np.imag(1.0 / (1.0 / y[i, j])) # 用支路电抗倒数近似 # 简化做法对无变压器网络bp -1/x 组装对角与非对角 bp np.zeros((n, n)) for i in range(n): for j in range(i 1, n): if abs(y[i, j]) 1e-12: xij -1.0 / np.imag(y[i, j]) bp[i, i] 1.0 / xij bp[j, j] 1.0 / xij bp[i, j] - 1.0 / xij bp[j, i] - 1.0 / xij th np.zeros(n) v np.ones(n) v[list(ref) list(pv)] [v_spec[i] for i in list(ref) list(pv)] p_idx [i for i in range(n) if i not in ref] q_idx list(pq) for it in range(max_iter): p, q calc_power(y, v, th) dp np.array([p_spec[i] - p[i] for i in p_idx]) dpv dp / v[p_idx] dth np.linalg.solve(bp[np.ix_(p_idx, p_idx)], dpv) th[p_idx] dth p, q calc_power(y, v, th) dq np.array([q_spec[i] - q[i] for i in q_idx]) dqv dq / v[q_idx] dv np.linalg.solve(bpp[np.ix_(q_idx, q_idx)], dqv) v[q_idx] dv if np.max(np.abs(dp)) eps and np.max(np.abs(dq)) eps: break return v, th这段代码先组装 B′阵再用np.ix_取出与未知量对应的子块。每轮迭代里先固定电压幅值解相角再固定相角解电压幅值两个子问题交替推进。dp / v[p_idx]和dq / v[q_idx]是 PQ 分解法的标准右端项处理方式把功率失配量除以当前电压幅值让线性化后的方程更接近常数矩阵形式。牛拉法是一轮迭代同时解出相角和电压幅值PQ 分解法则把一轮拆成两个半步。因此 PQ 分解法迭代次数多但内部矩阵不需要重新分解整体效率在规模足够大时优势明显。调试早期发现不收敛时可以把每轮的 dp、dq 最大分量打出来观察是停留在 P 侧还是 Q 侧判断瓶颈在哪个子问题。4.3 参数经验与无功越限处理参数建议值注意事项收敛容差1e-6 p.u.与牛拉法等精度标准一致最大迭代次数30~50PQ 分解法不应设 10 次以内的上限初始电压平启动 1.0∠0°重负荷网络可先用牛拉法得初值PV 无功越限转 PQ 并固定 Q注意避免频繁类型切换注意工程计算里 PV 节点的无功越限是高频坑。发电机无功达到上限后电压幅值已经无法维持设定值此时必须把节点类型从 PV 改成 PQ并以 Q_upper 作为给定无功。如果下一轮迭代又向限值内回落不要马上切回 PV否则会出现节点类型振荡迭代始终无法收敛。小系统里 PQ 分解法和牛拉法结果应完全一致因为收敛到的是同一个物理平衡点。若两个算法给出的支路功率差异超过 0.1%多半是 B′、B″ 矩阵与导纳矩阵的单位定义不一致先检查对地电纳是否被重复计入。5. 潮流计算不收敛时按这三个方向排查5.1 先查数据再查迭代迭代发散时第一件事不是调初值和收敛因子而是回到输入数据。导纳矩阵的对角元素是否包含所有与节点相连的支路是检查频率最高的错误点。孤岛节点在导纳矩阵里表现为整行几乎为零牛顿法会在该节点出现单步修正量巨大、残差反而增大的现象。另一个常见问题是负荷符号写反。负荷的 P、Q 应为负值如果把负荷写成注入整个断面会变成一台“虚拟发电机”潮流可能收敛到电压整体升高约 1.05 的结果物理上看起来像是系统在输出大量无功。接口程序里统一负荷方向并在计算前打印全网有功、无功加和作为粗校验能拦住大部分数据错误。5.2 三个必须做的结果校验收敛只代表残差为零不代表结果一定正确。第一项校验是功率平衡平衡节点计算得到的 P 应约等于全网负荷总 P 加上网损总 P偏差应小于收敛精度对应有名值。第二项校验是支路功率对称性对每条支路从两端分别计算的功率数值方向相反且两个流向的有功功率之和应为正这个正值就是支路损耗不可能为负。第三项校验是电压合理性PQ 节点电压落在 0.95~1.05 范围内是常规输电网的合理区间与这个区间偏离过大时即使迭代收敛也要复核无功补偿配置。5.3 一次收敛后中途振荡的处理经验现场分析中还会遇到一种情况前几次迭代都很顺利到第 8 轮左右残差在 1e-4 附近不再下降。这时优先怀疑 PV 节点无功越限的切换逻辑写在了迭代循环内部导致节点类型每两轮被重写一次。常见的修复是记录节点类型切换状态切换后的节点在后续至少 5 轮迭代中强制保持 PQ 类型不做反向切换给电压幅值留出足够的收敛空间。如果仍然振荡可以在 Q-V 迭代的修正量上加入阻尼系数把v[q_idx] dv改成v[q_idx] 0.8 * dv。牛拉法不需要这个手段但 PQ 分解法在两个子问题交替时容易出现过冲小幅阻尼是从工程经验里最容易学到的技巧。本文还有配套的精品资源点击获取
简介电力系统潮流计算是电力系统分析中的核心环节用于确定给定运行条件下各节点电压、支路电流及功率分布。这套PPT面向电气工程专业学生、电力系统分析初学者及需进行上机计算的工程师以计算机计算方法为主线依次讲解开式网、简单闭式网及复杂电力系统的潮流计算数学模型其中开式网可采用代数方法直接求解简单闭式网需考虑环路并运用回路电流法处理。资源压缩包内为1个PPT演示文稿大小约1.66MB内容围绕第11章“潮流计算计算机计算”展开涵盖节点电压方程、平衡节点/PV节点/PQ节点分类、功率平衡方程及电压、功率上下限等约束条件同时对高斯迭代法、高斯-赛得尔法和牛顿-拉夫逊法等求解非线性代数方程组的算法进行逐步推演并涉及PV节点无功越界处理、平衡节点不参与迭代等实现细节便于读者对照公式理解编程思路与收敛过程。已有184人浏览学习可作为课程复习、自学备考或编写潮流计算程序的入门参考。1. 电力系统潮流计算要解决的不只是解方程看到电力系统分析教材里那份“潮流计算.ppt”很多人的第一反应是又要面对一堆矩阵和迭代。实际做过电网断面计算后瓶颈通常不在解方程而在把节点类型、负荷方向和导纳数据组织成一套不自洽的输入再让迭代算法恰好收敛到物理上可行的解。潮流计算本身可以说得很简单给定发电机出力和负荷分布求出全网每个节点的电压幅值、相角进而算出每条线路和变压器的功率流动。它的结果服务于电压越限评估、网损统计、保护定值配合、新能源消纳分析这些具体业务。写代码、调参数、判断结果是否合理比背数学公式更能决定你能不能把一次潮流算明白。2. 潮流计算的数学模型从节点功率方程到雅可比矩阵2.1 三类节点决定方程个数和可解性潮流计算第一步不是列方程而是给每个节点定类型。电力系统分析里最常见的分类只有三种平衡节点、PV 节点、PQ 节点。平衡节点通常选在骨干电厂或者主网枢纽它的电压幅值和相角是已知的有功和无功是待求的PV 节点是发电机节点有功出力和电压幅值给定无功是待求的PQ 节点是负荷节点或无功可调但电压不控制的节点有功和无功给定电压幅值和相角待求。节点类型已知量待求量典型对象平衡节点VθV、θP、Q主网骨干电源PV 节点P、VQ、θ发电机组PQ 节点P、QV、θ负荷、新能源场站为什么必须有一个平衡节点核心原因是电力系统的网损在求解之前是未知的。所有发电机出力之和必须覆盖所有负荷加上网损而网损只有算出各支路电流之后才知道所以必须留一个节点不指定有功让它吸收整个系统的功率差额。如果把所有发电机有功都固定方程个数会超过未知量个数潮流无解。写成数学形式对每个非平衡节点列有功方程对每个 PQ 节点再列无功方程。系统有 n 个节点、m 个 PV 节点则未知量个数为 (n-1) 个相角加上 (n-1-m) 个 PQ 节点电压幅值方程式个数也正好是 (n-1)(n-1-m)方程闭合可以求解。2.2 极坐标下的功率失配方程与雅可比四个分块在极坐标下节点导纳矩阵写成 Y G jB节点 i 的有功和无功注入为P_i V_i Σ V_j (G_ij cosθ_ij B_ij sinθ_ij)Q_i V_i Σ V_j (G_ij sinθ_ij - B_ij cosθ_ij)其中 θ_ij θ_i - θ_j。潮流求解的实质是让残差趋近于零定义功率失配方程ΔP_i P_i_spec - P_i(V,θ) 0ΔQ_i Q_i_spec - Q_i(V,θ) 0对 PV 节点只列 ΔP对平衡节点两者都不列对 PQ 节点两者都列。牛顿-拉夫逊法的思路是把上述非线性方程在当前点做一阶泰勒展开用线性方程组解出修正量再更新电压和相角反复迭代直到残差满足精度。这里的雅可比矩阵按极坐标变量可以分为四个分块H 对应有功对相角的偏导N 对应有功对电压幅值的偏导K 对应无功对相角的偏导L 对应无功对电压幅值的偏导。非对角元的形式比较规整对角元需要额外叠加本节点的注入项。实际写程序时最容易错的就是对角元符号因为不同教材对 θ_ij 的正负定义不同照抄公式前务必和自己在潮流方程里用的夹角定义对齐。2.3 牛拉法和 PQ 分解法的选用边界牛顿-拉夫逊法具备二阶收敛速度正常情况下 3 到 5 次迭代就能把残差压到 1e-6 以下代价是每次迭代都要重新计算并分解雅可比矩阵。对于几百节点的小系统这个成本可以忽略对于上万节点的大电网稀疏分解的开销会明显影响计算速度。PQ 分解法是在牛拉法基础上做的近似。输电网支路的电阻远小于电抗节点电压幅值变化对有功潮流影响很小相角差变化对无功潮流影响也很小于是可以把雅可比矩阵固定成两个常数矩阵分别在 P-θ 迭代和 Q-V 迭代里使用。这样一来迭代次数会从几次增加到十几次甚至几十次但每轮迭代不再需要重新组装矩阵和重新分解总耗时往往比牛拉法低一个量级。算法单轮成本典型迭代次数适用场景牛顿-拉夫逊法高每轮重算并分解 J3~5 次任意 r/x 比值、中小规模PQ 分解法低只做三角回代10~30 次高压输电网、大规模在线分析选择哪种算法不是一个性能判断题而是收敛性和计算资源的权衡。下面两章分别给出两种算法的最小可运行实现。3. 用 Python 从零实现牛拉法潮流计算3.1 先建导纳矩阵支路参数表与 build_y潮流程序的第一份输入不是迭代算法而是支路参数。我用一个四节点系统作为教学算例包含 1 个平衡节点、1 个 PV 节点和 2 个 PQ 节点。所有数值都用标幺值100 MVA 基准下可以直接复现。支路首端末端rp.u.xp.u.对地电纳 b/2p.u.L1120.0050.100.02L2130.0100.150.03L3240.0150.200.04L4340.0080.100.02节点数据为节点 1 是平衡节点V1.06θ0节点 2 是 PV 节点P0.50V1.04节点 3 是 PQ 节点P-0.70Q-0.30节点 4 是 PQ 节点P-0.55Q-0.20。负荷的 P 和 Q 都带负号这是约定俗成的方向注入网络为正值消耗为负值。import numpy as np def build_y(n, branches): y np.zeros((n, n), dtypecomplex) for i, j, r, x, b_half in branches: z r 1j * x # 支路阻抗 ys 1.0 / z # 支路导纳 y[i, i] ys 1j * b_half # 首端自导纳加对地电纳 y[j, j] ys 1j * b_half # 末端自导纳加对地电纳 y[i, j] - ys # 互导纳为负 y[j, i] - ys return y这段代码把一个支路表变成节点导纳矩阵。自导纳由与该节点相连的所有支路导纳之和再加上该支路两端各一半的对地电纳互导纳等于支路导纳的负数。直接使用 dense 矩阵只适合教学和小系统验证工程实现应该用 scipy.sparse 的 CSC 格式保存后续稀疏 LU 分解速度会快几个数量级。3.2 牛拉法迭代核心数值雅可比版本我一般先把收敛流程跑通再用解析雅可比替换。下面的版本用中心差分组装雅可比矩阵逻辑最直观不容易被符号约定带偏。搞清楚流程后再按 2.2 节的四个分块改成解析实现。def calc_power(y, v, th): vc v * np.exp(1j * th) sc vc * np.conj(y vc) return sc.real, sc.imag def newton_lf(y, ref, pv, pq, p_spec, q_spec, v_spec, eps1e-6, max_iter15): n y.shape[0] v np.ones(n) th np.zeros(n) v[list(ref) list(pv)] [v_spec[i] for i in list(ref) list(pv)] theta_idx [i for i in range(n) if i not in ref] # 需要求相角的节点 v_idx list(pq) # 需要求电压幅值的节点 def mismatch(theta, volt): p, q calc_power(y, volt, theta) d [] for i in theta_idx: d.append(p_spec[i] - p[i]) for i in v_idx: d.append(q_spec[i] - q[i]) return np.array(d) x0 np.concatenate([th[theta_idx], v[v_idx]]) for it in range(max_iter): th_full np.zeros(n); th_full[theta_idx] x0[:len(theta_idx)] v_full np.ones(n); pv_and_ref set(pv) | set(ref); v_full[list(ref) list(pv)] [v_spec[i] for i in list(ref) list(pv)] v_full[v_idx] x0[len(theta_idx):] r mismatch(th_full, v_full) if np.max(np.abs(r)) eps: break J np.zeros((len(x0), len(x0))) h 1e-6 for k in range(len(x0)): xp x0.copy(); xp[k] h thp np.zeros(n); thp[theta_idx] xp[:len(theta_idx)] vp np.ones(n); vp[list(ref) list(pv)] [v_spec[i] for i in list(ref) list(pv)] vp[v_idx] xp[len(theta_idx):] xm x0.copy(); xm[k] - h thm np.zeros(n); thm[theta_idx] xm[:len(theta_idx)] vm np.ones(n); vm[list(ref) list(pv)] [v_spec[i] for i in list(ref) list(pv)] vm[v_idx] xm[len(theta_idx):] J[:, k] (mismatch(thp, vp) - mismatch(thm, vm)) / (2 * h) dx np.linalg.solve(J, -r) x0 x0 dx return v_full, th_full这里把整个状态向量拆成两块所有非平衡节点的相角放在前面所有 PQ 节点的电压幅值放在后面。每次计算残差时再从这两块拼回完整的电压向量。之所以不直接对导纳矩阵求导是为了先避开雅可比符号约定等程序能够稳定收敛后再把中心差分替换成解析式。中心差分每轮迭代要多算很多次功率只适合教学验证不能直接用在大规模算例上。np.linalg.solve(J, -r)这一步是解线性方程组-r 是牛顿法的修正方向。要注意残差量级和 eps 的单位都是标幺值100 MVA 基准下 1e-6 对应 0.1 kW已经远高于工程精度需求。3.3 收敛判据和四个关键参数参数常见取值说明收敛容差 eps1e-6 p.u.指最大功率失配不指电压变化量最大迭代次数10~15牛拉法正常 3~5 次收敛初始电压平启动 1.0∠0°适用于绝大多数输电网断面初值电压幅值PV 和平衡节点按给定值其他节点取 1.0提示判断收敛时不要同时看残差和电压变化量。潮流方程是功率残差方程只看 max(|ΔP|, |ΔQ|) 是否小于 eps不要在迭代早期混入有功功率不平衡判断否则会导致迭代提前终止在非收敛点。对小系统而言平启动初值几乎不会出问题。但系统重载、存在长距离重负荷线路时牛拉法也可能发散。这时候要检查的不只是初值还包括导纳矩阵里是否漏掉了对地电纳、并联电容器等无功补偿设备。4. PQ 分解法输电网潮流计算的快速通道4.1 两个近似和两个矩阵B′、B″PQ 分解法的基础是高压输电网中支路电抗远大于电阻电压幅值和相角两个问题近似解耦。有功功率主要和相角差有关无功功率主要和电压幅值有关。于是牛拉法里耦合的 2x2 分块雅可比被拆成两个独立的常系数矩阵。第一个矩阵 B′ 用于 P-θ 迭代由支路电抗倒数组成忽略电阻和对地电纳。第二个矩阵 B″ 用于 Q-V 迭代保留了对地电纳的影响通常直接取导纳矩阵虚部的负值。二者只需要在迭代开始前组装并分解一次后续每轮迭代只做两次线性回代计算量大减。这种近似不是对所有系统都成立。配电网 r/x 比值较大P-V 和 Q-θ 的耦合不可忽略PQ 分解法很容易发散或收敛到错误解。所以该方法只建议在高压输电网、主网分析这类场景使用。对比项牛拉法PQ 分解法雅可比矩阵每轮重算、重分解常数矩阵分解一次迭代次数3~5 次10~30 次有功与无功迭代联立求解交替求解对 r/x 比值的依赖不敏感较敏感4.2 PQ 分解法迭代代码与牛拉法的行为差异def pq_lf(y, ref, pv, pq, p_spec, q_spec, v_spec, eps1e-6, max_iter50): n y.shape[0] bpp -np.imag(y) # 取导纳矩阵虚部并取负作为 Q-V 矩阵 bp np.zeros((n, n)) for i in range(n): for j in range(n): if i ! j and abs(y[i, j]) 1e-12: bp[i, j] -np.imag(1.0 / (1.0 / y[i, j])) # 用支路电抗倒数近似 # 简化做法对无变压器网络bp -1/x 组装对角与非对角 bp np.zeros((n, n)) for i in range(n): for j in range(i 1, n): if abs(y[i, j]) 1e-12: xij -1.0 / np.imag(y[i, j]) bp[i, i] 1.0 / xij bp[j, j] 1.0 / xij bp[i, j] - 1.0 / xij bp[j, i] - 1.0 / xij th np.zeros(n) v np.ones(n) v[list(ref) list(pv)] [v_spec[i] for i in list(ref) list(pv)] p_idx [i for i in range(n) if i not in ref] q_idx list(pq) for it in range(max_iter): p, q calc_power(y, v, th) dp np.array([p_spec[i] - p[i] for i in p_idx]) dpv dp / v[p_idx] dth np.linalg.solve(bp[np.ix_(p_idx, p_idx)], dpv) th[p_idx] dth p, q calc_power(y, v, th) dq np.array([q_spec[i] - q[i] for i in q_idx]) dqv dq / v[q_idx] dv np.linalg.solve(bpp[np.ix_(q_idx, q_idx)], dqv) v[q_idx] dv if np.max(np.abs(dp)) eps and np.max(np.abs(dq)) eps: break return v, th这段代码先组装 B′阵再用np.ix_取出与未知量对应的子块。每轮迭代里先固定电压幅值解相角再固定相角解电压幅值两个子问题交替推进。dp / v[p_idx]和dq / v[q_idx]是 PQ 分解法的标准右端项处理方式把功率失配量除以当前电压幅值让线性化后的方程更接近常数矩阵形式。牛拉法是一轮迭代同时解出相角和电压幅值PQ 分解法则把一轮拆成两个半步。因此 PQ 分解法迭代次数多但内部矩阵不需要重新分解整体效率在规模足够大时优势明显。调试早期发现不收敛时可以把每轮的 dp、dq 最大分量打出来观察是停留在 P 侧还是 Q 侧判断瓶颈在哪个子问题。4.3 参数经验与无功越限处理参数建议值注意事项收敛容差1e-6 p.u.与牛拉法等精度标准一致最大迭代次数30~50PQ 分解法不应设 10 次以内的上限初始电压平启动 1.0∠0°重负荷网络可先用牛拉法得初值PV 无功越限转 PQ 并固定 Q注意避免频繁类型切换注意工程计算里 PV 节点的无功越限是高频坑。发电机无功达到上限后电压幅值已经无法维持设定值此时必须把节点类型从 PV 改成 PQ并以 Q_upper 作为给定无功。如果下一轮迭代又向限值内回落不要马上切回 PV否则会出现节点类型振荡迭代始终无法收敛。小系统里 PQ 分解法和牛拉法结果应完全一致因为收敛到的是同一个物理平衡点。若两个算法给出的支路功率差异超过 0.1%多半是 B′、B″ 矩阵与导纳矩阵的单位定义不一致先检查对地电纳是否被重复计入。5. 潮流计算不收敛时按这三个方向排查5.1 先查数据再查迭代迭代发散时第一件事不是调初值和收敛因子而是回到输入数据。导纳矩阵的对角元素是否包含所有与节点相连的支路是检查频率最高的错误点。孤岛节点在导纳矩阵里表现为整行几乎为零牛顿法会在该节点出现单步修正量巨大、残差反而增大的现象。另一个常见问题是负荷符号写反。负荷的 P、Q 应为负值如果把负荷写成注入整个断面会变成一台“虚拟发电机”潮流可能收敛到电压整体升高约 1.05 的结果物理上看起来像是系统在输出大量无功。接口程序里统一负荷方向并在计算前打印全网有功、无功加和作为粗校验能拦住大部分数据错误。5.2 三个必须做的结果校验收敛只代表残差为零不代表结果一定正确。第一项校验是功率平衡平衡节点计算得到的 P 应约等于全网负荷总 P 加上网损总 P偏差应小于收敛精度对应有名值。第二项校验是支路功率对称性对每条支路从两端分别计算的功率数值方向相反且两个流向的有功功率之和应为正这个正值就是支路损耗不可能为负。第三项校验是电压合理性PQ 节点电压落在 0.95~1.05 范围内是常规输电网的合理区间与这个区间偏离过大时即使迭代收敛也要复核无功补偿配置。5.3 一次收敛后中途振荡的处理经验现场分析中还会遇到一种情况前几次迭代都很顺利到第 8 轮左右残差在 1e-4 附近不再下降。这时优先怀疑 PV 节点无功越限的切换逻辑写在了迭代循环内部导致节点类型每两轮被重写一次。常见的修复是记录节点类型切换状态切换后的节点在后续至少 5 轮迭代中强制保持 PQ 类型不做反向切换给电压幅值留出足够的收敛空间。如果仍然振荡可以在 Q-V 迭代的修正量上加入阻尼系数把v[q_idx] dv改成v[q_idx] 0.8 * dv。牛拉法不需要这个手段但 PQ 分解法在两个子问题交替时容易出现过冲小幅阻尼是从工程经验里最容易学到的技巧。本文还有配套的精品资源点击获取
