CN104714928B - 一种基于对称稀疏矩阵技术的高斯消元法求取电力系统节点阻抗矩阵的方法 - Google Patents

Abstract

一种基于对称稀疏矩阵技术的高斯消元法求取电力系统节点阻抗矩阵的方法，属于电力系统分析计算领域。主要包括以下步骤：形成节点导纳矩阵Y；将Y阵与E_n阵形成增广阵B_n＝[YE_n]；根据对称稀疏性对B_n阵消元得B_n ^(n‑1)′＝[Y^(n‑1)′E_n ^(n‑1)′]；根据Y^(n‑1)′Z_n＝E_n ^(n‑1)′、稀疏性、对称性求Z_n阵对角元Z_nn以上和以左元素；根据Y^(n‑1)′阵得Y^(k‑1)′阵；根据Y^(k‑1)′Z_k＝E_k ^(k‑1)′、稀疏性、对称性求Z_k阵对角元Z_kk以上和以左元素。本发明方法利用对称稀疏性，避免了前代过程的所有无效计算，减少了约50％非零元素的计算；利用E阵元素结构的特点以及上三角元素的稀疏性按对称方式回代求取Z_k阵元素，大大加速了回代计算。用本发明方法对IEEE‑30、‑57、‑118节点等系统进行验算，与传统的高斯消元法和LDU三角分解法相比，对IEEE‑118节点系统计算速度可提高96～97％。

Description

一种基于对称稀疏矩阵技术的高斯消元法求取电力系统节点 阻抗矩阵的方法

技术领域

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

背景技术

电力系统中一般都用传统的LDU三角分解法、个别文献也介绍无规格化的高斯消元法求取节点阻抗矩阵Z。这2种方法求解Z阵时一般是将对n*n阶Z阵的求解转换成对n个Z_k阵(Z₁～Z_n)整列元素的求解，未用到Z阵元素的对称性进行计算，即未仅计算Z_k阵中对角元Z_kk及以上的元素Z_k-1,k～Z_1k，再根据Z_k-1,k～Z_1k直接获得对角元以左的元素Z_k,k-1～Z_k1。因此计算整列Z_k阵元素的方式比用对称性计算要多计算约50％的Z_k阵元素。

LDU三角分解法适合于求解常系数方程，但计算原理及计算过程均比高斯消元法复杂。由于LDU三角分解法的计算过程中含有规格化的因子矩阵，因此求解Z阵元素时其计算效率高于不含规格化的高斯消元法。实际上，含规格化的高斯消元法比不含规格化的高斯消元法的计算速度可高约30％。因此，如果用含规格化的高斯消元法将对一个n*n阶Z阵的求解转换成对n个Z_k阵的求解，则其计算速度应优于LDU三角分解法。

此外，由于节点导纳矩阵Y阵是对称稀疏矩阵，传统方法在消元过程中未有效地利用Y阵的对称性和稀疏性，在回代过程中未考虑E阵元素结构的特点和消元后矩阵上三角元素的稀疏性，从而使计算效率更是大大降低。

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

发明内容

为了克服上述现有技术的不足，本发明提供了一种基于对称稀疏矩阵技术的高斯消元法求取电力系统节点阻抗矩阵的方法。

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

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

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

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

步骤4：根据对称稀疏性对B_n阵进行n-1次含规格化的高斯消元得B_n ^(n-1)′＝[Y^(n-1)′E_n ^(n-1)′]，并记录Y^(n-1)′阵上三角中非零元素的位置；

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

(1)用对称稀疏矩阵技术快速消元，同时记录其上三角非零元素的位置。

消元前的Y阵如左下式，进行含规格化的消元后，Y阵变成Y^(n-1)′阵如右下式。

假设Y₃₁≠0，则对Y₃₁元素要进行消元。定义：Y₁₁元素为对角元素；Y₁₁右侧的所有元素Y_1j均为交叉元素；Y₃₁元素为消元元素；Y₃₁右侧的所有元素Y_3j为计算元素。

1)根据对角元以右非零的交叉元素的值和位置用对称性获得对角元以下非零的消元元素的值和位置，并仅计算其行列交互点上的计算元素。

Y^(n-1)′阵中，每行元素在消元前后、规格化之前，其对角元以右非零的交叉元素与对角元以下非零的消元元素数值相等、位置对称。而在规格化之后，对角元以右非零的交叉元素均除以了对角元素，此时与对角元以下非零的消元元素数值不等(只相差对角元的数值)、但位置仍然对称。

