CN109284476B - 非零元素随机存放和随机对称消元求取电力系统节点阻抗的方法 - Google Patents

Abstract

一种非零元素随机存放和随机对称消元求取电力系统节点阻抗的方法，打开Y阵数据文件Y(n,d)，将数据读入Y(n,d₁)数组；将Y(n,d₁)数组与E_n阵构成增广阵B_n＝[Y(n,d₁)E_n]；对B_n阵进行n‑1次含规格化的基于对称稀疏技术的高斯消元得

规定Z阵中Z_k阵的求取顺序为第n～1列、Z_k阵元素的求取顺序为z_kk～z_1k，再根据Y(n,d₁)^(k‑1)′Z_k＝E_k ^(k‑1)′分步回代求解Z_k阵中对角元z_kk及以上的元素，并根据对称性得z_kk以左的元素；求出Z阵并输出结果。用本申请分别对IEEE‑30、‑57、‑118、‑300系统的Y阵求解Z阵，与传统高斯消元法相比，不仅大大减少存贮单元量，且数据文件的读取及消元计算速度均得到了大幅提升。

Description

非零元素随机存放和随机对称消元求取电力系统节点阻抗的方法

技术领域

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

背景技术

大型电力系统节点导纳矩阵Y的形成、存贮、读写及消元过程中，如不考虑Y阵元素的稀疏性和对称性，会导致大量零元素和对称元素的存贮以及不必要的元素计算，从而使得形成Y阵所需时间较长、存贮空间极大、读写Y阵的数据文件时间较长、对Y阵的的前代和回代计算时间较长等。如Y阵的Y(n,2n)数组形式虽简单直观，便于数据处理，但大量零元素的存在使得Y的形成、存贮、数据文件的读写等过程效率低下。尽管其计算过程简单方便，但当利用稀疏矩阵技术时对其进行高斯消元时，虽然计算速度可大大提高，然而由于大量判断语句的使用，其计算速度的提高也并不理想。。

传统的考虑元素稀疏性Y阵元素的存贮方式有按坐标存贮、按顺序存贮、按链表存贮、三角存贮法、Ellpack-Itpack存贮法、CSR存贮法和超矩阵存贮法等等。尽管这些存贮方式可以省去不少存贮单元，但其结构复杂，且对角元素与非对角元素分开存贮也使得存取过程繁琐，也不能清晰地反映元素之间的对称或对应关系，不利于对Y阵的数据处理。这些存贮方式虽然可以省去大量零元素的存贮，但未利用Y阵元素的对称性和Y阵结构的特点，不但存贮单元的节省并未达到最佳效果，而且不便于快速形成Y阵，且数据的检索、修改、计算等极为不便，特别是无法直接对这些数据进行消元计算，因而其存贮效率无法达到最佳状态，计算效率也不高。

新型非零元素的Y(n,d)存贮方式虽然较好地解决了按坐标存贮、按顺序存贮、按链表存贮方式中存在的问题，存贮效率达到了最佳状态，且数据的检索、修改、计算等也较为方便。但在形成Y(n,d)的数组时，要求其上三角元素按顺序排列，即要求形成和读取支路数据(I、J、R、X、K)时要求按节点号i<j以及j₁<j₂<j₃<j₄<j₅<j₆的方式，从而使形成的上三角元素的列号也必须按j₁<j₂<j₃<j₄<j₅<j₆顺序存放。这个要求会大大增加形成Y阵时的判断和循环，从而大大影响形成Y阵的速度，而且与实际电力系统工程计算中数据的随机形成的情况不符。此外，顺序存贮的Y(n,d)数组利用稀疏矩阵技术对其直接进行高斯消元时的问题并未得到解决，而使其计算效率也无法达到最佳。

常用求解节点阻抗矩阵Z的方法有支路追加法和Y阵求逆法。传统方法中由于在求解Z阵时根据其系数矩阵的变化认为是求解常系数方程，因此在Y阵求逆法中最常用的是求解常系数方程的LDU三角分解法，但传统LDU三角分解法求解Z阵时是求解n个整列的Z_k阵，进而求得Z阵，即Z_k阵的计算顺序为：Z₁,Z₂,┄,Z_k,┄,Z_n-1,Z_n。尽管Z阵元素的对称性众所周知，这种计算方式却很难利用Z阵元素的对称性。这是用LDU三角分解法计算Z阵元素效率低的一个原因。此外，LDU三角分解法中对中间矩阵的求解由于没有利用E阵元素结构的特点，而使得其增加了大量不必要的计算，从而影响计算速度的提高。虽然几乎没有介绍用求解变系数方程的高斯消元法求解Z阵元素，但实际上，如果能充分考虑Z_k阵的计算顺序和Z_k阵元素的计算顺序，并结合E_k阵元素结构的特点，从而可利用求解变系数方程的方法用于求解常系数方程的Z阵。实际上，高斯消元法比传统LDU三角分解法的计算速度快约30％，从而大大提高求取Z阵的计算速度。

