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

(南昌大學(xué)信息工程學(xué)院，江西南昌 330031)

1 引言

高斯消元法[1-5](高斯法)一般對下三角元素消元，用上三角回代求解方程，而高斯-約當(dāng)消元法[6-10](約當(dāng)法)是同時(shí)對上下三角元素消元直接得到方程解。兩者計(jì)算原理相似，均用于求解變系數(shù)而不是常系數(shù)方程。由于約當(dāng)法上三角元素消元計(jì)算比高斯法上三角元素回代計(jì)算更為復(fù)雜，使其計(jì)算速度低于高斯法。此外，約當(dāng)法與高斯法問題類似，如將對n階方程的求解分解成對n個(gè)n*(n+1)階方程的求解，計(jì)算過程冗余且不易展現(xiàn)及利用相關(guān)矩陣元素計(jì)算規(guī)律；對對稱矩陣如節(jié)點(diǎn)導(dǎo)納矩陣Y下三角消元未利用元素對稱性；對Y陣上三角消元有大量的無效計(jì)算；對Y陣上下三角元素消元均須應(yīng)用計(jì)算公式，不利于過程理解和編程；未將取倒后的對角元素作為規(guī)格化因子，造成大量的除法計(jì)算；計(jì)算單位矩陣E元素時(shí)未考慮E陣元素結(jié)構(gòu)特點(diǎn)等，使約當(dāng)法應(yīng)用受到很大限制。

本文提出分段對稱反向約當(dāng)法(新方法)，包括：構(gòu)建特殊增廣陣；考慮Y陣元素結(jié)構(gòu)特點(diǎn)分段對其上下三角元素消元；對下三角元素采用正向消元及對稱計(jì)算，簡化所有下三角元素計(jì)算；對上三角元素采用反向消元，省略消元元素列以右元素的計(jì)算；將取倒的對角元素作為規(guī)格化因子以減少除法運(yùn)算；消元和規(guī)格化時(shí)僅計(jì)算E陣部分有效對角元素和下三角元素，省略所有上三角元素計(jì)算；計(jì)算Y、E陣元素均用四角規(guī)則而無需計(jì)算公式；用約當(dāng)法求解常系數(shù)節(jié)點(diǎn)阻抗矩陣Z，并根據(jù)E陣元素結(jié)構(gòu)特點(diǎn)完成Z陣元素對稱計(jì)算[13]。與高斯法、約當(dāng)法相比，新方法計(jì)算速度可大大提高。

2 特殊增廣陣的構(gòu)成及元素計(jì)算

2.1 傳統(tǒng)增廣陣的構(gòu)成

求解式(1)的主要問題是：分別求解n個(gè)n*(n+1)階方程計(jì)算過程冗余繁瑣，不利于其求解常系數(shù)方程，如分別求解n個(gè)Zk陣時(shí)也不利于Z陣元素對稱性應(yīng)用[1-3，11-12]；E陣元素結(jié)構(gòu)特點(diǎn)不明顯且未能利用該特點(diǎn)簡化E陣元素計(jì)算，使求解過程復(fù)雜；對式(1)從左至右、從上至下消元對對稱矩陣下三角元素難以實(shí)現(xiàn)對稱算法；對上三角元素消元時(shí)，消元元素右側(cè)大量元素被重復(fù)計(jì)算；未將取倒后的對角元素作為規(guī)格化因子等。導(dǎo)致計(jì)算過程復(fù)雜、大量冗余計(jì)算，使計(jì)算速度受到很大影響。此外，其上下三角元素消元均需參考或利用高斯法計(jì)算公式，不利計(jì)算過程理解和編程[6-10]。

(1)

2.2 特殊增廣陣的構(gòu)成

由于求解YZk=Ek方程時(shí)，對Y陣第n次約當(dāng)消元后其第1～k列計(jì)算結(jié)果與對Y陣第k次約當(dāng)消元后第1～k列計(jì)算結(jié)果完全相同，因此可用Y陣第n次約當(dāng)消元中第1～k列計(jì)算結(jié)果與各個(gè)Ek陣同時(shí)求解Zk陣，而不必反復(fù)對n個(gè)[YEk]陣進(jìn)行約當(dāng)消元。為此將Y、E陣構(gòu)成特殊增廣陣[YE]如式(2)，直接求取一個(gè)YZ=E方程，而不是分別求解n個(gè)方程YZk=Ek，從而大大提高計(jì)算效率。

(2)

2.3 特殊增廣陣的約當(dāng)法消元

對[YE]陣上下三角元素從左到右按列連續(xù)約當(dāng)消元，可得[Y(n)″E(n)″]陣如式(3)。元素上標(biāo)表示元素被計(jì)算次數(shù)，其Y(n)″陣元素上標(biāo)從左至右均為0～(n-1)，而E(n)″陣元素上標(biāo)均為n。

