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

Abstract

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

Description

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

技术领域

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

背景技术

电力系统中一般都用传统的LDU三角分解法求取节点阻抗矩阵Z。但实际上，传统方法在三角分解过程中均未考虑利用元素的对称性和稀疏性，在回代过程中未考虑E阵元素结构的特点和上三角元素的稀疏性，从而使计算效率大大降低。且由于计算过程和计算变量的不同，LDU三角分解法的计算效率低于CU三角分解法，因此用CU三角分解法求取Z阵元素是更好的选择。

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

传统的CU三角分解法在回代过程中由于是按整列求取Z_k阵元素，未能利用Z阵元素的对称性特点，也未能利用单位矩阵E元素结构的特点和U阵元素的稀疏性，其回代过程包括求解二个完整的方程CW_k＝E_k和UZ_k＝W_k。因此其回代过程存在大量不必要的计算，计算效率极低。

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

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

发明内容

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

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

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

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

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

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

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

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

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

根据式(1)可以看出：

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

2)第1行的u_i1元素均除以所在行的对角元素c₁₁，其余所有u_ij元素在对角元素c_ii左侧前一个的c_i,i-1元素被消元后，也均要除以所在行的对角元素c_ii，而其它元素则不必。

3)C、U阵各元素之间在任何情况下均有u_ij＝c_ji/c_ii的对应关系，这使得任一对角元以右的u_ij元素与对角元以下的c_ji元素虽数值不等，但非零元素的位置对称。利用此特点，计算对角元以下的非零元素c_ji，就可根据u_ij＝c_ji/c_ii关系得到对角元以右的非零元素u_ij，从而大大简化对非零元素u_ij的计算。

尽管3)中的计算过程较为直观。但按照习惯在三角分解过程中，先分步计算对角元c_ii以右的非零元素u′_ij，而对c_ii以左的元素均不计算，然后按对称性得到对角元以下的非零元素c_ji后，同时需记录非零元素u′_ij位置，再对对角元以右的非零元素u′_ij除以该行对角元c_ii得到u_ij。这种计算过程有利于理解CU三角分解过程对称性和稀疏性的应用。

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

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

1)先按“Г”方式确定并记录其对角元以右非零元素u′_ij的位置，再根据对称性由元素u′_ij得到对角元以下非零元素c_ji的数值和位置，再对对角元以右的u′_ij元素除以该行对角元c_ii得到u_ij。

2)然后按规格化的高斯消元方式计算“Г”所包含的整个矩形框中上三角中C阵的对角元素和U阵的非零元素，按对称性得到下三角中C阵的非零元素。

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

按1)～3)依次循环，计算C、U阵元素。

步骤4：根据CW_k＝E_k仅求取w_kk元素；

E_k阵元素的特点为：第k行的元素为1，其余元素全部为0，且本发明方法中对Z_k阵元素仅需求解其对角元Z_kk及以上的元素。因此，对方程CW_k＝E_k的求解可简化成仅求取W_k阵中的对角元素w_kk＝1/c_kk，而不需求解整个W_k阵，对应的也仅需使用E_k阵第k行对角元及以上的元素。此时，w_kk以上的元素仍然全部为0，w_kk以下的元素则不用考虑。

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

本发明方法对方程UZ_k＝W_k求解Z_k阵元素的计算顺序为：Z_n,…,Z_k,…,Z₁，且在计算各个Z_k阵过程中，仅计算Z_k阵对角元Z_kk及其以上的非对角元素，即Z_kk,Z_k,k-1,…,Z_k1。

由于U阵元素中存在大量的非零元素，利用E_k阵元素结构的特点以及U阵元素的稀疏性、对Z_k阵仅求取对角元及其以上的元素可大大加快回代计算速度。

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

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

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

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

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

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

(2)仅计算少量的U阵非零元素和C阵的对角元素c_ii，通过对称性得到C阵非零的非对角元素，并记录U阵非零元素位置以便回代过程中使用。

(3)回代过程利用单位矩阵E阵结构特点可将W_k阵的求解转换成仅对其对角元素w_kk＝1/c_kk的计算，大大减少回代过程元素的计算量。

(4)回代过程利用U阵元素的稀疏性、仅求取Z_k阵对角元Z_kk及以上元素、并根据对称性求Z_kk以左元素，也可大大减少回代过程元素的计算量。

附图说明

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

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

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

具体实施方式

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

实施例1。

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

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

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

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

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

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

实施例2。

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

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

根据表1可以看出：

(1)按传统的公式法形成CU的因子阵时，需求取C、U阵的全部元素，无法利用元素自身的对称性和稀疏性；本发明方法的过程法仅求取少量U阵的非零元素和c_ii元素，可利用对称性得C阵的非零元素，大大减少了三角分解过程中元素的计算量。

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

(3)本发明方法求解Z阵元素的过程仅需求解Z_k阵对角元Z_kk及其以上的元素。由于利用E_k阵元素的特点，对方程CW_k＝E_k的求解转换成计算对n个简单中间变量w_kk＝1/c_kk的计算。此外，还可利用U阵元素的稀疏性进行回代。

求解方程CW_k＝E_k所需计算元素个数为n，求解方程UZ_k＝W_k需计算变量个数为n(n+1)/2≈n²/2。

因此，上述比较可以发现本发明方法与传统方法相比，元素的计算量和计算过程大大简化。

实施例3。

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

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

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

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

T₂：本发明方法考虑对称稀疏性、利用C阵和U阵元素的关系、记录U阵非零元素位置；

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

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

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

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

T″₁：传统方法分解过程不考虑对称稀疏性、不利用C阵和U阵关系、计算所有元素；回代过程按整列回代；

T″₂：本发明方法分解过程考虑对称稀疏性、利用了C阵和U阵关系、记录U阵非零元素位置；回代过程按对称回代、考虑E阵特殊性、利用U阵元素稀疏性。

根据表2可以看出：

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

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

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

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

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

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

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

Claims (1)

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

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

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

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

步骤4：根据CW_k＝E_k仅求取w_kk元素；

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

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

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