CN104391823A - 一种基于a＝ldu三角分解法求解电力系统节点阻抗矩阵的方法 - Google Patents

Abstract

一种基于A＝LDU三角分解法求解电力系统节点阻抗矩阵的方法，属于电力系统分析计算领域。主要包括以下步骤：形成节点导纳矩阵Y，对Y阵进行LDU三角分解；对DH_k＝W_k(E_k)方程求取h_kk元素；对UZ_k＝H_k求取Z_k阵对角元Z_kk及以上的非对角元素；根据对称性求对角元Z_kk以左的非对角元素；写Z阵数据到数据文件。本发明方法主要利用了单位矩阵E阵中E_k阵的结构特点、Z_k阵元素的计算顺序以及Z阵元素的对称性，省去了W阵的计算，将对H阵的计算简化为仅对其对角元素h_kk的计算，对Z阵可省去50％非对角元素的计算，大幅度提高了Z阵元素的计算速度。用本发明方法对IEEE-57、-118、-300节点系统进行验算，与传统的LDU三角分解法相比，本发明方法计算速度可提高约35～45％。

Description

一种基于A＝LDU三角分解法求解电力系统节点阻抗矩阵的方法

技术领域

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

背景技术

节点阻抗矩阵Z在电力系统中应用十分广泛且具有重要的作用。传统的求解Z阵的方法有支路追加法、导纳矩阵Y消元求逆法、LDU三角分解法等。在传统方法中，LDU三角分解法相对而言计算速度最快，因而使用最多，其特点是利用了适于求解常系数线性方程的三角分解法，对Y阵进行LDU三角分解后，可将一个对n×n阶Z阵元素的求解分成对n个列矩阵Z_k元素的求解。

但是传统的LDU三角分解法在求解Z阵元素的过程中未考虑利用Z阵元素的对称性，因此要计算Z阵的全部元素，计算量大，计算时间长。此外，也未利用单位矩阵E中E_k阵的结构特点以及Z_k阵元素的计算顺序，对进行LDU三角分解后的三个方程都要进行求解，因此实际上其计算时间并不理想。

传统LDU三角分解法Z阵元素的计算顺序为：Z₁,…,Z_k,…,Z_n，求解过程如下：

发明内容

本发明的目的是克服现有方法的不足之处，提供一种基于A＝LDU三角分解法快速求解阻抗矩阵的方法。

本发明是通过以下技术方案实现的。

本发明包括以下步骤：

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

步骤1：读取各支路数据，形成节点导纳矩阵Y；

步骤2：对节点导纳矩阵Y进行LDU三角分解，得到L、D、U三个因子阵；

步骤2具体实施过程如下：

根据YZ_k＝E_k，令Y＝LDU，得LDUZ_k＝E_k。再将LDUZ_k＝E_k进一步分解为LW_k＝E_k，DH_k＝W_k，UZ_k＝H_k三个方程。W_k、H_k阵是仅用于中间计算的过渡变量矩阵。

步骤3：对DH_k＝W_k(E_k)方程，仅求取H_k阵中的对角元素h_kk，通过UZ_k＝H_k求取Z_k阵对角元Z_kk及以上的非对角元素，利用对称性得到对角元Z_kk以左的非对角元素；

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

(1)对方程LW_k＝E_k求解W_k阵的过程

本发明方法可省去对方程LW_k＝E_k的求解，而传统的LDU三角分解法要求解整个W_k阵。

本发明中Z_k阵元素的计算顺序为：Z_n,…,Z_k,…,Z₁，利用单位矩阵E中E_k阵元素结构的特点可以对Z_k阵仅求取其对角元Z_kk及以上的元素，因此对方程LW_k＝E_k也仅需求解W_k阵对角元及其以上的元素，对应的也仅需使用E_k阵第k行对角元及以上的元素。由于E_k阵元素结构的特点为：第k行的元素为1，其余元素全部为零。这种仅计算W_k阵对角元及其以上元素的算法，使得W″_k＝E_k成立，即所求得的W″_k阵对角元及其以上的元素与E_k阵对角元及其以上的元素完全相同。本发明方法可省去对W_k阵的计算。

(2)对方程DH_k＝W_k求解H_k阵的过程

本发明方法将对方程DH_k＝W_k求解简化成仅求取H_k阵中的对角元素h_kk，而传统的LDU三角分解法要求解整个H_k阵。

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

(3)对方程UZ_k＝H_k求解Z_k阵的过程

本发明方法仅计算Z_k阵对角元Z_kk及其以上的非对角元素，利用对称性得到对角元Z_kk以左的非对角元素，而传统的LDU三角分解法是求解每一列Z_k阵的全部元素。

本发明中利用Z阵元素的对称性，可减少50％非对角元素的计算。当然，如果需要，完全可以不求Z阵下三角的非对角元素，以进一步加快计算速度。

