电力网络方程求解技术

电力网络方程的求解技术

在电力系统的分析计算中，暂态分析一般关注电压和电流，电力网络模型常

为线性日勺节点电压方程；稳态分析一般关注功率和电压，其电力网络模型常为

非线性时尚方程，而非线性时尚方程也必须通过求解线性W、J修正方程才能得到

其解。因此，无论是电力系统稳态分析，还是暂态分析几乎都会波及线性方程组

的求解问题，并且线性方程组的求解往往是计算量最大的一部份工作。因此，研

究线性方程组的求解技术对电力系统分析计算有重要的意义。

线性方程组的解法可归纳为直接法和迭代法。从理论上来说，假定每一步

运算过程中没有舍入误差，直接法通过有限次运算，最终得到方程组的I解就是精

礁解。不过，这只是理想化时假定，在计算过程中，完全杜绝舍入误差是不也许

的，只能控制和约束由有限位算术运算带来H勺舍入误差的增长和危害，这样直接

法得到的解也不一定是绝对精确的。

迭代法就是用某种极限过程去逐渐迫近线性方程组精确解的措施。该措施

具有对计算机日勺存贮单元需求少，程序设计简朴、原始系数矩阵在计算过程中

不变等长处，是求解大型稀疏系数矩阵方程组的重要措施。迭代法不是用有限

步运算求精确解，而是通过迭代得到满足一定精度规定的方程组的近似解。

在数值计算历史上，直接解法和迭代解法交替生辉。一种解法日勺兴旺与计

算机的硬件环境和问题规模是亲密有关於J。一般说来，对同等规模的线性方程

组，直接法对计算机出J规定高于迭代法。对于中等规模的线性方程组，由于直

接法的精确性和可靠性高，一般都用直接法求解。对于高阶方程组和稀疏方程

组,用迭代法可防止直接法带来日勺高舍入误差。

第一节计算机在电力系统应用H勺初期，曾经由于内存容量的限制采用过迭代法

求解电力网络时线性方程式组。迭代法的致命缺陷是存在收敛性问题。自从稀

疏技术成功地在电力系统应用之后,迭代法几乎完全被所替代。但伴随电力系

统规模的迅速扩大，使得直接法很难满足在线应用的规定,规定采用并行计算技

术提高电力系统分析计算的速度。由于迭代法有很好的并行性，也许会再次得

到广泛W、J应用。

第二节由于电力网络构造出J特点,在以导纳矩阵表达日勺电力网络方程中系数矩

阵和常数矢量中非零元素非常少,这种状况下的矩阵和矢量是稀疏日勺。在与稀

疏矩阵和稀疏矢量有关的运算中,有零元素参与W、J运算是没有必要进行日勺,对零

元素W、J存储也是多出出J。因此，可以采用“排零存储”、“排零运算”W、J措施，

只存储稀疏矩阵和稀疏矢量中日勺非零元素及必要的检索信息，只取这些非零元

素来进行运算,省去对零元素日勺存储和与零元素进行日勺运算,这样可以大大减少

存储量，提高计算速度。这种作法用计算机程序来实现就是稀疏技术。它包括

了稀疏矩阵技术和稀疏矢量技术两方面。和不采用稀疏技术相比,采月稀疏技

术可以加紧计算速度几十甚至上百倍，并且对计算机的内存规定也可以大大减

少。电力系统规模越大，使用稀疏技术带来日勺效益就越明显。可以说，稀疏技术

的引入是对电力系记录算技术的J一次革命,使许多本来不能做出J电网计算可以

很轻易地实现。

第三节线性方程组的迭代解法

一、线性方程组的迭代解法的思绪

用迭代法求解线性方程组就是对方程组进行等价变换，构造同解方程组

,以此构造迭代关系式

Xa+,>=MX(A)4-g(2-1)

式中，称为迭代矩阵。任取初始矢量，代入式(2-1)中,经迭代计算得到解序

列。若解序列收敛，设的极限为，对迭代式两边取极限

limXa+,）=lim（MX（ft）4-g）

火一>00

即，是方程组日勺解，此时称迭代法收敛，否则称迭代法发散。迭代法日勺长

处是占有存储空间少，程序实现简朴，尤其合用于大型稀疏矩阵；不尽人意之处

