CN104391825A - 一种基于高斯消元法快速求解电力系统节点阻抗矩阵的方法 - Google Patents

Abstract

一种基于高斯消元法快速求解电力系统节点阻抗矩阵的方法，涉及电力系统分析计算领域，主要步骤如下：输入节点导纳矩阵Y阵数据；用Y阵和E阵的最后一列E_n阵形成增广矩阵B_n＝[YE_n]；对B_n阵进行n-1次含规格化的高斯消元得B_n ^(n-1)′＝[Y^(n-1)′E_n′]；根据Y^(n-1)′Z_n＝E_n′求解节点阻抗矩阵Z阵中的Z_n阵；根据对称性得Z_nn以左的所有元素；根据Y^(k-1)′Z_k＝E_k′等回代求解Z阵第n-1～1列的Z_k阵对角元及以上元素；根据对称性得Z_kk以左的所有元素。本发明原理简单、思路清晰，用本发明方法对IEEE-57、-118、-300节点系统的Y阵求解Z阵，与LDU三角分解法和不含规格化的高斯消元法相比，计算速度大大提高。

Description

一种基于高斯消元法快速求解电力系统节点阻抗矩阵的方法

技术领域

本发明涉及电力系统分析计算领域，主要涉及一种快速求取电力系统节点阻抗矩阵的计算方法。

背景技术

由于节点阻抗矩阵Z阵是满矩阵，其所含的信息比节点导纳矩阵Y阵阵多，因此Z阵与Y阵一样，在电力系统分析计算中意义重大、应用广泛，但求解Z阵远比求解Y阵复杂。常用求解Z阵的方法有支路追加法和Y阵求逆法。在Y阵求逆法中最常用的是LDU三角分解法，而不含规格化的高斯消元法由于计算速度较慢应用很少。

派生于高斯消元法的LDU三角分解法适合于求解常系数方程，但计算原理及计算过程均比高斯消元法复杂。用LDU三角分解法求解Z阵主要特点是利用该方法可将对n*n阶Z阵的求解转换成对n个Z_k阵的求解。由于LDU三角分解法的计算过程中含有规格化的因子矩阵，因此求解Z阵时其计算效率高于不含规格化的高斯消元法。但实际上经过详细地比较、分析、计算可以发现，在各种三角分解法中，由于LDU三角分解法形成因子阵时计算各个元素所需的元素个数及中间矩阵变量的个数均高于LR、CU三角分解法，因此LDU三角分解法的计算效率实际上最低。此外，许多文献用LDU三角分解法求解Z阵时是求解n个整列的Z_k阵，进而求得Z阵，即Z_k阵的计算顺序为：Z₁,Z₂,┄,Z_k,┄,Z_n-1,Z_n，完全不用Z阵元素的对称性。这也导致用LDU三角分解法计算Z阵元素效率的降低。LDU三角分解法求解Z阵的计算流程图见图2。

不含规格化的高斯消元法计算原理简单、直观、程序编写方便，但计算速度并不理想，因此应用很少。该方法一般采用直接对Y阵求逆矩阵Z阵，计算过程繁琐，计算量大，计算时间长。或者由于对Z阵或Z_k阵的求解不含规格化计算，从而导致其计算速度远不如用LDU三角分解法理想。不含规格化的高斯消元法求解Z阵的计算流程图见图1。

根据计算分析，含规格化的高斯消元法比不含规格化的高斯消元法的计算速度相比，计算速度可提高约30％。因此，如果用含规格化的高斯消元法将对一个n*n阶Z阵的求解转换成对n个Z_k阵的求解，并利用Z阵元素的对称性，则其计算速度一定优于用LDU三角分解法求解Z阵的速度。

发明内容

本发明的目的是提供一种基于高斯消元法快速求解电力系统节点阻抗矩阵的方法，能大大提高电力系统分析计算中求解Z阵的计算速度。

本发明是通过以下技术方案实现的，基本步骤如下：

步骤1：输入节点导纳矩阵Y阵数据；

步骤2：将Y阵与E阵最后一列的E_n阵形成增广阵B_n＝[YE_n]；

步骤3：对B_n阵进行n-1次含规格化的高斯消元得B^(n-1)′＝[Y^(n-1)′E_n′]；

步骤4：根据Y^(n-1)′Z_n＝E_n′求解Z阵中的Z_n阵；

步骤5：根据对称性得Z_nn以左的所有元素；

步骤6：根据Y^(k-1)′Z_k＝E_k′等回代求解Z阵第n-1～1列的Z_k阵中对角元Z_kk及其以上的元素；

步骤7：根据对称性得Z_kk以左的所有元素；

步骤8：求出Z阵并输出结果。

步骤3中，本发明增加了规格化运算，因此可大大提高计算速度。但由于E阵结构的特殊性，E_n阵的对角元以上的元素不发生变化，仅是其对角元素从1变成Y^(n-1)′阵中第n行实际对角元素的倒数

步骤5中，根据步骤4得到Z_n阵以及Z阵元素的对称性可直接得到对角元素Z_nn以左的所有元素，即Z阵的最后一行的元素。这种计算方式可减少50％非对角元素的计算，因此可进一步提高计算速度。步骤8中也同样如此。

步骤6中，由于Y^(k-1)′阵第k行及以上的元素与Y^(n-1)′阵第k行及以上的元素完全相同，因此可直接用Y^(n-1)′阵第k行及以上的元素代替Y^(k-1)′阵，并与E_k′阵一起求解Z_k阵，而此时E_k′阵与E_k阵对角元及其以上元素唯一的不同就是其对角元素从1变成Y^(n-1)′阵中第k行实际对角元素的倒数且在求解Z_k阵时，其对角元以下的元素前期计算已经得出，无需计算。所以用E_k′阵求解Z_k阵时，Z_k阵与E_k′阵对角元以下的元素根本无需考虑。由于本发明依次从右向左、从下向上逐列、逐行地求取Z阵元素，其回代求解Z_k阵元素的个数将逐步减少，到第一列时只需求解Z₁₁一个元素。

