电力系统潮流计算的计算机算法（四）——PQ快速解耦潮流算法

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

本篇为这一章的第四篇笔记。上一篇传送门。

潮流计算的快速解耦法

牛顿-拉夫逊法潮流计算，主要的工作量在于形成雅可比矩阵和求解修正方程。由于雅可比矩阵的阶数为n+m-1，约为节点总数的两倍，非对称矩阵，且在迭代过程中需要不断的变化，所以大规模的电力系统中应用该算法很费时费力。

PQ分解法潮流计算的修正方程

有一种快速解耦法，是最有效且应用最广。电力系统中线路的电阻远小于电抗，即$R\ll X$，换成导纳表示即$G\ll B$，说明可以忽略所有的G，令$G=0$。前面一章已经推理过，从物理意义上，有功功率P主要和电压相角$\theta$相关，无功功率主要和电压有效值U相关。

节点的分类保持不变。修正方程式如下：
$$\left[\begin{array}{c}\Delta \mathbf{P}\\ \Delta \mathbf{Q}\end{array}\right]=-\left[\begin{array}{cc}\mathbf{H}& \mathbf{N}\\ \mathbf{M}& \mathbf{L}\end{array}\right]\left[\begin{array}{c}\Delta \mathbf{\theta }\\ \Delta \tilde{\mathbf{U}}\end{array}\right]$$

根据有功功率P主要和电压相角$\theta$相关，无功功率主要和电压有效值U相关，可以将$\mathbf{N},\mathbf{M}$置为零，从而将修正方程化为：
$$\left[\begin{array}{c}\Delta \mathbf{P}\\ \Delta \mathbf{Q}\end{array}\right]=-\left[\begin{array}{cc}\mathbf{H}& 0\\ 0& \mathbf{L}\end{array}\right]\left[\begin{array}{c}\Delta \mathbf{\theta }\\ \Delta \tilde{\mathbf{U}}\end{array}\right]$$

所以得到两个独立的方程，第一个式子只和相角和有功功率有关，第二个式子只和电压和无功功率有关：
$$\begin{gathered} -\mathbf{H}\Delta \mathbf{\theta }=\Delta \mathbf{P} \\ -\mathbf{L}\Delta \tilde{\mathbf{U}}=\Delta \mathbf{Q} \end{gathered}$$

另外，各个节点之间的相角差很小，即$\theta \approx 0°$，再考虑到$G\ll B$，所以有$\cos \theta \approx 1,G\cos \theta \ll B$。

所以，雅可比矩阵的各个元素可以化简为：
$$\begin{array}{rl}& H_{ij}=U_i{U}_j{B}_{ij}\\ & H_{ii}=-U_i\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^n{U}_j{B}_{ij}=-U_i\sum _{j=1}^n{U}_j{B}_{ij}+U_i^2{B}_{ii}=Q_i+U_i^2{B}_{ii}\\ & L_{ij}=U_i{U}_j{B}_{ij}\\ & L_{ii}=U_i\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^n{U}_j{B}_{ij}+2U_i^2{B}_{ii}=U_i\sum _{j=1}^n{U}_j{B}_{ij}+U_i^2{B}_{ii}=-Q_i+U_i^2{B}_{ii}\end{array}$$

按照导纳的定义，上面两个式子中$U_i^2{B}_{ii}$应是各元件电抗远大于电阻的前提下，除了节点i之外的其他节点都接地的时候，由节点i注入的无功功率。它必然远大于正常运行时节点i的注入无功功率，即$U_i^2{B}_{ii}\gg Q_i$，则可以化简：
$$\begin{gathered} H_{ii}=U_i^2{B}_{ii} \\ {L}_{ii}=U_i^2{B}_{ii} \end{gathered}$$

这样，矩阵$\mathbf{H},\mathbf{L}$具有相同的表达式，但是他们的阶数不同，前者为(n-1)阶，后者为(m-1)阶。

矩阵$\mathbf{H}$可以表示为：
$$boldsymbolH=\left[\begin{array}{cccc}{U}_1^2{B}_{11}& U_1{U}_2{B}_{12}& \cdots & U_1{U}_{n-1}{B}_{1,n-1}\\ U_2{U}_1{B}_{21}& U_2^2{B}_{22}& \cdots & U_2{U}_{n-1}{B}_{2,n-1}\\ ⋮& ⋮& \ddots & ⋮\\ U_{n-1}{U}_1{B}_{n-1,1}& U_{n-1}{U}_2{B}_{n-1,2}& \cdots & U_{n-1}^2{B}_{n-1,n-1}\end{array}\right]=\left[\begin{array}{cccc}{U}_1& 0& \cdots & 0\\ 0& U_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & U_{n-1}\end{array}\right]\left[\begin{array}{cccc}{B}_{11}& B_{12}& \cdots & B_{1.n-1}\\ B_{21}& B_{22}& \cdots & ⋮\\ ⋮& ⋮& \ddots & ⋮\\ B_{n-1,1}& B_{n-1,2}& \cdots & B_{n-1,n-1}\end{array}\right]\left[\begin{array}{cccc}{U}_1& 0& \cdots & 0\\ 0& U_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & U_{n-1}\end{array}\right]$$