是要面对判断迭代与否收敛和收敛速度日勺问题。

迭代式（2・1）收敛与否完全决定于迭代矩降的性质，与迭代初始值的选用无

关。可以证明迭代矩阵的谱半径

p（M）=rnax|2z|<1（2-2）

是迭代收敛口勺充足必要条件，其中是矩阵的特性根。

因此,称谱半径不不小于1的矩阵为收敛矩阵。计算矩阵的谱半径,需规定

解矩阵H勺特性值,一般这是较为繁重H勺工作。可以通过计算矩阵的范数等措施

简化判断收敛日勺工作,其中,计算矩阵的1范数和范数的措施比较简朴。（向

量的1范数等于向量元素绝对值之和，向量的范数范数等于向量元素绝对值的

最大值。矩阵日勺1范数等干矩阵列向量日勺1范数的最大值：矩阵的范数等于

矩阵行向量的1范数H勺最大值。）

网『嚅白闻（2-3）

J/=1

1网8=..力叫（2-4）

XI

式（2-3）、式（2-4）分别是矩阵1范数和范数的计算公式。可以证明,

只要迭代矩阵满足或，就可以判断迭代序列是收敛的。但要注意的是，当

或时,可以有，因此不能判断迭代序列发散。

二、在计算中当相邻两次日勺向量误差的某种范数不不小于给定精度时，则

停止迭代计算，并视为方程组口勺近似解。

三、雅可比（Jacobi）迭代法

a

0\2...

falu

。仇2…b[〃ll\\

...这

力2|0…b2n=0

B=%(2-8)

・・・・•*••

••••

?«1%…。一_SL

ann

(2-9)

得雅克比迭代式口勺矩阵形式

X(z)=BX")+g(2-10)

2、式中，为雅克比迭代矩阵。

3、雅可比迭代收敛条件

当方程组的系数矩阵具有某些性质时，可直接鉴定由它生成日勺雅可比迭

代矩阵是收敛日勺。

定理：若方程组口勺系数矩阵，满足下列条件之一，则其雅可比迭代法是收

敛口勺。

（1）为行对角占优阵，即

（2）为列对角占优阵，即

4、定理：若方程组H勺系数矩阵为对称正定阵，并且也为对称正定，则

雅可比迭代收敛。（为日勺对角线元素构成的对角线矩阵）

5、雅可比迭代算法

loopi=Vn〃形成雅克比迭代矩阵不嗒数项

loopj=1,Z-1

%=-4/%

endloop

loopj=Z+1,7z

琢=-%/%

endloop

endloop

xl={0,0，…⑼〃解矢量中间值，存放迭代前的值

工2={1,1,--,1}〃解矢量，存放迭代后的直，开始时存放解的神直

while||xl-x2||。£〃迭代循环

x\=x2

loopi=\,n

乌=&

loopj-1,/-1

x2t=x2,+%*x\.

endloop

loopj=z+l,7z

x2t=吗+%

endloop

endloop

endwhile

四、高斯•塞德尔(Gauss-Seidel)迭代法

1、高斯・塞德尔迭代格式

在雅可比迭代中，的计算公式是

靖力声产+&(2-U)

j=\；=;+1

实际上，在计算前，已经得到时值，不妨将已算出的分量直接代入迭代式

中，及时使用最新计算出的分量值。因此日勺计算公式可改为:

广)二尸)+£/"7+&(2-12)

式(2・12)称为高斯一塞德尔迭代式。

2、高斯一塞德尔迭代的收敛条件

定理:若方程组系数矩阵A为列或行对角优时,则高斯塞德尔迭代收敛。

3、定理:若方程组系数矩阵A为对称正定阵,则高斯塞德尔迭代收敛。

4、对于方程组，假如由它构造高斯-塞德尔迭代和雅可比迭代都收敛，那

么，多数状况下高斯一塞德尔迭代比雅可比迭代日勺收敛效果要好，不过状况并非

总是如此。

5、高斯一塞德尔迭代算法

loopi=\,n〃形成迭代矩阵和常数项

loopj=i,z-i

%=-%/%

endloop

loopj=i+\,n

endloop

endloop

3={0,0,...,0}//解矢量，存放迭代前的直

x2={l,l,…1}//解矢量，存放迭代后的直

while1kl-九2口〃>£//迭，弋循环

