CN104572584A - 一种基于对称稀疏矩阵技术的改进lr三角分解求取电力系统节点阻抗矩阵的方法 - Google Patents

Abstract

一种基于对称稀疏矩阵技术的改进LR三角分解求取电力系统节点阻抗矩阵Z的方法，属于电力系统分析计算领域。包括以下步骤：读取数据文件；形成节点导纳矩阵Y；根据稀疏性和对称性对Y阵进行LR三角分解求R阵元素；根据R阵元素稀疏性求Z_k阵对角元Z_kk及以上元素；按对称性求Z_kk以左元素；写Z阵数据到数据文件。本发明方法根据稀疏性和对称性按过程法对Y阵进行LR三角分解，仅求上三角的R阵元素，大幅提高三角分解速度；利用单位矩阵E阵结构特点省略W阵元素的计算，利用R阵元素的稀疏性直接解方程RZ_k＝E_k，大幅提高回代求解速度。用本发明方法对IEEE-30、-57、-118节点系统进行验算，与传统的LR三角分解法相比，计算速度可提高约83～98％。

Description

一种基于对称稀疏矩阵技术的改进LR三角分解求取电力系统节点阻抗矩阵的方法

技术领域

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

背景技术

在电力系统中用三角分解法求取节点阻抗矩阵Z时一般都使用LDU三角分解法，但实际上，由于计算过程和计算变量的不同，LDU三角分解法的计算效率低于LR三角分解法，因此用LR三角分解法求取Z阵元素是更好的选择。

传统的LR三角分解法对Y阵进行三角分解时，L、R二个因子阵都要形成。且形成过程中各个元素均按“┘”(反L)方式或按“行”方式一个个用计算公式一步形成(简称公式法)。因此，在其形成因子阵的过程中根本无法利用L阵、R阵元素的稀疏性和对称性特点，从而导致对大量零元素和部分非零元素不必要的计算，使计算效率大大降低，三角分解过程的速度极慢。

传统的LR三角分解法在回代过程中由于是按整列求取Z_k阵元素，未利用Z阵元素的对称性特点，也未利用单位矩阵E元素结构的特点，其回代过程包括求解方程LW_k＝E_k和RZ_k＝W_k，且对方程LW_k＝E_k要求解整个W_k阵。因此其回代过程存在大量不必要的计算，计算效率极低。此外，传统的LR三角分解法在回代过程也未利用稀疏性，从而进一步导致计算效率的降低。

其它三角分解法在三角分解和回代过程中也存在类似的问题。

电力系统计算中稀疏矩阵技术运用很广，主要为省去大量零元素的存贮及计算，加快高斯消元法的计算速度。矩阵元素的存贮方案也很多，如按坐标存贮、按顺序存贮、按链表存贮等等。尽管这些存贮方式可以省去不少存贮单元，但计算速度并没有达到最优效果，而且这些存贮方式结构复杂，且对角元素与非对角元素分开存贮也使得存取过程繁琐，特别不利于对称矩阵中的数据处理。实际上，这些存贮方式主要为减少存贮单元，对存贮过程的简化或存贮速度的提高并没有特别优势。而且这些存贮方式主要用于高斯消元法中，很难用于公式法的三角分解法中。且由于传统的稀疏矩阵技术一般不考虑矩阵元素结构的特点对非零元素进行存贮，因此其存贮方式在用公式法进行LR三角分解时无法利用L、R因子阵元素的稀疏性、对称性及其相互间的关系等特点。其它三角分解法也有类似问题。

发明内容

为了克服上述现有技术的不足，本发明提供了一种基于对称稀疏矩阵技术的改进LR三角分解快速求取电力系统节点阻抗矩阵的方法。

本发明是通过以下技术方案实现的，主要包括以下步骤：

步骤1：读入n节点系统各线路支路数据文件；

步骤2：形成节点导纳矩阵Y；

