快速LDU三角分解法的研究

席小青，陆节涣，庄广宇，陈 恳

（南昌大学信息工程学院，南昌 330031）

快速LDU三角分解法的研究

席小青，陆节涣，庄广宇，陈 恳

（南昌大学信息工程学院，南昌 330031）

由于传统LDU三角分解法中各因子阵元素之间关系不够清晰，导致计算过程复杂、不易理解，且存贮单元及计算量较大。为此，本文提出快速LDU三角分解法，引入的合成阵既可体现L、D、U元素关系又能大大减少存贮单元；综合应用“逐行规格化，按列消元”和四角规则方式，可无需依赖计算公式直接完成三角分解；在计算过程中改变元素的计算过程，可大大减少计算所需元素的总数。对各种IEEE节点系统编程计算，证明了本文所提方法的高效可行。该方法可用于电力系统等各个工程领域对常系数线性方程组的快速求解。

线性方程；LDU三角分解法；高斯消元；规格化；四角规则；电力系统

三角分解法主要用于对常系数线性方程组的系数矩阵A进行三角分解，其条件是A阵的前n-1阶主子式不等于零。常用的三角分解法有LR、LDU、CU3种，其中LR三角分解法介绍较多[1-14]，LDU三角分解法其次[6-15]，CU三角分解法较少[10-14]，而实际应用较多的是LDU三角分解法（LDU法）。三角分解法在电力系统阻抗矩阵的求解[16-20]、动态无功优化[21]、潮流计算[22-23]等电力系统计算方面应用广泛，因线性方程组的求解只是其中一小部分。然而，传统LDU法中L、D、U因子阵的形成虽有多种方式，但由于其计算过程依赖计算公式（公式法），均较为复杂，特别不利于编程计算。且其L、D、U因子阵元素单独存放需要更多的存贮数组，各因子阵元素只能用计算公式一次一个地按顺序求出，无法利用因子阵元素之间的关系，计算效率不高。

本文为解决上述问题，提出快速LDU法。由于快速LDU法在计算过程不用计算公式（过程法），因此简单直观。快速LDU法引入了合成阵的概念，使L、D、U因子阵元素关系一目了然并大大节省存贮单元；“逐行规格化，按列消元”使计算过程清晰明了；用四角规则分步计算L、D、U因子阵元素可大大简化计算过程并提高计算效率；根据L、D、U因子阵元素结构的特点进行计算，可减少约40%的元素计算个数。

1 传统LDU三角分解法

求解n阶线性方程组X，其系数矩阵为A，线性方程的矩阵形式为

A阵可分解为单位下三角矩阵L、对角阵D和单位上三角矩阵U的乘积，即

展开式（2）可分别得L、D、U3个因子阵，其元素计算公式为

传统LDU法其形成过程复杂繁琐、计算量大、不易理解，式（3）～式（5）用于编程计算也比较麻烦。由于传统LDU法中分别存放A、L、D、U因子阵元素需要4个数组，计算过程中各元素之间相应的关系无法体现和利用，所有的lij、dii、uij元素完全依赖式（3）～式（5）每次一个计算完成，其计算过程很难简化或扩展，可将其称为公式法。

2 快速LDU三角分解法

以四阶矩阵为例，将A阵分解为L、D、U因子阵的乘积，找出其元素之间对应的关系，即

将L、D、U因子阵相乘并与A阵各元素相比较，转换后可得L、D、U因子阵元素与A阵元素的关系如图1中LDU合成阵所示。

图1 L、D、U因子阵元素与A阵元素的关系Fig.1 Relationships among the elements in L，D and U and those in A

图1中虽然L、U因子阵的对角元素lii（=1）、uii（=1）与dii元素重合，但并不影响前代和回代计算。图1的合成阵就是快速LDU法的基础。

根据图1，可以按列消元方式对A阵进行LDU三角分解后合成阵中各元素的计算，而这些计算公式与式（3）～式（5）完全相同，这说明图1的合成阵体现了传统LDU法的全部内容，但合成阵的应用使得快速LDU法比传统LDU法的计算过程更简单易解。然而，此时的快速LDU法计算所有元素所要求的元素个数与传统LDU法相同，因此其计算速度还不存在差异。

根据LDU合成阵，可归纳快速LDU法有以下几个特点。

（1）可以按列消元法方式对A阵进行LDU三角分解，在A阵的基础上直接形成LDU合成阵，只需1个数组。