本发明的方法可以根据Y阵快速求解Z阵，且原理及过程简单，程序编写方便，计算速度快。本发明方法的计算流程图见图3。用本发明对IEEE-57、-118、-300节点系统的Y阵求解Z阵，与传统的LDU三角分解法和不含规格化的高斯消元法方法相比，计算速度可分别提高约40％～60％(见实施例1)。

附图说明

图1为不含规格化的高斯消元法对Y阵求取Z阵的计算流程图。

图2为LDU三角分解法对Y阵求取Z阵的计算流程图。

图3为本发明对Y阵求取Z阵的计算流程图。

具体实施方式

本发明将通过以下实施例作进一步说明。

本方法由于从Z阵的最后一列开始求解，可将Y阵和E阵的最后一列E_n阵形成增广阵B_n＝[YE_n]。对B_n阵进行含规格化的高斯消元后得阵，展开如下。

由于E阵结构的特殊性，E_n阵对角元以上的元素并不变化，仅是其对角元素从1变成Y^(n-1)′阵中第n行实际对角元素的倒数

根据式(1)构成的方程Y^(n-1)′Z_n＝E_n′可求解Z阵的倒数第一列Z_n阵，再根据矩阵的对称性，得到Z阵的倒数第一行元素。然后继续依次求解Z阵的倒数第二列和第二行元素，直至最后。

将Y阵与E阵的第k列E_k阵形成增广阵B_k＝[YE_k]，并进行k-1次含规格化的高斯消元后可得阵，展开如下。

此时E_k阵对角元以上的元素也不变化，仅是其对角元素从1变成Y^(k-1)′阵中第k行实际对角元素的倒数比较式(1)和式(2)可以发现，Y^(k-1)′阵第k行及以上的元素与Y^(n-1)′阵第k行及以上的元素完全相同。因此后续计算中均可用Y^(n-1)′阵第k行及以上的元素以及相应的E_k′阵求解相应的Z_k阵元素。

式(2)对应的方程组为

根据式(3)的第k行可得式(4)，第j行可得式(5)。

$Z_{\mathrm{kk}}=\frac{1}{Y_{\mathrm{kk}}^{(k-1)}}-\underset{m=k+1}{\overset{n}{\mathrm{\Σ}}}{Y}_{\mathrm{km}}^{\left(k\right)}{Z}_{\mathrm{mk}},(k=n,n-1,n-2,......,\mathrm{2,1})---\left(4\right)$

$Z_{\mathrm{jk}}=-\underset{m=j+1}{\overset{n}{\mathrm{\Σ}}}{Y}_{\mathrm{jm}}^{\left(j\right)}{Z}_{\mathrm{mk}},(j=k-1,k-2,......,\mathrm{2,1})---\left(5\right)$

本发明的计算规律是依次从右向左、从下向上逐列、逐行地求取Z阵元素，对任一Z_k阵先求对角元Z_kk及其以上的所有元素，再根据对称性得Z_kk以左的所有元素，而对角元Z_kk以下的元素是先前的计算结果。

实施例1。

对IEEE-57、-118、-300节点系统在不考虑元素稀疏性的情况下求取其Y阵的逆矩阵Z阵，并比较不含规格化的高斯消元法、LDU三角分解法和本发明计算方法的计算时间。计算结果比较如下表所示。

T₁：不含规格化高斯消元法计算时间；T₂：LDU三角分解法计算时间；T₃：本方法计算时间。

T₂/T₁：；LDU三角分解法与不含规格化高斯消元法计算时间的百分比；

T₃/T₁：本方法与颔规格化高斯消元法计算时间的百分比；

T₃/T₂：本方法与LDU三角分解法计算时间的百分比。

从表1可以看出，LDU三角分解法与不含规格化高斯消元法相比，计算速度快约25～35％；本方法与不含规格化高斯消元法相比，计算速度可提高约60％；本方法与LDU三角分解法相比，计算速度可提高约40～50％。上述计算结果足以表明本发明方法的优势所在。

由于不含规格化的高斯消元法、LDU三角分解法在求解系统Z阵时都求取了整个Z阵元素，所以这里三种方法都按求取整个Z阵元素考虑。如果考虑Z阵元素的对称性仅求取Z阵的上三角元素，则本发明方法的计算速度还有一定程度的增加。

本发明可以采用任何一种编程语言和编程环境实现，这里采用C++编程语言，开发环境是Visual C++。

Claims (1)

1.一种基于高斯消元法快速求解电力系统节点阻抗矩阵的方法，其特征是包括以下步骤：

步骤1：输入节点导纳矩阵Y阵数据；

步骤2：Y阵与E阵最后一列的E_n阵形成增广阵B_n＝[YE_n]；

步骤3：对B_n阵进行n-1次含规格化的高斯消元得B^(n-1)′＝[Y^(n-1)′E_n′]；

步骤4：根据Y^(n-1)′Z_n＝E_n′求解Z阵中的Z_n阵；

步骤5：根据对称性得Z_nn以左的所有元素；

步骤6：根据Y^(k-1)′Z_k＝E_k′等回代求解Z阵第n-1～1列的Z_k阵中对角元Z_kk及其以上的元素；

步骤7：根据对称性得Z_kk以左的所有元素；

步骤8：求出Z阵并输出结果。