CN109241492A - 快速求取电力系统节点阻抗矩阵的高斯-约当消元新算法 - Google Patents

Abstract

本发明提供一种快速求取电力系统节点阻抗矩阵Z的高斯‑约当消元新算法，涉及电力系统分析计算领域，主要包括以下步骤：读取节点导纳矩阵Ｙ数据；将n阶Ｙ阵与n阶单位矩阵E构成n*2n阶特殊增广阵B＝YE；对B阵中Ｙ阵第1～n‑1列下三角元素消元，用四角规则计算Ｙ阵对角元素和上三角元素以及E阵对角元素和下三角元素，而通过每行规格化前的上三角元素赋值得到Ｙ阵的下三角元素，得B^(n‑1)′＝Y^(n‑1)′E^(n‑1)′；用四角规则对B^(n‑1)′＝Y^(n‑1)′E^(n‑1)′中Y^(n‑1)′阵的上三角元素从第n～2列消元，仅计算E^(n‑1)′阵的对角元素和下三角元素，得B^2(n‑1)′＝Y^2(n‑1)′E^2(n‑1)′；根据E^2(n‑1)′阵可直接得到Z阵的对角元素及下三角元素；根据对称性得Z阵的上三角元素；将Z阵写入数据文件。用本申请方法求取IEEE‑57～‑300节点系统Z阵，比传统高斯‑约当消元法计算速度高约65％。

Description

快速求取电力系统节点阻抗矩阵的高斯-约当消元新算法

技术领域

本发明涉及到一种快速求取电力系统节点阻抗矩阵的高斯-约当消元新算法，可用于电力系统等各个工程领域的相关计算。

背景技术

在电力系统中，节点阻抗矩阵Z是满矩阵，所含电力系统节点信息比节点导纳矩阵Y更多，因此求解Z阵的过程要比求解Y阵更为复杂。常见的求解n*n阶Z阵的方法有支路追加法、LDU三角分解法(LDU法)、高斯消元法(高斯法)等。其中支路追加法不但计算过程繁琐，且对已经形成Y阵的系统实际意义也不大。而传统的LDU法、高斯法需将对方程YZ＝E中Z阵的求解转换成对n个方程YZ_k＝E_k中Z_k阵的求解，即将Y阵与E_k阵分别构成n个n+1阶的增广阵，进而求解n个Z_k阵的列向量，从而得到整个Z阵。虽然改进的LDU法和高斯法可仅求解Z_k阵的对角元及以上元素，利用对称性求解完整的Z阵，但由于LDU法和高斯法均含前代和回代计算，且LDU法中还要求解n个中间矩阵W_k、H_k，因此其计算过程繁琐而影响计算速度。

派生于高斯法的高斯-约当消元法(约当法)，计算原理与高斯法类似，只是高斯法一般是对下三角消元，而约当法是同时对上下三角消元，没有回代过程，因此其计算过程简单。但用传统约当法求解变系数方程时也存在几个问题：①目前文献中没有用约当法求解Z阵的内容介绍；②将求解一个n阶的方程分解成求解n个n*(n+1)阶的方程，一个个求解，计算过程繁琐复杂；③用约当法对上下三角进行消元时均按从左到右的顺序，计算效率低下；④求解对称方程时未能利用元素的对称性进行求解等等。

发明内容

为了克服现有技术方法计算效率不佳的情况，本申请提出了快速求取电力系统节点阻抗矩阵的高斯-约当消元新算法。本申请是通过以下技术方案实现，主要包括以下步骤：

步骤1：读取Y阵数据；

步骤2：将n阶Y阵与n阶单位矩阵E构成n*2n阶特殊增广阵B＝YE；

传统方法中均是将n阶Y阵与n阶E阵中第k列的E_k阵(E阵中的1列)构成n个n*(n+1)阶的增广阵B_k＝YE_k一列列地求解。本申请方法是将整个Y阵与整个E阵构成n*2n阶的特殊增广阵B＝YE一次性求解，其展开后如式(1)。

