电力系统潮流计算的计算机算法（三）——牛顿-拉夫逊潮流算法

本篇为本科课程《电力系统稳态分析》的笔记。

本篇为这一章的第三篇笔记。上一篇传送门，下一篇传送门。

极坐标形式的牛顿-拉夫逊潮流算法

潮流方程式

一般选取平衡节点的电压相角为参考相角。

对于有n个节点的电力系统，设有一个平衡节点，m个PQ节点，n-m-1个PV节点。编号顺序为PQ节点、PV节点、平衡节点，即：

节点编号	种类	个数
$1,2,\cdots ,m$	PQ节点	$m$
$m+1,m+2,\cdots ,n-1$	PV节点	$n-m-1$
$n$	平衡节点	$1$

1. PQ节点的功率方程
   $$\begin{gathered} \Delta P_i=P_{Gi}-P_{Li}-U_i\sum _{j=1}^n{U}_j\left(G_{ij}\cos {\theta }_{ij}+B_{ij}\sin {\theta }_{ij}\right) \\ \Delta Q_i=Q_{Gi}-Q_{Li}-U_i\sum _{j=1}^n{U}_j\left(G_{ij}\sin {\theta }_{ij}-B_{ij}\cos {\theta }_{ij}\right) \\ (i=1,\cdots ,m) \end{gathered}$$
   其中，$P_i^{{}^{\prime }}=P_{Gi}-P_{Li},Q_i^{{}^{\prime }}=Q_{Gi}-Q_{Li}$是给定条件，而$U_i,{\theta }_i$为待求量。
2. PV节点的功率方程
   $$\begin{gathered} \Delta P_i=P_{Gi}-P_{Li}-U_i\sum _{j=1}^n{U}_j\left(G_{ij}\cos {\theta }_{ij}+B_{ij}\sin {\theta }_{ij}\right) \\ \Delta Q_i=Q_{Gi}-Q_{Li}-U_i\sum _{j=1}^n{U}_j\left(G_{ij}\sin {\theta }_{ij}-B_{ij}\cos {\theta }_{ij}\right) \\ (i=m+1,\cdots ,n-1) \end{gathered}$$
   其中，$P_i^{{}^{\prime }},U_i$是给定条件，而$Q_i^{{}^{\prime }},{\theta }_i$为待求量。
3. 平衡节点的功率方程
   $$\begin{gathered} \Delta P_n=P_{Gn}-P_{Ln}-U_n\sum _{j=1}^n{U}_j\left(G_{nj}\cos {\theta }_{nj}+B_{nj}\sin {\theta }_{nj}\right) \\ \Delta Q_n=Q_{Gn}-Q_{Ln}-U_n\sum _{j=1}^n{U}_j\left(G_{nj}\sin {\theta }_{nj}-B_{nj}\cos {\theta }_{nj}\right) \end{gathered}$$
   其中，$U_n,{\theta }_n$是给定条件，且${\theta }_n=0$，而$P_n^{{}^{\prime }},Q_n^{{}^{\prime }}$为待求量。

综合上面三类节点的功率方程，$\Delta P_i$的方程在PQ节点和PV节点存在，$\Delta Q_i$的方程在PQ节点存在，则实际需要求解的方程组为：
$$\begin{gathered} \Delta P_i=P_{Gi}-P_{Li}-U_i\sum _{j=1}^n{U}_j\left(G_{ij}\cos {\theta }_{ij}+B_{ij}\sin {\theta }_{ij}\right)(i=1,\cdots ,n-1) \\ \Delta Q_i=Q_{Gi}-Q_{Li}-U_i\sum _{j=1}^n{U}_j\left(G_{ij}\sin {\theta }_{ij}-B_{ij}\cos {\theta }_{ij}\right)(i=1,\cdots ,m) \end{gathered}$$

