电力网络方程求解核心技术

第二章电力网络方程求解技术

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

为线性节点电压方程；稳态分析普通关注功率和电压，其电力网络模型常为非

线性潮流方程，而非线性潮流方程也必要通过求解线性修正方程才干得到其

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

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

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

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

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

解。但是，这只是抱负化假定，在计算过程中，完全杜绝舍入误差是不也许，

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

解也不一定是绝对精准。

迭代法就是用某种极限过程去逐渐逼近线性方程组精准解办法。该办法具

备对计算机存贮单元需求少，程序设计简朴、原始系数矩阵在计算过程中不变

等长处，是求解大型稀疏系数矩阵方程组重要办法。迭代法不是用有限步运算

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

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

机硬件环境和问题规模是密切有关。普通说来，对同等规模线性方程组，直接

法对计算机规定高于迭代法。对于中档规模线性方程组(n200)，由于直接法

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

迭代法可避免直接法带来高舍入误差。

计算机在电力系统应用初期，曾经由于内存容量限制采用过迭代法求解电

力网络线性方程式组。迭代法致命缺陷是存在收敛性问题。自从稀疏技术成功

地在电力系统应用之后，迭代法儿乎完全被所代替。但随着电力系统规模迅速

扩大，使得直接法很难满足在线应用规定，规定采用并行计算技术提高电力系

统分析计算速度。由于迭代法有较好并行性，也许会再次得到广泛应用。

由于电力网络构造特点，在以导纳矩阵表达电力网络方程中系数矩阵和常

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

疏矢量有关运算中，有零元素参加运算是没有必要进行，对零元素存储也是多

余。因此，可以采用“挂零存储”、“排零运算”办法，只存储稀疏矩阵和稀

疏矢量中非零元素及必要检索信息，只取这些非零元素来进行运算，省去对零

元素存储和与零元素进行运算，这样可以大大减少存储量，提高计算速度。这

种作法用计算机程序来实现就是稀疏技术。它涉及了稀疏矩阵技术和稀疏矢量

技术两方面。和不采用稀疏技术相比，采用稀疏技术可以加快计算速度几十甚

至上百倍，并且对计算机内存规定也可以大大减少。电力系统规模越大，使用

稀疏技术带来效益就越明显。可以说，稀疏技术引入是对电力系记录算技术一

次革命，使许多本来不能做电网计算可以很容易地实现。

第一节线性方程组迭代解法

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

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

同解方程组X=MX+g，以此构造迭代关系式

Xs+i)=MXd)+g(2・1)