而矩阵$\mathbf{L}$可以表示为：
$$\mathbf{L}=\left[\begin{array}{cccc}{U}_1^2{B}_{11}& U_1{U}_2{B}_{12}& \cdots & U_1{U}_m{B}_{1m}\\ U_2{U}_1{B}_{21}& U_2^2{B}_{22}& \cdots & U_2{U}_m{B}_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ U_m{U}_1{B}_{m1}& U_m{U}_2{B}_{m2}& \cdots & U_m^2{B}_{mm}\end{array}\right]=\left[\begin{array}{cccc}{U}_1& 0& \cdots & 0\\ 0& U_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & U_m\end{array}\right]\left[\begin{array}{cccc}{B}_{11}& B_{12}& \cdots & B_{1m}\\ B_{21}& B_{22}& \cdots & B_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ B_{m1}& B_{m2}& \cdots & B_{mm}\end{array}\right]\left[\begin{array}{cccc}{U}_1& 0& \cdots & 0\\ 0& U_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & U_m\end{array}\right]$$

然后将变化过的矩阵$\mathbf{H},\mathbf{L}$带入上面说过的两个独立方程式内，可得新的修正方程，第一个式子为：
$$-\left[\begin{array}{cccc}{U}_1& 0& \cdots & 0\\ 0& U_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & U_{n-1}\end{array}\right]\left[\begin{array}{cccc}{B}_{11}& B_{12}& \cdots & B_{1,n-1}\\ B_{21}& B_{22}& \cdots & B_{2,n-1}\\ ⋮& ⋮& \ddots & ⋮\\ B_{n-1,1}& B_{n-1,2}& \cdots & B_{n-1,n-1}\end{array}\right]\left[\begin{array}{c}{U}_1\Delta {\theta }_1\\ U_2\Delta {\theta }_2\\ ⋮\\ U_{n-1}\Delta {\theta }_{n-1}\end{array}\right]=\left[\begin{array}{c}\Delta P_1\\ \Delta P_2\\ ⋮\\ \Delta P_{n-1}\end{array}\right]$$
进一步化简为：
$$-\left[\begin{array}{cccc}{B}_{11}& B_{12}& \cdots & B_{1,n-1}\\ B_{21}& B_{22}& \cdots & B_{2,n-1}\\ ⋮& ⋮& \ddots & ⋮\\ B_{n-1,1}& B_{n-1,2}& \cdots & B_{n-1,n-1}\end{array}\right]\left[\begin{array}{c}{U}_1\Delta {\theta }_1\\ U_2\Delta {\theta }_2\\ ⋮\\ U_{n-1}\Delta {\theta }_{n-1}\end{array}\right]=\left[\begin{array}{c}\Delta P_1/U_1\\ \Delta P_2/U_2\\ ⋮\\ \Delta P_{n-1}/U_{n-1}\end{array}\right]$$

第二个式子为：
$$-\left[\begin{array}{cccc}{U}_1& 0& \cdots & 0\\ 0& U_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & U_m\end{array}\right]\left[\begin{array}{cccc}{B}_{11}& B_{12}& \cdots & B_{1m}\\ B_{21}& B_{22}& \cdots & B_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ B_{m1}& B_{m2}& \cdots & B_{mm}\end{array}\right]\left[\begin{array}{c}\Delta U_1\\ \Delta U_2\\ ⋮\\ \Delta U_m\end{array}\right]=\left[\begin{array}{c}\Delta Q_1\\ \Delta Q_2\\ ⋮\\ \Delta Q_m\end{array}\right]$$
进一步化简为：
$$-\left[\begin{array}{cccc}{B}_{11}& B_{12}& \cdots & B_{1m}\\ B_{21}& B_{22}& \cdots & B_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ B_{m1}& B_{m2}& \cdots & B_{mm}\end{array}\right]\left[\begin{array}{c}\Delta U_1\\ \Delta U_2\\ ⋮\\ \Delta U_m\end{array}\right]=\left[\begin{array}{c}\Delta Q_1/U_1\\ \Delta Q_2/U_2\\ ⋮\\ \Delta Q_m/U_m\end{array}\right]$$