这里是一个n+m-1阶的适定非线性方程组，包括n-1个有功功率方程和m个无功功率方程。共计n+m-1个未知量，包括${\theta }_1,{\theta }_2,\cdots ,{\theta }_{n-1}$和$U_1,U_2,\cdots ,U_{n-1}$，而且待求相角个数等于有功功率方程的个数，待求电压值的个数等于无功功率方程的个数。

修正方程式及其求解

将电压相角和电压有效值写成向量形式。
$$\begin{gathered} \mathbf{\theta }=[{\theta }_1,{\theta }_2,\cdots ,{\theta }_{n-1}{]}^T \\ \mathbf{U}=[U_1,U_2,\cdots ,U_m{]}^T \end{gathered}$$

所以变量可以写为：
$$\mathbf{x}=\left[\begin{array}{c}\mathbf{\theta }\\ \mathbf{U}\end{array}\right]$$

有功功率和无功功率误差的求法：
$$\begin{gathered} \Delta P_i(\mathbf{x})=P_{Gi}-P_{Li}-U_i\sum _{j=1}^n{U}_j\left(G_{ij}\cos {\theta }_{ij}+B_{ij}\sin {\theta }_{ij}\right)=P_i^{{}^{\prime }}-P_i(\mathbf{x}) \\ \Delta Q_i(\mathbf{x})=Q_{Gi}-Q_{Li}-U_i\sum _{j=1}^n{U}_j\left(G_{ij}\sin {\theta }_{ij}-B_{ij}\cos {\theta }_{ij}\right)=Q_i^{{}^{\prime }}-Q_i(\mathbf{x}) \end{gathered}$$

然后写出修正方程，下面全部略去了表示迭代次数的上标k：
$$\left[\begin{array}{c}\Delta P_1\left(\mathbf{x}\right)\\ \Delta P_2\left(\mathbf{x}\right)\\ ⋮\\ \Delta P_{n-1}\left(\mathbf{x}\right)\\ \Delta Q_1\left(\mathbf{x}\right)\\ \Delta Q_2\left(\mathbf{x}\right)\\ ⋮\\ \Delta Q_m\left(\mathbf{x}\right)\end{array}\right]=\left[\begin{array}{cccccccc}{H}_{11}& H_{12}& \cdots & H_{1,m-1}& N_{11}& N_{12}& \cdots & N_{1,m}\\ H_{21}& H_{22}& \cdots & H_{2,m-1}& N_{21}& N_{22}& \cdots & N_{2,m}\\ ⋮& ⋮& ⋮& \cdots & ⋮& ⋮& ⋮& \cdots \\ H_{n-1,1}& H_{n-1,2}& \cdots & H_{n-1,m-1}& N_{n-1,1}& N_{n-1,2}& \cdots & N_{n-1,m}\\ M_{11}& M_{12}& \cdots & M_{1,m-1}& L_{11}& L_{12}& \cdots & L_{1,m}\\ M_{21}& M_{22}& \cdots & M_{2,m-1}& L_{21}& L_{22}& \cdots & L_{2,m}\\ ⋮& ⋮& \cdots & ⋮& ⋮& ⋮& \cdots & ⋮\\ M_{m,1}& M_{m,2}& \cdots & M_{m,m-1}& L_{m,1}& L_{m,2}& \cdots & L_{m,m}\end{array}\right]\left[\begin{array}{c}\Delta {\theta }_1\\ \Delta {\theta }_2\\ ⋮\\ \Delta {\theta }_{n-1}\\ \Delta U_1/U_1\\ \Delta U_2/U_2\\ ⋮\\ \Delta U_m/U_m\end{array}\right]$$

矩阵中，使用的是$\Delta U_i/U_i$而不是$\Delta U_i$，只是为了让表达式在形式上保持一致。