此外，无论用传统的LDU三角分解法还是用传统的高斯消元法求解Z阵时其稀疏矩阵技术的应用并不方便。对LDU三角分解法还必须改变其元素前代计算过程，即将传统方法中利用元素公式每次计算一个完整元素的计算方式改成类似消元过程的元素的分步计算，并需要大量的判断语句和对上三角非零元素的记录；而在其回代过程中必须根据上三角非零元素按顺序取用相应Z_k阵中的元素，这些都容易导致其计算速度的提高很难达到最佳状态，因此稀疏矩阵技术的应用并未达到最佳。对高斯消元法无论在其前代过程中还是在其回代过程中，稀疏矩阵技术的应用也有类似的问题。而且无论在进行LDU三角分解法的过程、高斯消元的前代过程，还是各种方法的回代过程，其计算必须按顺序完成，从而使得计算过程的简化受到一定的限制。

发明内容

为了克服上述现有技术的不足，本发明提出一种非零元素随机存放和随机对称消元求取电力系统节点阻抗的方法，可直接根据非零元素的Y(n,d₁)随机存贮方式直接快速进行高斯消元求解Z阵，较好地解决了形成非零元素Y(n,d)存贮方式数据的随机读取，Y(n,d₁)数组的随机形成、随机存贮和随机消元等问题，以便在同一数组中完成“排零存贮”和“排零运算”的稀疏矩阵技术。较好地解决了高斯消元求解Z阵时元素的对称性应用问题，并省去了原稀疏矩阵技术中记录上三角非零元素的环节，可直接按随机顺序在回代过程完成元素的计算。

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

本发明所述的一种非零元素随机存放和随机对称消元求取电力系统节点阻抗的方法，步骤如下：

步骤1：打开Y阵数据文件Y(n,d)，将数据读入Y(n,d₁)数组；

(1)打开Y(n,d)数据文件之前，已经按随机顺序读取支路数据(I、J、R、X、K)，并按随机顺序形成了仅含对角元素和非零的上三角元素的Y阵数据文件Y(n,d)。Y(n,d)数据文件中除要求行号i小于列号j外，并不要求各行中的列号按j₁<j₂<j₃<j₄<j₅<j₆顺序排列，即其列号可以随机排列。Y(n,d)数据文件分为3组，第1组仅1列，存贮与对角元连接的静态的非零的非对角元数之和S_i，(不包括对角元)，S_i值由程序自动累加以保证快速地读写对应的参数，以免对Y(n,d)数据文件中多余存贮单元的读取；第2组共3列为对角元组，存贮对角元的行号和参数g_ii、b_ii；第3组共3l_max列为非零的非对角元组，其中l_max为系统中各节点上三角的最大连接支路数，而各节点上三角的连接支路数为l_i，存贮与该对角元连接的所有非零的非对角元的列号j及参数g_ij、b_ij。Y(n,d)数据文件的列数为d＝3l_max+4，其静态结构如表1所示。

表1Y(n,d)数据文件的静态结构

注：Y(n,d)数据文件元素的静态参数用g_ij、b_ij表示。

(2)打开Y(n,d)数据文件之后，将Y(n,d)数据文件的数据读入到Y(n,d₁)数组。Y(n,d₁)数组的动态结构和Y(n,d)数据文件的静态结构有以下几点不同：

1)Y(n,d)数据文件中的非零元素数S_i是Y阵中的静态非零元素数，不包括消元过程中所产生的新的非零元素数；Y(n,d₁)数组(见表2)中的动态非零元素数S′_i是Y阵消元过程中的动态非零元素数，包括了消元过程中所产生的新的非零元素数。

2)Y(n,d)数据文件中S_i仅与各节点上三角连接的支路数l_i有关，对角元组和非零的非对角元组均为3列；Y(n,d₁)数组中对角元组仍为3列，而非零的非对角元组均为5列。增加的第4～5列用以存放消元计算过程中规格化后的数据。与不增加2列而反复乘除对角元的方法相比，可使消元计算时间减少约10％。

3)Y(n,d)数据文件的列数为d＝3l_max+4，Y(n,d₁)数组的列数为d₁＝5S′_max+4，其中S′_max为Y阵在消元过程中上三角所产生的最大的动态非零元素数，且S′_max>>l_max。

4)Y(n,d)数组的静态结构不能直接用于高斯消元，而Y(n,d₁)数组的动态结构可直接用于高斯消元。将Y(n,d)数据文件的数据读入到Y(n,d₁)数组第1～3列相应的单元中，Y(n,d)数据文件的S_i也被读入到Y(n,d₁)数组中，此时有S_i＝S′_i，但在消元过程中S′_i的数值将不断变化，而Y(n,d)数据文件的S_i的数值不变。消元过程中所产生的新的非零元素被连续、随机地放到相应行紧邻的新的非零元素组的第2～3列，而规格化后的元素放在同组元素的第4～5列。Y(n,d₁)数组的动态结构见表2。

表2Y(n,d₁)数组的动态结构结构

注：Y(n,d₁)数组规格化前的动态参数用g′_ij、b′_ij表示，规格化后的用g″_ij、b″_ij表示。

