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

Abstract

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

Description

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

技术领域

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

背景技术

电力系统中一般都用传统的LDU三角分解法求取节点阻抗矩阵Z。但实际上，传统方法在三角分解过程中未考虑利用L、U阵元素的对称性和稀疏性，在回代过程中未考虑利用E阵元素结构的特点和U阵元素的稀疏性，从而使计算效率大大降低。

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

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

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

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

发明内容

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

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

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

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

步骤3：根据对称性和稀疏性对Y阵进行LDU三角分解，并记录U阵非零元素位置；

步骤3中具体实施过程如下：

(1)本发明方法采用过程法对L、D、U阵元素分步形成，并利用L、U阵元素的对应关系，仅计算U阵元素就可按对称性获得L阵元素。

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

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

根据式(1)可以看出：

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

2)第i行对角元d_ii元素以右的u_ij元素和第i列对角元以下的l_ji元素有u_ij＝l_ji的关系，这表明u_ij元素和l_ji元素不但非零元素的位置对称，而且数值也完全相等。因此可根据计算出的非零元素u_ij，按对称性直接得到相对应的非零元素l_ji。

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

(4)“Г”方式过程法的计算过程如下：

1)先按“Г”方式确定其对角元以右非零元素u_ij的数值及位置，再根据对称性确定对角元以下非零元素素l_ji的数值及位置。

2)然后按类似高斯消元的方式分步计算“Г”所包含的整个矩形框中相应的L、D、U阵元素。

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

按1)～3)依次循环，计算所有D、U阵元素，并分别按对称性得L阵元素。

步骤4：对方程DH_k＝W_k＝E_k仅求取h_kk元素；

步骤4中具体实施过程如下：

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

本发明方法对Z_k阵仅需求取其对角元Z_kk及以上的元素，因此对方程LW_k＝E_k也仅需求取W_k阵对角元及其以上的元素，对应的也仅需使用E_k阵第k行对角元及以上的元素。由于E_k阵元素除第k行为1外，其余全部为0。所以除去W_k阵对角元以下元素时有W_k＝E_k成立，即所求得的W_k阵对角元及其以上的元素与E_k阵对角元及其以上的元素完全相同，因此本发明方法中对方程LW_k＝E_k的求解可完全省略。

(2)对方程DH_k＝W_k＝E_k仅求取h_kk元素。

本发明方法中同样利用E_k阵元素结构的特点以及Z_k阵元素的计算顺序，且由于W_k＝E_k，可将对方程DH_k＝W_k的求解转换成对方程DH_k＝E_k的求解，而将对H_k阵的求解简化成仅求取H_k阵中的对角元素h_kk＝1/d_kk的计算。此时，h_kk以上的元素仍然全部为0，h_kk以下的元素不用考虑。

根据上述分析可知，由于E_k阵元素结构的特点，对方程LW_k＝E_k的求解可完全省略，对方程DH_k＝E_k的求解可简化为仅求取h_kk元素。

步骤5：根据UZ_k＝H_k应用U阵元素的稀疏性回代求取Z_k阵对角元Z_kk及以上的非对角元素；

本发明方法求解Z_k阵的计算顺序为：Z_n,…,Z_k,…,Z₁，且在计算各个Z_k阵中，仅计算对角元Z_kk及其以上的元素。由于U阵元素中存在大量的非零元素，因此利用U阵元素的稀疏性可大大加快计算速度。

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

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

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

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

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

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

(2)仅计算少量的U阵非零元素，通过对称性得到L阵非零元素，大大减少了三角分解过程元素的计算，记录U阵非零元素位置以便回代过程中使用。

(3)利用单位矩阵E阵结构特点，可省略对方程LW_k＝E_k的求解；对方程DH_k＝E_k仅求取h_kk元素。

(4)对方程UZ_k＝H_k，利用U阵元素的稀疏性求取Z_k阵对角元Z_kk及以上元素，并根据对称性求Z_kk以左元素，同样大大减少回代过程元素的计算量。

附图说明

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

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

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

具体实施方式

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

实施例1。

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

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

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

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

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

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

实施例2。

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

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

根据表1可以看出：

(1)按公式法形成传统的LDU的因子阵时，需求解L、D、U阵的全部元素，无法利用元素自身的对称性和稀疏性；而本发明方法利用对称稀疏矩阵技术，只需求解少量U阵的非零元素和D阵元素，少量的L阵非零元素可根据对称性由U阵元素得到，大大减少了三角分解过程中元素的计算量。

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

(3)本发明方法求解Z阵元素的过程仅需求解Z_k阵对角元Z_kk及其以上的元素。由于利用E_k阵元素的特点，仅计算W_k阵对角元及其以上元素的算法使得求出的W_k阵与E_k阵对角元及其以上的元素完全相同，即W_k＝E_k成立，所以本发明中对方程LW_k＝E_k的求解完全可省去。

此时对方程DH_k＝W_k的求解直接转化成对方程DH_k＝E_k的求解。同样由于E_k阵元素的特点，对方程DH_k＝E_k的求解可直接转化成对对角元素h_kk＝1/d_kk简单变量的求解。

求解方程DH_k＝E_k过程所需计算元素个数为n，求解方程UZ_k＝H_k过程所需计算元素个数为n(n+1)/2≈n²/2，中间矩阵变量仅1个。

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

实施例3。

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

表2传统方法与本发明方法在“分解”、“回代”和“分解+回代”过程时间比较

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

T₁：传统方法不考虑对称稀疏性、不利用u_ij＝l_ji特性、计算所有元素；

T₂：本发明方法考虑对称稀疏性、仅计算D阵元素和少量U阵非零元素、并记录其位置、利用u_ij＝l_ji特性得L阵非零元素；

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

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

T′₂：本发明方法按对称回代、考虑E阵结构特殊性、考虑U阵元素稀疏性；

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

T″₁：传统方法分解过程不考虑对称稀疏性、不利用u_ij＝l_ji特性、计算所有元素；回代过程按整列回代；

T″₂：本发明方法分解过程考虑对称稀疏性、仅计算计算D阵元素和少量U阵非零元素并记录其位置、利用u_ij＝l_ji的特性得L阵非零元素；回代过程考虑E阵结构特殊性、考虑U阵元素稀疏性、按对称回代。

根据表2可以看出：

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

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

(1)“分解”过程中，本发明方法仅占传统方法时间的1.12％。

(2)“回代”过程中，本发明方法仅占传统方法时间的2.12％。

(3)“分解+回代”过程中，本发明方法仅占传统方法时间的1.89％。

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

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

Claims (1)

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

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

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

步骤3：根据对称性和稀疏性对Y阵进行LDU三角分解，并记录U阵非零元素位置；

步骤4：对方程DH_k＝W_k＝E_k仅求取h_kk元素；

步骤5：根据UZ_k＝H_k应用U阵元素的稀疏性回代求取Z_k阵对角元Z_kk及以上的非对角元素；

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

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