（2）lji、dii和uij元素在同一个合成阵中，虽然lii、uii与dii元素重合，但并不影响前代和回代计算。前代过程由lji和dii元素完成，回代过程由uij元素完成。

（3）合成阵中lji、uij元素在任何情况下均有规律性的对应关系：如lji=uij，且其非零元素的位置对称。因此在三角分解过程中，计算对角元以右非零的uij元素，可直接得到对角元以下非零的lji元素，从而省去对非零的lji元素的计算。

（4）以按列消元方式分步计算合成阵中lji、dii、uij元素，即对合成阵中同一行lji、dii、uij元素的计算类似于对对角元左侧的lji元素逐个进行消元。由于所有元素均是分步计算得到，因此可称其为过程法。过程法无需使用任何计算公式、简单易解，特别方便程序编写。但更重要的是在过程法中，可适当改变计算顺序而获得更高的计算效率。

根据合成阵中lji、dii、uij元素结构的特点，为进一步提高计算效率，规定计算过程如下。

（1）以“逐行规格化，按列消元”方式，先对第i行对角元以右的元素规格化得uij元素，但uij元素并不马上赋值给对角元素以下相应的lji元素，或者lji′元素（也就是元素）暂不除该列的对角元dii，即此时的lji元素只是lji′元素。

（2）对第i列的lji′元素进行消元，得其右侧的所有计算元素

（3）消元完成后再将第i行的uij元素赋值给第i列的lji元素，或用对角元dii除lji′元素得lji元素。

表1 公式法与过程法计算元素所需元素个数的比较Tab.1 Comparison of the number of required elements between the calculation formula and the proposed method

通过公式法与过程法计算元素所需元素个数的比较，在对lji′右侧的元素进行消元计算时所用的元素数量方面，上述过程可省去约40%所需元素个数。大大提高计算效率。这种计算方式只有在快速LDU法中容易实现，而在公式法中不易。

3 四角规则

尽管快速LDU法的应用非常简单，但在实际编程过程中应用式（3）～式（5）并不方便。故提出直接用LDU合成阵元素进行计算的四角规则。采用四角规则对图1中LDU合成阵各元素的计算过程进行分析。

1）第1次计算

第1次计算是对图2中阴影部分进行计算，过程如下：

（1）对第1行元素规格化，得到第1行的u1j元素，而对角元以下的元素仍为lj1′；

图2 第1次计算Fig.2 First calculation

（2）对第1列的lj1′元素消元，得其右侧所有计算元素，但除d22外，其余元素均未完成计算，即还未转换成相应的lji、dii、uij元素。

2）第2次计算

第2次计算是对图3中阴影部分进行计算，过程如下：

（1）将第1行的u1j元素赋值给第1列相应的lj1元素，或用d11除lj1′元素得相应的lj1元素；

（2）对第2行元素规格化，得到第2行的u2j元素，而对角元以下的元素仍为lj2′；

（3）对第2列的lj2′元素消元，得到其右侧所有计算元素，但除d33外，其余元素均未完成计算。

图3 第2次计算Fig.3 Second calculation

3）第3次计算

第3次计算对图4中阴影部分进行计算，过程如下：

（1）将第2行的u2j元素赋值给第2列相应的lj2元素，或用d22除lj2′元素得相应的lj2元素；

（2）对第3行元素规格化，得到第3行的u34元素，而对角元以下的元素仍为l43′；

（3）对第3列的l43′元素消元，得到d44。

图4 第3次计算Fig.4 Third calculation

4）第4次计算

第4次计算将第3行的u34元素赋值给第3列的l43元素，或用d33除l43′元素得l43元素，完成计算。

图5 第4次计算Fig.5 Fourth calculation

根据图2～图5，合成阵各元素的计算过程可简化如下。

步骤1对第i行元素规格化得uij元素。

步骤2对第i列的lji′元素消元，得其右侧的所有计算元素

步骤3将第i行的uij元素赋值给第i列的lji元素，或用对角元dii除lji′元素得lji元素。

步骤4依次循环，完成所有计算。

步骤2中参加计算的元素在合成阵中的位置始终如图6所示。图中dii为对角元素；uik为交叉元素，位于对角元素同行以右；lji′（lji）为消元元素，位于对角元素同列以下，括号外（内）为赋值前（后）的值；为计算元素，位于消元元素所在行右侧与交叉元素所在列的交互点上，分步完成计算后根据其所在位置转化成相应的dii、lji′（lji）、uik元素，括号外（内）为计算前（后）的值。