本发明方法中Z阵元素的计算顺序为：Z_n,…,Z_k,…,Z₁，求解过程如下，其中带上标“′”的元素表示根据对称性得到的非对角元素。

步骤4：将Z阵写入数据文件以备后续程序使用及结果验证。

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

本发明方法由于利用了E_k阵的结构特点、Z_k阵元素的计算顺序以及Z阵元素的对称性，与传统的LDU三角分解法相比主要具有以下3个优势：

(1)对第一个方程LW_k＝E_k，W阵无需计算，从而省去对该方程的求解。

(2)对第二个方程DH_k＝W_k，可将对方程DH_k＝W_k的求解简化成直接对对角元素h_kk＝1/d_kk的计算，大大简化了对该方程的计算。

(3)对第三个方程UZ_k＝H_k，仅计算Z_k阵对角元Z_kk及其以上的非对角元素，而利用对称性得到对角元Z_kk以左的非对角元素，减少了该方程50％非对角元素的计算。如果需要，也可以不求Z阵下三角的非对角元素，以进一步加快计算速度。

(4)综上所述，传统LDU三角分解法求取Z_k阵元素有3个方程组和2个中间矩阵变量W_k及H_k，每个方程的计算元素个数均为n×n。因此其计算元素的总数为3n²。

而本发明方法虽有2个方程组，但仅有1个中间矩阵变量H_k，且其计算元素个数仅为n，而最后一个方程组实际计算的元素个数为n(n+1)/2，即使加上根据对称性得到的元素个数n(n-1)/2。因此其计算元素的总数为n²+n。

因此，本发明方法的计算元素总数约为传统LDU三角分解法的1/3，大大减少了元素计算量，从而大幅提高计算速度。

附图说明

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

图2为传统法的LDU三角分解法求取Z阵元素的计算流程图

具体实施方式

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

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

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

根据表1可以看出：

(1)传统的LDU三角分解法求解Z阵元素的过程：将一列列Z_k阵元素全部求解出来。整个过程可表示为：根据YZ_k＝E_k，令Y＝LDU，得LDUZ_k＝E_k。LDUZ_k＝E_k可进一步分解为：LW_k＝E_k，DH_k＝W_k，UZ_k＝H_k三个方程组。这三个方程组各有n个方程，每个方程均要求解n个变量。因此其计算元素总数为3n²个，而中间矩阵变量有2个。

(2)本发明方法求解Z阵元素的过程：对方程LW_k＝E_k，可省去整个W阵元素的计算；对方程DH_k＝W_k，方程数为n，每次仅需求解1个简单变量；对方程UZ_k＝H_k，所需计算元素个数为n(n+1)/2，如需要Z阵下三角的非对角元素，可按对称性直接得到，元素个数为n(n-1)/2。

因此，本发明方法计算元素的总数仅占传统的LDU三角分解法计算元素总数的1/3左右，而且少求解2个方程组和2个中间矩阵变量，只要求n个简单中间变量。

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

表2  两种方法求解IEEE各个系统的Z阵元素在“分解”、“回代”及“分解+回代”过程中平均计算时间的比较

T₁、T₂：传统LDU三角分解法、本发明方法在“回代”过程的平均计算时间；

T′₁、T′₂：传统LDU三角分解法、本发明方法在“分解+回代”过程的平均计算时间。

T₂/T₁：本发明方法与传统LDU三角分解法“回代”过程计算时间的比值。

T′₂/T′₁：本发明方法与传统LDU三角分解法“分解+回代”过程计算时间的比值。

从表2可以看出，本发明方法与传统LDU三角分解法相比有以下几点不同：

(1)本发明方法的“回代”过程计算时间减少了约65％；

(2)本发明方法的“分解+回代”过程计算时间减少了约50％；

(3)本发明方法的优势随着系统节点数的增加而略有增加。

以上几点足以表明本发明方法在求取Z阵元素时与传统LDU三角分解法相比，在计算速度上具有极大的优势。

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

Claims (1)

1.一种基于A＝LDU三角分解法求解电力系统节点阻抗矩阵的方法，其特征是包括以下步骤：

步骤1：读取各支路数据，形成节点导纳矩阵Y；

步骤2：对节点导纳矩阵Y进行LDU三角分解，得到L、D、U三个因子阵：

根据YZ_k＝E_k，令Y＝LDU，得LDUZ_k＝E_k，再将LDUZ_k＝E_k进一步分解为LW_k＝E_k，DH_k＝W_k，UZ_k＝H_k三个方程；

步骤3：对DH_k＝W_k(E_k)方程，仅求取H_k阵中的对角元素h_kk，通过UZ_k＝H_k求取Z_k阵对角元Z_kk及以上的非对角元素，利用对称性得到对角元Z_kk以左的非对角元素；

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