尧图网页设计 
简介基于Python的牛顿-拉夫逊法电力系统潮流计算项目面向电力系统初学者、课程设计与毕设人群解决传统手算潮流迭代步骤繁琐的问题实现从节点数据输入到迭代求解的完整流程。资源共7个文件核心为4个Python脚本分别承担数据预处理、迭代计算、主程序运行与功能测试另附README说明文档和开源许可证压缩包仅11KB轻量易用便于快速阅读与二次开发。该资源已有211人学习浏览适合作为课程作业或入门实践的起步模板。借助代码中的节点导纳矩阵构建、雅可比矩阵修正与收敛判据等关键模块读者可直观理解牛顿-拉夫逊法的数值迭代思想并在此基础上扩展为更大规模电力网络的潮流分析工具。1. 为什么潮流计算绕不开牛顿-拉夫逊法电力系统潮流计算的本质是求一组满足节点功率平衡方程的电压幅值和相角。电网规模一大方程就是高维非线性的没法直接求解析解。业界最常见的做法是牛顿-拉夫逊法把非线性方程组在初值处做泰勒展开保留一阶项反复迭代修正。它的收敛速度是二次的通常迭代 3 到 5 次就能把偏差压到 1e-6 以下这是高斯-赛德尔法做不到的。用 Python 实现牛顿-拉夫逊法好处在于 NumPy 把矩阵运算和线性方程组求解都封装好了核心的雅可比矩阵迭代更新不到 100 行代码。相比 MATLAB 的传统教学实现Python 在数据预处理、结果可视化和与调度系统对接上都更顺手。适合的人群是电力系统方向的学生、刚接触潮流计算的工程师以及想在 Python 里复现教材算法的人。2. 从节点导纳矩阵到功率方程牛顿-拉夫逊的建模基础2.1 节点导纳矩阵怎么构建潮流计算的第一步是建立节点导纳矩阵 Y。它的规则很简单对角线元素 Yii 是连接节点 i 的所有支路导纳之和非对角线元素 Yij 是节点 i 和 j 之间支路导纳的负值。这里的导纳是复数单位是西门子 S。用 Python 构建时最直接的做法是先用零矩阵初始化然后遍历支路参数填充。实际中我会把支路数据定义成一个列表每个元素包含首端节点、末端节点、电阻、电抗和电纳代码清晰也好调试。import numpy as np def build_ybus(branch_data, n): # branch_data: [from_bus, to_bus, r, x, b]b 为对地电纳的一半 # n: 节点总数 Y np.zeros((n, n), dtypecomplex) for f, t, r, x, b in branch_data: z r 1j * x y 1 / z Y[f, f] y 1j * b Y[t, t] y 1j * b Y[f, t] - y Y[t, f] - y return Y构建导纳矩阵的时候有几个容易踩的坑。第一所有参数必须是标幺值不是有名值否则矩阵数值会差很多个数量级第二变压器支路要考虑变比非 1:1 变比时非对角元素要乘以变比系数第三线路对地电纳通常以总电纳的一半并入两端节点写数据时别搞混。Ybus 构建完成后它是整个牛顿-拉夫逊迭代里唯一不变的量每次迭代都用它计算注入功率和雅可比矩阵所以建议用函数封装起来避免重复计算。2.2 节点分类和功率方程在潮流计算里节点分三类。PQ 节点已知有功 P 和无功 Q待求电压幅值 V 和相角 θ普通负荷节点大多是这一类PV 节点已知有功 P 和电压幅值 V待求无功 Q 和相角 θ发电机组一般按 PV 节点处理平衡节点电压幅值 V 和相角 θ 固定通常取 1.0∠0°待求全网的功率平衡一般选容量最大的发电厂节点。每个节点的功率方程有两个有功和无功分别写成Pi Vi * Σ(Yij * Vj * cos(θi - θj - φij)) Qi Vi * Σ(Yij * Vj * sin(θi - θj - φij))其中 φij 是导纳矩阵元素 Yij 的相角。这两个方程是所有后续计算的起点牛顿-拉夫逊的核心就是让每个节点的 P 和 Q 计算值逐步逼近给定值。在代码实现层面更常见的做法是用复数形式计算注入功率NumPy 的向量化能省掉一层循环def calc_power_injection(V_mag, V_ang, Ybus): V V_mag * np.exp(1j * V_ang) S V * np.conj(Ybus V) return S.real, S.imag # 返回有功和无功向量这种写法简洁但初学者容易忽略np.conj这一步少了共轭算出来的功率方向是反的迭代大概率发散。2.3 直角坐标和极坐标选哪种牛顿-拉夫逊法有两种坐标形式。直角坐标下电压表示为 e jf方程是二次的雅可比矩阵结构固定迭代计算量稍小极坐标下电压表示为 V∠θ方程带三角函数雅可比矩阵维度比直角坐标少一半因为 PV 节点的 Q 方程不参与迭代。教材里最常见的推导用的是极坐标原因是物理意义直观电压幅值和相角可以直接和运行约束挂钩。工程实现上极坐标雅可比矩阵的四个分块H、N、J、L各有明确的物理含义调试时也好定位问题。所以我后面的实现统一用极坐标形式。3. 核心迭代雅可比矩阵的构造和求解3.1 有功和无功的不平衡量牛顿法的本质是对不平衡量做线性化修正。PQ 节点的有功和无功、PV 节点的有功无功不参与都会和给定值产生偏差这个偏差量就是迭代的驱动力。对于节点 iΔPi P_spec_i - P_calc_i ΔQi Q_spec_i - Q_calc_i当所有 ΔP 和 ΔQ 的绝对值都小于收敛阈值比如 1e-6时潮流计算就算收敛了。这一步计算量不大但要注意给定值 P_spec 里要扣除节点上并联补偿电容或电抗器的功率。3.2 雅可比矩阵各分块的公式极坐标下将功率方程分别对 θ 和 V 求偏导得到 4 个分块矩阵。i ≠ j 时公式相对简单i j 时需要在求和项基础上加上自身节点的项。具体来说H 是 ΔP 对 θ 的偏导N 是 ΔP 对 V 的偏导J 是 ΔQ 对 θ 的偏导L 是 ΔQ 对 V 的偏导。i ≠ j: Hij Vi * Vj * (Gij * sin(θi - θj) - Bij * cos(θi - θj)) Nij Vi * (Gij * cos(θi - θj) Bij * sin(θi - θj)) Jij -Vi * Vj * (Gij * cos(θi - θj) Bij * sin(θi - θj)) Lij Vi * (Gij * sin(θi - θj) - Bij * cos(θi - θj))这里的 Gij 和 Bij 分别是 Ybus 元素的实部和虚部。对角元素的公式比非对角项多出自己的注入功率项实现时最容易出错的地方就在这里。3.3 代码实现迭代主循环把公式落成 Python 代码时我会做一个完整的迭代函数。代码里要注意平衡节点不参与修正方程PV 节点的无功方程不参与修正所以需要把对应的行列剔除。def newton_raphson_power_flow(Ybus, S_spec, V_init, pv_nodes, pq_nodes, tol1e-6, max_iter10): n len(V_init) V_mag np.abs(V_init) V_ang np.angle(V_init) # 参与迭代的节点PQ PV不含平衡节点 pvpq pv_nodes pq_nodes pqpq pq_nodes # 只有 PQ 节点参与无功修正 # 雅可比矩阵维度 n_pvpq len(pvpq) n_pq len(pqpq) # 给定功率标幺值 P_spec S_spec.real Q_spec S_spec.imag for iteration in range(max_iter): # 计算当前电压下的注入功率 V V_mag * np.exp(1j * V_ang) S_calc V * np.conj(Ybus V) P_calc S_calc.real Q_calc S_calc.imag # 计算不平衡量 dP P_spec[pvpq] - P_calc[pvpq] dQ Q_spec[pqpq] - Q_calc[pqpq] # 检查收敛 if np.max(np.abs(np.r_[dP, dQ])) tol: break # 构造雅可比矩阵详细函数见下文 J build_jacobian(Ybus, V_mag, V_ang, pvpq, pqpq, P_calc, Q_calc) # 求解修正方程 dX np.linalg.solve(J, np.r_[dP, dQ]) # 分离角度和幅值的修正量并更新 dV_ang dX[:n_pvpq] dV_mag dX[n_pvpq:n_pvpq n_pq] V_ang[pvpq] dV_ang V_mag[pqpq] dV_mag return V_mag, V_ang, iteration 1np.linalg.solve是核心求解函数它内部调用 LAPACK 的 dgesv 例程做 LU 分解对小规模电网几百个节点以内效率足够。值得注意的是np.r_是将两个向量拼接成一个方便统一求解。3.4 雅可比矩阵构造的完整实现上面主循环里调用的build_jacobian是牛顿-拉夫逊法的关键部分值得单独拆开看。它要按 H、N、J、L 四个分块拼接成一个大矩阵然后删掉平衡节点对应行列。def build_jacobian(Ybus, V_mag, V_ang, pvpq, pqpq, P_calc, Q_calc): # 极坐标下的雅可比矩阵构建 n len(V_mag) G Ybus.real B Ybus.imag # 预分配内存 H np.zeros_like(Ybus, dtypefloat) N np.zeros_like(Ybus, dtypefloat) J np.zeros_like(Ybus, dtypefloat) L np.zeros_like(Ybus, dtypefloat) # 计算角度差矩阵 d_ang V_ang[:, None] - V_ang[None, :] # 非对角项 H V_mag[:, None] * V_mag[None, :] * (G * np.sin(d_ang) - B * np.cos(d_ang)) N V_mag[:, None] * (G * np.cos(d_ang) B * np.sin(d_ang)) J -V_mag[:, None] * V_mag[None, :] * (G * np.cos(d_ang) B * np.sin(d_ang)) L V_mag[:, None] * (G * np.sin(d_ang) - B * np.cos(d_ang)) # 对角项修正 for i in range(n): # 矩阵对角线需要叠加 Vi^2 * (Gi jBi) 项 H[i, i] Q_calc[i] - V_mag[i]**2 * B[i, i] N[i, i] P_calc[i] / V_mag[i] V_mag[i] * G[i, i] J[i, i] -P_calc[i] V_mag[i]**2 * G[i, i] L[i, i] Q_calc[i] / V_mag[i] - V_mag[i] * B[i, i] # 拼接雅可比矩阵第一行行是 [H, N]第二行是 [J, L] J_full np.vstack([ np.hstack([H, N]), np.hstack([J, L]) ]) # 剔除不需要的行列 pvpq_idx [i for i in pvpq] pqpq_idx [i for i in pqpq] keep pvpq_idx [n i for i in pqpq_idx] # 同时剔除平衡节点对应行和列 bal_nodes [i for i in range(n) if i not in pvpq and i not in pqpq] drop bal_nodes [n i for i in bal_nodes] keep [i for i in keep if i not in drop] return J_full[np.ix_(keep, keep)]这段代码里有一个关键的细节雅可比矩阵表达式的下标关系。很多教材用 i 表示节点行j 表示节点列但组装成大矩阵时$\Delta P$ 部分在上、$\Delta Q$ 部分在下对应的列分别是 $\Delta \theta$ 和 $\Delta V$顺序不能颠倒否则求解出来的修正量会和变量对不上号。3.5 和不动点迭代的差异初学者经常把牛顿-拉夫逊和简单迭代混为一谈。不动点迭代是直接代入计算期望值比如 $X^{(k1)} f(X^{(k)})$形式简单但是收敛性差很多时候线性的收敛速度意味着要几十次甚至上百次迭代。牛顿-拉夫逊每轮迭代都要重新计算雅可比矩阵并解一次线性方程组单轮计算量大但收敛极快。在高电压输电网络中牛顿法的收敛半径也比较大初值取平启动所有 PQ 节点 V1.0, θ0通常就能收敛。相比之下快速分解法PQ 分解法虽然单轮更快但适用范围有限遇到高 R/X 比的配电网会发散。所以做输电网潮流牛顿-拉夫逊还是首选。4. 用 Python 跑通一个完整算例从数据到结果4.1 算例数据和预处理我选一个经典的三节点两电源系统来跑通流程节点 1 是平衡节点节点 2 是 PV 节点节点 3 是 PQ 节点。支路参数和节点注入功率用标幺值表示基准容量取 100 MVA。# 支路数据: [首端, 末端, 电阻(pu), 电抗(pu), 对地电纳一半(pu)] branch_data [ [0, 1, 0.01, 0.05, 0.02], [1, 2, 0.015, 0.06, 0.02], [0, 2, 0.012, 0.055, 0.015] ] # 节点功率标幺值: P jQ注入为正负荷为负 S_spec np.array([ complex(0.5, 0.1), # 节点1: 平衡节点发电机 complex(0.3, 0.2), # 节点2: PV节点发电机 complex(-0.4, -0.15) # 节点3: PQ节点负荷 ], dtypecomplex) pv_nodes [1] pq_nodes [2] V_init np.array([1.0, 1.0, 1.0], dtypefloat) # 电压幅值初值 Ybus build_ybus(branch_data, 3) V_mag, V_ang, iters newton_raphson_power_flow( Ybus, S_spec, V_init, pv_nodes, pq_nodes ) print(迭代次数:, iters) for i in range(3): print(f节点{i1}: V{V_mag[i]:.4f}, θ{np.degrees(V_ang[i]):.2f}°)运行这个脚本迭代收敛后输出结果。可以看到平衡节点的注入功率会和我们给定的 0.50.1j 不同这是正常的因为平衡节点的功率本来就是待定量用于补偿全网的损耗和功率不平衡。4.2 迭代过程可视化观察收敛速度牛顿-拉夫逊的优势在于收敛速度快这里我加一段代码记录每轮的最大不平衡量观察它是怎么下降的。常见的表现是第一轮残差在 0.1 量级第二轮降到 1e-3第三轮到 1e-6呈平方收敛趋势。如果在你的算例里残差呈线性缓慢下降多半是雅可比矩阵构造有误或迭代步长不对。def newton_raphson_with_trace(Ybus, S_spec, V_init, pv_nodes, pq_nodes, tol1e-6, max_iter10): # 和标准版本相同但记录每轮的残差 ... residuals [] for iteration in range(max_iter): ... residual np.max(np.abs(np.r_[dP, dQ])) residuals.append(residual) ... return V_mag, V_ang, iteration 1, residuals然后把residuals画成半对数坐标图Y 轴取 log。平方收敛的特征非常明显图线在最初的几轮里每个迭代步落差都增大而不是均匀下降。4.3 参数设置的经验值和调试建议几个关键参数的选择直接决定迭代是否收敛。收敛阈值 tol 一般取 1e-6对大多数输电网分析足够再严也没有太大实际意义因为模型本身的误差远大于这个量级。最大迭代次数在常规电网取 10 次基本够用但要小心病态网络重负荷和电压偏低时可能需要 15 到 20 次。初值 V_init 的作用常被忽视。平启动时所有节点取 1.0∠0° 是一种标准做法但在某些重负荷条件下直接把初值改成 0.95 反而更容易收敛。如果遇到不收敛的情况我一般的调试顺序是先查 Ybus 是否正确然后查功率方向约定是否有误再检查雅可比矩阵对角项是否漏了项最后才是调初值。4.4 常见报错和处理方法迭代发散时最常碰到的报错是np.linalg.solve抛出LinAlgError比如矩阵奇异或近奇异。这说明雅可比矩阵中出现了线性相关的行或列常见原因包括PV 节点电压幅值设置得太低导致无功越限或者导纳矩阵数值差异过大有的元素是 10^-4有的是 10^2造成矩阵条件数过差。解决方法是先打印雅可比矩阵的条件数——np.linalg.cond(J)——如果大于 1e12就要做数值尺度归一化。还有一种情况是虽然没报错但残差迭代 50 次也不收敛。这时候多半是负荷节点的无功给定值超出了实际能力在物理上就不存在解。比如重负荷下无功需求 1.0 pu而该节点电压只有 0.7此时应当减小负荷或增加无功补偿而不是硬调算法参数。5. 提升收敛鲁棒性的两个实用手段5.1 最优乘子法下垂系统在重负荷场景下牛顿-拉夫逊法可能出现迭代过程中的逐次振荡也就是残差大小反复跳动不下降。一个有效的解决方法是采用最优乘子 $\mu$把每一步的修正量乘以一个系数再进行更新x(k1) x(k) μ * Δx$\mu$ 的选取方式是让下一步的残差最小化可以用一维搜索或二次插值近似。在 Python 里的实现很简单就是在外层循环再加一个内层搜索def line_search_update(J, dX, V_mag, V_ang, Ybus, S_spec, pvpq, pqpq): # 用二次插值近似最优步长 mu 1.0 best_residual calc_residual(V_mag mu * dV_mag, V_ang mu * dV_ang) for _ in range(5): mu_new mu * 0.5 residual_new calc_residual(...) if residual_new best_residual: best_residual residual_new mu mu_new return mu加了最优乘子后收敛域明显扩大代价是每轮迭代多算几次功率对中小型网络完全可接受。5.2 稀疏矩阵和大型系统扩展上面实现的雅可比矩阵是稠密的节点数少时没有问题但一旦超过 1000 个节点稠密 LU 分解的时间和内存都会呈立方增长。此时应该切换到 SciPy 的稀疏矩阵工具核心改动只有两处用scipy.sparse构建 Ybus 和雅可比矩阵用splu或bicgstab求解修正方程。from scipy.sparse import csr_matrix, csc_matrix from scipy.sparse.linalg import splu # Ybus 构建时改用 lil_matrix 逐步填充 # 最后转成 csr_matrix 格式加速求解 Y_sparse csr_matrix(Y) ... J_sparse csc_matrix(J_full) # LU分解要求CSC格式 lu splu(J_sparse) dX lu.solve(np.r_[dP, dQ])预处理技术比如 LU 分解的节点重排序能显著降低填充量SCIPY 内置的splu已经自动做了列排序不需要手动干预。实际工程中2000 节点以下的电网用稀疏牛顿法在个人电脑上几秒内就能完成一次潮流计算这也是很多开源库如 pandapower 的实现思路。最后建议在项目里保留一个残差追踪函数把每次迭代的 $|\Delta P|\infty$ 和 $|\Delta Q|\infty$ 写进日志。一旦出现不收敛日志能直接指出是哪一个节点、哪一个物理量的残差消不下去这比盯着叠代次数空想要快得多。本文还有配套的精品资源点击获取 用Python实现牛顿-拉夫逊法潮流计算：原理、代码与实战
简介基于Python的牛顿-拉夫逊法电力系统潮流计算项目面向电力系统初学者、课程设计与毕设人群解决传统手算潮流迭代步骤繁琐的问题实现从节点数据输入到迭代求解的完整流程。资源共7个文件核心为4个Python脚本分别承担数据预处理、迭代计算、主程序运行与功能测试另附README说明文档和开源许可证压缩包仅11KB轻量易用便于快速阅读与二次开发。该资源已有211人学习浏览适合作为课程作业或入门实践的起步模板。借助代码中的节点导纳矩阵构建、雅可比矩阵修正与收敛判据等关键模块读者可直观理解牛顿-拉夫逊法的数值迭代思想并在此基础上扩展为更大规模电力网络的潮流分析工具。1. 为什么潮流计算绕不开牛顿-拉夫逊法电力系统潮流计算的本质是求一组满足节点功率平衡方程的电压幅值和相角。电网规模一大方程就是高维非线性的没法直接求解析解。业界最常见的做法是牛顿-拉夫逊法把非线性方程组在初值处做泰勒展开保留一阶项反复迭代修正。它的收敛速度是二次的通常迭代 3 到 5 次就能把偏差压到 1e-6 以下这是高斯-赛德尔法做不到的。用 Python 实现牛顿-拉夫逊法好处在于 NumPy 把矩阵运算和线性方程组求解都封装好了核心的雅可比矩阵迭代更新不到 100 行代码。相比 MATLAB 的传统教学实现Python 在数据预处理、结果可视化和与调度系统对接上都更顺手。适合的人群是电力系统方向的学生、刚接触潮流计算的工程师以及想在 Python 里复现教材算法的人。2. 从节点导纳矩阵到功率方程牛顿-拉夫逊的建模基础2.1 节点导纳矩阵怎么构建潮流计算的第一步是建立节点导纳矩阵 Y。它的规则很简单对角线元素 Yii 是连接节点 i 的所有支路导纳之和非对角线元素 Yij 是节点 i 和 j 之间支路导纳的负值。这里的导纳是复数单位是西门子 S。用 Python 构建时最直接的做法是先用零矩阵初始化然后遍历支路参数填充。实际中我会把支路数据定义成一个列表每个元素包含首端节点、末端节点、电阻、电抗和电纳代码清晰也好调试。import numpy as np def build_ybus(branch_data, n): # branch_data: [from_bus, to_bus, r, x, b]b 为对地电纳的一半 # n: 节点总数 Y np.zeros((n, n), dtypecomplex) for f, t, r, x, b in branch_data: z r 1j * x y 1 / z Y[f, f] y 1j * b Y[t, t] y 1j * b Y[f, t] - y Y[t, f] - y return Y构建导纳矩阵的时候有几个容易踩的坑。第一所有参数必须是标幺值不是有名值否则矩阵数值会差很多个数量级第二变压器支路要考虑变比非 1:1 变比时非对角元素要乘以变比系数第三线路对地电纳通常以总电纳的一半并入两端节点写数据时别搞混。Ybus 构建完成后它是整个牛顿-拉夫逊迭代里唯一不变的量每次迭代都用它计算注入功率和雅可比矩阵所以建议用函数封装起来避免重复计算。2.2 节点分类和功率方程在潮流计算里节点分三类。PQ 节点已知有功 P 和无功 Q待求电压幅值 V 和相角 θ普通负荷节点大多是这一类PV 节点已知有功 P 和电压幅值 V待求无功 Q 和相角 θ发电机组一般按 PV 节点处理平衡节点电压幅值 V 和相角 θ 固定通常取 1.0∠0°待求全网的功率平衡一般选容量最大的发电厂节点。每个节点的功率方程有两个有功和无功分别写成Pi Vi * Σ(Yij * Vj * cos(θi - θj - φij)) Qi Vi * Σ(Yij * Vj * sin(θi - θj - φij))其中 φij 是导纳矩阵元素 Yij 的相角。这两个方程是所有后续计算的起点牛顿-拉夫逊的核心就是让每个节点的 P 和 Q 计算值逐步逼近给定值。在代码实现层面更常见的做法是用复数形式计算注入功率NumPy 的向量化能省掉一层循环def calc_power_injection(V_mag, V_ang, Ybus): V V_mag * np.exp(1j * V_ang) S V * np.conj(Ybus V) return S.real, S.imag # 返回有功和无功向量这种写法简洁但初学者容易忽略np.conj这一步少了共轭算出来的功率方向是反的迭代大概率发散。2.3 直角坐标和极坐标选哪种牛顿-拉夫逊法有两种坐标形式。直角坐标下电压表示为 e jf方程是二次的雅可比矩阵结构固定迭代计算量稍小极坐标下电压表示为 V∠θ方程带三角函数雅可比矩阵维度比直角坐标少一半因为 PV 节点的 Q 方程不参与迭代。教材里最常见的推导用的是极坐标原因是物理意义直观电压幅值和相角可以直接和运行约束挂钩。工程实现上极坐标雅可比矩阵的四个分块H、N、J、L各有明确的物理含义调试时也好定位问题。所以我后面的实现统一用极坐标形式。3. 核心迭代雅可比矩阵的构造和求解3.1 有功和无功的不平衡量牛顿法的本质是对不平衡量做线性化修正。PQ 节点的有功和无功、PV 节点的有功无功不参与都会和给定值产生偏差这个偏差量就是迭代的驱动力。对于节点 iΔPi P_spec_i - P_calc_i ΔQi Q_spec_i - Q_calc_i当所有 ΔP 和 ΔQ 的绝对值都小于收敛阈值比如 1e-6时潮流计算就算收敛了。这一步计算量不大但要注意给定值 P_spec 里要扣除节点上并联补偿电容或电抗器的功率。3.2 雅可比矩阵各分块的公式极坐标下将功率方程分别对 θ 和 V 求偏导得到 4 个分块矩阵。i ≠ j 时公式相对简单i j 时需要在求和项基础上加上自身节点的项。具体来说H 是 ΔP 对 θ 的偏导N 是 ΔP 对 V 的偏导J 是 ΔQ 对 θ 的偏导L 是 ΔQ 对 V 的偏导。i ≠ j: Hij Vi * Vj * (Gij * sin(θi - θj) - Bij * cos(θi - θj)) Nij Vi * (Gij * cos(θi - θj) Bij * sin(θi - θj)) Jij -Vi * Vj * (Gij * cos(θi - θj) Bij * sin(θi - θj)) Lij Vi * (Gij * sin(θi - θj) - Bij * cos(θi - θj))这里的 Gij 和 Bij 分别是 Ybus 元素的实部和虚部。对角元素的公式比非对角项多出自己的注入功率项实现时最容易出错的地方就在这里。3.3 代码实现迭代主循环把公式落成 Python 代码时我会做一个完整的迭代函数。代码里要注意平衡节点不参与修正方程PV 节点的无功方程不参与修正所以需要把对应的行列剔除。def newton_raphson_power_flow(Ybus, S_spec, V_init, pv_nodes, pq_nodes, tol1e-6, max_iter10): n len(V_init) V_mag np.abs(V_init) V_ang np.angle(V_init) # 参与迭代的节点PQ PV不含平衡节点 pvpq pv_nodes pq_nodes pqpq pq_nodes # 只有 PQ 节点参与无功修正 # 雅可比矩阵维度 n_pvpq len(pvpq) n_pq len(pqpq) # 给定功率标幺值 P_spec S_spec.real Q_spec S_spec.imag for iteration in range(max_iter): # 计算当前电压下的注入功率 V V_mag * np.exp(1j * V_ang) S_calc V * np.conj(Ybus V) P_calc S_calc.real Q_calc S_calc.imag # 计算不平衡量 dP P_spec[pvpq] - P_calc[pvpq] dQ Q_spec[pqpq] - Q_calc[pqpq] # 检查收敛 if np.max(np.abs(np.r_[dP, dQ])) tol: break # 构造雅可比矩阵详细函数见下文 J build_jacobian(Ybus, V_mag, V_ang, pvpq, pqpq, P_calc, Q_calc) # 求解修正方程 dX np.linalg.solve(J, np.r_[dP, dQ]) # 分离角度和幅值的修正量并更新 dV_ang dX[:n_pvpq] dV_mag dX[n_pvpq:n_pvpq n_pq] V_ang[pvpq] dV_ang V_mag[pqpq] dV_mag return V_mag, V_ang, iteration 1np.linalg.solve是核心求解函数它内部调用 LAPACK 的 dgesv 例程做 LU 分解对小规模电网几百个节点以内效率足够。值得注意的是np.r_是将两个向量拼接成一个方便统一求解。3.4 雅可比矩阵构造的完整实现上面主循环里调用的build_jacobian是牛顿-拉夫逊法的关键部分值得单独拆开看。它要按 H、N、J、L 四个分块拼接成一个大矩阵然后删掉平衡节点对应行列。def build_jacobian(Ybus, V_mag, V_ang, pvpq, pqpq, P_calc, Q_calc): # 极坐标下的雅可比矩阵构建 n len(V_mag) G Ybus.real B Ybus.imag # 预分配内存 H np.zeros_like(Ybus, dtypefloat) N np.zeros_like(Ybus, dtypefloat) J np.zeros_like(Ybus, dtypefloat) L np.zeros_like(Ybus, dtypefloat) # 计算角度差矩阵 d_ang V_ang[:, None] - V_ang[None, :] # 非对角项 H V_mag[:, None] * V_mag[None, :] * (G * np.sin(d_ang) - B * np.cos(d_ang)) N V_mag[:, None] * (G * np.cos(d_ang) B * np.sin(d_ang)) J -V_mag[:, None] * V_mag[None, :] * (G * np.cos(d_ang) B * np.sin(d_ang)) L V_mag[:, None] * (G * np.sin(d_ang) - B * np.cos(d_ang)) # 对角项修正 for i in range(n): # 矩阵对角线需要叠加 Vi^2 * (Gi jBi) 项 H[i, i] Q_calc[i] - V_mag[i]**2 * B[i, i] N[i, i] P_calc[i] / V_mag[i] V_mag[i] * G[i, i] J[i, i] -P_calc[i] V_mag[i]**2 * G[i, i] L[i, i] Q_calc[i] / V_mag[i] - V_mag[i] * B[i, i] # 拼接雅可比矩阵第一行行是 [H, N]第二行是 [J, L] J_full np.vstack([ np.hstack([H, N]), np.hstack([J, L]) ]) # 剔除不需要的行列 pvpq_idx [i for i in pvpq] pqpq_idx [i for i in pqpq] keep pvpq_idx [n i for i in pqpq_idx] # 同时剔除平衡节点对应行和列 bal_nodes [i for i in range(n) if i not in pvpq and i not in pqpq] drop bal_nodes [n i for i in bal_nodes] keep [i for i in keep if i not in drop] return J_full[np.ix_(keep, keep)]这段代码里有一个关键的细节雅可比矩阵表达式的下标关系。很多教材用 i 表示节点行j 表示节点列但组装成大矩阵时$\Delta P$ 部分在上、$\Delta Q$ 部分在下对应的列分别是 $\Delta \theta$ 和 $\Delta V$顺序不能颠倒否则求解出来的修正量会和变量对不上号。3.5 和不动点迭代的差异初学者经常把牛顿-拉夫逊和简单迭代混为一谈。不动点迭代是直接代入计算期望值比如 $X^{(k1)} f(X^{(k)})$形式简单但是收敛性差很多时候线性的收敛速度意味着要几十次甚至上百次迭代。牛顿-拉夫逊每轮迭代都要重新计算雅可比矩阵并解一次线性方程组单轮计算量大但收敛极快。在高电压输电网络中牛顿法的收敛半径也比较大初值取平启动所有 PQ 节点 V1.0, θ0通常就能收敛。相比之下快速分解法PQ 分解法虽然单轮更快但适用范围有限遇到高 R/X 比的配电网会发散。所以做输电网潮流牛顿-拉夫逊还是首选。4. 用 Python 跑通一个完整算例从数据到结果4.1 算例数据和预处理我选一个经典的三节点两电源系统来跑通流程节点 1 是平衡节点节点 2 是 PV 节点节点 3 是 PQ 节点。支路参数和节点注入功率用标幺值表示基准容量取 100 MVA。# 支路数据: [首端, 末端, 电阻(pu), 电抗(pu), 对地电纳一半(pu)] branch_data [ [0, 1, 0.01, 0.05, 0.02], [1, 2, 0.015, 0.06, 0.02], [0, 2, 0.012, 0.055, 0.015] ] # 节点功率标幺值: P jQ注入为正负荷为负 S_spec np.array([ complex(0.5, 0.1), # 节点1: 平衡节点发电机 complex(0.3, 0.2), # 节点2: PV节点发电机 complex(-0.4, -0.15) # 节点3: PQ节点负荷 ], dtypecomplex) pv_nodes [1] pq_nodes [2] V_init np.array([1.0, 1.0, 1.0], dtypefloat) # 电压幅值初值 Ybus build_ybus(branch_data, 3) V_mag, V_ang, iters newton_raphson_power_flow( Ybus, S_spec, V_init, pv_nodes, pq_nodes ) print(迭代次数:, iters) for i in range(3): print(f节点{i1}: V{V_mag[i]:.4f}, θ{np.degrees(V_ang[i]):.2f}°)运行这个脚本迭代收敛后输出结果。可以看到平衡节点的注入功率会和我们给定的 0.50.1j 不同这是正常的因为平衡节点的功率本来就是待定量用于补偿全网的损耗和功率不平衡。4.2 迭代过程可视化观察收敛速度牛顿-拉夫逊的优势在于收敛速度快这里我加一段代码记录每轮的最大不平衡量观察它是怎么下降的。常见的表现是第一轮残差在 0.1 量级第二轮降到 1e-3第三轮到 1e-6呈平方收敛趋势。如果在你的算例里残差呈线性缓慢下降多半是雅可比矩阵构造有误或迭代步长不对。def newton_raphson_with_trace(Ybus, S_spec, V_init, pv_nodes, pq_nodes, tol1e-6, max_iter10): # 和标准版本相同但记录每轮的残差 ... residuals [] for iteration in range(max_iter): ... residual np.max(np.abs(np.r_[dP, dQ])) residuals.append(residual) ... return V_mag, V_ang, iteration 1, residuals然后把residuals画成半对数坐标图Y 轴取 log。平方收敛的特征非常明显图线在最初的几轮里每个迭代步落差都增大而不是均匀下降。4.3 参数设置的经验值和调试建议几个关键参数的选择直接决定迭代是否收敛。收敛阈值 tol 一般取 1e-6对大多数输电网分析足够再严也没有太大实际意义因为模型本身的误差远大于这个量级。最大迭代次数在常规电网取 10 次基本够用但要小心病态网络重负荷和电压偏低时可能需要 15 到 20 次。初值 V_init 的作用常被忽视。平启动时所有节点取 1.0∠0° 是一种标准做法但在某些重负荷条件下直接把初值改成 0.95 反而更容易收敛。如果遇到不收敛的情况我一般的调试顺序是先查 Ybus 是否正确然后查功率方向约定是否有误再检查雅可比矩阵对角项是否漏了项最后才是调初值。4.4 常见报错和处理方法迭代发散时最常碰到的报错是np.linalg.solve抛出LinAlgError比如矩阵奇异或近奇异。这说明雅可比矩阵中出现了线性相关的行或列常见原因包括PV 节点电压幅值设置得太低导致无功越限或者导纳矩阵数值差异过大有的元素是 10^-4有的是 10^2造成矩阵条件数过差。解决方法是先打印雅可比矩阵的条件数——np.linalg.cond(J)——如果大于 1e12就要做数值尺度归一化。还有一种情况是虽然没报错但残差迭代 50 次也不收敛。这时候多半是负荷节点的无功给定值超出了实际能力在物理上就不存在解。比如重负荷下无功需求 1.0 pu而该节点电压只有 0.7此时应当减小负荷或增加无功补偿而不是硬调算法参数。5. 提升收敛鲁棒性的两个实用手段5.1 最优乘子法下垂系统在重负荷场景下牛顿-拉夫逊法可能出现迭代过程中的逐次振荡也就是残差大小反复跳动不下降。一个有效的解决方法是采用最优乘子 $\mu$把每一步的修正量乘以一个系数再进行更新x(k1) x(k) μ * Δx$\mu$ 的选取方式是让下一步的残差最小化可以用一维搜索或二次插值近似。在 Python 里的实现很简单就是在外层循环再加一个内层搜索def line_search_update(J, dX, V_mag, V_ang, Ybus, S_spec, pvpq, pqpq): # 用二次插值近似最优步长 mu 1.0 best_residual calc_residual(V_mag mu * dV_mag, V_ang mu * dV_ang) for _ in range(5): mu_new mu * 0.5 residual_new calc_residual(...) if residual_new best_residual: best_residual residual_new mu mu_new return mu加了最优乘子后收敛域明显扩大代价是每轮迭代多算几次功率对中小型网络完全可接受。5.2 稀疏矩阵和大型系统扩展上面实现的雅可比矩阵是稠密的节点数少时没有问题但一旦超过 1000 个节点稠密 LU 分解的时间和内存都会呈立方增长。此时应该切换到 SciPy 的稀疏矩阵工具核心改动只有两处用scipy.sparse构建 Ybus 和雅可比矩阵用splu或bicgstab求解修正方程。from scipy.sparse import csr_matrix, csc_matrix from scipy.sparse.linalg import splu # Ybus 构建时改用 lil_matrix 逐步填充 # 最后转成 csr_matrix 格式加速求解 Y_sparse csr_matrix(Y) ... J_sparse csc_matrix(J_full) # LU分解要求CSC格式 lu splu(J_sparse) dX lu.solve(np.r_[dP, dQ])预处理技术比如 LU 分解的节点重排序能显著降低填充量SCIPY 内置的splu已经自动做了列排序不需要手动干预。实际工程中2000 节点以下的电网用稀疏牛顿法在个人电脑上几秒内就能完成一次潮流计算这也是很多开源库如 pandapower 的实现思路。最后建议在项目里保留一个残差追踪函数把每次迭代的 $|\Delta P|\infty$ 和 $|\Delta Q|\infty$ 写进日志。一旦出现不收敛日志能直接指出是哪一个节点、哪一个物理量的残差消不下去这比盯着叠代次数空想要快得多。本文还有配套的精品资源点击获取