图6 参加计算的元素在合成阵中的位置Fig.6 Positions of the calculated elements in the synthetic matrix

根据图6可以看出，在合成阵元素的计算过程中，每次与计算有关的4个元素均在不同矩形的4个角上，因此将这种计算方法称为四角规则。由于计算等式为即实际计算中只用除对角元外的3个元素，所以可将该计算过程看成是规格化后的消元计算。

由于四角规则的计算结果与式（3）～式（5）的结果完全一致，因此形成LDU合成阵的过程根本无需使用式（3）～式（5），而用四角规则以“逐行规格化，按列消元”方式可更快、更简单地计算出LDU合成阵的所有元素，并使三角分解过程计算程序的编写变得极为简单直观。

4 计算举例

表2给出公式法与过程法对IEEE-30、IEEE-57、IEEE-118节点系统的复数导纳矩阵Y在不考虑元素稀疏性时进行三角分解所需时间的比较。表中t1为公式法对A阵进行三角分解所需时间；t2为过程法对A阵进行三角分解所需时间，消元完成后对lji′元素进行lji=lji′/dii计算；t3为过程法对A阵进行三角分解所需时间，消元完成后直接将uij元素赋值给lji元素。

从表2可以看出，对IEEE-30、IEEE-57、IEEE-118节点系统，过程法均比公式法快约55%，但消元计算后直接将uij元素赋值给lji元素所需时间比对每个lji′元素进行lij=lji′/dii计算所需时间还要少2%～4%，主要由于除法计算比赋值语句需要更多的计算时间。

表2 公式法与过程法对三角分解法所需时间的比较Tab.2 Comparison of the required time for the triangular factorization algorithm between calculation formula and process method

5 结 语

本文针对传统LDU法中的诸多不便，提出快速LDU法。在快速LDU法中，引入了合成阵的概念，将L、D、U因子阵元素直接计算后存放在原A阵数组中。合成阵清晰地表明了lji、dii和uij元素关系，应用“逐行规格化，按列消元”方式，可简单地用四角规则分步计算L、D、U因子阵的元素，而无需使用繁琐的计算公式。改变元素的计算过程，可大大减少计算所需元素的总数，从而大大提高编程及计算效率。对IEEE-30、IEEE-57、IEEE-118节点系统编程计算，过程法均比公式法快约55%，且快速LDU法为进一步简化各种三角分解法奠定了良好的基础。

[1]关根泰次.电力系统潮流计算[M].日本：电气书院，1971．

[2]陈珩.电力系统稳态分析[M].3版.北京：中国电力出版社，2007．

[3]王祖佑.电力系统稳态运行计算机分析[M].上海：水利电力出版社，1987．

[4]王艳天（Wang Yantian）.解线性方程组的LU分解法（LU decomposition method for solving linear equations）[J].科技创新导报（Science and Technology Innovation Herald），2009（4）：245-245．

[5]Hadi Saadat.Power System Analysis[M].New York：Mcgraw-Hill International Editions，1999.

[6]西安交通大学，清华大学，浙江大学等.电力系统计算[M].北京：水利电力出版社，1978．

[7]曾祥金，吴华安.矩阵分析及其应用[M].武汉：武汉大学出版社，2007．

[8]周全仁，张清益.电子计算机在电网计算机中的应用—电网计算与常用程序[M].湖南：湖南科学技术出版社，1983．

[9]陈永琳.电力系统继电保护的计算机整定计算[M].北京：水利电力出版社，1994．

[10]钱焕延，赵晓彬.计算方法[M].西安：西安电子科技大学出版社，1990．

[11]何仰赞，温增银.电力系统分析[M].3版.武汉：华中科技大学出版社，2002．

[12]周杰.矩阵分析及其应用[M].四川：四川大学出版社，2008．

[13]李庆杨，王能超，易大义.数值分析[M].5版.北京：清华大学出版社，2008．

[14]邱晓燕，刘天琪.电力系统分析的计算机算法[M].北京：中国电力出版社，2009．

[15]张伯明.高等电力网络分析[M].2版.北京：清华大学出版社，2007．