步骤3：根据稀疏性和对称性对Y阵进行LR三角分解，但仅求R阵元素，并记录其非零元素位置；

本发明所述的步骤3中具体实施过程如下：

(1)本发明方法中找出了L、R阵元素的对应关系，在三角分解时仅需形成用于回代过程的R阵，将L阵作为R阵的映射矩阵处理。

下述分析中认为L阵及其元素只是形式上存在，而实际上并未形成。对计算中所需使用的L阵元素通过其与R阵元素的对应关系获取。

(2)以4阶的L、R阵构成的合成阵为例，其各元素与Y阵元素的关系如式(1)。

$\left[\begin{array}{cccc}{r}_{11}=Y_{11}& r_{12}=Y_{12}& r_{13}=Y_{13}& r_{14}=Y_{14}\\ l_{21}=Y_{21}/r_{11}& r_{22}=Y_{22}-l_{21}{r}_{12}& r_{23}=Y_{23}-l_{21}{r}_{13}& r_{24}=Y_{24}-l_{21}{r}_{14}\\ l_{31}=& l_{32}=& r_{33}=& r_{34}=\\ Y_{31}/r_{11}& (Y_{32}-l_{31}{r}_{12})/r_{22}& Y_{33}-(l_{31}{r}_{13}+l_{32}{r}_{23})& Y_{34}-(l_{31}{r}_{14}+l_{32}{r}_{24})\\ l_{41}=& l_{42}=& l_{43}=& r_{44}=\\ Y_{41}/r_{11}& (Y_{42}-l_{41}{r}_{12})/r_{22}& [Y_{43}-(l_{41}+r_{13}+l_{42}{r}_{23})]/r_{33}& Y_{44}=(l_{41}{r}_{14}+l_{42}{r}_{24}+l_{43}{r}_{34})\end{array}\right]---\left(1\right)$

根据式(1)可以看出：

1)同一行的l_ij、r_ii、r_ij元素之间的关系类似于对对角元左侧的l_ij元素逐个消元。因此，本发明方法为了有效地利用元素的稀疏性和对称性特点，对L、R阵元素的计算不是一次完成，而是多次分步完成，简称为过程法。

2)第1列的l_i1元素由直接除以所在列的对角元r₁₁得到。所有l_ij元素在其左侧的前一个l_i,j-1元素被消元后，均要除以所在列的对角元r_jj，而其它元素则不必。

3)L、R阵各元素之间在任何情况下均有l_ji＝r_ij/r_ii的对应关系，这使得任一对角元以右的r_ij元素与对角元以下的l_ji元素虽数值不等，但非零元素的位置对称。

4)根据3)，可以用对角元以右非零的r_ij元素的位置来确定对角元以下非零的l_ji元素的位置，从而有利于LR三角分解过程稀疏性和对称性的利用。

5)根据3)，在三角分解过程中，只需计算对角元以右非零的r_ij元素，就可根据l_ji＝r_ij/r_ii关系得到对角元以下非零的l_ji元素，但同时需记录非零r_ij元素的位置。

(3)根据式(1)以及l_ji＝r_ij/r_ii的关系，将R阵中l_ji的元素用r_ij元素替代，而L阵中l_ji的元素作为映射矩阵元素存在，实际上l_ji不用计算。这样可得R阵各元素与Y阵元素的关系如式(2)。

$\left[\begin{array}{cccc}{r}_{11}=Y_{11}& r_{12}=Y_{12}& r_{13}=Y_{13}& r_{14}=Y_{14}\\ l_{21}& r_{22}=Y_{22}-r_{12}{r}_{12}/r_{11}& r_{23}=Y_{23}-r_{12}{r}_{13}/r_{11}& r_{23}=Y_{24}-r_{12}{r}_{14}/r_{11}\\ l_{31}& l_{32}& \begin{array}{c}{r}_{33}=\\ Y_{33}-(r_{13}{r}_{13}/r_{11}+r_{23}{r}_{23}/r_{22})\end{array}& \begin{array}{c}{r}_{34}=\\ Y_{34}-(r_{13}{r}_{14}/r_{11}+r_{23}{r}_{24}/r_{22})\end{array}\\ l_{41}& l_{42}& l_{43}& \begin{array}{c}{r}_{44}=\\ Y_{44}-(r_{14}{r}_{14}/r_{11}+r_{24}{r}_{24}/r_{22}+r_{34}{r}_{34}/r_{33})\end{array}\end{array}\right]---\left(2\right)$