步骤2：将Y(n,d₁)数组与E_n阵构成增广阵B_n＝[Y(n,d₁)E_n]；

步骤3：对B_n阵进行n-1次含规格化的基于对称稀疏技术的高斯消元得

(1)用四角规则完成结构对称的Y(n,2n)数组的消元计算

1)四角规则完成Y(n,2n)数组的消元计算。

传统高斯消元法中相应元素的计算需依赖计算公式，因此计算过程复杂繁琐，编程困难。如果找出这些元素的计算规律，就可将其变成极为简单的运算。

Y(n,2n)数组中，规则排列的第k行的对角元素、规格化前(后)的交叉元素、第k列的消元元素、及需计算的相应计算元素在计算前(后)的简化矩阵如式1所示。

式1为四角规则、非零元素判断、计算元素确定、元素对称算法原理图。

各元素定义及在简化矩阵中的位置如下：

对角元素(终值，参考元素)；

交叉元素规格化前(后，终值)，对角元素同行以右；

消元元素(终值)，对角元素同列以下；

计算元素前值和新值(可能完成或未完成计算的元素)，位于与非零的消元元素同行、与非零的交叉元素同列的交互点上。

式1表明消元计算过程中各元素的关系为：以对角元素为参考元素，“计算元素的新值＝其原值-非零的消元元素*非零的交叉元素”，即只需计算对角元以下非零的消元元素所在行和对角元以右非零的交叉元素所在列交互点上的元素。因此，根据消元过程中元素的位置，可无需消元计算公式直接用四角规则写出以下计算式。

2)Y(n,2n)数组中非零元素的快速判断和计算元素的快速确定。

对称矩阵以按列消元方式进行高斯消元计算时，其对角元以右非零元素和对角元以下非零元素的位置始终对称。因此如式1中若对角元素以右的交叉元素

则可得对角元素以下的消元元素

而根据对角元以右非零元素所在列和对角元以下非零元素所在行的交互点就可确定有效的计算元素为

四个。

3)Y(n,2n)数组中计算元素的对称算法。

对于对称矩阵，在高斯消元法用元素的对称算法可进一步减少下三角元素的计算。如对含规格化高斯消元法，式1中对角元素以右规格化前的

元素与对角元素以下的元素

完全对称，而规格化后的

元素与

是不完全对称。利用该特性，可在对第k列的

元素按列消元的过程中，仅计算其所在行对角元及其以右的元素，即仅计算式1中的

三个元素，而不用计算其所在行对角元以左的元素如

元素。在对第j-1列的元素消元完成后和第j行的元素规格化前，先将第j行的非零元素

赋给第j列

元素，然后对第j行的

元素规格化，再继续对第j列的

元素消元，依此循环。这种计算方式使得所有对角元以下的元素均通过赋值、而不是通过计算得到，从而大大简化对角元以下元素计算。

4)Y(n,2n)数组中上三角非零元素的记录。

根据系统的不同，可在Y(n,2n)数组左侧增加约n/4～n/2的列数，记录第一次消元过程中每行上三角非零元素的个数S′_i和列号，以便在后续的前代过程和对X阵的按行回代过程无需再做非零判断和记录直接应用第一次消元过程中的非零判断。S′_i是由程序自动累加，可进一步提高对称稀疏矩阵技术的应用效率。

(2)用四角规则完成结构不对称的Y(n,d₁)数组的消元计算

1)直接用四角规则完成Y(n,d₁)数组中随机存放元素的消元计算。

Y(n,2n)数组中用四角规则并考虑元素的对称算法后的要点是“以对角元作为参考元素，仅计算对角元以右非零的交叉元素和对角元以下非零的消元元素交互点上的元素”。虽然Y(n,d₁)数组中每行元素的存放顺序随机，但Y(n,d₁)数组中非零的非对角元素的存贮方式是以组为单位，每组第1列存放列号；第2～3列存放是与下三角元素完全相等的规格化前的参数；第4～5列是存放规格化后的上三角元素的参数。因此对Y(n,d₁)数组中非零的非对角元进行消元，就是对非零的非对角元组中第2～3列存放的规格化前的参数进行消元。所以，可以对角元为参考元素，将第k行规格化前每个参数的列号i与行号k分别互换，即将

并将各个y_ik分别作为第k列、第i行的消元元素，而将第k行各个规格化后的参数作为交叉元素，就可按以对角元素为参考元素，“计算元素的新值＝其原值-消元元素*交叉元素”的四角规则直接完成Y(n,d₁)数组中随机存放元素的消元计算。

2)Y(n,d₁)数组中非零元素的快速判断和计算元素的快速确定。

A.Y(n,d₁)数组中非零元素的存放方式决定了在Y(n,d₁)数组中无需进行非零元素的判断，而对于计算元素的快速确定，只要将第k行各个规格化前每个参数的列号i与行号k分别互换，即将

并分别计算第k列、第i行的各个y_ik与第k行各个规格化后参数的列号在交互点上的元素(如第k行的列号分别为i、j、p、m，包括y_ki的列号)y′_ii、y′_ij、y′_ip、y′_im。