[16]冯天民，刘宝柱，鲍海（Feng Tianmin，Liu Baozhu，Bao Hai）.一类含CCCS网络形成节点阻抗矩阵的新算法（A novel algorithm for building Z-matrix of electric power network including CCCS）[J].电工技术学报（Transactions of China Electrotechnical Society），2009，24（2）：139-144．

[17]杨美佳，刘宝柱（Yang Meijia，Liu Baozhu）.针对大量接地支路电网形成节点阻抗矩阵的改进算法（An improved algorithm for forming Z-matrix of power grid containing large amount of grounded-branches）[J].电力系统保护与控制（Power System Protection and Control），2010，38（22）：161-165．

[18]朱永利，宋少群，冯建衡（Zhu Yongli，Song Shaoqun，Feng Jianheng）.互联电网节点阻抗阵实时修改与边界等值化简的并行计算方法（A parallel algorithm of realtime modification and boundary equivalents of impedance matrices of interconnected networks）[J].电工技术学报（Transactions of China Electrotechnical Society），2007，22（9）：167-173．

[19]任丽娜，陈军（Ren Lina，Chen Jun）.双复数阻抗矩阵节点导纳法实现低压短路电流计算（Double complex impedance matrix node admittance method to achieve low voltage circuit current calculation）[J].电气应用（Electrotechnical Application），2009，28（1）：56-61．

[20]乐全明，郁惟镛，杜俊红（Yue Quanming，Yu Weiyong，Du Junhong）.一种形成节点阻抗矩阵的改进算法（An improved novel algorithm for building Z-matrix）[J].中国电机工程学报（Proceedings of the CSEE），2005，25（2）：34-39．

[21]赖永生，刘明波，陈燕梅（Lai Yongsheng，Liu Mingbo，Chen Yanmei）.大规模电网的动态无功优化算法（Algorithm for dynamic reactive power optimization problem in large power grid）[J].电力系统及其自动化学报（Proceedings of the CSU-EPSA），2012，24（5）：7-12．

[22]聂永辉，肖白，刘凤兰（Nie Yonghui，Xiao Bai，Liu Fenglan）.电力系统最优潮流新模型及其内点法实现（New optimal power flow model and its solution by using nonlinear interior point method）[J].电力系统及其自动化学报（Proceedings of the CSU-EPSA），2014，26（11）：53-57．

[23]初壮，于群英，李笑薇（Chu Zhuang，Yu Qunying，Li Xiaowei）.求解含小阻抗支路配电网潮流的牛顿法（Newton method for solving power flow of distribution networks with small impedance branches）[J].电力系统及其自动化学报（Proceedings of the CSU-EPSA），2016，28（9）：36-41.

Study on Fast LDU Triangular Factorization Algorithm

XI Xiaoqing，LU Jiehuan，ZHUANG Guangyu，CHEN Ken
（School of Information Engineering，Nanchang University，Nanchang 330031，China）

Since the relationship among the elements in factor matricesL，DandUis not clear，the calculation process of the conventionalLDUtriangular factorization algorithm is complicated and hard to understand.Besides，the number of storage units and the corresponding calculation load become larger.To solve the above problems，a fastLDUtriangular factorization algorithm is presented in this paper.The introduced synthetic matrix can reflect the relationship among the elements inL，DandU，and reduce the number of storage units obviously.By applying the principle of“Normalization by rows and elimination by columns”as well as the Four-angle Rule，triangular factorization can be completed without relying on the calculation formulas.Moreover，with a new calculation process，the total number of required elements can be reduced obviously.The proposed algorithm is proved to be efficient and practical by programming and calculation on different IEEE node systems.The proposed algorithm can be used to fast solve the constant coefficient linear equations in power systems and other engineering fields.

linear equation；LDUtriangular factorization algorithm；Gaussian elimination；normalization；four-angle rule；power system

TM315

A

1003-8930（2017）10-0118-05

10.3969/j.issn.1003-8930.2017.10.020

2016-03-09；

2017-06-29

陈 恳(1956—)，男，博士，教授，研究方向为电力系统稳态分析及优化。E-mail：2494337178@qq.com

席小青（1991—），女，硕士研究生，研究方向为电力系统稳态分析及优化。Email：1164565914@qq.com

陆节涣（1991—），男，硕士研究生，研究方向为电力系统稳态分析及优化。Email：357447051@qq.com

庄广宇（1991—），男，硕士研究生，研究方向为电力系统稳态分析及优化。Email：1615206374@qq.com