對稱CU三角分解法及其應用

陳 懇，魏藝君，程思潔，魏倩文，丁 戈

（1.南昌大學信息工程學院，江西 南昌 330031；2.國網江西省電力公司上饒供電公司，江西 上饒 334000）

LR、LDU、CU三角分解法［1-6，15-17］與因子表法的計算原理類似，均用于常系數線性方程AX=F中的求解。由于因子表法和各種分解法計算過程不同，其計算效率和計算速度有所不同，然而因子表法和各種分解法的計算過程或計算方式均存在不合理情況。如因子表法中元素計算一般采用按行消元、逐行規格化［1-3］，而各種分解法是按從左上角到右下角計算因子陣元素或以按行方式從左到右計算各行元素［1-4］，這些計算方式均無法應用對稱矩陣元素本身的特點；因子表法中是形成因子表后再將對角元素取倒，可提高回代過程計算速度，但對其前代計

算速度沒有幫助，各種分解法中未考慮對角元素取倒，前代和回代計算速度均受影響［14-15］；因子表法中所有元素存放在同一個數組，便于應用元素之間關系，而各種分解法中由于各因子陣元素分別存放在不同數組中，根本無法利用和展示元素間的對應關系；因子表法無中間矩陣，可直接完成前代和回代計算，各種分解法需通過中間矩陣W或W和H完成前代和回代計算；因子表法和各種分解法均需依賴消元計算公式或元素計算公式完成計算等。上述問題均不利于計算過程優化和程序編寫，更不利于計算速度提高。

CU分解法比LR、LDU分解法的文獻介紹和應用均較少［1-6］。但實際上，CU分解法的計算過程類似于按行消元、逐行規格化的高斯消元法［4］，且其計算速度與LR、LDU分解法相比有明顯優勢［16］。然而，由于CU分解法中C、U陣元素分開存放、元素計算需依靠計算公式、沒有考慮對稱矩陣的特點及對角元取倒，其計算效率遠未達到最佳。文獻［17］的快速CU分解法雖計算速度有一定程度提高，但仍然遠未達到其應有的計算效率，而文獻［18］只是擴展了CU分解法的應用范圍，對計算速度提高無幫助。

本文針對CU分解法中存在的問題及其計算過程的特點，提出對稱CU分解法，其中包括：用可展示各元素之間關系的CU合成陣同時存放c、u元素；用四角規則按行分步計算u元素而無需應用元素計算公式，直接根據對稱性按列賦值得到c元素而無需計算；對對角元素提前取倒以減少除法計算等。將因子表法、CU分解法、對稱CU分解法用于計算極坐標PQ分解法潮流，并比較其計算速度。

1 CU分解法及其存在的問題

CU分解法將系數矩陣A分解為下三角矩陣C和單位上三角矩陣U的乘積，即A=CU。需三個數組分別存放A、C、U陣元素，元素之間關系不清。其元素計算順序為c元素按行、u元素按列或均按行方式根據式（1）～（2）每次計算一個完整的c或u元素。

式（1）～式（2）表明，c、u元素計算過程復雜，不適合編程，且無論采用c元素按行、u元素按列或均按行方式都無法實現c、u元素對稱計算。此外，式（2）中用cii元素做除數，勢必影響C、U因子陣前代過程計算速度。

對方程AX=F進行CU分解時需將對原方程的求解轉換分別對CW=F和UX=W二個方程求解，W陣前代和X陣回代過程分別展開如式（3）～（4）。

式（3）中求取W陣也始終用cii元素做除數，同樣影響形成C、U因子陣后對后續F陣前代過程的計算速度。

2 合成陣及元素的拆分和分步計算

要解決CU分解法中上述問題，首先要建立可清晰地體現c、u元素關系的CU合成陣。以四階A44陣為例，由于A44=C44U44，根據式（1）～（2）可將C44、U44陣與A44陣元素關系均在CU合成陣中展示如式（5）。如式（5）。由于C陣中對角元素為cii，U陣中對角元素為1，而式（5）人為地將C、U陣合在一起，沒有U陣對角元素，因此CU合成陣與A陣之間并不滿足數學關系，但并不影響前代和回代計算。此外CU合成陣可直接在A陣數組中形成，而無需再占用其它數組。

