分段对称反向高斯-约当消元法及其应用

陈 恳，熊哲浩，魏艺君，廖嘉文

(南昌大学信息工程学院，江西南昌 330031)

1 引言

高斯消元法[1-5](高斯法)一般对下三角元素消元，用上三角回代求解方程，而高斯-约当消元法[6-10](约当法)是同时对上下三角元素消元直接得到方程解。两者计算原理相似，均用于求解变系数而不是常系数方程。由于约当法上三角元素消元计算比高斯法上三角元素回代计算更为复杂，使其计算速度低于高斯法。此外，约当法与高斯法问题类似，如将对n阶方程的求解分解成对n个n*(n+1)阶方程的求解，计算过程冗余且不易展现及利用相关矩阵元素计算规律；对对称矩阵如节点导纳矩阵Y下三角消元未利用元素对称性；对Y阵上三角消元有大量的无效计算；对Y阵上下三角元素消元均须应用计算公式，不利于过程理解和编程；未将取倒后的对角元素作为规格化因子，造成大量的除法计算；计算单位矩阵E元素时未考虑E阵元素结构特点等，使约当法应用受到很大限制。

本文提出分段对称反向约当法(新方法)，包括：构建特殊增广阵；考虑Y阵元素结构特点分段对其上下三角元素消元；对下三角元素采用正向消元及对称计算，简化所有下三角元素计算；对上三角元素采用反向消元，省略消元元素列以右元素的计算；将取倒的对角元素作为规格化因子以减少除法运算；消元和规格化时仅计算E阵部分有效对角元素和下三角元素，省略所有上三角元素计算；计算Y、E阵元素均用四角规则而无需计算公式；用约当法求解常系数节点阻抗矩阵Z，并根据E阵元素结构特点完成Z阵元素对称计算[13]。与高斯法、约当法相比，新方法计算速度可大大提高。

2 特殊增广阵的构成及元素计算

2.1 传统增广阵的构成

求解式(1)的主要问题是：分别求解n个n*(n+1)阶方程计算过程冗余繁琐，不利于其求解常系数方程，如分别求解n个Zk阵时也不利于Z阵元素对称性应用[1-3，11-12]；E阵元素结构特点不明显且未能利用该特点简化E阵元素计算，使求解过程复杂；对式(1)从左至右、从上至下消元对对称矩阵下三角元素难以实现对称算法；对上三角元素消元时，消元元素右侧大量元素被重复计算；未将取倒后的对角元素作为规格化因子等。导致计算过程复杂、大量冗余计算，使计算速度受到很大影响。此外，其上下三角元素消元均需参考或利用高斯法计算公式，不利计算过程理解和编程[6-10]。

(1)

2.2 特殊增广阵的构成

由于求解YZk=Ek方程时，对Y阵第n次约当消元后其第1～k列计算结果与对Y阵第k次约当消元后第1～k列计算结果完全相同，因此可用Y阵第n次约当消元中第1～k列计算结果与各个Ek阵同时求解Zk阵，而不必反复对n个[YEk]阵进行约当消元。为此将Y、E阵构成特殊增广阵[YE]如式(2)，直接求取一个YZ=E方程，而不是分别求解n个方程YZk=Ek，从而大大提高计算效率。

(2)

2.3 特殊增广阵的约当法消元

对[YE]阵上下三角元素从左到右按列连续约当消元，可得[Y(n)″E(n)″]阵如式(3)。元素上标表示元素被计算次数，其Y(n)″阵元素上标从左至右均为0～(n-1)，而E(n)″阵元素上标均为n。

和约当法一样，直接对[YE]阵下三角消元难以利用元素对称性[15-16]，对上三角消元有大量的无效计算；若消元过程中对角元素不取倒或按因子表法形成因子表后将对角元素取倒[1-4]，则规格化或回代计算中均有大量除法运算；若不考虑E阵元素结构特点，E阵元素的计算量也大大增加等。此时[YE]阵并未改变约当法计算速度不佳的状况。

2.4 特殊增广阵的分段对称反向约当消元

因此可对[YE]阵进行分段对称反向约当消元，即先对Y阵下三角元素按列从左至右正向对称消元，并将取倒后的对角元素作为规格化运算因子，得[Y(n-1)′E(n-1)′]阵；再对Y(n-1)′阵上三角元素按列从右至左反向消元，不计Y(n-1)′阵中消元列以右元素得[Y(n)″E(n)″]阵如式(4)。式(4)中Y(n)″阵对角元素上标为k-1，同行以右上三角元素的上标为k；下三角元素只有一次赋值，其上标为“1”；E(n)″阵所有上三角元素均未被计算，下三角元素上标均从n开始按“-1”递减到对角元素为(n-k+1)。比较式(4)与式(3)可知，新方法中Y(n)″、E(n)″阵元素计算及除法计算大大减少，因此计算速度大幅度提升。

(3)

(4)

(5)

2.5 约当法的四角规则

约当法对上下三角元素消元均可参考高斯法计算公式，但应用和编程均极为不便。

为此，给出约当法对Y阵上下三角第k列元素消元前后简化的[Y(k)″E(k)″]阵中元素相应状态如式(5)，各元素定义如下：

由于参加计算的元素在矩形的四个角上，故称四角规则[15-16]。因此根据四角规则，即元素在矩阵中位置可直接写出相应计算式而无需计算公式，也无需考虑元素上下标。

2.6 分段对称反向约当法中元素的计算