x\=x2

loopi=\,n

x2i=8i

loopj=U-l

x2f.=x2j+bg*x2-

endloop

loopj=i+\,n

x2j=x2j+"Jx2j

endloop

endloop

enclwhile

五、逐次超松弛(SOR)迭代法

1、逐次超松弛迭代格式

方程组的雅可比迭代形式，记其中是下三角矩阵，是上三角矩阵。

得高斯・塞德尔的I迭代形式:

Xa+,)=£x(z)+UXa)+g(2-13)

记，有

X"+D=X")+AX"')(2-14)

这样可视为的修正量,假如将改为加上修正量乘一种因子，迭代格改

为：

X"+i)=X")+3AX")

X'"D=X(力+o(X("+D-X'幻)

X(i+,)=X®+G(tx"卬+Cx⑶+g-X"))

整顿得

Xa+,)=(1-MX")+6y(LXa+1)+UX(A)+g)(2-15)

2、这里为修正因子,称为松弛因子，而式(2-15)称为松弛迭代c

3、逐次超松弛迭代收敛的条件

定理:逐次超松弛迭代收敛的必要条件。

定理:若为正定矩阵，则当忖、逐次超松弛迭代恒收敛。

4、以上定理给出了逐次超松弛迭代因子口勺范围。对于每个给定的方程组,

确定究竟取多少为最佳，这是比较困难的问题,

5、一般，把W、J迭代称为亚松弛迭代，把的迭代称为高斯・塞德尔迭代，而

把的迭代称为松弛迭代。

6、逐次超松弛迭代算法

loopi=l,〃〃形成迭代矩阵和常数项

gj二)"因

loopj-1,M

endloop

loopj=i+\,n

%=_%/%

endloop

endloop

xl={(),()/..,(“//解矢量，存放迭代前的值

x2={l,l,…,1}〃解矢量,存放迭代后的值

6=l-co

while卜1-x211P>£〃迭弋循环

x\=x2

loopi=l,/i

loopj=I,/-1

x2f=x2j+①*bi)*x2j

endloop

loopy=Z+1,7?

<=9+0*bi)*x2j

endloop

endloop

endwhile

第四节线性方程组的直接解法

线性方程组可以用消去法直接求解，虽然是很古老的J措施,不过计算实践表明，对电力系统来说是很有

效的。这是由于电力系统中常见的大型线性方程组H勺系数矩阵，如导纳矩阵是十分稀疏H勺，因此当充足运

用柜阵的稀疏性时，直接解法的计算速度很快。与上节简介的迭代法相比，虽然直接解法占用计算机口勺内

存量要大某些,不过它没有收敛性问题。

一、本节对消去法进行一般日勺数学描述,给出合用于电子数字计算机的体

现式,并简介它日勺常用变态形式一因子表解法和三角分解解法。

二、高斯消去法

1、高斯次序消去法的基本思想是：对线性代数方程组所对应的增广矩阵

(A|B)进行一系列“把某一行元素倍加到另一行上”的初等变换，使得

(A|B)中A日勺对角线如下日勺元素消去为0,从而使原方程组等价Ef、J转化为轻易

求解的上三角形线性代数方程组，再通过回代得到上三角形线性代数方程组的

解，即可求得原方程组日勺解。高斯消去法求解线性方程组分为两个环节：消去运

算(前代运算)和回代运算。

2、消去(前代)运算

设有线性方程组:

AX=B(2-16)

将系数矩阵和常数向量合并写成增广矩阵:

a\\a\2a\n瓦

A=[AB]=；a；b：(2-17)

消去运算有两种措施，按列消去和按行消去。首先讨论按列消去日勺过程，其

环节是：

第一步，消去第一列。

首先，把增广矩阵的第一行规格化为

1%...〃l.n⑴〃⑴1(2-18)

式中:

~(j=2,3,…#)

<%(2-19)

然后，用式(2.18)所示日勺行消去日勺第一列对角线如下各元素，成果使日勺第

2力行其他元素化为

%。)-(j-2,3,...,〃工一2,3,...,〃)

(2-20)

=b「aM”(i=2,3,...〃)

式中：上标（1）表达该元素第一次运算的成果。这时矩阵变为：

1%⑴…V4⑴

(1)(1)>(1)

a22…a2n%

