开断潮流：从牛顿拉夫逊到支路潮流修正

开断潮流基于牛顿拉夫逊潮流计算结果引入灵敏度矩阵和雅可比矩阵计算支路功率对故障点注入功率的灵敏度进而计算故障后所有支路潮流的修正量在电力系统分析领域潮流计算是一项基础且关键的任务。牛顿拉夫逊潮流计算是常用方法之一而在此基础上进行扩展能进一步挖掘系统故障后潮流变化的奥秘比如引入灵敏度矩阵和雅可比矩阵计算支路功率对故障点注入功率的灵敏度以此算出故障后所有支路潮流的修正量可谓是开启了潮流分析的新潮流。牛顿拉夫逊潮流计算基础牛顿拉夫逊法是基于非线性方程的迭代求解方法。在潮流计算中电力系统的潮流方程是非线性的一般形式可表示为\[ \begin{cases}Pi Vi \sum{j1}^{n} Vj (G{ij} \cos \theta{ij} B{ij} \sin \theta{ij}) \\Qi Vi \sum{j1}^{n} Vj (G{ij} \sin \theta{ij} - B{ij} \cos \theta{ij})\end{cases} \]其中 \(Pi\) 和 \(Qi\) 分别是节点 \(i\) 的有功和无功功率注入\(Vi\) 和 \(Vj\) 是节点电压幅值\(G{ij}\) 和 \(B{ij}\) 是节点导纳矩阵元素\(\theta_{ij}\) 是节点 \(i\) 和 \(j\) 电压的相角差。在Python中实现简单的牛顿拉夫逊潮流计算核心代码片段如下为简化只展示关键迭代部分# 假设已有节点导纳矩阵Ybus初始电压幅值V和相角theta max_iter 100 tol 1e-6 for k in range(max_iter): P [] Q [] for i in range(len(V)): Pi 0 Qi 0 for j in range(len(V)): Pi V[i] * V[j] * (Ybus[i, j].real * np.cos(theta[i] - theta[j]) Ybus[i, j].imag * np.sin(theta[i] - theta[j])) Qi V[i] * V[j] * (Ybus[i, j].real * np.sin(theta[i] - theta[j]) - Ybus[i, j].imag * np.cos(theta[i] - theta[j])) P.append(Pi) Q.append(Qi) # 计算不平衡功率 dP [P[i] - P_spec[i] for i in range(len(P))] dQ [Q[i] - Q_spec[i] for i in range(len(Q))] # 检查收敛 if np.linalg.norm(dP) tol and np.linalg.norm(dQ) tol: break # 计算雅可比矩阵 # 此处省略复杂的雅可比矩阵计算代码 # 更新电压幅值和相角 # 此处省略更新代码这段代码首先在循环中计算每个节点的有功和无功功率然后算出功率不平衡量。当不平衡量小于设定的容差时认为计算收敛。引入灵敏度矩阵和雅可比矩阵在得到牛顿拉夫逊潮流计算结果后我们开始引入灵敏度矩阵 \(S\) 和雅可比矩阵 \(J\)。灵敏度矩阵 \(S\) 描述了支路功率 \(P{branch}\) 对故障点注入功率 \(P{fault}\) 的变化率即 \(S\frac{\partial P{branch}}{\partial P{fault}}\)。开断潮流基于牛顿拉夫逊潮流计算结果引入灵敏度矩阵和雅可比矩阵计算支路功率对故障点注入功率的灵敏度进而计算故障后所有支路潮流的修正量雅可比矩阵在潮流计算中用于线性化非线性的潮流方程在我们扩展的计算中同样重要。以极坐标形式的潮流方程为例雅可比矩阵的元素有如下形式部分元素\[ J{1,1} \frac{\partial Pi}{\partial \thetaj} -Vi Vj (G{ij} \sin \theta{ij} - B{ij} \cos \theta_{ij}) \]\[ J{1,2} \frac{\partial Pi}{\partial Vj} \frac{Vi}{Vj} (G{ij} \cos \theta{ij} B{ij} \sin \theta_{ij}) \]在Python中计算雅可比矩阵部分代码如下假设已有节点导纳矩阵Ybus电压幅值V和相角thetaJ np.zeros((2 * len(V), 2 * len(V))) for i in range(len(V)): for j in range(len(V)): if i j: J[2 * i, 2 * i] -V[i] * sum([V[k] * (Ybus[i, k].real * np.sin(theta[i] - theta[k]) - Ybus[i, k].imag * np.cos(theta[i] - theta[k])) for k in range(len(V)) if k! i]) J[2 * i 1, 2 * i 1] -V[i] * sum([V[k] * (Ybus[i, k].real * np.cos(theta[i] - theta[k]) Ybus[i, k].imag * np.sin(theta[i] - theta[k])) for k in range(len(V)) if k! i]) else: J[2 * i, 2 * j] -V[i] * V[j] * (Ybus[i, j].real * np.sin(theta[i] - theta[j]) - Ybus[i, j].imag * np.cos(theta[i] - theta[j])) J[2 * i, 2 * j 1] V[i] * (Ybus[i, j].real * np.cos(theta[i] - theta[j]) Ybus[i, j].imag * np.sin(theta[i] - theta[j])) J[2 * i 1, 2 * j] V[i] * V[j] * (Ybus[i, j].real * np.cos(theta[i] - theta[j]) Ybus[i, j].imag * np.sin(theta[i] - theta[j])) J[2 * i 1, 2 * j 1] -V[i] * (Ybus[i, j].real * np.sin(theta[i] - theta[j]) - Ybus[i, j].imag * np.cos(theta[i] - theta[j]))这段代码根据雅可比矩阵元素的公式通过双重循环填充雅可比矩阵。计算支路潮流修正量有了灵敏度矩阵和雅可比矩阵后我们就可以计算故障后支路潮流的修正量。假设故障点注入功率有一个增量 \(\Delta P{fault}\)那么支路潮流的修正量 \(\Delta P{branch}\) 可以通过灵敏度矩阵计算\[ \Delta P{branch} S \cdot \Delta P{fault} \]用Python代码表示如下假设已有灵敏度矩阵S和故障点注入功率增量deltaPfaultdelta_P_branch np.dot(S, delta_P_fault)这行代码简单地通过矩阵乘法得到支路潮流的修正量。通过以上步骤我们基于牛顿拉夫逊潮流计算结果成功引入灵敏度矩阵和雅可比矩阵计算出了故障后所有支路潮流的修正量为电力系统故障分析和运行优化提供了有力的工具。