雅可比矩阵各个非对角元素的计算为：
$$\begin{array}{c}{H}_{ij}=\frac{\partial \Delta P_i}{\partial {\theta }_j}=-U_i{U}_j\left(G_{ij}\sin {\theta }_{ij}-B_{ij}\cos {\theta }_{ij}\right)\\ N_{ij}=\frac{\partial \Delta P_i}{\partial U_j}{U}_j=-U_i{U}_j\left(G_{ij}\cos {\theta }_{ij}+B_{ij}\sin {\theta }_{ij}\right)\\ M_{ij}=\frac{\partial \Delta Q_i}{\partial {\theta }_j}=U_i{U}_j\left(G_{ij}\cos {\theta }_{ij}+B_{ij}\sin {\theta }_{ij}\right)\\ L_{ij}=\frac{\partial \Delta Q_i}{\partial U_j}{U}_j=-U_i{U}_j\left(G_{ij}\sin {\theta }_{ij}-B_{ij}\cos {\theta }_{ij}\right)\end{array}$$

而且有以下的关系：
$$\begin{gathered} H_{ij}=L_{ij} \\ {N}_{ij}=-M_{ij} \end{gathered}$$

各个对角元素的计算为：
$$\begin{array}{rl}& H_{ii}=\frac{\partial \Delta P_i}{\partial {\theta }_i}=U_i\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^n{U}_j\left(G_y\sin {\theta }_y-B_y\cos {\theta }_y\right)=U_i^2{B}_{ii}+Q_i\left(x\right)\\ & N_{ii}=\frac{\partial \Delta P_i}{\partial U_i}{U}_i=-U_i\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^n{U}_j\left(G_y\cos {\theta }_y+B_y\sin {\theta }_y\right)-2U_i^2{G}_{ii}\\ & =-U_i^2{G}_{ii}-P_i\left(\mathbf{x}\right)\\ & M_{ii}=\frac{\partial \Delta Q_i}{\partial {\theta }_i}=-U_i\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^n{U}_j\left(G_{ij}\cos {\theta }_{ij}+B_{ij}\sin {\theta }_{ij}\right)=U_i^2{G}_{ii}-P_i\left(x\right)\\ & L_{ii}=\frac{\partial \Delta Q_i}{\partial U_i}{U}_i=-U_i\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^n{U}_j\left(G_{ij}\sin {\theta }_{ij}-B_{ij}\cos {\theta }_{ij}\right)+2U_i^2{B}_{ii}\\ & =U_i^2{B}_{ij}-Q_i\left(\mathbf{x}\right)\end{array}$$

修正方程可以简写为：
$$\left[\begin{array}{c}\Delta {\mathbf{P}}^{(k)}\\ \Delta {\mathbf{Q}}^{(k)}\end{array}\right]=-\left[\begin{array}{cc}{\mathbf{H}}^{(k)}& {\mathbf{N}}^{(k)}\\ {\mathbf{M}}^{(k)}& {\mathbf{L}}^{(k)}\end{array}\right]\left[\begin{array}{c}\Delta {\mathbf{\theta }}^{(k)}\\ \Delta {\tilde{\mathbf{U}}}^{(k)}\end{array}\right]=-{\mathbf{J}}^{(k)}\left[\begin{array}{c}\Delta {\mathbf{\theta }}^{(k)}\\ \Delta {\tilde{\mathbf{U}}}^{(k)}\end{array}\right]$$

其中，

