对称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分解法潮流，无论是分解过程还是整个潮流计算时间，新方法的计算速度远远快于后者。新方法同样可用于各工程领域常系数方程的快速求解。