A[=[AIB1]=

/2⑴…心⑴

与之对应的方程组是，它与同解。矩阵未标出日勺元素为零，下同。

第二步，消去第二列。

首先，把增广矩阵的第二行规格化为

01的3⑵…％/2)为⑵(2-21)

式中:

q⑴

旬⑵=朱(；=3,4,...,/?)

«22（2-22）

声)_嫂.

然后，用式（2-21）所示的行消去的第二列对角线如下各元素，，,成果

使日勺第3〜n行其他元素化为

4•⑵=47一项)%⑵(，=3,4,...〃;i=3,4,...M(2.23)

b："=〃⑴-a；；%?⑵(/=3,4,...,〃)

式中：上标（2）表达该元素第二次运算的成果。这时矩阵变为：

一般地,在消去第k列时要做如下日勺运算:

。=々+1,…⑼

(2-24)

aJ'=尾1)—或一%J"(J=左十1,;i=〃十1,…川

(2-25)

b‘>=b>7-娟'葭)(i=k+1,...〃)

式(2.24)称为规格化运算，式(2-25)称为消去运算。

通过对矩阵的n次消去运算，即k从1依次取到n按式(2-24),(2・25)

运算,使矩阵A对角线如下的元素所有化为零,从而得到增广矩阵

1牝⑴心⑴…《〃⑴々⑴

1%⑵…小⑵打⑵

A?=fAnBJ=1…⑸⑶打⑶(2-26)

.1b：匚

与之对应的方程组是，即

2+《2⑴％2+43ax+…%=4⑴

々+。23⑵七+…+。2”⑵％=%⑵

“七十•••+"3〃")x"=4"(2-27)

A=b;)

它与原方程组AX=B同解。

以上算法首先消去中第一列对角线下W、J元素，然后消去第二列对角线下日勺

元素，依次直到对角线下日勺所有元素都被消去为止。这种消去算法称为按列消

去法。

下面简介按行消去法，它首先消去中第二行对角线左边的元素,然后消去第

三行对角线左边日勺元素,依次直到对角线左边R勺所有元素都被消去为止。其环

节是:

第一步,首先对增广矩阵口勺第一行做规格化运算,成果为：

⑴⑴

1《2%...4I⑴〃⑴I(2-28)

式中:

(2-29)

甲=2

即

第二步，首先用式(2.28)所示的行消去的第二行对角线左边的元素，成果使

口勺第2行其他元素化为

，⑴・、、