和約當(dāng)法一樣，直接對[YE]陣下三角消元難以利用元素對稱性[15-16]，對上三角消元有大量的無效計(jì)算；若消元過程中對角元素不取倒或按因子表法形成因子表后將對角元素取倒[1-4]，則規(guī)格化或回代計(jì)算中均有大量除法運(yùn)算；若不考慮E陣元素結(jié)構(gòu)特點(diǎn)，E陣元素的計(jì)算量也大大增加等。此時(shí)[YE]陣并未改變約當(dāng)法計(jì)算速度不佳的狀況。

2.4 特殊增廣陣的分段對稱反向約當(dāng)消元

因此可對[YE]陣進(jìn)行分段對稱反向約當(dāng)消元，即先對Y陣下三角元素按列從左至右正向?qū)ΨQ消元，并將取倒后的對角元素作為規(guī)格化運(yùn)算因子，得[Y(n-1)′E(n-1)′]陣；再對Y(n-1)′陣上三角元素按列從右至左反向消元，不計(jì)Y(n-1)′陣中消元列以右元素得[Y(n)″E(n)″]陣如式(4)。式(4)中Y(n)″陣對角元素上標(biāo)為k-1，同行以右上三角元素的上標(biāo)為k；下三角元素只有一次賦值，其上標(biāo)為“1”；E(n)″陣所有上三角元素均未被計(jì)算，下三角元素上標(biāo)均從n開始按“-1”遞減到對角元素為(n-k+1)。比較式(4)與式(3)可知，新方法中Y(n)″、E(n)″陣元素計(jì)算及除法計(jì)算大大減少，因此計(jì)算速度大幅度提升。

(3)

(4)

(5)

2.5 約當(dāng)法的四角規(guī)則

約當(dāng)法對上下三角元素消元均可參考高斯法計(jì)算公式，但應(yīng)用和編程均極為不便。

為此，給出約當(dāng)法對Y陣上下三角第k列元素消元前后簡化的[Y(k)″E(k)″]陣中元素相應(yīng)狀態(tài)如式(5)，各元素定義如下：

由于參加計(jì)算的元素在矩形的四個(gè)角上，故稱四角規(guī)則[15-16]。因此根據(jù)四角規(guī)則，即元素在矩陣中位置可直接寫出相應(yīng)計(jì)算式而無需計(jì)算公式，也無需考慮元素上下標(biāo)。

2.6 分段對稱反向約當(dāng)法中元素的計(jì)算

用四角規(guī)則根據(jù)元素在矩陣中的位置可直接寫出相應(yīng)計(jì)算式而無需計(jì)算公式。但此時(shí)約當(dāng)法中對下三角消元未利用元素對稱性；對上三角消元計(jì)算效率極低；未將取倒后的對角元素作為規(guī)格化因子；對E陣元素規(guī)格化和消元均未利用E陣元素結(jié)構(gòu)特點(diǎn)，約當(dāng)法計(jì)算速度仍然不如高斯法。

分段對稱約當(dāng)法利用了Y陣元素對稱性和E陣元素結(jié)構(gòu)特點(diǎn)，分段對[YE]陣上下三角元素消元；將取倒后的對角元素作為規(guī)格化因子以減少除法計(jì)算；對Y陣下三角元素按第1～n-1列、從上至下順序正向消元，并采用對稱算法，簡化下三角元素計(jì)算，得[Y(n-1)′E(n-1)′]陣；再對[Y(n-1)′E(n-1)′]陣中Y(n-1)′陣上三角元素按第n～2列、從下至上順序反向消元，不計(jì)算Y(n-1)′陣消元元素以右元素，僅計(jì)算E(n-1)′陣元素，得[Y(n)″E(n)″]陣如式(4)。

2.6.1 下三角元素正向?qū)ΨQ消元

對Y陣第k列下三角元素消元，利用Y陣元素對稱性，僅計(jì)算Y陣中相應(yīng)的對角元素及上三角元素，根據(jù)對稱性得到下三角元素。計(jì)算后的簡化矩陣如式(6)。

對第k列元素進(jìn)行含規(guī)格化消元前，其第k列以右和第k行以下所有上三角元素上標(biāo)均為k-1；下三角元素均未被計(jì)算，上標(biāo)均為0。Y、E陣元素計(jì)算過程及特點(diǎn)如下：

1)第k行規(guī)格化前交叉元素對稱賦值給第k列消元元素，則第k列消元元素上標(biāo)為1。

2)第k行對角元素取倒并對第k行元素規(guī)格化，第k行交叉元素上標(biāo)為k。

3)對第k列元素消元僅計(jì)算第k列消元元素以右和第k行交叉元素以下所有的對角元素和上三角元素，計(jì)算后元素上標(biāo)均為k。下三角元素均不計(jì)算，上標(biāo)仍均為0。

4)對E陣第k行元素規(guī)格化，僅規(guī)格化其第1～k列元素，而不規(guī)格化第k+1～n列元素。

5)對Y陣第k列下三角元素消元，利用E陣元素結(jié)構(gòu)特點(diǎn)，僅計(jì)算E陣第k行以下、第1～k列所有對角元素和下三角元素，而第k列以右元素和所有上三角元素均不計(jì)算。