根據式（5）可看出以下幾點：①計算c、u元素可看成是分步對其所在行左側的c元素進行類似含規格化的按列消元，計算完成后u元素均要除以所在行的對角元素cii。②對稱矩陣中c、u元素是按cii比例對稱。③假設分步計算中未完成計算的元素為計算第i行uij元素的前一步元素，其除以對角元素cii后得uij元素，此時的元素就是與其對應的cji元素。因此可在得到元素后，將其直接賦值給對應的cji元素，從而省去所有下三角cji元素的計算。④對cji元素賦值后，將元素除以對角元素cii可得uij元素。因此利用式（5）不但可拆分、并分步計算c、u元素，還可利用c、u元素的比例對稱關系而大大簡化計算，從而提高計算速度。

由于計算uij元素的最后一步均為用元素除以所在行的對角元素cii，導致程序中大量除法運算而影響分解速度。因此可在最后計算第i行的uij元素前先將cii元素取倒，再將元素賦值給對應的cji元素，然后用（1／cii）乘以元素得到uij元素，從而幾乎可免去除法運算。

因此可將式（5）的合成陣寫成式（6）。與因子表法中形成因子表后再將對角元取倒、僅省略其回代計算過程的除法運算方式相比，本方法可幾乎完全省略前代和回代計算過程的全部除法運算，而式（2）～（3）的計算方式很難消除除法運算。

將式（1）～ （2）的計算式拆分為式（7）～ （10）。

此時式（7）～（10）已經將式（1）～（2）每次計算一個完整的c或u元素的方式拆分成僅對對角元素和上三角元素cii元素的分步計算，計算出元素直接賦值得到cji元素，再對元素除以所在行的對角元素cii得到uij元素。根據式（7）～（10）中c、u元素的拆分，和式（6）元素的分步計算，不但可實現c、u元素的對稱計算，還可省略大量的除法運算，式（3）中的除法運算也可免去。

3 四角規則及對稱CU分解法

3.1 CU分解法中應用四角規則

先不考慮矩陣對稱性進行傳統CU分解。假設c、u元素分步計算過程為：┄c?、u′?→c″、u″→c′、u′→c、u，當對n階對稱的CU合成陣完成了其第k行以上和第k列及以左c、u元素的分步計算，將繼續對第k行及以下和第k列以右的c、u元素繼續分步計算時所對應的CU合成陣的簡化矩陣如式（11）。為簡化分析，同時給出傳統法中式（11）和對稱法中式（12）相關元素定義如下：

cii（1／cii）：對角元素（對稱法取倒所得）；

u′i，i+1、u′ip、u′in（ui，i+1、uip、uin）：交叉元素（傳統法除以對角元素所得，對稱法乘以對角元素倒數所得），對角元素同行以右；

ci+1，i、cpi、cni（ai+1，i、api、ani）：消元元素（對稱法賦值前），對角元同列以下；

c?i+1，i+1～u?i+1，n、c?p，i+1～u?pn、c?n，i+1～u?nn（c″i+1，i+1～u″i+1，n、c″p，i+1～u″pn、c″n，i+1～u″nn）：計算元素的前值或新值（中間值或終值，僅給出傳統法的狀態），消元元素所在行與交叉元素ui，i+1、uip、uin所在列的交互點。

如將式（1）～（2）改寫成分步計算，則其文字表述為：“以對角元素為參考元素，計算元素的新值c″ij（u″ij）=其前值c?ij（u?ij）-消元元素c＊i×交叉元素ui＊”。

因此在計算計算元素的新值前，先將所有交叉元素u′ij元素除以其所在行對角元素cii得uij元素。再參照式（11）根據該文字表述可直接寫出相應計算元素c″i+1，i+1～u″i+1，n、c″p，i+1～u″pn、c″n，i+1～u″nn該步的計算結果。

上述計算結果表明用合成陣計算c、u元素時，每次與計算有關的四個元素均在矩形的四個角上，因此將其稱為四角規則。四角規則的計算結果與式（1）～（2）完全相同，因此根據CU合成陣可無需式（1）～（2）而直接用四角規則以類似“逐行規格化，按列消元”的方式計算CU合成陣的所有元素。這種計算方式使CU分解過程極為簡單直觀，特別利于程序編寫。