(Ci(1)=a,-a,a(/j=2,C3),…如

2j2j2iiJ(2-30)

也⑴=%-生自⑴

这时矩阵变为：

4”⑴4⑴

6⑴*

A[=[A[BJ=

ah

a,a22nnn

然后，对增广矩阵的第二行做规格化运算,变为:

八(2)

2)

o123...aJb『(2-31)

式中:

//⑵=~7iT(J=3,4,...〃)

，增鹏与(2-32)

。22

这时矩阵变为：

-1%⑴…4〃⑴々⑴.

A；=[A2B2]=：…生〃A

■・•・・・••・••••

_/14〃2­­■%〃bn.

显然，与之对应口勺方程组是，它与同解。

第三步，首先用式(2-28)所示的行消去向第三行对角线左边的第一元素

;然后用式(2-31)所示的行消去的第三行对角线左边的第二元素；最终

用第三行对角线元素对第三行做规格化运算。

一般地,在消去第k行时要做如下的运算：

=唠T一说丁％/")(，〃=12…水一1；J=4,3,…,〃)(233)

'W»akb,「)("1,2,…,”1)-

加-1)

a「=飞(万Z+l,•••,〃)

^%(2-34)

,g)

3、通过n次消去运算，得到与式（2・26）相似的增广矩阵，和与式（2-26）

相似的同解方程。

4、回代运算

将方程组展开为:

(2-35)

可见，经消去运算后系数矩阵变为一上三角矩阵了。回代运算就是由式（2-

35）求出原方程解日勺过程。可以采用按行回代算法或按列回代算法。

按行回代以行自下而上的次序进行。其过程为：首

先由第个方程得到解：

…?（2-36）

然后，将代入第个方程，解出:

（2-37）

再将将，代入第个方程可解出：

——詈e/z-4肾./（2-38）

如此，如已得到的解分量，，…，，得出求解分量的算法

王=/?.*'-£1a(i=〃，〃一(2-39)

j=i+\

式（2-39）就是按行回代的一般公式。

按列回代的计算公式是:

XJ_牙）一斓）'〃⑴]

X

=—Un-\~

Z.2百0x2一…一n-2,n-\%(2-40)

■Vi00

5__b?一000

三、也是按，，，W、J次序依次求各位置数。

四、因子表法

1、因子表

从上节高斯消去过程可以看出①线性方程组常数项不影响系数矩阵的消去

成果；②常数项的消去运算只与系数矩阵的下三角中即将被化为1或0的元素

有关；③回代求方程的解只与消去运算完毕后的上三角元素有关。

在实际计算中，常常碰到方程需多次求解,每次仅变化常数项,而系数矩阵保

持不变。在这种状况下，如每次求解都做一次对系数矩阵和常数项完整消去运

算，很明显将会有大量反复的I运算。假如常数项变化时，只需对常数项做消去运

算就行了，这就是因子表法。

因子表就是记录线性方程组求解过程中消去（前代）和回代运算所需数据

口勺一种表格。回代运算所需数据由对系数矩阵消去运算所得口勺上三角矩阵元素

确定。为了对常数项进行消去（前代）运算，还必须记录消去（前代）过程中

所需日勺运算数据。消去（前代）过程又分为规格化和消去运行，以按列消去为

例，由式（2-24）和（2-25）可知，消去（前代）过程对常数项第个元素H勺运

算包括

/y=45_或7也⑻（k=1,2,…,"1）（2-41）

申=翼（2-42）

4T

可见，常数项的消去运算只与系数矩阵的下三角部分和常数项有关。将式

（2-31）和（2-42）中的运算因子按它们下标指定的I位置放在下三角部分，和消

去运算得到的上三角矩阵放在一起,就得到了因子表

如a⑴4?*…嫌)

422喈媚）...喈

喈⑶

“32“34…“3”(2-43)

乙⑴,，⑵,,(4)

a

4\a42...Cl4ll

••••••••••••••••••

a⑴

an\“"3...苏〃川

式（2-43）因子表中对角线元素为消去过程中规格化为1之前的对角元

素；下三角元素为消去过程中消去为0之前的下三角元素。式中，对角线及下

三角部分用于对常数项的消去（前代）运算,上三角元素用来进行回代元算。

因子表也常写成如下的形式:

八“12

，2142

hi,32(2-44)

,41,42

42

"1nn

式中

d产娟)(/=1,2,...,/?)

囱=a-p(i<j,i=1,2,…=i+

/..=")(j<3=1,2,・・・,〃，,/=1,2,...,/-1)

2、使用因子表解线性方程组

（1）形成因子表（按列消去）

loopk=\,n〃七为消去的计数

loopj=k-\-\,n

akj=akj/akk〃规格化第攵行第/列元算

loopi=A+1/〃对第Z+1行到〃行的第/列元素做消去运算

%=%-4人%

endloop

endloop

endloop

(2)对常数项做消去运算

LoopZ=l,/7〃i正在消去常数项的序号

Loopk=〃用下三角元素进行消去运算

b.=2一6也

EndLoop

bi=bj%//规格化运算

EndLoop

(3)回代运算得方程的解

Loopi=n,l〃从后往前回代

Xj=a〃赋解的初值

Loopj=i+l,n〃用已求得的解回代

xi=xi-aijxj

EndLoop

EndLoop

五、三角分解法

1、矩阵的三角分解

设方阵A可用一种下三角矩阵L与一种单位

上三角矩阵U日勺乘积表达，即:

A=LU(2-45)

展开为:

flll<712"13a\nh1%"13…

hl，22

“21a22。23a2n1423"2"

一

“31”32"33〃3rr，32，331…%”(2-46)

…/期_

_°n\an2%3…〃，工Jnl1

将式(2-46)右边两矩阵相乘，其元素应与左边矩阵对应元素H勺值相等。比

较第一行元素，得:

aw=%〃对角线元素

%2

W空=*)

a\l~，1必2I2=-j-=

Zll

_"13

M

%3=A113=

"13=一。13(2-47)

Ml%

aalj

\j=111lj-=="=喈(/=2,3,...川

h%