步骤3：对B阵中Y阵第1～n-1列下三角元素消元，用四角规则计算Y阵对角元素和上三角元素以及E阵对角元素和下三角元素，而通过每行规格化前的上三角元素赋值得到Y阵的下三角元素，得B^(n-1)′＝Y^(n-1)′E^(n-1)′；

(1)按第1～n-1列正向顺序对Y阵的下三角元素消元；

(2)用四角规则完成对Y阵下三角元素消元的计算；

传统方法中必须利用消元计算公式，计算过程繁琐复杂、不利于编程计算。本申请方法中用四角规则可直接完成消元计算，无需计算公式。对第k列元素消元前后简化矩阵中四角规则应用如式(2)，设对角元素已取倒。

其中，为取倒后的对角元素；为规格化前(后)的交叉元素；为消元元素； 为计算元素的前值和新值。由于四个元素在矩形的四个角上，所以称为四角规则。Y、E阵元素计算方法一样，因此以Y阵元素计算为例。

对式(2)中元素消元，传统方法需用消元计算公式分别计算四个元素。根据式(2)，可不用消元计算公式而直接用四角规则写出计算元素的表达式如下。

四角规则也可用文字表示为“以对角元素为参考元素，计算元素的新值＝其前值-消元元素*交叉元素”，从而大大简化编程计算。

(3)对Y阵下三角元素消元，仅计算Y阵对角元素和上三角元素；

传统的消元计算没有利用对称矩阵元素结构的特点，需计算第k列消元元素以右和第k行交叉元素以下的所有元素，如对式(2)中元素消元，需计算这些元素，计算量大。

本申请方法则利用了对称矩阵元素结构的特点，对第k列元素消元时，仅计算第k列消元元素以右和第k行交叉元素以下所有的对角元素和上三角元素，即仅计算所有消元元素所在行、其对角元素及其以右的元素(上三角元素)，而对消元元素所在行、对角元以左的元素(下三角元素)可通过规格化前的上三角元素赋值得到，从而省略Y阵中所有下三角非对角元素的计算。因而仅需计算这些元素，而元素是通过上三角元素赋值得到。

(4)对Y阵下三角元素消元，仅计算E阵对角元素和下三角元素；

传统方法对Y阵下三角元素消元时，没有利用E阵元素结构的特点，因此需计算E阵中的所有元素，不但导致E阵中大量的无效计算，而且使得Z阵元素的对称性难以应用。

本申请方法则利用E阵元素结构的特点，对Y阵第k列元素消元时，仅计算第k列消元元素以右和E阵第k行以下、其第1～k列的所有元素，而E阵第k行以下、其第k列以右的元素均不计算，即仅计算E阵中部分对角元素及部分下三角元素，且不计算E阵中所有上三角的非对角元素。这种利用E阵元素结构的特点计算的方式，使得后续计算中Z阵元素的对称性可以得以应用。将式(2)右侧E阵部分元素写出可得式(3)。根据式(3)可以看出，由于E阵元素结构的特点，E阵中第k行以下、第k列以右的所有元素都不计算。

(5)规格化前将对角元素取倒，减少程序中的除法运算；

传统消元计算中，未考虑规格化前将对角元素取倒，如在形成因子表的过程中，对角元素并未取倒，而是在形成因子表后将对角元素取倒，因此规格化计算中有不少的除法运算。此外，对第k行元素规格化时，计算了第k行E阵的全部元素，也导致不少无效计算。

本申请方法则在规格化前将对角元素取倒，从而将后续规格化计算中的除法运算全部转化为乘法运算以提高计算速度。此外，对第k行E阵元素规格化时，仅规格化第k行E阵中的第1～k列的元素，其第k列以右的元素不规格化，大大减少无效的规格化计算。

(6)根据Y阵第k+1行规格化前元素得到Y阵第k+1列的下三角元素，继续对第k+1列消元，直至完成第n-1列消元，得B^(n-1)′阵；