式（11）計算過程表明，由于未利用c、u元素的比例對稱關系uik=cki／cii分步計算c、u元素，因此必須計算所有c、u元素。這與式（1）～（2）的計算過程幾乎完全相同，只是將式（1）～（2）每次計算一個完整c或u元素變成分步計算，對計算速度提升影響不大。

3.2 對稱CU分解法中應用四角規則

考慮矩陣對稱性進行對稱CU分解的CU合成陣簡化矩陣如式（12），與式（11）相比，所有下三角元素ai+1，i、api～ap，p-1、ani～an，n-1等均未進行任何計算，仍為A陣值。式（12）的計算過程可分為二步：

（1）將第i行元素u′i，i+1、u′ip、u′in對應賦值給第i列的ai+1，i、api、ani得消元元素ci+1，i、cpi、cni；將對角元素cii取倒；分別乘以u′i，i+1、u′ip、u′in得ui，i+1、uip、uin；對第i列ci+1，i、cpi、cni消元，僅計算消元元素所在行的對角元素和上三角元素c″i+1，i+1～u″i+1，n、c″pp～u″pn、c″nn等，此時所計算的元素遠比式（11）中少。計算完成后，除c″i+1，i+1=ci+1，i+1為終值外，其余元素仍然均為未完成計算的中間值，而所有的下三角元素ap，i+1～an，i+1、anp等仍為A陣值。

（2）將第i+1行元素u″i+1，p、u″i+1，n對應賦值給第i+1列ap，i+1、an，i+1得cp，i+1、cn，i+1；將ci+1，i+1取倒；分別乘以u″i+1，p、u″i+1，n得ui+1，p、ui+1，n。然后對第i+1列的cp，i+1、cn，i+1消元，并依此循環。

式（12）計算過程表明，如果利用c、u元素比例對稱關系uik=cki／cii分步計算c、u元素，則只須計算相應行的對角元素ckk和上三角ukp元素，所有下三角元素cpk是逐列根據未乘以對角元倒數的上三角ukp元素對應賦值得到。因此對稱CU分解法可省略所有下三角元素計算，可大大提高計算速度。

4 算例分析

分別用因子表法、CU分解法、對稱CU分解法求解IEEE-30～-118系統極坐標PQ分解法潮流，收斂判據為ε≤10-5。其中，形成B′陣時考慮參數r影響，忽略c、k、bT影響，形成B″陣時考慮參數c、k、bT影響，忽略r影響［8-13］。計算時間比較如表1所示。

表1 極坐標PQ分解法潮流計算時間比較

t1、t2、t3：因子表法、CU分解法、對稱CU分解法用于極坐標PQ分解法潮流的分解過程計算時間。

t11、t22、t33：因子表法、CU分解法、對稱CU分解法用于極坐標PQ分解法潮流的總計算時間。

根據表1中可以看出：

（1）CU分解法與因子表法相比，前者分解過程計算時間和總潮流計算時間均慢于后者約1%和4%～7%，說明因子表法求解常系數方程更具優勢。

（2）對稱CU分解法與因子表法相比，前者分解過程計算時間和總潮流計算時間均快于因子表法約23～50%和2%～19%，說明對稱CU分解法遠優于因子表法。

（3）對稱CU分解法與CU分解法相比，前者分解過程計算時間和總潮流計算時間均快于因子表法約24～50%和6%～25%，同樣說明對稱CU分解法遠優于因子表法。

5 結論

通過針對CU三角分解法中元素對應關系不清、對角元素運算處理不當、元素計算過程固化且需用計算公式、程序編寫效率低、計算速度不理想等問題，提出對稱CU三角分解法。新方法引入CU合成陣；拆分c、u元素計算過程，并用四角規則分步計算對角元素和上三角元素而無需計算公式，下三角元素直接賦值而無需計算；改變對角元素的運算方式以免去除法運算。新方法大大簡化了CU三角分解法的計算過程、幾乎免去了程序中的除法計算，并可大大提高程序編寫效率。將新方法、因子表法和CU三角分解法用于求取IEEE-30～-118節點系統的極坐標PQ分解法潮流，無論是分解過程還是整個潮流計算時間，新方法的計算速度遠遠快于后者。新方法同樣可用于各工程領域常系數方程的快速求解。