可见,L中第一列出J对角线元素与A中第一列出J对角线元素相似,U中第

一行非对角线元素等于A中第一行非对角线元素用对角线元素规格化后来日勺

值。比较第一列元素，得:

(2-48)

(，=2,3,...4)

可见,L中第一列非对角线元素与A中第一列非对角线元素相似。比较第二行

()，得:

22=，2必2十，22122=a22_'2必2=Cl22~Cl2\~L~=^22

%

_a23_/2必3___21噌_感_(2)/QJQ\

'23=，2I〃13+，22〃23"G-嫂F”(2-49)

'.=/w.+/w=2J2…=2j…=磊=喈(J=2,3,…M

2y21ly222y

l224242

可见,L中第二列H勺对角线元素等于A中第二列对角线元素经第一次消去后的

值。U中第二行非对角线元素等于A中第二行上三角非对角线元素经第二次

消去后的值。比较第一列()

。32=，31%2+,32，32=。32一，3l"l2=。32—a3\an~^32

々42=^41^12+,42,42=〃42—，4|"12=〃42—〃41〃12=。42

(2-50)

《2=4防2+li2li2=ai2-4M2=ai2-孙娟)=娟)(l=3,4,…J2)

可见,L中第二列非对角线元素等于A中第二列下三角非对角线元素经第一

次消去后时值。

以上分析得到的结论是：三角分解的L矩阵中时下三角元素与式(2・44)

因子表的下三角部分(包括对角线元素)相等；U矩阵中的上三角元素与式

(2-34)因子表的上三角部分(不包括对角线元素)相等。L矩阵、U矩阵称

为因子表矩阵,可由高斯消去法得到。

另一种三角分解KJ措施是将A分解为单位下三角矩阵L、对角线矩阵D和

单位上三角矩阵日勺乘积

A=LDU(2-51)

展开为

式中，如取U与式(2・45)中时U相似；D为角线矩阵,

由式(2-45)中L矩阵的对角线元素构成。则式(2-51)中L口勺非对角元素为

式(2-45)中L中的非对角元除以所在列的对角线元素而得。当A为对称矩阵

时,有

A=Ar=(LDU)T=UTDLT=LDU

L=UT

U=LT

2、即，当A为对称矩阵时，式(2-51)中日勺L和U互为转置。

3、用三角分解法解线性方程组

设线性方程组:

Ax=b(2-53)

系数矩阵分解为，引入中间矢量和，并令:

Ux=y(2-54)

Dy=z(2-55)

则方程组(2-53)变为

Lz=b(2-56)

这样,解线性方程组(2-53)可由如下几种环节完毕：

1)由式(2.56)求出，称为前代运算。

2)由式(2-55)求出，称为除法运算。

3)由式(2-54)求出，称为回代运算。

(1)前代运算

将L分解为一种单位矩阵与一种严格下三角矩阵的和,

式(2-56)变为:

(I+L)z=b=>z+Lz=b=>z=b-Lz(2-57)

将式(2-57)展开，得

(2-58)

按式(2-58)写出求解的计算流程如式(2.

59），称为按行前代。

z<-b

Loopi=2,n

Loopj=1,z-1

(2-59)

Zi=Z「hjZj

EndLoop

EndLoop

式(2-57)展开也可写成如下口勺形式:

(2-60)

按式（2-60）写出求解的另一计算流程为如式

（2-61），称为按列前代。

z<—b

Loopz=1,«-1

Loopj=i+\,n

(2-61)

Zj=Zj-lj,

EndLoop

EndLoop

前代运算从下标小口勺元素开始，由前向后进行。实际上就是用对常数项

的消去运算。

(2)除法运算

按式(2-55)

y,=zjd.(z=1,2,...,n)(2-62)

(3)回代运算

将U分解为一种单位矩阵与一种严格上三角矩陈的和,

式（2-54）式变为:

(1+0)x=ynx+Ux=ynx=y-Ux(2-63)

将式（2-52）两边展开

W

0W]2•••|(W-1)

••••••0・・・・・・•••・.•

=—W(w-2)(n-l)“(”-2)〃*2(2-64)

)'〃一1°"s-n”

%.0J.

按式（2-64）写出求解日勺计算流程如式（2-65）,

称为按行回代。

x<-y

Loopi=n-1,1

Loopj=i+\,n