B.如果消元计算过程中没有产生新的非零计算元素，则将计算元素的新值直接替换原第2～3列中计算元素的值；如果产生新的非零计算元素，则在对应行紧邻右侧的新的非零元素组中第1～3列分别加入新产生的非零元素的列号和参数。

在Y(n,2n)数组中，假设第k行规格化前的元素为

规格化后的元素为

其对应第k列的非零的消元元素为

由于只要计算

所在行与

元素所在列交互点上的元素，而

元素的行号实际上就是Y(n,d₁)数组中规格化前的元素

的列号。因此在Y(n,d₁)数组中就是将每行每个规格化前元素的列号j、p分别与行号k互换，分别计算j、p行与第k行每个规格化后的元素的列号在交互点上的元素，就可同时完成非零元素的快速判断和计算元素的快速确定。

3)Y(n,d₁)数组中计算元素的对称算法。

根据Y阵的对称性可知，Y(n,2n)数组中其对角元以右非零的非对角元素规格化前与对角元以下非零的非对角元素完全相等，因此在消元过程中，可完全不用计算各行对角元以左(或对角元以下)的元素，而可以在对某行元素规格化之前通过对称性先得到对角元以下的元素，再对该行元素规格化，再对该列元素消元计算。

由于在Y(n,d₁)数组中将每个非零的非对角元组分为5列，第1列为其列号，第2～3列为规格化前的参数，第4～5列为规格化后的参数。因此可分别将第k行各个规格化前元素的列号与行号互换，再将其列号)分别与第k行各个规格化后元素的列号进行比较，并按下述三种情况处理：

A.如果规格化前元素的列号i(→行号i)大于规格化后元素的列号j，则说明该元素为下三角非零的非对角元素，完全可不用计算；

B.如果规格化前元素的列号i(→行号i)等于规格化后元素的列号j，则说明该元素为对角元素，需将其计算后的值y′_ii＝g′_ii+jb′_ii直接替换第i行原对角元的值y_ii＝g_ii+jb_ii；

C.如果规格化前元素的列号i(→行号i)小于规格化后元素的列号j，则说明该元素为上三角元素，需将其列号j分别与第i行元素的各个列号进行比较，如果有相同的列号，则将其计算后的值y′_ij＝g′_ij+jb′_ij直接替换第i行该组元素的值y_ij＝g_ij+jb_ij；如果没有相同的列号，则将其列号j及计算后的值y′_ij＝g′_ij+jb′_ij直接放在紧邻的新的非对角元组的第1～3列(列号及参数)。

Y(n,2n)数组中尽管也不用计算下三角的元素，但必须在对对应的元素规格化之前将其赋值给相应的下三角元素，或者要反复进行对角元的乘除计算。而Y(n,d₁)数组中只是将规格化前的参数和规格化后的参数分开存放，无需赋值可直接完成消元计算。

4)Y(n,d₁)数组中计算元素的存放方式。

Y(n,2n)数组中，在计算非零元素交互点上的元素时，其交互点上无论非零元素的新值或新产生的非零元素均存放在其交互点上。而在Y(n,d₁)数组中，在计算非零元素交互点上的元素时，非零元素的新值或新产生的非零元素均存放在非零元素组的第2～3列(包括对角元)，在对该行元素进行规格化后，其第4～5列才会出现相应的元素。

步骤4：规定Z阵中Z_k阵的求取顺序为第n～1列以及Z_k阵元素的求取顺序为z_kk～z_1k，考虑E_k阵元素结构的特点，再根据Y(n,d₁)^(k-1)′Z_k＝E_k ^(k-1)′分步回代求解Z_k阵中对角元z_kk及以上的元素，并根据对称性得z_kk以左的元素；

(1)传统方法中求取Z_k阵的不足之处

1)由Y(n,2n)数组表示的Y阵分别与E_k阵构成增广阵B_k＝[Y E_k]。由于传统方法中没有利用Y^(k-1)′阵第k行及其以上的元素与Y^(n-1)′阵第k行及其以上的元素完全相同的特点，因此需对B_k阵分别进行k-1次含规格化的基于对称稀疏技术的高斯消元得B_k ^(k-1)′＝[Y^(k-1)′E_k ^(k-1)′]，分别记录其Y^(k-1)′阵上三角所有非零元素的位置；且由于在根据E_k阵求取E_k ^(k-1)′阵的过程中未考虑E_k阵元素结构的特点，因此需求取E_k ^(k-1)′阵的全部元素。

2)根据Y^(k-1)′Z_k＝E_k ^(k-1)′求解Z_k阵时，Z_k阵的求取顺序从第1～n列，且未考虑E_k阵元素结构的特点，Z_k阵元素的求取顺序为z_1k～z_nk，即求取了Z_k阵的所有元素，因此很难应用Z阵元素的对称性进行求解。

3)Z_k阵中对角元z_kk和非对角元z_jk的回代计算式分别如下：

