1、版权所有,二 线性稀疏矩阵方程的求解,版权所有,主要内容,高斯消去法,三角分解法,节点编号顺序的优化,版权所有,直接解法,线性方程组可以用直接解法，计算实践表明，对电力系统来说它很有效。 电力系统中常见的大型线性方程组的系数矩阵十分稀疏，直接解法的计算速度很快。 与迭代法相比，没有收敛性问题。,版权所有,一、高斯消去法 1、按列消元按行回代的算法,设有n阶线性方程组a11x1+a12x2+a1nxn=b1 a21x1+a22x2+a2nxn=b2. (21) an1x1+an2x2+annxn=bn,或缩记为：AX=B (22),版权所有,求解的具体步骤如下：,（1）若a110，由（21）第1
2、式解出x1=b1-(a12x2+a1nxn)/a11,a11x1+ a12x2 + a1nxn =b1a22(1)x2+a2n(1)xn=b2(1). (23)an2(1)x2+ann(1)xn=bn(1),代入第2至第n式消去x1 ,有,式中，,i=2,3,n;j=i,i+1,n,版权所有,(2)若a22(1)0，由（23）第2式解出,代入第3至第n式以消去,便得:,式中，,(i=3,4,n；j=i,i+1,n),版权所有,由此只要 ,消元过程就可继续，在作完第k步消元后，原方程组将变为,(24),（i=k+1,. ,n； j=i,i+1, n+1）,式中 ,(25),版权所有,如果不满足，
3、可将尚待继续消元的那部分方程式重新排列次序，使第k个方程中 的系数不为零即可。,由于 与 的算法相同，若把 记为 ，在利用公式（25）时可把列标 j 一直取到n+1。以求取,只要： ，消元就可以继续,版权所有,(26),式中 ， （i=1,2, ,n ; j=i,i+1,n+1）,(27),经过n-1次消元，最后得到的方程为。,版权所有,消元的结果是把原方程组（21）演化成系数矩阵呈上三角形的方程组 （26）。这两组方程组有同解 。,回代过程求解X,版权所有,(28),回代过程:,按列消元对应于网络变换时的消去节点，物理概念清晰,版权所有,例21 按列消元按行回代的高斯消去法,按列消去,按行回
4、代: ( 1)由第三式解得x3=4 (2) 由此及第二式解得x2=2 (3) 由x3、x2此第一式解得x1=1,版权所有,电力系统计算中，高斯消去法的另一种常用计算格式是按行消元逐行规格化的算法。 具体做法如下：,（1）若 则以 乘方程组（21）中的第1式，使之规格化，得到 （29）式中, （j=2,3,n+1）,2按行消元并进行规格化的算法,版权所有,（2）对方程组（21）中的第2式作运算。首先进行消元，用-a21乘（29）全式，再同式（21）的第2式相加，便得到 ：,式中，,假定 ，用 去乘上式规格化，便得,这样，得到了经过消元和规格化处理的第2个方程式。,（j=3，4， ，n+1）,式中
5、,版权所有,消元过程逐行进行。对原方程组（21）中的第i个方程式的演算包括，先作i-1次消元，利用已完成消元和规格化处理的i-1个方程式依次消去 ,然后作一次规格化计算，使 的系数变为1。以k代表消元次数，逐次消元计算通式为 ：,（210）,作完i-1次消元后，规格化计算公式为,（j=i+1,i+2,.,n+1） （211）,版权所有,一般，经过i-1步按行消去运算，增广矩阵变成,版权所有,按照上述步骤，对方程组（21）的全部方程式作完消元和规格化演算，便得到了以下的方程组,（212）,式中的系数表达式为：,（i=1,2,n-1;j=i+1,i+2, ,n+1）,（213）,版权所有,利用方程
6、组 （212）,通过回代计算，即可求得全部的未知变量，其计算通式为,（i=n,n-1,.,1 ) (2-14),注意（212），（213），（214）与式（26），（27），（28）本质上一致。,版权所有,例22 按行消元逐行规格化的高斯消去法,S1. 规格化第一行,S2. 一、二行相消,S3. 规格化第二行,S4.一、三行相消,S5. 二、三行相消,S6. 第三行规格化,版权所有,最后得到：,其中，依次取1/2，3，2/5，5，-23/2，-1/12为运算因子。 由后向前取虚线上三角中元素进行回代运算,S1. 取1,S2. 取1/2,S3. 取3/2,版权所有,因子表,版权所有,高斯消去法的
7、每一步演算都相当于进行矩阵的初等变换。以按列消元的算法为例，第一步消元时所用的初等矩阵为,其中,（i=2，3，,n）,二、三角分解法 1、将非奇异方阵A分解为单位下三角矩阵L和上三角矩阵R的乘积,版权所有,用 左乘式（22）的两端，将结果展开便得到方程组 (23） 以后的每一步消元都是对上次变换所得的结果再作一次初等变换。所用的变换矩阵都是单列单位下三角矩阵，第k步消元所用的矩阵为,（215）,（i=k+1, ,n）,（216）,版权所有,依次作完n-1次变换后，便得到方程组 (26)。若将方程组 (26）的系数矩阵记为R，则其元素为,(i=1,2, ,n;j=i,i+1, ,n),从演算过程
8、可知,（217）,因为初等矩阵非奇，故有,根据单列单位下三角矩阵的性质可知,版权所有,这样，便得到,A=LR,（218）,非奇方阵A被表示为矩阵L和R的乘积，这两个三角矩阵称为A的因子矩阵.,版权所有,从计算公式可见，要使分解得以进行下去，必须有 。为满足这个条件，要求矩阵A的各阶主子式都不等于零。如果矩阵A非奇异，通过对它的行（或列）的次序的适当调整，这个条件是能满足的。,（219）,利用公式 (25）和 (27），可确定两个因子矩阵的元素计算公式：,版权所有,LF=B RX=F,或者展开写成,（220）,（221）,将A=LR代入线性方程组(22），便得LRX=B.这个方程又可以分解为以下
9、两个方程,版权所有,这两组方程式的系数矩阵都是三角形矩阵，求解极为方便。先由方程组（220）自上而下地依次算出 其计算通式为,方程组（221）与式（26）一致，求解属于回代过程，类似式（28）。,（222）,版权所有,例23 三角分解法(LR),利用DOLITTLE分解,版权所有,三角分解后，求解过程如下：,LF=B RX=F,版权所有,2、将非奇方阵A分解为单位下三角矩阵L、对角线矩阵D和单位上三角矩阵U的乘积,如果A非奇，则上三角矩阵R的对角线元素都不等于零。矩阵R又可分解为对角线矩阵D和单位上三角矩阵U的乘积，即R=DU，或展开写成,比较两方的对应元素可得,由此可知,（223）,（224
10、）,版权所有,这样便得A=LDU,（225）,这种分解称为方阵A的一种LDU分解。若A的各阶主子式均不为零，则这种分解是唯一的。,版权所有,（226）,利用公式（219），计及式（223），可得因子矩阵的元素表达式如下,版权所有,例24 三角分解法(LDU),利用DOLITTLE分解,利用LDU分解,版权所有,将式（225）代入式（22），可得LDUX=B,（227）,这个方程又可分解为以下三个方程组,即（220）,版权所有,根据式（222）和（213）可知 。因此，求解这组方程相当于对经消元变换后的右端常数作一次规格化演算。,（229）,由此可得,方程组DHF可展开为,（228）,版权所有,
11、当A的各阶主子式均不为零时，根据分解的唯一性，应有,因此,（230）,方程组UX=H展开后即是式（212） 若A为对称矩阵，则应有：,版权所有,利用公式（226），计及 ，便得各因子矩阵的元素表达式为,由于三角矩阵U和L互为转置，只需算出其中的一个即可。,（231）,版权所有,若令LD=C，则矩阵C仍为下三角矩阵，其元素为,这样便得: A=CU,这种分解亦称为Crout分解。,（232）,（233）,3、将非奇方阵A分解为下三角矩阵C和单位上三角矩阵U的乘积,版权所有,利用公式（226），计及式（232），可得因子矩阵的元素表达式如下,（234）,不难验证，用Crout分解求解线性方程组的算法
12、，相当于按行消元逐行规格化的高斯消去法,版权所有,例25 三角分解法(CROUT),利用CROUT分解(2-34),版权所有,对比按行消元逐行规格化的高斯消去法因子表：,版权所有,三角分解后，求解过程如下：,CH=B UX=H,版权所有,4、因子表及其应用,网络方程需要求解多次，每次只是改变方程右端的常数向量B，使用的系数矩阵A相同 对线性方程组的系数矩阵A进行三角分解，所得的下三角矩阵用于消元运算，而上三角因子矩阵则用于回代运算。 对于需要多次求解的方程组，可以把三角形因子矩阵的元素以适当的形式贮存起来以备反复应用。 不同的形成因子表方法，对应求解过程（公式）有所不同。,版权所有,作LR分解
13、时，可把因子矩阵L和R的元素排列成：,作LDU分解时，把各因子矩阵的元素排列成：,矩阵L和U的对角元素都是1，不必存放。 上述因子矩阵的元素正好占据矩阵A的对应元素位置。因此，以上几种排列格式都可以称为矩阵A的因子表。,作Crout分解时，把因子矩阵C和U的元素排列成：,对矩阵A，作三种因子分解时的因子矩阵元素,版权所有,以按行消元逐行规格化的算法为例，这种算法需要保留矩阵C和U的元素。由于对角线元素Cii在计算过程中都作为除数出现，在计算机中乘法要比除法节省时间。因此，在实际使用的因子表中，对角线位置都是存放Cii的倒数1/Cii。由于对称矩阵的因子矩阵L和U互为转置矩阵，在因子表中保留上三
14、角部分（或下三角部分），而对角线位置则存放矩阵D的对应元素的倒数。(常用),版权所有,利用因子表(CROUT分解)：,消元回代：,（i=n,n-1,.,1 ) (2-14),（222）,版权所有,利用因子表(LDU分解)：,消元回代：,（i=n,n-1,.,1 ) (2-14),（222）,版权所有,右端的常数向量分别取为：,解：用CROUT分解形成因子表如下：,例26 用因子表求解方程组AX=B。,版权所有,先做消元运算,再做回代,版权所有,先做消元运算,再做回代,版权所有,设四节点网络的节点编号分别如图（a）(b)所示。,三、节点编号顺序的优化 1、节点编号顺序与稀疏度的关系,（a）,（b
15、）,版权所有,对应这两种编号方案的节点方程分别为：,（235）,分别进行三次、一次消元运算消去系数矩阵中第一列后，这两个系数矩阵中非零元素的分布将如下式所示。 “*”表示原非零元：“”表示消元后新出现的非零元，称注入元。,（236）,版权所有,再分别进行三次、两次消元运算，消去其中第二、第三列，得上三角矩阵中的非零元素分布如下式所：,（237）,按方案一编号时，需经六次消元进入回代； 按方案二编号时，仅需三次消元就可进入回代。 方案二回代过程也较简单。 差别关键在于消元过程中是否会出现注入元，取决于网络节点编号的顺序为保持节点导纳矩阵从而保持因子表的稀疏度，降低对存贮空间的需求、减少运算量，必
16、须尽可能优化节点编号的顺序。,版权所有,2高斯消元与消去节点的关系,先将式(236)以及相应的右端项部分展开如下：,（238a）,（238b）,版权所有,上两式中虚线以下部分其实就是由式(235)式第一行解出U1并用以消去其它各行中U1后的所得。 这一消去U1 的过程，对应用负荷移置消去节点1，进行网络简化的过程 不同点仅仅在于，采用负荷移置法时，是以全网节点电压都等于额定值的假设条件下，不以电流而以功率表示节点的注入。以图（a）及式（238a）说明上述网络简化论点。 第一步 由图(a)套用节点电流移置公式,版权所有,从而节点2、3、4注入电流分别改变为：,这就是（238a）式的右端列向量中各
17、元素,第一步电流得证。,版权所有,第二步 由图(a)套用节点导纳移置公式，消去节点1后，节点2与3、4、0之间导纳变化量分别为：,从而节点2的自导纳及互导纳分别为：（对应（238a）中第二行）,相似地，可以得系数矩阵第三、四行各元素 第二步 导纳得证。,版权所有,几个结论 以高斯消元法逐列消元，对应于以消去节点法逐个消去节点 消元过程中的注入元，在物理意义上对应于由于消去某节点而出现新的互联支路导纳。 就形成因子表而言，三角分解法与高斯消元法完全等效，而以高斯消元法逐列消元又对应于以消去节点法逐个消去节点，因此可通过考察消去节点以考察因子表的形成 基于如上关系，高斯消元后如出现注入元，该注入元
18、也将出现在三角分解后所得的上、下三角矩阵中，并将出现在所形成的因子表中。 因子表中是否会出现注入元等价于网络消去节点后是否会出现新的互联支路。,版权所有,3节点编号的最优顺序,节点编号改变会影响注入元素数目 例如，对简单五节点网络节点编号有三种方案：,方案1 方案2 方案3,分析: 以高斯消去法建立因子表 以三角分解法建立因子表,版权所有,（1）以高斯消去法建立因子表,消去节点1后,消去节点2后,消去节点3后,版权所有,版权所有,以高斯消去法建立的因子表中非零因子的分布,版权所有,利用下式后，可得因子表中非零因子，与高斯法完全一致,（2）以三角分解法建立因子表,版权所有,三种节点优化编号方法,
19、静态优化法按静态联结支路数的多少编号半动态优化法按动态联结支路数的多少编号最常用 动态优化法按动态增加支路数的多少编号,三种节点优化编号方法编号结果不同,版权所有,静态优化法按静态联结支路数的多少编号,如果编号为,消去1，2，3后，如果继续消去4，则出现新支路BD，出现注入元：,缺点是在连接支路一样时，有编号顺序问题。,如果编号如右则无问题：,版权所有,半动态优化法按各节点动态联结支路数的多少编号,先编号F，然后消去F,取E，并消去：,缺点是对环网，可能出现注入元,对放射型网络，保证不会产生注入元。,再取A，然后消去A,重复。,版权所有,例题：如下图,按静态法编号为,1，2消去后，不产生注入元
20、 3消去后，产生注入元（B，F）（6，4） 4消去后，产生注入元（B，D）（6，7） 产生两个注入元,版权所有,如果按半动态编号：,先取A（1），再取G（2）,再取C，可见仍然是两个注入元,取F(#3)，产生注入(E,D),取B(#4)，产生注入(E,C),版权所有,动态优化法按动态增加支路数的多少编号,其中：Jn表示与节点n相连的节点数。Dn表示与节点n相连的Jn个节点间原来有的支路数,消去某节点时，可能增加的支路数 等于：,如消去f时，增加的支路为1028个：,版权所有,四、稀疏向量法,参考文献： W.F.Tinney, V.Brandwajn, S.M.Chen, Sparse Vector Methods, IEEE Trans., Vol. PAS-104, No.2,1985.,版权所有,结 束,待续,