式(2)表明，不用L阵的元素也可直接计算出R阵的所有元素。

(4)利用Y阵与L、R阵元素的相应关系和L、R阵各元素之间的对称关系，按“Г”(倒L)方式分步形成R阵各个元素，而消元计算中需用到的L阵各元素通过其与R阵各个元素的关系获得。用“Г”方式的过程法获得R阵各个元素的过程中可充分利用L阵、R阵元素的稀疏性和对称性，从而大大加快三角分解的速度。

实施过程如下：

(1)先按“Г”方式确定其对角元以右非零元素r_ij的位置，再根据对称性确定对角元以下非零元素l_ji的位置(只是认为l_ji元素存在，下同)。

(2)然后按类似高斯消元的方式分步计算“Г”所包含的R阵中的相应元素，L阵的所有元素均不计算。计算R阵元素中所需用到的L阵的元素通过l_ji＝r_ij/r_ii获得。

(3)采用对称稀疏性技术计算R阵元素，即以对角元r_ii元素为界，对r_ii左侧不为零的的l_ij元素逐个进行消元，分别计算l_ij元素所在行与不为零的r_jk元素所在列相交的、对角元r_ii及其右侧的r_ik元素。

按(1)～(3)依次循环，计算所有R阵元素。

步骤4：根据RZ_k＝E_k应用R阵元素的稀疏性回代求Z_k阵对角元Z_kk及以上的非对角元素；

(1)本发明方法对方程LW_k＝E_k的求解可完全省略。

本发明方法对Z_k阵仅需求取其对角元Z_kk及以上的元素，因此对方程LW_k＝E_k也仅需求取W_k阵对角元及其以上的元素，对应的也仅需使用E_k阵第k行对角元及以上的元素。由于E_k阵元素除第k行为1外，其余全部为零。假设W_k阵对角元及其以上的元素用W″_k阵表示，此时有W″_k＝E_k成立，即所求得的W″_k阵对角元及其以上的元素与E_k阵对角元及其以上的元素完全相同，此时方程LW_k＝E_k可变换成方程LW″_k＝E_k。

由于W″_k＝E_k成立，因此本发明方法中对方程LW_k＝E_k的求解可完全省略。

(2)利用R阵元素的稀疏性，对方程RZ_k＝W_k仅求解对角元及其以上的元素。

本发明方法求解Z_k阵的计算顺序为：Z_n,…,Z_k,…,Z₁，且在计算各个Z_k阵中，仅计算对角元Z_kk及其以上的元素。由于W″_k＝E_k成立，因此对方程RZ_k＝W_k的求解可直接变换成对方程RZ_k＝E_k对角元及其以上的元素的求取。由于R阵元素中存在大量的非零元素，因此利用R阵元素的稀疏性可大大加快计算速度。

根据上述分析可知，由于E_k阵元素结构的特点以及本发明方法中对Z_k阵仅求取对角元及其以上的元素，对方程LW_k＝E_k的求解可完全省略。因此本发明方法的回代过程用R阵元素就可完成，完全可省略L阵元素的计算。

步骤5：根据对称性求对角元Z_kk以左的非对角元素；

由于Z阵元素对称，Z_k阵对角元以左的元素均可按对称性获得。因此仅求取对角元及其以上元素的计算方式具有较高的计算效率，这样可减少Z阵50％非对角元素的计算。

步骤6：将Z阵写入数据文件以备后续程序使用。