因此，在求解Y^(k-1)′Z_k＝E_k ^(k-1)′方程时要引入稀疏矩阵技术，就要根据前代过程所“记录的Y^(k-1)′阵上三角所有非零元素的位置”，而传统方法中记录非零元素位置的方法并不理想，如对角元和非对角元分开记录、没有每行的非零元素个数等等，使得稀疏矩阵技术的应用效果也不理想。

(2)本申请方法中具体实施过程如下：

1)对Y(n,d₁)数组与E_k阵构成增广阵B_k＝[Y(n,d₁)E_k]进行k-1次含规格化的基于对称稀疏技术的高斯消元得B_k ^(k-1)′＝[Y(n,d₁)^(k-1)′E_k ^(k-1)′]。由于Y(n,d₁)^(k-1)′数组结构仅存放非零元素，新产生的非零元素只是按顺序随机排放在相邻的新的非对角元组中，因此其前代和回代过程中均不用再判断并记录Y(n,d₁)^(k-1)′数组中上三角非零元素的位置；由于Y(n,d₁)^(k-1)′数组第k行及其以上的元素与Y(n,d₁)^(n-1)′数组中第k行及其以上的元素完全相同，因此在后续的前代计算中可直接用Y(n,d₁)^(n-1)′阵第k行及以上的元素代替Y(n,d₁)^(k-1)′阵与E_k ^(k-1)′阵求解Z_k阵；由于求取Z_k阵元素时将仅求取对角元及以上的元素z_kk～z_1k，根据E_k阵元素结构的特点，在用E_k ^(k-1)′阵求解Z_k阵时，E_k ^(k-1)′阵与E_k阵对角元以上元素完全相同，而Z_k阵对角元以下元素已经根据对称性得到，因此E_k ^(k-1)′阵对角元以下的元素没有作用，从而无需求解。因此对E_k阵求取E_k ^(k-1)′阵也仅需求取E_k ^(k-1)′阵的对角元素，此时E_k阵的对角元素从e_kk＝1变成E_k ^(k-1)′阵的对角元素

2)由于Y^(k-1)′Z_k＝E_k ^(k-1)′求解Z_k阵时，Z_k阵的求取顺序从第n～1列、Z_k阵元素的求取顺序从z_kk～z_1k，仅求取Z_k阵对角元及以上元素，而对角元以下元素可根据Z阵元素的对称性直接得到。

3)Z_k阵中对角元z_kk和非对角元z_jk的回代计算式分别如下：

其中，l_ki和l_ji分别为第k行和第k-1～1行存放的非零的非对角元素列号。

因此，在求解Y(n,d₁)^(k-1)′Z_k＝E_k ^(k-1)′方程时可直接根据Y(n,d₁)^(k-1)′数组结构应用稀疏矩阵技术，而其对角元与非对角元的同数组存放以及每行非零的非对角元个数的存在使得稀疏矩阵技术的应用效果极为理想。注意：Y(n,d₁)^(k-1)′数组中仅非对角元组的第4～5列元素可用于对方程Y(n,d₁)^(k-1)′Z_k＝E_k ^(k-1)′进行回代。

4)根据对称性得z_kk以左的元素。

步骤5：求出Z阵并输出结果。

本申请方法提出一种基于对称稀疏矩阵的非零元素随机存放和随机对称消元求取节点阻抗矩阵的方法，主要特点如下：

(1)建立了按随机顺序读取支路数据、按随机顺序存放Y阵上三角元素的列号及参数的Y(n,d)的数据文件，用S_i控制其上三角静态的非零的非对角元数，以保证快速准确地读取相关数据。

(2)对应Y(n,d)数据文件，建立了可随机进行对称性高斯消元的Y(n,d₁)数组，用S′_i控制Y阵消元过程中其上三角变化的及新产生的动态非零元素数，以保证快速消元。

(3)直接用四角规则在结构不对称的Y(n,d₁)数组中完成了基于稀疏矩阵技术的对称高斯消元计算，其中包括非零元素的快速判断、计算元素的快速确定、计算元素的对称算法等。

(4)规定Z阵中Z_k阵的求取顺序为第n～1列以及Z_k阵元素的求取顺序为z_kk～z_1k，并考虑E_k阵元素结构的特点，再根据Y(n,d₁)^(k-1)′Z_k＝E_k ^(k-1)′可直接利用Y(n,d₁)^(n-1)′阵第k行及以上的元素代替Y(n,d₁)^(k-1)′阵，与E_k ^(k-1)′阵一起分步回代求解所有Z_k阵中对角元z_kk及以上的元素，并将对E_k ^(k-1)′阵的求解简化成仅仅对其对角元素

的求解，再分步根据对称性得z_kk以左的元素。

用本申请方法分别对IEEE-30、-57、-118、-300系统的Y阵求解Z阵，与传统高斯消元法相比，不仅大大减少了所需的存贮单元量，且在数据文件的读取及计算速度方面均得到了大幅的提升。

附图说明

图1为不考虑稀疏和对称性的传统高斯消元法求取Z阵计算流程图。

图2为本申请方法求取Z阵计算流程图。

具体实施方式

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

实施例1。