式中，M称为迭代矩阵。任取初始矢量X(o)=Lo)^0)...”(J,代入式

12

(2-1)中，经迭代计算得到解序列X⑴,X⑵,…。若解序列收敛，设X,极限为

X*，对迭代式两边取极限

limX(*+o=lim(MX(*)+g)

即X“=MX—g，X*是方程组AX二y解，此时称迭代法收敛，否则称迭代法发

散。迭代法长处是占有存储空间少，程序实现简朴，特别合用于大型稀疏矩

阵；不尽人意之处是要面对判断迭代与否收敛和收敛速度问题。

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

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

p（M）=maxk|<1（2-2）

是迭代收敛充分必要条件，其中入是矩阵M特性根。

/

因而，称谱半径不大于1矩阵为收敛矩阵。计算矩阵谱半径，需规定解矩

阵特性值，普通这是较为繁重工作。可以通过计算矩阵范数等办法简化判断收

敛工作，其中，计算矩除1范数和a范数办法比较简朴°（向量1范数等于向

量元素绝对值之和，向量喻数范数等于向量元素绝对值最大值。矩阵1范数

等于矩阵列向量1范数最大值；矩阵8范数等于矩阵行向量1范数最大值。）

Mil=max2Tm（2-3）

1以“日"

Ml=max^\m|（2-4）

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

迭代矩阵M满足同|<1或心寸<1，就可以判断迭代序列XkD=MX-+g是收

100

敛。但要注意是，当||M>1或卜H>1时，可以有P（M）<1，因而不能判断迭代序

18

列发散。

在计算中当相邻两次向量误差某种范数Ik-x,||不大于给定精度时，则

P

停止迭代计算，并视x5为方程组AX=y近似解。

二、雅可比（Jacobi）迭代法

1、雅可比迭代格式

设M元线性方程组

(2-9)

得雅克比迭代式矩阵形式

X(A-+1)=BX(A)+g(2-10)

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

2、雅可比迭代收敛条件

当方程组—=丫系数矩阵A具备某些性质时，可直接鉴定由它生成雅可比

迭代矩阵是收敛。

定理：若方程组。=丫系数矩阵A，满足下列条件之一，则其雅可比迭代

法是收敛。

（1）A为行对角占优阵，即幺|>女"=12…,〃

⑵A为列对角占优阵，即/|>XhI.

「/「vM=11,…,〃

定理：若方程组。=丫系数矩阵A为对称正定阵，并且2D-A也为对称正

定，则雅可比迭代收敛。（D为A对角线元素构成对角线矩阵）

3、雅可比迭代算法

loopi=tn//形成雅克比迭代矩阵和常数项

gI=y[IaH..

loop/=1,/-1

b=-aIa

yVii

endloop

loopj=i+1,〃

b=-aIa

yij"

endloop

endloop

A1={0,0,…,o}〃解矢量中间值,存放迭代前的值

工2=%,»,1}〃解矢量，存放迭代后的值，开始时存放解的初值

while||x1-x2||>£//迭'弋循环

p

x\=x2

loopZ=1,7?

x2I=gI.

loopj=l,z-l

x2=x2+)*xl

i>>Jj

endloop

loopj=/+l,n

x2=x2+"xl

••ijj

endloop

endloop

endwhile

三、高斯•塞德尔(GaussSeideD迭代法

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

在雅可比迭代中，工优用计算公式是

i

*4+i>_w+Z〃w)+g(2-11)

1~0>jjijji

j=\/=i+l

事实上，在计算X(N)前，已经得到WZ),,双旬值，不妨将已算出分量直接

i1i-l

代入迭代式中，及时使用最新计算出分量值。因而松川计算公式可改为：

i

刈+1)*&+i)+w)+g(2-12)

i~Uijjijji

y=lj=Ml

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

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

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

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

对于方程组AX=y，如果由它构造高斯・塞德尔迭代和雅可比迭代都收敛，

那么，多数状况下高斯一塞德尔迭代比雅可比迭代收敛效果要好，但是状况并

非总是如此。

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

loopi=tnII形成迭代矩阵和常数项

gI=..IH

loopJ=1,z-1

b=-aIa

VV>>

endloop

loopj=i

b=-aIa

UtJH

endloop

endloop

,6={0,0,10}//解矢量，存放迭代前的值

工2=%,h-,1}〃解矢量，存放迭代后的值

while||¥1-A-2||>£//迭代循环

p

x1=x2

loopi=1,A?

x2I=gI.

loopj=1,z-1

x2=x2+b*x2

i•Vj

endloop

loop7=z+1,n

x2=x2+0*x2

，•UJ

endloop

endloop

endwhile

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

1、逐次超松弛迭代格式

