分段對稱反向高斯-約當消元法及其應用

陳 懇，熊哲浩，魏藝君，廖嘉文

(南昌大學信息工程學院，江西南昌 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陣元素的對稱計算等等。大大簡化了傳統約當法計算過程，與高斯法、約當法相比，計算速度大大提高，并顯著提高編程效率。