以5阶矩阵为例，对Y(n,2n)和Y(n,d₁)数组高斯消元步骤进行比较，表3为该5阶矩阵Y(n,2n)和Y(n,d₁)数组的初始状态。

(1)Y(n,2n)和Y(n,d₁)数组上三角的初始状态

Y(n,2n)数组中是顺序排列，而Y(n,d₁)数组中是随机排列(下同)；Y(n,d₁)数组中第一行非零的非对角元数为S′₁＝3，相关的节点的列号按生成顺序为l₁＝[4,5,2]；第二行S′₂＝2，l₂＝[5,3]；第三行S′₃＝1，l₃＝[4]；第四行S′₄＝1，l₄＝[5]；第五行S′₅＝0。为简化起见，令y_ij＝g_ij+jb_ij；Y(n,d₁)数组中各组的第2列为规格化前的元素y′_ij，第3列为规格化后的元素y″_ij。Y(n,2n)和Y(n,d₁)数组元素的状态见表3。

表3 5阶矩阵Y(n,2n)和Y(n,d₁)数组元素的状态

(2)对第1行规格化和第1列消元后，Y(n,2n)⁽¹⁾′、Y(n,d₁)⁽¹⁾′数组的状态

1)对第1行元素规格化，Y(n,d₁)数组中实际规格化元素的顺序为y₁₄、y₁₅、y₁₂，所得元素y⁽¹⁾ ₁₄、y⁽¹⁾ ₁₅、y⁽¹⁾ ₁₂分别放在y₁₄、y₁₅、y₁₂元素的右侧；对第1列元素消元，实际消元顺序为第1行规格化前的非对角元素y₁₄、y₁₅、y₁₂的对称元素y₄₁、y₅₁、y₂₁。因此，对第1列元素消元实际上是对第1行规格化前的非对角元素进行消元。见表4。

2)对y₁₄(y₄₁)消元，如考虑元素的对称算法，将其列号4与该行所有的规格化后元素的列号4、5、2进行比较，仅需计算规格化后元素的列号大于和等于列号4的对应元素，并将列号4作为行号，而规格化后元素的列号作为列号。因此对y₁₄消元，仅需计算y⁽¹⁾ ₄₄、y⁽¹⁾ ₄₅，即对角元及以右的元素；而无需计算规格化后元素的列号小于列号4的对应元素，即无需计算y⁽¹⁾ ₄₂，也就是对角元以左的元素。同理，对y₁₅(y₅₁)消元，仅需计算规格化后元素的列号大于和等于列号5的对应元素y⁽¹⁾ ₅₅；无需计算规格化后元素的列号小于列号5的对应元素y⁽¹⁾ ₅₂、y⁽¹⁾ ₅₄。对y₂₁消元，仅需计算规格化后元素的列号大于和等于列号2的对应元素，y⁽¹⁾ ₂₂、y⁽¹⁾ ₂₄、y⁽¹⁾ ₂₅。计算过程见表4。

上述计算过程表明，在计算非零元素交互点上的元素时，无论是非零元素的新值还是新产生的非零元素均存放在非零元素组的第2～3列，而不会放在第4～5列。仅是规格化后的元素才会放在第4～5列。

表4第1行规格化、第1列消元后的Y(n,2n)⁽¹⁾′和Y(n,d₁)⁽¹⁾′数组的状态

(3)对第2行规格化、第2列消元后的Y(n,2n)⁽²⁾′和Y(n,d₁)⁽²⁾′数组的状态

1)对第2行元素规格化，Y(n,d₁)数组中实际规格化元素的顺序为y⁽¹⁾ ₂₅、y₂₃、y⁽¹⁾ ₂₄，所得元素y⁽²⁾ ₂₅、y⁽¹⁾ ₂₃、y⁽²⁾ ₂₄分别放在y⁽¹⁾ ₂₅、y₂₃、y⁽¹⁾ ₂₄元素的右侧；对第2列元素消元，实际消元顺序为第2行规格化前的元素y⁽¹⁾ ₂₅、y₂₃、y⁽¹⁾ ₂₄的对称元素y⁽¹⁾ ₅₂、y₃₂、y⁽¹⁾ ₄₂。因此对第2列元素消元也就是对第2行规格化前的非对角元素进行消元。见表5。

2)对y⁽¹⁾ ₂₅(y⁽¹⁾ ₅₂)消元，仅需计算y⁽²⁾ ₅₅，无需计算y⁽¹⁾ ₅₃、y⁽²⁾ ₅₄，此时Y(n,2n)数组中y⁽¹⁾ ₅₂的位置上的实际元素虽仍然是y₅₂，但在对第2行元素规格化之前，需先将第2行的非零元素赋值给第二列对应元素。而在Y(n,d₁)数组中是直接对y⁽¹⁾ ₂₅进行消元，无需赋值；对y₂₃(y₃₂)消元，需计算y⁽¹⁾ ₃₃、y⁽¹⁾ ₃₄、y⁽¹⁾ ₃₅；对y⁽¹⁾ ₂₄(y⁽¹⁾ ₄₂)消元，仅需计算y⁽²⁾ ₄₄、y⁽²⁾ ₄₅，无需计算y⁽¹⁾ ₄₃。见表5。

