电力系统潮流计算：Newton-Raphson方法及其收敛性分析

引言

潮流计算（Load Flow / Power Flow）是电力系统分析的基础工具，用于确定电网各节点电压幅值、相角及线路功率分布。Newton-Raphson（NR）法因其快速收敛（二次收敛速度）成为工业标准算法。本文推导NR潮流方程的雅可比矩阵结构，并讨论其收敛特性。

理论基础

节点功率方程

节点 $i$ 的注入视在功率：

$$S_i = P_i + j Q_i = V_i \left( \sum_{k=1}^{N} Y_{ik} V_k \right)^*$$

其中 $Y_{ik} = G_{ik} + j B_{ik}$ 为节点导纳矩阵元素。写成极坐标 $V_i = |V_i| e^{j\theta_i}$ 形式：

$$P_i = |V_i| \sum_{k=1}^{N} |V_k| \left( G_{ik} \cos\theta_{ik} + B_{ik} \sin\theta_{ik} \right)$$

$$Q_i = |V_i| \sum_{k=1}^{N} |V_k| \left( G_{ik} \sin\theta_{ik} - B_{ik} \cos\theta_{ik} \right)$$

其中 $\theta_{ik} = \theta_i - \theta_k$。

Newton-Raphson 迭代格式

令状态向量 $x = [\theta , |V|]^T$（PV节点仅 $\theta$），失配量 $\Delta F(x) = F^{\text{spec}} - F(x)$：

$$J \Delta x = \Delta F$$

雅可比矩阵分块结构：

$$J = begin{bmatrix} frac{partial P}{partial theta} & frac{partial P}{partial |V|} \ frac{partial Q}{partial theta} & frac{partial Q}{partial |V|} end{bmatrix}$$

对角块计算示例（$i = j$）：

$$\frac{\partial P_i}{\partial \theta_i} = -|V_i| \sum_{k \neq i} |V_k| \left( G_{ik} \sin\theta_{ik} - B_{ik} \cos\theta_{ik} \right)$$

工程数据：IEEE 9节点算例

节点	类型	$P_G$ (MW)	$Q_G$ (Mvar)	$P_L$ (MW)	$Q_L$ (Mvar)
1	平衡	—	—	0	0
2	PV	163	—	0	0
3	PV	85	—	0	0
5	PQ	0	0	125	50
6	PQ	0	0	90	30
8	PQ	0	0	100	35

NR 迭代通常 3-5 步收敛至 $10^{-6}$ pu。快速解耦法（Fast Decoupled Load Flow, FDLF）利用 $P-\theta$ 与 $Q-|V|$ 的弱耦合关系近似，每次迭代仅需计算一次雅可比，计算量约为 NR 的 20%。

数值实现

以下用 Newton-Raphson 法求解一个两节点系统（平衡节点 + PQ 负荷节点）：

import numpy as np

# 线路导纳 y = 1/(0.02 + j0.1) pu，节点1为平衡节点 V1=1.0∠0°
y = 1.0 / (0.02 + 0.1j)
G, B = y.real, y.imag
P_spec, Q_spec = -1.0, -0.5   # 节点2注入功率 (pu)

def mismatch(V2, th2):
    P2 = V2*(V2*G) + V2*1.0*(G*np.cos(th2) + B*np.sin(th2))
    Q2 = V2*(-V2*B) + V2*1.0*(G*np.sin(th2) - B*np.cos(th2))
    return P_spec - P2, Q_spec - Q2

V2, th2 = 1.0, 0.0            # 平坦启动
for it in range(10):
    dP, dQ = mismatch(V2, th2)
    if max(abs(dP), abs(dQ)) < 1e-8:
        break
    # 数值雅可比（工程实现用解析式）
    eps = 1e-6
    J = np.zeros((2, 2))
    for k, (dv, dt) in enumerate([(eps, 0), (0, eps)]):
        dP2, dQ2 = mismatch(V2 + dv, th2 + dt)
        J[0, k] = (dP2 - dP) / eps * (-1)
        J[1, k] = (dQ2 - dQ) / eps * (-1)
    dV, dth = np.linalg.solve(J, [dP, dQ])
    V2, th2 = V2 + dV, th2 + dth
    print(f"iter {it}: V2={V2:.6f} pu, th2={np.degrees(th2):+.4f}°")

收敛性分析

NR 的收敛特性确保：

对足够接近解的初值（平坦启动 $V=1.0+\theta=0$），二次收敛
重载或病态系统可采用最优乘子法（Optimal Multiplier）增强全局收敛性
连续潮流（Continuation Power Flow, CPF）用于 PV 曲线追踪与电压稳定裕度评估

结论

Newton-Raphson 法是潮流计算的黄金标准，其雅可比矩阵的稀疏性（典型电网仅 2-3% 非零）使其可高效求解上万节点系统。结合稀疏线性方程求解器（如 SuperLU, KLU），是现代 EMS/SCADA 系统的数学核心。