考虑到程序的结构化，所形成的Z阵数据文件可由下一个程序调用执行。

本发明方法具有以下优点。

(1)采用过程法计算三角分解中各个元素，可根据元素本身的对称性和稀疏性进行计算，大大减少三角分解过程元素的计算量。

(2)仅求上R阵元素，不计算L阵元素，又大大减少了三角分解过程元素的计算量，记录三角分解过程中R阵非零元素位置以便回代过程中使用。

(3)回代过程利用单位矩阵E阵结构特点省略W阵元素的计算，可直接根据方程RZ_k＝E_k，求取Z_k阵对角元Z_kk及以上元素，并根据对称性求Z_kk以左元素。

(4)回代过程利用R阵元素的稀疏性同样可大大减少回代过程元素的计算量。

附图说明

图1为传统的不考虑元素稀疏性和对称性的LR三角分解法求取Z阵元素流程图。

图2为本发明方法求取Z阵元素流程图。

图3为本发明方法求取Z阵元素计算流程图。

具体实施方式

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

实施例1。

以式2为例，说明本发明方法中元素定义及其对称稀疏性技术的应用。

假设l₃₁≠0，则对r₁₁元素之下的l₃₁元素要进行消元。先定义r₁₁元素为对角元素，r₁₁右侧的所有r_1j元素均定义为交叉元素，l₃₁元素定义为消元元素，则l₃₁右侧的所有元素均定义为计算元素。

(1)消元过程中若不考虑元素稀疏性，则需分步计算l₃₁元素右侧的所有计算元素l₃₂、r₃₃、r₃₄。

(2)若采用稀疏性技术，则仅分步计算l₃₁元素右侧其所在行与不为零的r_1j元素所在列相交的所有计算元素l₃₂、r₃₃、r₃₄(需考虑对角元以左的l₃₂元素)。

(3)若采用对称稀疏性技术，则仅需分步计算l₃₁元素所在行与不为零的r_1j元素所在列相交的、对角元及以右的所有计算元素r₃₃、r₃₄(无需考虑对角元以左的l₃₂元素)。

根据实施例1可以看出：采用稀疏性技术与传统的LR三角分解法相比，分解过程中元素的计算量大大降低；若采用对称稀疏性技术，则还能在采用稀疏性技术LR三角分解法的基础上减少50％非零元素的计算。且对称稀疏性技术，特别适合于仅计算R阵元素的三角分解过程。

实施例2。

以n×n阶节点系统为例，分别比较传统的LR三角分解法和本发明方法求解Z阵元素过程的不同。比较结果如表1所示。

表1传统的LR三角分解法和本发明方法求解Z阵元素过程的比较

(1)传统LR三角分解法按公式法求解L、R阵的全部元素，无法利用对称稀疏矩阵技术；本发明方法按过程法仅求解R阵的非零元素，利用对称稀疏矩阵技术得L阵的非零元素，大大减少了三角分解过程中元素的计算量。

(2)传统LR三角分解法求取Z阵元素的过程是将一列列Z_k阵元素全部求解出来。因此需求解各有n个方程的2个方程组，每个方程均要求解n个变量，其变量元素总数为2n²个，中间矩阵变量1个。

(3)本发明方法求解Z阵元素的过程仅需求解Z_k阵对角元Z_kk及其以上的元素。由于利用E_k阵元素的特点，因此仅计算W_k阵对角元及其以上元素的算法使得求出的W″_k阵与E_k阵对角元及其以上的元素完全相同，即W″_k＝E_k成立。因此对方程RZ_k＝W″_k的求解可直接转化成对方程RZ_k＝E_k中对角元及其以上的元素的求解。

所以本发明中对方程LW_k＝E_k的求解完全可省去，该过程所需计算元素个数为0。求解方程RZ_k＝E_k过程所需计算元素个数为n(n+1)/2≈n²/2。即使根据对称性得到Z阵的下三角元素，本发明方法计算变量的总数为n²。