• $\Delta {\mathbf{\theta }}^{(k)}=[\Delta {\theta }_1^{(k)}\Delta {\theta }_2^{(k)}\cdots \Delta {\theta }_{n-1}^{(k)}{]}^T$，$\Delta {\tilde{\mathbf{U}}}^{(k)}=[\Delta U_1^{(k)}\Delta U_2^{(k)}\cdots \Delta U_{n-1}^{(k)}{]}^T$。
• 而雅可比矩阵${\mathbf{J}}^{(k)}$中的分块矩阵${\mathbf{H}}^{(k)},{\mathbf{N}}^{(k)},{\mathbf{M}}^{(k)},{\mathbf{L}}^{(k)}$是在未知量的取值为${\mathbf{\theta }}^{(k)},{\tilde{\mathbf{U}}}^{(k)}$时，带入上面元素的求取公式而得到的。修正方程已知雅可比矩阵和误差的矩阵，则可以求出来变量修正量矩阵，即$\Delta {\mathbf{\theta }}^{(k)},\Delta {\tilde{\mathbf{U}}}^{(k)}$，从而可以得出下一代的变量取值：
  $$\begin{gathered} {\theta }_i^{(k+1)}={\theta }_i^{(k)}+\Delta {\theta }_i^{(k)}(i=1,\cdots ,n-1) \\ {U}_i^{(k+1)}=U_i^{(k)}+\Delta U_i^{(k)}\times {\tilde{U}}_i^{(k)}(i=1,\cdots ,m) \end{gathered}$$

以上就是计算潮流的迭代格式，可见形成雅可比矩阵和求解修正方程式牛顿潮流计算中的主要步骤。

雅可比矩阵具有以下的特点：

1. 矩阵的阶数为n+m-1.
2. 如果互导纳$Y_{ij}=0$，那么其中的分块矩阵$\mathbf{H},\mathbf{N},\mathbf{M},\mathbf{L}$的相应元素也为零。所以雅可比矩阵是稀疏矩阵。
3. 矩阵的两个对角子矩阵$\mathbf{H},\mathbf{L}$是不对称的，因为：
   $$\frac{\partial \Delta P_i}{\partial {\theta }_j}\ne \frac{\partial \Delta P_j}{\partial {\theta }_i};\frac{\partial \Delta Q_i}{\partial U_j}\ne \frac{\partial \Delta Q_j}{\partial U_i}$$
4. 矩阵的各个元素都是节点电压的有效值和相角的函数，所以在整个迭代过程中，所有元素都会随着节点电压相量的逐次修正而不断变化。所以牛顿潮流计算中的雅可比矩阵要在每一次迭代进行计算。

牛顿潮流计算收敛判据一般为：
$$||\Delta {\mathbf{P}}^{(k)}|{|}_{\infty }=\underset{i}{\max }|\Delta P_i^{(k)}|⩽\epsilon 且||\Delta {\mathbf{Q}}^{(k)}|{|}_{\infty }=\underset{i}{\max }|\Delta Q_i^{(k)}|⩽\epsilon$$

初值的给定

一般取电压有效值为1，即相等于额定电压，电压相角为0。

牛顿法潮流计算的步骤

1. 输入系统的原始数据。
2. 形成节点导纳矩阵。
3. 给定功率误差的容许值$\epsilon$，给定平衡节点电压相角为0，是参考相角，然后给定节点电压有效值和相角初值：$U_i^0=1(i=1,\cdots ,m),{\theta }_i^0=1(i=1,\cdots ,n-1)$，将各个PQ节点电压有效值的初值组成相量${\mathbf{U}}^{(0)}$，各PQ节点和PV节点电压相角的初值组成向量${\mathbf{\theta }}^{(0)}$。
4. 设置迭代次数k=0。
5. 应用${\mathbf{\theta }}^{(k)},{\mathbf{U}}^{(k)}$，各个PV节点电压有效值和平衡节点电压有效值和相角，按照有功功率和无功功率的误差计算公式计算各个PQ、PV节点的有功功率误差$\Delta P_i^{(k)}$，还有各个PQ节点的无功功率误差$\Delta Q_i^{(k)}$，并组成功率误差向量$\Delta {\mathbf{P}}^{(k)},\Delta {\mathbf{Q}}^{(k)}$。
6. 按照收敛判据，判断迭代是否收敛。如果收敛，则转向第11步，否则进行下一步。
7. 应用${\mathbf{\theta }}^{(k)},{\mathbf{U}}^{(k)}$，计算出雅可比矩阵的各个元素，这样形成了修正方程中的雅可比矩阵${\mathbf{J}}^{(k)}$。
8. 根据求出来的雅可比矩阵${\mathbf{J}}^{(k)}$和第5步计算出来的误差向量$\Delta {\mathbf{P}}^{(k)},\Delta {\mathbf{Q}}^{(k)}$，按照修正方程式计算出$\Delta {\mathbf{\theta }}^{(k)},\Delta {\tilde{\mathbf{U}}}^{(k)}$。
9. 然后根据修正量$\Delta {\mathbf{\theta }}^{(k)},\Delta {\tilde{\mathbf{U}}}^{(k)}$，求出下一步的电压有效值和相角。
10. 令k=k+1，返回第5步进入下一轮的迭代。
11. 计算平衡节点发出的有功功率和无功功率，计算PV节点的无功功率等。