因此可根据对角元以右非零的交叉元素的值和位置用对称性获得对角元以下非零的消元元素的值和位置，并利用非零的交叉元素所在行以及非零的消元元素所在列的交互点来确定需计算的计算元素。这种用非零的交叉元素根据对称性来确定非零的消元元素、再用它们行与列的交互点来确定所需计算的计算元素的方式可大大减少对元素的无效减少，加快高斯消元的计算速度。

2)按对称稀疏性，仅计算非零的消元元素所在行其对角元及其以右和非零的交叉元素所在列交互点上的计算元素。

假设Y₃₁≠0。

不考虑元素的稀疏性：要计算Y₃₁元素右侧所有的交叉元素Y₃₂、Y₃₃、…；

考虑元素的稀疏性，仅计算Y₃₁元素所在行之右、与非零的交叉元素Y_1j所在列相交的所有计算元素Y_3j(1<j<n)，这可大大减少零元素的计算；

考虑元素的对称稀疏性：仅仅计算Y₃₁元素所在行的对角元素Y₃₃及其以右、与非零的交叉元素Y_1j所在列相交的所有计算元素Y_3j(3≤j<n)，这可再减少约50％非零元素的计算。因为对角元Y₃₃以下的非零元素均可根据对称性从Y₃₃以右的非零元素中获得。

3)将规格化前对角元以右的非零元素赋值给对角元以下的相应元素

在Y₃₃的前一个元素Y₃₂消元完成后，先将对角元Y₃₃以右的非零元素Y_3j的值按对称性赋给Y₃₃以下的非零元素Y_j3，然后再对Y₃₃以右的非零元素Y_3j规格化。

4)记录上三角矩阵中非零元素的位置以便回代过程中应用上三角元素的稀疏性。

步骤5：根据Y^(n-1)′Z_n＝E_n ^(n-1)′以及上三角元素的稀疏性求Z_n阵元素；

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

步骤7：根据Y^(n-1)′阵得Y^(k-1)′阵，并根据Y^(k-1)′Z_k＝E_k ^(k-1)′以及上三角元素的稀疏性等求第n-1～1列Z_k阵的对角元Z_kk及以上的非对角元素；

由于E_k阵元素结构的特点，在对第k行元素规格化后，E_k阵的对角元从e_kk＝1变化 为

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

本发明方法对Y^(k-1)′Z_k＝E_k ^(k-1)′求解Z_k阵元素的计算顺序为：Z_n,…,Z_k,…,Z₁，且在计算各个Z_k阵过程中，仅计算Z_k阵对角元Z_kk及其以上的非对角元素，即计算Z_kk,Z_k-1,k,Z_k-2,k,…,Z_1k，再根据对称性得Z_kk以左的所有元素Z_k,k-1,Z_k,k-2,…,Z_k1。这种计算方式具有较高的计算效率。

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

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

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

(1)可根据元素本身的对称性和稀疏性进行计算，不但大大减少消元过程中对元素无效计算，而且可减少约50％非零元素的计算，从而大幅提高前代计算速度。

(2)消元过程中下三角的非零元素可通过与之对应的、规格化前的上三角的非零元素获得。

(3)回代过程利用单位矩阵E元素结构特点，可将E_k ^(k-1)′阵的求解转换成对其对角元素的计算，使得用对角元素和上三角元素就可完成相应计算，大大减少回代过程元素的计算量。

(4)回代过程利用E阵元素结构的特点以及上三角元素的稀疏性按对称方式回代求取Z_k阵对角元Z_kk及以上元素、根据对称性求Z_kk以左元素，也可大幅提高回代计算速度。

本发明对传统的含规格化的高斯消元法在前代过程中利用了Y^(n-1)′阵元素的对称稀疏性、E阵元素结构的特点，回代过程中改变了Z_k阵的求解顺序和求解方式，利用了Y^(n-1)′阵上三角元素的稀疏性等等，可大大提高计算速度。

附图说明

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

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

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

具体实施方式

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

实施例1。

不考虑元素的稀疏性、考虑元素的稀疏性以及考虑元素的对称稀疏性计算过程的比较。进行含规格化的消元后，Y阵变成如下Y^(n-1)′阵。

假设Y₃₁≠0，则对Y₃₁元素要进行消元。定义：Y₁₁元素为对角元素；Y₁₁右侧的所有元素Y_1j均为交叉元素；Y₃₁元素为消元元素；Y₃₁右侧的所有元素Y_3j为计算元素。