因此，比较三角分解过程和回代过程可以发现本发明方法的元素计算量和计算过程大大简化。

实施例3。

分别用传统的LR三角分解法(图1)以及本发明方法1、2(图2、3)对IEEE-30、-57、-118节点系统的Y阵求其Z阵元素，并比较其“分解”和“分解+回代”过程的平均计算时间。计算结果如表2所示。

表2稀疏性及对称性对LR三角分解法“分解”和“分解+回代”过程时间的影响

(1)“分解”过程平均迭代时间：

T₁：传统方法不考虑稀疏性且不利用L阵和R阵元素的关系特性、计算所有元素；

T₂：本发明方法1考虑对称稀疏性、仅计算R阵元素；

T₃：本发明方法2考虑对称稀疏性、仅计算R阵元素、全部记录R阵非零元素位置。

(2)“回代”过程平均迭代时间：

T′₁：传统方法按整列回代；

T′₂：本发明方法1按对称回代(求Z_k阵对角元Z_kk及以上元素、按对称性求对角元Z_kk以左元素)、考虑E阵特殊性；

T′₃：本发明方法2按对称回代、考虑E阵特殊性、考虑R阵元素稀疏性。

(3)“分解+回代”过程的平均迭代时间：

T″₁：传统方法分解过程不考虑对称稀疏性、不利用L阵和R阵元素的关系特性；回代过程按整列回代；

T″₂：本发明方法1分解过程考虑对称稀疏性、仅计算R阵元素；回代过程按对称回代、考虑E阵特殊性；

T″₃：本发明方法2分解过程考虑对称稀疏性、仅计算R阵元素、全部记录R阵非零元素位置；回代过程按对称回代、考虑E阵特殊性、考虑R阵元素稀疏性。

根据表2可以看出：

无论是“分解”过程还是“分解+回代”过程，本发明方法的计算速度大大优于传统LR三角分解法的计算速度。

以IEEE-118节点系统为例，结果比较如下：

(1)“分解”过程中，本发明方法考虑对称稀疏性、仅计算R阵元素的平均计算时间仅占传统方法时间的0.98％；若考虑对称稀疏性、仅计算R阵元素、全部记录R阵非零元素位置仅占传统方法时间的1.07％。

(2)“回代”过程中，本发明方法按对称回代、考虑E阵结构特点的平均计算时间仅占传统方法时间的28.25％；若按对称回代、考虑E阵特殊性、考虑R阵元素稀疏性的平均计算时间仅占传统方法时间的2.38％。

(3)“分解+回代”过程中，本发明方法分解过程考虑对称稀疏性、仅计算R阵元素；回代过程按对称回代、考虑E阵特殊性的平均计算时间仅占传统方法时间的22.18％；分解过程考虑对称稀疏性、仅计算R阵元素、全部记录R阵非零元素位置；回代过程按对称回代、考虑E阵特殊性、考虑R阵元素稀疏性的平均计算时间仅占传统方法时间的2.08％。

上述计算结果表明，本发明提出的基于对称稀疏矩阵技术的改进LR三角分解法与传统的LR三角分解法相比，无论在三角分解过程还是回代过程均可大大加快计算速度，从而大大加快求取电力系统Z阵元素的速度。

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

Claims (1)

1.一种基于对称稀疏矩阵技术的改进LR三角分解求取电力系统节点阻抗矩阵的方法，其特征包括以下步骤：

步骤1：读入n节点系统各线路支路数据文件；

步骤2：形成节点导纳矩阵Y；

步骤3：根据稀疏性和对称性对Y阵进行LR三角分解，但仅求R阵元素，并记录其非零元素位置；

步骤4：根据RZ_k＝E_k应用R阵元素的稀疏性回代求Z_k阵对角元Z_kk及以上的非对角元素；

步骤5：根据对称性求对角元Z_kk以左的非对角元素；

步骤6：将Z阵写入数据文件。