表5第2行规格化、第2列消元后的Y(n,2n)⁽²⁾′和Y(n,d₁)⁽²⁾′数组的状态

(4)对第3行规格化、第3列消元后的Y(n,2n)⁽³⁾′和Y(n,d₁)⁽³⁾′数组的状态

1)对第3行元素规格化，Y(n,d₁)数组中规格化的顺序为y⁽¹⁾ ₃₄、y⁽¹⁾ ₃₅，所得元素y⁽²⁾ ₃₄、y⁽²⁾ ₃₅分别放在y⁽¹⁾ ₃₄、y⁽¹⁾ ₃₅元素的右侧；对第3列元素消元，实际消元顺序为第3行规格化前的元素y⁽¹⁾ ₃₄、y⁽¹⁾ ₃₅的对称元素y⁽¹⁾ ₄₃、y⁽¹⁾ ₅₃。因此对第3列元素消元实际上也就是对第3行规格化前的元素进行消元。见表6。

2)对y⁽¹⁾ ₃₄(y⁽¹⁾ ₄₃)消元，需计算y⁽²⁾ ₄₄、y⁽²⁾ ₄₅，此时Y(n,2n)数组中y⁽¹⁾ ₄₃的位置上的实际元素虽仍然是y₄₃，但在对第3行元素规格化之前，需先将第3行的非零元素赋值给第3列对应元素。而在Y(n,d₁)数组中是直接对对y⁽¹⁾ ₃₄进行消元，而无需赋值；对y⁽¹⁾ ₃₅(y⁽¹⁾ ₅₃)消元，仅需计算y⁽³⁾ ₅₅，无需计算y⁽³⁾ ₅₄。见表6。

表6第3行规格化、第3列消元后的Y(n,2n)⁽³⁾′和Y(n,d₁)⁽³⁾′数组的状态

(5)对第4行规格化、第4列消元后的Y(n,2n)⁽⁴⁾′和Y(n,d₁)⁽⁴⁾′数组的状态

1)对第4行元素规格化，Y(n,d₁)数组中仅规格化元素y⁽³⁾ ₄₅，所得元素y⁽⁴⁾ ₄₅放在y⁽³⁾ ₄₅元素的右侧；对第4列元素消元，就是第4行规格化前的元素y⁽³⁾ ₄₅的对称元素y⁽³⁾ ₅₄进行消元。见表7。

2)对y⁽³⁾ ₄₅(y⁽³⁾ ₅₄)消元仅需计算y⁽⁴⁾ ₅₅。见表7。

表7第4行规格化、第4列消元后的Y(n,2n)⁽⁴⁾′和Y(n,d₁)⁽⁴⁾′数组的状态

(6)由于对角元不参加规格化计算，因此对第5行规格化后的Y(n,2n)⁽⁵⁾′和Y(n,d₁)⁽⁵⁾′数组的状态仍然如表7所示。

实施例2。对IEEE-30、-57、-118、-300系统分别用传统高斯消元法和本申请方法求取Z阵，其计算时间的比较结果如表8所示。

表8传统方法和本申请方法求解Z阵计算时间的比较

t_1r：读取Y(n,2n)结构数据文件时间；

t_2r：读取Y(n,d)结构数据文件时间；

t_2r/t_1r：读取Y(n,d)与读取Y(n,2n)时间的百分比；

t_1f：传统高斯消元法前代过程计算时间；

t_2f：本申请方法前代过程计算时间；

t_2f/t_1f：本申请方法与传统高斯消元法前代过程计算时间的百分比；

t_1fb：传统高斯消元法前代+回代过程计算时间；

t_2fb：本申请方法前代+回代过程计算时间；

t_1fb/t_2fb：本方法与传统高斯消元法前代+回代过程计算时间的百分比。

从表8可以看出，对IEEE-30、-57、-118、-300系统，读取Y(n,d)结构存贮的数据所需的时间约占读取Y(n,2n)结构存贮的数据所需时间的14.5％～1.5％。本申请方法的前代过程与传统高斯消元法相比，计算速度可提高约70％～99％，而本申请方法的前代+回代过程，计算速度可提高约78％～98％。上述计算结果足以表明本申请方法不仅仅可大大减少存贮单元，还能大幅提高数据文件的读取速度以及消元计算速度方面。

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

Claims (1)

1.一种非零元素随机存放和随机对称消元求取电力系统节点阻抗的方法，其特征是包括以下步骤：

步骤1：打开按随机顺序形成的仅含对角元素和非零的上三角元素的节点导纳矩阵Y阵数据文件Y(n,d)，将数据读入Y(n,d₁)数组；

所述的Y(n,d)数据文件分为3组：第1组为静态非零元素计数组，存贮Y阵中与对角元连接的非零的非对角元数S_i，仅1列；第2组为对角元组，共3列，存贮对角元的行号和参数g_ii、b_ii；第3组为非零的非对角元组，共3l_max列，其中l_max为系统中各节点上三角的最大连接支路数，而各节点上三角的连接支路数为l_i，且有S_i＝l_i，存贮与对角元连接的所有非零的非对角元的列号j及参数g_ij、b_ij；Y(n,d)数据文件的列数为d＝3l_max+4；