2.6.2 上三角元素反向消元

對Y(n-1)′陣上三角元素按第2～n列正向消元，需計(jì)算Y(n-1)′陣中消元元素以右所有上三角元素；若未利用E陣元素結(jié)構(gòu)特點(diǎn)，還需計(jì)算E(n-1)′陣中所有元素，計(jì)算量巨大。

根據(jù)式(6)，繼續(xù)對Y(n-1)′陣上三角元素按第n～2列從下至上反向消元，利用消元后交叉元素為零特點(diǎn)，可不計(jì)算Y(n-1)′陣中任何元素；利用E陣元素結(jié)構(gòu)特點(diǎn)，可僅計(jì)算E(n-1)′陣中部分對角元素和下三角元素，大大減少[Y(n)″E(n)″]陣冗余計(jì)算。對Y(n-1)′陣第k列元素反向消元前后簡化矩陣如式(7)。

(6)

(7)

對Y(n-1)′陣第k列元素反向消元，Y、E陣元素計(jì)算過程及特點(diǎn)如下：

2)僅計(jì)算E(n-1)′陣中相應(yīng)行中的第1個(gè)元素到對角元素，對角元素以右元素均不計(jì)算，即僅計(jì)算E(n-1)′陣中第1～k列的下三角元素和對角元素，不計(jì)算第k+1～n列的下三角元素及任何上三角元素。

對[Y(n-1)′E(n-1)′]陣上三角元素分段對稱約當(dāng)消元后得式(4)，即[Y(n)″E(n)″]。

3 求取節(jié)點(diǎn)阻抗矩陣

常用求取Z陣的方法有求逆法、支路追加法、LDU三角分解法[1-3](LDU法)等。當(dāng)系統(tǒng)節(jié)點(diǎn)數(shù)較大時(shí)，求逆法不適用于求取Z陣，而支路追加法計(jì)算過程繁瑣[13]，對形成Y陣的系統(tǒng)意義不大。求解常系數(shù)方程的LDU法是將對方程YZ=E中Z陣的求解轉(zhuǎn)換為分別對n個(gè)Zk陣求解，含前代和回代計(jì)算；要求解2n個(gè)中間矩陣Wk、Hk；且無法使用Z陣元素對稱性，必須求解整個(gè)Z陣元素，計(jì)算速度受到很大影響。

至今沒有文獻(xiàn)用求解變系數(shù)方程的約當(dāng)法求解常系數(shù)方程YZ=E，但若適當(dāng)應(yīng)用約當(dāng)法的特點(diǎn)，則可實(shí)現(xiàn)該功能。

如前所述，對[YE]陣分段對稱約當(dāng)消元后得[Y(n)″E(n)″]陣，此時(shí)Y(n)″=E，E(n)″陣就是Z陣的對角元素及下三角元素，按對稱性可得Z陣上三角元素[14]。

4 例分析

分別用高斯法、約當(dāng)法、新方法求取IEEE-57、-118、-300節(jié)點(diǎn)系統(tǒng)Z陣，計(jì)算時(shí)間比較如表1。

根據(jù)表1可以看出，高斯法比約當(dāng)法快約5～13%，隨著系統(tǒng)規(guī)模加大逐漸減小，但約當(dāng)法始終沒有速度優(yōu)勢；新方法比約當(dāng)法快約66%，與系統(tǒng)規(guī)模關(guān)系不大；新方法比高斯法快約62～64%，也隨著系統(tǒng)規(guī)模加大而逐漸增加。因此與高斯法、約當(dāng)法相比，新方法可大幅度提高求解電力系統(tǒng)線性方程組或Z陣的計(jì)算速度。

表1 求取Z陣計(jì)算時(shí)間比較

t1、t2、t3：約當(dāng)法、高斯法、新方法的計(jì)算時(shí)間。

t2/t1、t3/t1、t3/t2：高斯法與約當(dāng)法、新方法與約當(dāng)法、新方法與高斯法計(jì)算時(shí)間的百分比。

5 結(jié)語

本文提出分段對稱反向約當(dāng)法，并將其用于求解常系數(shù)方程的Z陣，包括構(gòu)建特殊增廣陣；對Y陣下三角元素采用正向消元及對稱計(jì)算；對Y陣上三角元素采用反向消元；將取倒的對角元素作為規(guī)格化因子；規(guī)格化或?qū)ι舷氯窃叵獣r(shí)，E陣元素均僅計(jì)算其部分有效對角元素和下三角元素；Y、E陣上下三角元素均用四角規(guī)則計(jì)算而無需計(jì)算公式；綜合利用E陣元素結(jié)構(gòu)特點(diǎn)完成Z陣元素的對稱計(jì)算等等。大大簡化了傳統(tǒng)約當(dāng)法計(jì)算過程，與高斯法、約當(dāng)法相比，計(jì)算速度大大提高，并顯著提高編程效率。