方程组AX=y雅可比迭代形式Xg=BX⑹+g，记口=[+守其中[是下三角

矩阵，。是上三角矩阵。得高斯-塞德尔迭代形式：

X伏+I)=LX(A-+D+UX⑻+g(2-13)

t己AX伏）=X（#+i）—x（外,有

X伙+i）=X（^）4-AX（A）(2-14)

这样AX⑻可视为X⑻修正量，如果将XE改为X（灯加上修正量AX⑻乘一种因子

①，迭代格改为：

X（i+1）=X（&）+wAX（Jl）

X（*+i）=Xot）+a）（X伏+1）—X伏））

X优+1）=X伏）+co（LX优+1）+UX（K）+g-X（A-））

整顿得

XQ+1）=（1-co）X（A）+（D（LX（＜+1）+UX伙）+g）（2・15）

这里①为修正因子，称为松弛因子，而式（2・15）称为松弛迭代。

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

定理：逐次超松弛迭代收敛必要条件0〈卬＜2。

定理：若妫正定矩阵，则当0＜。＜2时，逐次超松弛迭代恒收敛。

以上定理给出了逐次超松弛迭代因子范畴。对于每个给定方程组，拟定力

究竟取多少为最佳，这是比较困难问题。

普通，把。〈。〈1迭代称为亚松弛迭代，把。=1迭代称为高斯・塞德尔迭

代，而把。〈田＜2迭代称为松弛迭代。

3、逐次超松弛迭代算法

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

loopy=1,M

b=-aIa

•Jtju

endloop

loopj=i4-1,/?

b=-aIa

•jij"

endloop

endloop

N={0,0广・,0}//解矢量，存放迭代前的值

不2=%,1,・，1}〃解矢量，存放迭代后的值

8=1-co

while|卜1-x211>£//迭代循环

p

x\=x2

loopi=\ji

x2=3*元1+co*g

loopj=l,z-l

x2=x2+（o*h*x2

••>Jj

endloop

loopj=i+\,n

x2=X2*x2

•>VJ

endloop

endloop

endwhile

第二节线性方程组直接解法

线性方程组可以用消去法直接求解，虽然是很古老办法，但是计算实践表

白，对电力系统来说是很有效。这是由于电力系统中常用大型线性方程组系数

矩阵，如导纳矩阵是十分稀疏，因此当充分运用矩阵稀疏性时，直接解法计算

速度不久。与上节简介迭代法相比，虽然直接解法占用计算机内存量要大某

些，但是它没有收敛性问题。

本节对消去法进行普通数学描述，给出合用于电子数字计算机表达式，并

简介它惯用变态形式一因子表解法和三角分解解法。

一、高斯消去法

高斯顺序消去法基本思想是：对线性代数方程组所相应增广矩阵（A|B）

进行一系列“把某一行元素倍加到另一行上”初等变换，使得（A|B）中A对

角线如下元素消去为0,从而使原方程组等价转化为容易求解上三角形线性代

数方程组，再通过代得到上三角形线性代数方程组解，即可求得原方程组

解。高斯消去法求解线性方程组分为两个环节：消去运算（前代运算）和回代

运算。

1、消去（前代）运算

设有线性方程组:

(2-16)

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

A=AB

消去运算有两种办法，按列消去和按行消去。一方面讨论按列消去过程,

其环节是：

第一步，消去第一列。

一方面，把增广矩阵飞第一行规格化为

ID...(1)〃⑴(2-18)

式中:

(7=2,3,...,〃)

(2-19)

然后，用式（2.18）所示行消去瓦第一列对角线如下各元素,成果

使A第2〜n行其她元素化为

々⑴=〃—aa(i)(7=2,31..,〃；/=2,3"..,力)

9•}“'7(2-20)

〃⑴=〃-〃〃(])(7=2,3,...,//)

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

1a(1)

1F=[AB]=

与之相应方程组是AX=B，它与AX=B同解。矩阵未标出元素为零，下同。

11

第二步，消去第二列。

一方面，把增广矩阵鼠第二行规格化为

I

01々3⑵…〃⑵【）3（2-21）

2,n2

式中：

0(1)

a(2)。=3,4,…，力

<2/嗤(2-22)

1内=如

I120(.)

、22

然后，用式(2.21)所示行消去A第二列对角线如下各元素……⑴，，

13242

〃⑴，成果使五第3〜n行其她元素化为

〃21

〃⑵)=〃⑴-〃(1)4(2)(7=3,4,3,4,…

4步i/i22j(2-23)

〃⑵=。⑴一a(\)b<2)(/=3,4,...,〃)

///22

式中：卜标(2)表达该元素第二次运算成果「这时矩阵T变为

!2

Ia（1）a⑴—a(1)

12131n1

1a⑵…a⑵〃（2）

7T=|AB]=232n2