直角坐标形式的牛顿-拉夫逊潮流算法

仍然假设$1,2,\cdots ,m$是PQ节点，$m+1,\cdots ,n-1$是PV节点，n是平衡节点。

有功功率和无功功率的误差求法为
$$\begin{gathered} \Delta P_i=P_{Gi}-P_{Li}-e_i\sum _{j=1}^n\left(G_{ij}{e}_j-B_{ij}{f}_j\right)-f_i\sum _{j=1}^n\left(G_{ij}{f}_j+B_{ij}{e}_j\right)\left(i=1,2,\cdots ,n-1\right) \\ \Delta Q_i=Q_{Gi}-Q_{Li}-f_i\sum _{j=1}^n\left(G_{ij}{e}_j-B_{ij}{f}_j\right)+e_i\sum _{j=1}^n\left(G_{ij}{f}_j+B_{ij}{e}_j\right)\left(i=1,2,\cdots ,m\right) \\ \Delta U_i^2=U_S^2-\left(e_i^2+f_i^2\right)=0\left(i=m+1,m+2,\cdots ,n-1\right) \end{gathered}$$
其中，$U_{Si}$是PV节点i的电压给定值，平衡节点的电压$e_n=U_{Sn},f_n=0$

修正方程式为：
$$\left[\begin{array}{c}\Delta P_1\\ \Delta Q_1\\ ⋮\\ \Delta P_m\\ \Delta Q_m\\ \cdots \\ \Delta P_{m+1}\\ \Delta U_{m+1}^2\\ ⋮\\ \Delta P_{m-1}\\ \Delta U_{n-1}^2\end{array}\right]=-\left[\begin{array}{ccccccccccc}{H}_{1,1}& N_{1,1}& \cdots & H_{1,m}& N_{1,m}& ⋮& H_{1,m+1}& N_{1,m+1}& \cdots & H_{1,m-1}& N_{1,m-1}\\ M_{1,1}& L_{1,1}& \cdots & M_{1,m}& L_{1,m}& ⋮& M_{1,m+1}& L_{1,m+1}& \cdots & M_{1,m-1}& L_{1,m-1}\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ H_{m,1}& N_{m,1}& \cdots & H_{m,m}& N_{m,m}& ⋮& H_{m,m+1}& N_{m,m+1}& \cdots & H_{m,m+1}& N_{m,m+1}\\ M_{m,1}& L_{m,1}& \cdots & M_{m,m}& L_{m,m}& ⋮& M_{m,m+1}& L_{m,m+1}& \cdots & M_{m,m+1}& L_{m,m-1}\\ \cdots & \cdots & \cdots & \cdots & \cdots & ⋮& \cdots & \cdots & \cdots & \cdots & \cdots \\ H_{m+1,1}& N_{m+1,1}& \cdots & H_{m+1,m}& N_{m+1,m}& ⋮& H_{m+1,m+1}& N_{m+1,m+1}& \cdots & H_{m+1,n-1}& N_{m+1,n-1}\\ R_{m+1,1}& S_{m+1,1}& \cdots & R_{m+1,m}& S_{m+1,m}& ⋮& R_{m+1,m+1}& S_{m+1,m+1}& \cdots & R_{m+1,n-1}& S_{m+1,n-1}\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ H_{n-1,1}& N_{n-1,1}& \cdots & H_{n-1,m}& N_{n-1,m}& ⋮& H_{n-1,m+1}& N_{n-1,m+1}& \cdots & H_{n-1,n-1}& N_{n-1,n-1}\\ R_{n-1,1}& S_{n-1,1}& \cdots & R_{n-1,m}& S_{n-1,m}& ⋮& R_{n-1,m+1}& S_{n-1,m+1}& \cdots & R_{n-1,n-1}& S_{n-1,n-1}\end{array}\right]\left[\begin{array}{c}\Delta f_1\\ \Delta e_1\\ ⋮\\ \Delta f_m\\ \Delta e_m\\ \cdots \\ \Delta f_{m+1}\\ \Delta e_{m+1}\\ ⋮\\ \Delta f_{n-1}\\ \Delta e_{n-1}\end{array}\right]$$