对B阵下三角进行n-1次约当消元后得B^(n-1)′阵，其展开后如式(4)。

步骤4：用四角规则对B^(n-1)′＝Y^(n-1)′E^(n-1)′中Y^(n-1)′阵的上三角元素从第n～2列消元，仅计算E^(n-1)′阵的对角元素和下三角元素，得B^2(n-1)′＝Y²⁽ⁿ-¹⁾′E^2(n-1)′；

(1)按第n～2列反向顺序对Y^(n-1)′阵的上三角元素消元；

(2)用四角规则完成对Y^(n-1)′阵上三角元素消元的计算；

传统方法对Y^(n-1)′阵上三角元素消元时也必须用消元计算公式，同样不直观、不适于编程。而本申请方法中用四角规则也直接完成对Y^(n-1)′阵上三角元素的消元计算，只是此时各元素的位置与对下三角元素消元时有所不同。对上三角第k列元素消元前后简化矩阵中四角规则的应用如式(5)。

式(5)与式(2)相比，其对角元素、交叉元素、计算元素的位置不同，但相应的四个元素仍然均在矩形的四个角上，所以根据式(5)中元素结构，也可用四角规则直接写出其计算元素的表达式，从而大大简化编程计算。

(3)对Y^(n-1)′阵上三角元素消元，仅计算E^(n-1)′阵中的对角元素和下三角元素，不计算E^(n-1)′阵的上三角元素和Y^(n-1)′阵中所有元素；

传统方法对Y^(n-1)′阵上三角元素消元时按第2～n列的正向顺序，要计算Y^(n-1)′阵中消元元素以右的所有上三角元素，且由于没有利用E阵元素结构特点，还需计算E^(n-1)′阵中所有元素，计算量巨大。

1)本申请方法对Y^(n-1)′阵上三角元素消元按第n～2列的反向顺序，与对Y阵第1～n-1列下三角元素消元顺序相反，同时利用了交叉元素为零无需消元计算的基本技巧和E阵元素结构的特点，因此无需计算Y^(n-1)′阵中上三角消元元素以右的Y^(n-1)′阵中的任何上三角元素，

2)对Y^(n-1)′阵上三角第k列元素消元，仅需计算消元元素以右E^(n-1)′阵中的对角元素和下三角元素。如对式(5)中元素消元，仅计算

其中对元素消元，仅计算其中

上述计算方式可省略对Y^(n-1)′阵中上三角消元时，Y^(n-1)′、E^(n-1)′阵中所有上三角元素计算。因此在式(5)中只有对e元素计算，没有对y元素的计算。其计算过程特点如下：

1)不计算Y^(n-1)′阵的任何元素，仅计算E^(n-1)′阵对角元素和下三角元素。

2)仅计算E^(n-1)′阵中消元元素所在行的第1个元素～对角元素为止，对角元素以右的元素均不计算。

3)由于规格化计算在对下三角元素消元时已经完成，因此对上三角元素消元时不用再考虑。

(4)继续对第k-1列元素消元，直至完成第2列消元，得B^2(n-1)′阵；

对B^(n-1)′阵进行n-1次的上三角约当消元后得B^2(n-1)′阵，其展开后为B^2(n-1)′＝Y^{2(n -1)}′E^2(n-1)′＝EE²⁽ⁿ-¹⁾′如式(6)，此时Y^2(n-1)′阵实际上已经变为E阵，而E^2(n-1)′阵实际上就是Z阵的对角元素及下三角元素。

式(6)Y^2(n-1)阵与式(4)Y^(n-1)′阵中上、下三角元素完全一致是因为对Y^2(n-1)，阵从右往左消元过程不用计算Y^(n-1)′阵中的任何上三角元素。因为根据四角规则可看出，从右往左消元时，Y^(n-1)′阵中消元元素以右上三角的交叉元素由于已经被消元为零，因此对Y^(n-1)′阵中其余上三角元素的计算可省略。