所述的Y(n,d₁)数组也分为3组，第1组为动态非零元素计数组S′_i，存贮Y阵中与对角元连接的、消元过程中上三角动态的非零的非对角元数，共1列，包括消元过程中上三角新产生的非零的非对角元；第2组为对角元组，共3列，存贮对角元的行号和消元过程中变化的参数g′_ii、b′_ii；第3组为非零的非对角元，共5S′_max列，S′_max为Y阵在消元过程中上三角所产生的最大的动态非零元素数，而各节点在消元过程中上三角所产生的动态非零元素数为S′_i，且S′_max>>l_max，存贮与该对角元连接的所有非零的非对角元的列号j及规格化前的参数g′_ij、b′_ij和规格化后的参数g″_ij、b″_ij，Y(n,d₁)数组的列数为d₁＝5S′_max+4；

步骤2：将Y(n,d₁)数组与E_n阵构成增广阵B_n＝[Y(n,d₁)E_n]；

步骤3：对B_n阵进行n-1次含规格化的基于对称稀疏技术的高斯消元得

(1)直接用四角规则完成Y(n,d₁)数组中随机存放元素的消元计算；

以对角元素为参考元素，将第k行规格化前每个参数的列号i与行号k分别互换：

并将各个y_ik分别作为第i行的消元元素，而将第k行各个规格化后的参数作为交叉元素，按四角规则直接完成Y(n,d₁)数组中随机存放元素的消元计算；所述的四角规则为：计算元素的新值＝其原值-消元元素*交叉元素；

(2)Y(n,d₁)数组中非零元素的快速判断和计算元素的快速确定；

1)将第k行各个规格化前每个参数的列号i与行号k分别互换：

并分别计算第k列、第i行各个y_ik与第k行各个规格化后参数的列号在交互点上的元素y′_ii、y′_ij、y′_ip、y′_im；

2)如果消元计算过程中没有产生新的非零计算元素，则将计算元素的新值直接替换原第2～3列中计算元素的原值；如果产生新的非零计算元素，则在对应行紧邻右侧的新的非零元素组中第1～3列分别加入新产生的非零元素的列号和参数；

(3)Y(n,d₁)数组中计算元素的对称算法

分别将第k行各个规格化前元素的列号与行号互换，分别与第k行各个规格化后元素的列号进行比较，并按下述三种情况处理：

1)如果第k行规格化前元素的列号i大于规格化后元素的列号j，则该元素为下三角非零的非对角元素，完全可不用计算；

2)如果第k行规格化前元素的列号i等于规格化后元素的列号j，则该元素为对角元素，需将其计算后的值y′_ii＝g′_ii+jb′_ii直接替换第i行原对角元的值y_ii＝g_ii+jb_ii；

3)如果第k行规格化前元素的列号i小于规格化后元素的列号j，则该元素为上三角元素，需将其列号j分别与第i行元素的各个列号进行比较，如果有相同的列号，则将其计算后的值y′i_j＝g′i_j+jb′_ij直接替换第i行原该组元素的值y_ij＝g_ij+jb_ij；如果没有相同的列号，则将其列号j及计算后的值y′_ij＝g′_ij+jb′_ij直接放在紧邻的新的非对角元组的第1～3列；

(4)Y(n,d₁)数组中计算元素的存放方式

计算得到的非零元素新值或新产生的非零元素均存放在非零元素组的第2～3列，对其规格化后的元素存放在第4～5列；

步骤4：规定Z阵中Z_k阵的求取顺序为第n～1列以及Z_k阵元素的求取顺序为z_kk～z_1k，再根据Y(n,d₁)^(k-1)′Z_k＝E_k ^(k-1)′分步回代求解Z_k阵中对角元z_kk及以上的元素，并根据对称性得z_kk以左的元素；

(1)Y(n,d₁)^(k-1)′数组第k行及其以上的元素与Y(n,d₁)^(n-1)′数组中第k行及其以上的元素完全相同，因此用Y(n,d₁)^(n-1)′阵第k行及以上的元素代替Y(n,d₁)^(k-1)′阵与E_k ^(k-1)′阵求解Z_k阵；

(2)在Y(n,d₁)^(k-1)′数组的消元过程中，原有非零元素数值的变化直接用该非零元素的新值替换，而新产生的非零元素仍然随机地按顺序排放在紧邻的新的非零元素组中；

(3)用E_k ^(k-1)′阵求解Z_k阵对角元及以上的元素z_kk～z_1k时，对E_k ^(k-1)′阵仅需求取其对角元素

(4)用Y(n,d₁)^(k-1)′数组中规格化后的元素对方程Y(n,d₁)^(k-1)′Z_k＝E_k ^(k-1)′进行回代求取Z_k阵中对角元z_kk和及以上的元素z_jk；

(5)根据对称性得z_kk以左的元素；

步骤5：求出Z阵并输出结果。

Images