a⑵...a（2）〃⑶

222333〃3

•••・・・...…

b（«）

一〃⑵­••a⑵

“3nnfl

普通地，在消去第k列时要做如下运算:

4U-I)

a(A)=——(八上+1,…，力

穹原九一1)

核（2-24）

必

户)=-1)

a(k-\)

kk

a（A）=a（k-\）—（/）{j=k+1,...,n\i=k+1,...,/z）

ijifik电（2-25）

〃⑷=〃（i）-au-i）b<*）（，=/+1//）

/i洸/

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

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

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

「1

a(Da⑴・・・a(1)。⑴

12131

1a(2)•••a(2)。（2）

232〃2(2-26)

AB]=1a⑶〃⑶

nnn3〃3

1〃(")

一//

与之相应方程组是AX=B，即

nn

⑴X+〃⑴X+•••+〃(Dx=。⑴

11221331"H1

X+〃(2))+...+〃(2)X=〃⑵

2232/.//2(2-27)

〈X+…+〃(3汝二〃⑶

33“,”3

X=b㈤

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

以上算法一方面消去A中第一列对角线下元素，然后消去第二列对角线下

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

法。

下面简介按行消去法，它一方面消去A中第二行对角线左边元素，然后消

去第二行对角线左边元素，依次直到对角线左边所有元素都被消去为止-其环

节是：

第一步，一方面对增广矩阵T第一行做规格化运算，成果为：

1彳2⑴。3⑴…〃⑴。⑴（2-28）

1

式中：

a⑴=—LL（/=2,3,

.17/n（2-29）

"1）=.

11广

11

第二步，一方面用式（2・28）所示行消去A第二行对角线左边元素，，成

21

果使了第2行其她元素化为

4⑴（尸23...,〃）（2.30）

。⑴=b-a〃⑴

22211

这时矩阵A变为了:

1

1a(i)...a(i)b⑴

12In|

I

^7=[ABJ=a(i)...Q⑴b⑴

1；222n2

aa...ab

nln2nnn

然后，对增广矩阵5■第二行做规格化运算，变为：

I

U

01°23a（?）b⑵(2-31)

2m2

式中：

(2-32)

这时矩阵T变为

1

WB/=：

a

nl

显然，与之相应方程组是AX=B，它与AX=B同解。

22

第三步，一方面用式（2-28）所示行消去A第三行对角线左边第一元素

2

Q；然后用式（2-31）所示行消去A第三行对角线左边第二元素Q⑴；最后用

31232

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

普通地，在消去第k行时要做如下运算:

[a(m)=a(m-l)-Q(m-l)a(m)(m=12…，A-1;/=k,3,…，〃)

Jkj为kmmj(2-33)

(m)=}("T)-5)(m=,/c-1)

kkkmm

Q(hl)

a(k>=3_(j=k+l,…

IgQ("l)

«kk(2-34)

[bik)=_k__

I1kQ(k-\)

ikk

通过n次消去运算，得到与式(2-26)相似增广矩阵，和与式(2-26)相似

同解方程。

2、回代运算

将方程组AX=B展开为:

1爪1)47(1)...如)「印

「t一111।

1213141〃XJ

1a(2)0(2)...a(2)42)I

21

23242〃I2

43)

1〃⑶...43)I£

343〃=3(2-35)

I

L

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

(2-35)求出原方程解过程。可以采用按行回代算法或按列回代算法。

按行回代以行自下而上顺序进行。其过程为：一方面由第〃个方程得到解

x=1M(2-36)

然后，将X代入第〃一1个方程，解出X：

nv-1

x=伙，T—x(2-37)

//—1(〃-1)〃n

再将将X，X代入第〃-2个方程可解出X：

n/A—1n—2

X=4〃-2)—4(〃-2)X-x(2-38)

,L2n-2(L2)(L1)("-2)"〃

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

对i+2”i

X=帅一£a(i)X(z=〃，/7—1,...2,1)(2-39)

•iyJ

尸外1

式(2・39)就是按行回代普通公式。

按列回代计算公式是:

I

"1)〃⑴

12

11“11,"

4"-2)I

X=0八Q-■•■-X—加7)I.X(2-40)

2ZT-2,W-1I

n-2少2n-2,n|n

X4呷)00"(I)

।Iw-1,»

IX施)_0_0_0