则第一个方程式可以简写为：
$${\mathbf{B}}^{\prime }\Delta \mathbf{\theta }\circ \mathbf{U}=\Delta \mathbf{P}/\mathbf{U}$$

其中，$\circ$表示两个矩阵逐元素相乘，$/$表示两个矩阵逐元素相除。矩阵${\mathbf{B}}^{\prime }$是$-\mathbf{B}$的前n-1行前n-1列，是对称常定矩阵，$\Delta \mathbf{\theta },\Delta \mathbf{P},\mathbf{U}$均为n-1行的列向量。

第二个方程式可以简写为：
$${\mathbf{B}}^{\prime \prime }\Delta \mathbf{U}=\Delta \mathbf{Q}/\mathbf{U}$$

其中，$/$表示两个矩阵逐元素相除。矩阵${\mathbf{B}}^{\prime \prime }$是$-\mathbf{B}$的前m行前m列，是对称常定矩阵，$\Delta \mathbf{\theta },\mathbf{U},\Delta \mathbf{P}$均为m行的列向量。

将上述方程化简，可得：
$$\begin{gathered} {\mathbf{B}}^{\prime }\Delta \tilde{\mathbf{\theta }}=\Delta \tilde{\mathbf{P}} \\ {\mathbf{B}}^{\prime \prime }\Delta \mathbf{U}=\Delta \tilde{\mathbf{Q}} \end{gathered}$$

其中，$\Delta \tilde{\mathbf{\theta }}=[U_1\Delta {\theta }_1{U}_2\Delta {\theta }_2\cdots U_{n-1}\Delta {\theta }_{n-1}{]}^T$，$\Delta \tilde{\mathbf{P}}=[\Delta P_1/U_1\Delta P_2/U_2\cdots \Delta P_{n-1}/U_{n-1}{]}^T$，\Delta\widetilde{\boldsymbol{Q}}=[\Delta Q_1/U_1\quad \Delta Q_2/U_2\quad \cdots \quad \Delta Q_{m}/U_{m}]^T。第一个方程式为相位修正方程，第二个方程为幅值修正方程。

PQ分解法潮流计算的步骤

计算方法的前三步和牛顿法相同，而后面的迭代过程为：

1. 形成两个系数矩阵${\mathbf{B}}^{\prime },{\mathbf{B}}^{\prime \prime }$，并对他们做三角分解，即Cholesky分解，如下所示：
   $$\begin{gathered} {\mathbf{B}}^{\prime }={\mathbf{L}}^{\prime }{\mathbf{D}}^{\prime }({\mathbf{L}}^{\prime }{)}^T \\ {\mathbf{B}}^{\prime \prime }={\mathbf{L}}^{\prime \prime }{\mathbf{D}}^{\prime \prime }({\mathbf{L}}^{\prime \prime }{)}^T \end{gathered}$$
2. 设置迭代次数为k=0。
3. 用${\mathbf{\theta }}^{(k)},{\mathbf{U}}^{(k)}$、各PV节点的电压有效值和平衡节点的电压有效值以及相角，按照有功功率误差的计算公式计算出各个PQ节点和PV节点的$\Delta P_i^{(k)}$，然后再进一步计算出$\Delta {\tilde{P}}_i^{(k)}$。这一步不要计算无功功率误差。
4. 通过前代和回代求解出修正方程的第一个式子，得到各节点的相角修正值$\Delta {\theta }_i^{(k)}(i=1,2,\cdots ,n-1)$。
5. 修正各个节点的相角值，即${\theta }_i^{(k+1)}={\theta }_i^{(k)}+\Delta {\theta }_i^{(k)}(i=1,2,\cdots ,n-1)$。
6. 利用${\mathbf{\theta }}^{(k+1)},{\mathbf{U}}^{(k)}$、各PV节点的电压有效值和平衡节点的电压有效值以及相角，按照无功功率误差的计算公式计算出各个PQ节点、的$\Delta Q_i^{(k)}$，然后再进一步计算出$\Delta {\tilde{Q}}_i^{(k)}$。
7. 通过前代和回代求解出修正方程的第二个式子，得到各节点的电压有效值的修正值$\Delta U_i^{(k)}(i=1,2,\cdots ,m)$。
8. 修正各个节点的电压有效值，即$U_i^{(k+1)}=U_i^{(k)}+\Delta U_i^{(k)}(i=1,2,\cdots ,m)$。
9. 使用收敛判据判断是否收敛，若收敛，则执行牛顿法潮流计算步骤中的第11步，然后结束，若不收敛，则设置k=k+1，然后返回第2步，进行下一轮的迭代。