但如果按传统方法对Y^(n-1)′阵的上三角元素采用从第2～n列，即从左往右正向消元时，就要对Y^(n-1)′阵中消元元素以右的所有上三角元素进行计算元素，从而大大增加计算量。

此外，式(6)为Y^2(n-1)′阵元素的实际状态，而其理论状态为Y^2(n-1)′＝E。按式(6)的计算方式编程可省去大量不必要的计算，大大提高了计算效率。

步骤5：根据E^2(n-1)′阵可直接得到Z阵的对角元素及下三角元素；

步骤6：根据对称性得Z阵的上三角元素；

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

本申请方法与传统方法相比，具有几个优点：

(1)用约当法一次性求解n阶Z阵的对角元素及下三角元素，利用Z阵元素的对称性求取Z阵的上三角元素，而不是将Z阵分成n个Z_k阵求解Z_k阵的所有元素。

(2)用四角规则直接完成消元计算，无需消元计算公式。

(3)对Ｙ阵下三角元素按从左往右正向顺序消元，仅计算Ｙ阵的对角元素和上三角元素，根据对称性得到Ｙ阵的下三角元素；仅计算E阵中对角元素和下三角元素。

(4)对Y^(n-1)′阵上三角元素按从右往左反向顺序消元，仅计算E^(n-1)′阵中对角元素和下三角元素，不计算E^(n-1)′阵的上三角元素和Y^(n-1)′阵中所有元素。

(5)规格化前将对角元素取倒，减少程序中的除法运算。

(6)计算原理及过程更加简单直观，计算速度快。

附图说明

图1传统高斯消元法求解节点阻抗矩阵Z的计算流程图

图2传统高斯-约当消元法求解节点阻抗矩阵Z的计算流程图

图3本申请方法3求解节点阻抗矩阵Z的计算流程图

具体实施方式

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

不考虑应用稀疏矩阵技术，分别用传统高斯法、传统约当法、改进高斯法、以及为说明本申请方法中各种快速算法的影响，将本申请方法分为本申请方法1(上下三角均从左往右消元，同时应用Ｙ、Z阵对称性)、本申请方法2(下三角从左往右、上三角从右往左，仅应用Z阵对称性，未应用Ｙ阵对称性)和本申请方法3(下三角从左往右、上三角从右往左消元，同时应用Ｙ、Z阵对称性)，对IEEE-57、-118、-300节点系统求取Z阵，各种方法计算时间的比较如表1所示。

由于传统约当法中常将n阶方程分成n个列矩阵分别进行求解计算效率较低，因此本实施例作在用传统约当法计算时，也将Ｙ、E阵构成特殊增广阵B＝YE求解Z阵，而不是按增广阵B＝YE_k分别求解Z_k阵进而求解Z阵，但对高斯法和改进高斯法仍按B＝YE_k分别求解Z_k阵进而求解Z阵。

此处传统方法计算Z阵时，均未用Ｙ、Z阵对称性；改进高斯法未用Ｙ阵对称性，但用Z阵对称性；本申请方法1和3均用Ｙ、Z阵对称性；本申请方法2仅用Z阵对称性，未用Ｙ阵对称性。

其中，t₁：传统约当法计算时间；t₂：传统高斯法计算时间；t₃：改进高斯法计算时间(未用Y阵、但用Z阵对称性)；t₄：本申请方法1计算时间(上下三角均从左往右消元，用Y、Z阵对称性)；t₅：本申请方法2计算时间(下三角从左往右、上三角从右往左，未用Y阵、但用Z阵对称性)；t₆：本申请方法3计算时间(下三角从左往右、上三角从右往左，用Y、Z阵对称性)。

表1各种方法求取Z阵计算时间的比较

t₂/t₁：传统高斯法计算时间与传统约当法计算时间百分比；

t₃/t₂：改进高斯法与传统高斯法计算时间百分比；

t₄/t₁：本申请方法1与传统约当法计算时间百分比；