//■J〃■\

也是按“，X，，X顺序依次求各位置数。

n1

二、因子表法

1、因子表

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

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

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

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

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

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

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

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

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

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

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

由式（2・24）和（2・25）可知，消去（前代）过程对常数项第i个元素b运算涉

/

及

b（k）=b（*-D-a（*-i）b（*）（k=1,2,,/-1）（241）

7iikk

（2-42）

it

可见，常数项消去运算只与系数矩阵下三角某些和常数项关于。将式（2-

31）和（2.42）中运算因子按它们下标指定位置放在下三角某些，和消去运算

得到上三角矩阵放在一起，就得到了因子表

a。⑴...a⑴1

I11121314In

4

。⑴々⑵•••a⑵

2223242〃

⑶

a。⑴a⑵a⑶•••a

313233343"(2-43)

a。⑴a⑵67(3)•••6/(4)

414243444”

aa⑴4(2)a⑶・•・a5-1)

L

n1”2”3"4nn」

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

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

些用于对常数项消去（前代）运算，上三角元素用来进行回代元算。因子表也

常写成如下形式:

dUUU...U

.11.1213141M

IdUU•••的“

721产7324u(2-44)

Cpi产,33产...3〃

/Ia

•424344.4〃

***・•

•*•*

////...d

n1n2n3nn

式中

d—a(i)(z=1,2,...,n)

iiU

li=4⑺(/</,/=1/=Z+1,...,72)

iJiJ

I=4((J<i,i=1,2,...,n,j=1,2,...,/-1)

iJU

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

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

loopk=1,77〃人为消去的计数

loopy=Z:+1,7?

a=Qla〃规格化第Z行第/•列元算

勺kjkk

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

a=a-aa

ijij・kj

endloop

endloop

endloop

（2）对常数项做消去运算

\Mopi=\,〃//z.正在消去常数项的序号

Loopk=1,/-1〃用下三角元素进行消去运算

h=b—ab

iiikk.

UndL^)op

b=b/a〃规格化运算

//li

End\j)op

(3)回代运算得方程解

Loop/=77,1//从后往前回代

〃赋解的初值

lj)op//用已求得的解回代

x=x—ax

ii引j

EndIJ)op

EndLop

三、三角分解法

1、矩阵三角分解

设方阵A可用一种下三角矩阵L与一种单位上三角矩阵U乘积表达，即:

A=LU(2-45)

展开为:

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

较第一行元素，得:

〃对角线元素

itiiaa

a=IM11=_12_=_LL=^1)

12111212/a12

a-H〃===cK\>(2-47)

—W--H-

13II13\31a13

iiii

aa

a=/〃〃=T*=-A=〃（i）（7=2,3,...,〃）

Viiiy•y/au

IIII

可见，L中第一列对角线元素/与A中第一列对角线元素彳］相似，U中第一

11

行非对角线元素等于A中第一行非对角线元素用对角线元素规格化后来值。比

较第一列元素，得:

'=a

2121

'=a

3131(2-48)

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

(z=2,y>2)，得:

a

+/'=4-///=〃a—H-=.7(1)

22222221122221a22

II

a一/M-a

=/〃+/////=232|v-----H—H-(2-49)

232113222323/)1)

2222

a-///a—a〃⑴〃⑴/.no\

a=/〃+/〃力露I—-------ut,.—=4⑵1J=J,…，")

2j21\J222J。⑴2/

222222

可见，L中第二列对角线元素/等于A中第二列对角线元素经第一次消去后值

22

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

消去后值。比较第二列（?>2,尸2）

a=/z/+/a一/U=a-a〃⑴：=a（\）

323112323232311232311232

a=/u+/%2—/Ua(

42411242424112=%2—〃4必，二=li