其中，非对角元素$(i\ne j)$计算方法为：
$$\begin{array}{rl}& H_{ij}=\frac{\partial \Delta P_i}{\partial f_j}=B_{ij}{e}_i-G_{ij}{f}_i{N}_{ij}=\frac{\partial \Delta P_i}{\partial e_j}=-G_{ij}{e}_i-B_{ij}{f}_i\\ & M_{ij}=\frac{\partial \Delta Q_i}{\partial f_j}=B_{ij}{f}_i+G_{ij}{e}_i=-N_{ij}{L}_{ij}=\frac{\partial \Delta Q_i}{\partial e_j}=-G_{ij}{f}_i+B_{ij}{e}_i=H_{ij}\\ & R_{ij}=\frac{\partial \Delta U_i^2}{\partial f_j}=0S_{ij}=\frac{\partial \Delta U_i^2}{\partial e_j}=0\end{array}$$

对角元素$(i=j)$的计算方法：
$$\begin{gathered} H_{ii}=\frac{\partial \Delta P_i}{\partial f_i}=-\sum _{j=1}^n\left(G_{ij}{f}_j+B_{ij}{e}_j\right)-G_{ii}{f}_i+B_{ii}{e}_i \\ {N}_{ii}=\frac{\partial \Delta P_i}{\partial e_i}=-\sum _{j=1}^n\left(G_{ij}{e}_j-B_{ij}{f}_j\right)-G_{ii}{e}_i-B_{ii}{f}_i \\ {J}_{ii}=\frac{\partial \Delta Q_i}{\partial f_i}=-\sum _{j=1}^n\left(G_{ij}{e}_j-B_{ij}{f}_j\right)+G_{ii}{e}_i+B_{ii}{f}_i \\ {L}_{ii}=\frac{\partial \Delta Q_i}{\partial e_i}=\sum _{j=1}^n\left(G_{ij}{f}_j+B_{ij}{e}_j\right)-G_{ii}{f}_i+B_{ii}{e}_i \\ {R}_{ii}=\frac{\partial \Delta U_i^2}{\partial f_i}=-2f_i \\ {S}_{ii}=\frac{\partial \Delta U_i^2}{\partial e_i}=-2e_i \end{gathered}$$

将修正方程写成下面的缩写形式：
$$\left[\begin{array}{c}\Delta \mathbf{P}\\ \Delta \mathbf{Q}\\ \Delta {\mathbf{U}}^2\end{array}\right]=-\left[\begin{array}{cc}\mathbf{H}& \mathbf{N}\\ \mathbf{M}& \mathbf{L}\\ \mathbf{R}& \mathbf{S}\end{array}\right]\left[\begin{array}{c}\Delta \mathbf{f}\\ \Delta \mathbf{e}\end{array}\right]=-\mathbf{J}\left[\begin{array}{c}\Delta \mathbf{f}\\ \Delta \mathbf{e}\end{array}\right]$$