t₅/t₁：本申请方法2与传统约当法计算时间百分比；

t₅/t₃：本申请方法2与改进高斯法计算时间百分比；

t₆/t₁：本申请方法3与传统约当法计算时间百分比；

t₆/t₃：本申请方法3与改进高斯法计算时间百分比；

t₆/t₄：本申请方法3与本申请方法1计算时间百分比；

t₆/t₅：本申请方法3与本申请方法2计算时间百分比。

根据表1可以看出，本申请方法3利用了本申请提出的所有计算技巧，对IEEE各节点系统的计算时间最少、计算效率最高。传统高斯法与传统约当法相比(t₂/t₁)，计算速度约高5～13％，且随着系统的加大两者计算速度趋近相同；本申请方法1与传统约当法相比(t₄/t₁)，计算速度可提高约45％；本申请方法2与改进高斯法相比(t₅/t₃)，计算速度可提高约6～13％，且随着系统的加大优势更加明显，使得本申请方法2的计算效率高于改进高斯法；本申请方法2与传统约当法相比(t_s/t₁)，计算速度可提高约60％；本申请方法3与本申请方法1相比(t₆/t₄)，计算速度提高约40％；本申请方法3与本申请方法2相比(t₆/t_５)，计算速度可提高约24％；本申请方法3与传统约当法相比(t₆/t₁)，计算速度可提高约66％。上述数据比较说明了本申请方法计算效率上的优势，可大大提高求解Z阵的计算速度。本申请方法也可用于电力系统分析计算中各种线性方程组的快速求解，并大大提高计算效率。

本申请方法可采用任何编程语言和编程环境实现。本申请中采用C++编译语言，开发环境是Visual C++。

Claims (1)

1.快速求取电力系统节点阻抗矩阵的高斯-约当消元新算法，其特征在于，包括以下步骤：

步骤1：读取节点导纳矩阵Y数据；

步骤2：将n阶Y阵与n阶单位矩阵E构成n*2n阶特殊增广阵B＝YE；

步骤3：对B阵中Ｙ阵第1～n-1列下三角元素消元，用四角规则计算Ｙ阵对角元素和上三角元素以及E阵对角元素和下三角元素，而通过每行规格化前的上三角元素赋值得到Ｙ阵的下三角元素，得B^(n-1)′＝Y^(n-1)′E^(n-1)′；

(1)按第1～n-1列正向顺序对Ｙ阵的下三角元素消元；

(2)用四角规则完成对Y阵下三角元素消元的计算；

(3)对Y阵下三角元素消元，仅计算Ｙ阵对角元素和上三角元素；

(4)对Ｙ阵下三角元素消元，仅计算E阵对角元素和下三角元素；

(5)规格化前将对角元素取倒，减少程序中的除法运算；

(6)根据Ｙ阵第k+1行规格化前上三角元素得到Y阵第k+1列的下三角元素，继续对第k+1列消元，直至完成第n-1列消元，得B^(n-1)′阵；

步骤4：用四角规则对B^(n-1)′＝Y^(n-1)′E^(n-1)′中Y^(n-1)′阵的上三角元素从第n～2列消元，仅计算E^(n-1)′阵的对角元素和下三角元素，得B^2(n-1)′＝Y^2(n-1)′E^2(n-1)′；

(1)按第n～2列反向顺序对Y^(n-1)′阵的上三角元素消元；

(2)用四角规则完成对Y^(n-1)′阵上三角元素消元的计算；

(3)对Y^(n-1)′阵上三角元素消元，仅计算E^(n-1)′阵中的对角元素和下三角元素，不计算E^(n-1)′阵的上三角元素和Y^(n-1)′阵中所有元素；

(4)继续对第k-1列元素消元，直至完成第2列消元，得B^2(n-1)′阵；

步骤5：根据E^2(n-1)′阵可直接得到Z阵的对角元素及下三角元素；

步骤6：根据对称性得Z阵的上三角元素；

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