用四角规则根据元素在矩阵中的位置可直接写出相应计算式而无需计算公式。但此时约当法中对下三角消元未利用元素对称性；对上三角消元计算效率极低；未将取倒后的对角元素作为规格化因子；对E阵元素规格化和消元均未利用E阵元素结构特点，约当法计算速度仍然不如高斯法。

分段对称约当法利用了Y阵元素对称性和E阵元素结构特点，分段对[YE]阵上下三角元素消元；将取倒后的对角元素作为规格化因子以减少除法计算；对Y阵下三角元素按第1～n-1列、从上至下顺序正向消元，并采用对称算法，简化下三角元素计算，得[Y(n-1)′E(n-1)′]阵；再对[Y(n-1)′E(n-1)′]阵中Y(n-1)′阵上三角元素按第n～2列、从下至上顺序反向消元，不计算Y(n-1)′阵消元元素以右元素，仅计算E(n-1)′阵元素，得[Y(n)″E(n)″]阵如式(4)。

2.6.1 下三角元素正向对称消元

对Y阵第k列下三角元素消元，利用Y阵元素对称性，仅计算Y阵中相应的对角元素及上三角元素，根据对称性得到下三角元素。计算后的简化矩阵如式(6)。

对第k列元素进行含规格化消元前，其第k列以右和第k行以下所有上三角元素上标均为k-1；下三角元素均未被计算，上标均为0。Y、E阵元素计算过程及特点如下：

1)第k行规格化前交叉元素对称赋值给第k列消元元素，则第k列消元元素上标为1。

2)第k行对角元素取倒并对第k行元素规格化，第k行交叉元素上标为k。

3)对第k列元素消元仅计算第k列消元元素以右和第k行交叉元素以下所有的对角元素和上三角元素，计算后元素上标均为k。下三角元素均不计算，上标仍均为0。

4)对E阵第k行元素规格化，仅规格化其第1～k列元素，而不规格化第k+1～n列元素。

5)对Y阵第k列下三角元素消元，利用E阵元素结构特点，仅计算E阵第k行以下、第1～k列所有对角元素和下三角元素，而第k列以右元素和所有上三角元素均不计算。

2.6.2 上三角元素反向消元

对Y(n-1)′阵上三角元素按第2～n列正向消元，需计算Y(n-1)′阵中消元元素以右所有上三角元素；若未利用E阵元素结构特点，还需计算E(n-1)′阵中所有元素，计算量巨大。

根据式(6)，继续对Y(n-1)′阵上三角元素按第n～2列从下至上反向消元，利用消元后交叉元素为零特点，可不计算Y(n-1)′阵中任何元素；利用E阵元素结构特点，可仅计算E(n-1)′阵中部分对角元素和下三角元素，大大减少[Y(n)″E(n)″]阵冗余计算。对Y(n-1)′阵第k列元素反向消元前后简化矩阵如式(7)。

(6)

(7)

对Y(n-1)′阵第k列元素反向消元，Y、E阵元素计算过程及特点如下：

2)仅计算E(n-1)′阵中相应行中的第1个元素到对角元素，对角元素以右元素均不计算，即仅计算E(n-1)′阵中第1～k列的下三角元素和对角元素，不计算第k+1～n列的下三角元素及任何上三角元素。

对[Y(n-1)′E(n-1)′]阵上三角元素分段对称约当消元后得式(4)，即[Y(n)″E(n)″]。

3 求取节点阻抗矩阵

常用求取Z阵的方法有求逆法、支路追加法、LDU三角分解法[1-3](LDU法)等。当系统节点数较大时，求逆法不适用于求取Z阵，而支路追加法计算过程繁琐[13]，对形成Y阵的系统意义不大。求解常系数方程的LDU法是将对方程YZ=E中Z阵的求解转换为分别对n个Zk阵求解，含前代和回代计算；要求解2n个中间矩阵Wk、Hk；且无法使用Z阵元素对称性，必须求解整个Z阵元素，计算速度受到很大影响。

至今没有文献用求解变系数方程的约当法求解常系数方程YZ=E，但若适当应用约当法的特点，则可实现该功能。

如前所述，对[YE]阵分段对称约当消元后得[Y(n)″E(n)″]阵，此时Y(n)″=E，E(n)″阵就是Z阵的对角元素及下三角元素，按对称性可得Z阵上三角元素[14]。

4 例分析

分别用高斯法、约当法、新方法求取IEEE-57、-118、-300节点系统Z阵，计算时间比较如表1。

根据表1可以看出，高斯法比约当法快约5～13%，随着系统规模加大逐渐减小，但约当法始终没有速度优势；新方法比约当法快约66%，与系统规模关系不大；新方法比高斯法快约62～64%，也随着系统规模加大而逐渐增加。因此与高斯法、约当法相比，新方法可大幅度提高求解电力系统线性方程组或Z阵的计算速度。

表1 求取Z阵计算时间比较

t1、t2、t3：约当法、高斯法、新方法的计算时间。

t2/t1、t3/t1、t3/t2：高斯法与约当法、新方法与约当法、新方法与高斯法计算时间的百分比。

5 结语

本文提出分段对称反向约当法，并将其用于求解常系数方程的Z阵，包括构建特殊增广阵；对Y阵下三角元素采用正向消元及对称计算；对Y阵上三角元素采用反向消元；将取倒的对角元素作为规格化因子；规格化或对上下三角元素消元时，E阵元素均仅计算其部分有效对角元素和下三角元素；Y、E阵上下三角元素均用四角规则计算而无需计算公式；综合利用E阵元素结构特点完成Z阵元素的对称计算等等。大大简化了传统约当法计算过程，与高斯法、约当法相比，计算速度大大提高，并显著提高编程效率。