(1)不考虑元素的稀疏性，要计算Y₃₁元素右侧所有的计算元素Y₃₂、Y₃₃、…。

(2)考虑元素的稀疏性，仅计算Y₃₁元素所在行之右、与非零的交叉Y_1j元素所在列交互的所有计算元素Y_3j(1<j<n)。与(1)相比，大大减少了零元素计算。

(3)考虑元素的对称稀疏性，仅仅计算Y₃₁元素所在行的对角元素Y₃₃及其以右、与非零的交叉Y_1j元素所在列交互的所有计算元素Y_3j(3<j<n)。与(2)相比，再减少50％非零元素的计算。对Y₃₃以左的元素无需进行消元计算，只是在对Y₃₃左侧元素Y₃₂消元完成后，需将Y₃₃以右的非零元素按对应关系赋值给Y₃₃以下的非零元素。本发明方法按(3)进行。

实施例2。

分别用传统的无规格化的高斯消元法(图1)、LDU三角分解法(图2)以及本发明方法(图3)对IEEE-30、-57、-118节点系统的Y阵求其Z阵元素，并比较其平均计算时间。计算结果如表1所示。

表1 无规格化的高斯消元法、LDU三角分解法与本发明求解Z阵计算时间的比较

T₁：无规格化高斯消元法平均计算时间；

T₂：LDU三角分解法平均计算时间；

T₃：本发明方法平均计算时间；

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

T₃/T₁：本发明方法与无规格化高斯消元法平均计算时间百分比；

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

根据表1可以看出：

LDU三角分解法的计算时间比无规格化的高斯消元法快30％～40％，但随着节点数越多，两者计算时间的差距在减小。

而本发明方法比无规格化的高斯消元法快80％～97％，比LDU三角分解法快65％～96％，且随着节点数越多，本发明方法计算时间的优势越明显。如对IEEE-118节点系统，与传统的高斯消元法和LDU三角分解法相比，本发明方法计算速度可提高约96～97％。

上述计算结果表明，本发明提出的基于对称稀疏矩阵技术的高斯消元法求取电力系统节点阻抗矩阵的方法与传统的无规格化的高斯消元法和LDU三角分解法相比，计算速度大大加快。

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

Claims (1)

1.一种基于对称稀疏矩阵技术的高斯消元法求取电力系统节点阻抗矩阵的方法，其特征包括以下步骤：

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

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

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

步骤4：根据对称稀疏性对B_n阵进行n-1次含规格化的高斯消元得B_n ^(n-1)′＝[Y^(n-1)′ E_n ^{(n -1)}′]，并记录Y^(n-1)′阵上三角中非零元素的位置；

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

(1)判断对角元以右非零的交叉元素用对称性确定对角元以下非零的消元元素；

(2)以按列消元方式，计算非零的交叉元素和非零的消元元素交互点上、并位于各行对角元以右的计算元素；

(3)将规格化前的对角元以右的非零元素赋值给对角元以下的相应元素；

(4)记录第一次消元过程中上三角非零元素的位置以便在后续的反复回代和前代过程中应用；

步骤5：根据Y^(n-1)′Z_n＝E_n ^(n-1)′以及上三角元素的稀疏性求Z_n阵元素；

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

(1)规定Z_k阵的求取顺序按从第n～1列、Z_k阵元素的求取顺序从Z_kk～Z_1k；

(2)对方程Y^(n-1)′Z_n＝E_n ^(n-1)′中仅求取E_n阵中的对角元素

(3)根据步骤4记录的Y^(n-1)′阵上三角元素的稀疏性，对方程Y^(n-1)′Z_n＝E_n ^(n-1)′的求解，仅计算Z_n阵对角元Z_nn及其以上的非对角元素；

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

步骤7：根据Y^(n-1)′阵得Y^(k-1)′阵，并根据Y^(k-1)′Z_k＝E_k ^(k-1)′以及上三角元素的稀疏性等求第n-1～1列Z_k阵的对角元Z_kk及以上的非对角元素；

(1)根据Y^(n-1)′阵第k行及以上的元素得Y^(k-1)′阵；

(2)对方程YZ_k＝E_k中仅求取E_k阵中的对角元素

(3)根据步骤4记录的Y^(n-1)′阵第k行及以上上三角元素的稀疏性，对方程Y^(k-1)′Z_k＝E_k ^(k-1)′的求解，仅计算Z_k阵对角元Z_kk及其以上的非对角元素；

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

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