第七章-电力系统小干扰稳定分析

电力系统在运行过程中无时不遭受到一些小的干扰，例如负荷的随机变化及随后的发电机组调节；因风吹引起架空线路线间距离变化从而导致线路等值电抗的变化，等等.这些现象随时都在发生。和第6章所述的大干扰不同，小干扰的发生一般不会引起系统结构的变化。电力系统小干扰稳定分析研究遭受小干扰后电力系统的稳定性。

系统在小干扰作用下所产生的振荡如果能够被抑制，以至于在相当长的时间以后，系统状态的偏移足够小，则系统是稳定的。相反，如果振荡的幅值不断增大或无限地维持下去,则系统是不稳定的。遭受小干扰后的系统是否稳定与很多因素有关，主要包括:初始运行状态，输电系统中各元件联系的紧密程度，以及各种控制装置的特性等等。由于电力系统运行过程中难以避免小干扰的存在，一个小干扰不稳定的系统在实际中难以正常运行.换言之，正常运行的电力系统首先应该是小干扰稳定的。因此，进行电力系统的小干扰稳定分析,判断系统在指定运行方式下是否稳定,也是电力系统分析中最基本和最重要的任务。 虽然我们可以用第6章介绍的方法分析系统在遭受小干扰后的动态响应，进而判断系统的稳定性，然而利用这种方法进行电力系统的小干扰稳定分析，除了计算速度慢之外，最大的缺点是当得出系统不稳定的结论后，不能对系统不稳定的现象和原因进行深入的分析.李雅普诺夫线性化方法为分析遭受小干扰后系统的稳定性提供了更为有力的工具。借助于线性系统特征分析的丰富成果，李雅普诺夫线性化方法在电力系统小干扰稳定分析中获得了广泛的应用。 下面我们首先介绍电力系统小干扰稳定分析的数学基础。

李雅普诺夫线性化方法与非线性系统的局部稳定性有关。从直观上来理解，非线性系统在小范围内运动时应当与它的线性化近似具有相似的特性。

将式(6—290)所描述的非线性系统在原点泰勒展开,得

fxexfx式中：A如果hx在邻域内是x的高阶无穷小

xxx0xxe量，则往往可以用线性系统

的稳定性来研究式(6-288)所描述的非线性系统在点xe的稳定性[1]：

(1）如果线性化后的系统渐近稳定，即当A的所有特征值的实部均为负，那么实际的非线性系统在平衡点是渐近稳定的。

(2)如果线性化后的系统不稳定，即当A的所有特征值中至少有一个实部为正，那么实际的非线性系统在平衡点是不稳定的。

（3）如果线性化后的系统临界稳定,即当A的所有特征值中无实部为正的特征值,但至少有一个实部为零的特征值,那么不能从线性近似中得出关于实际非线性系统稳定性的任何结论。

显然，李雅普诺夫线性化方法的基本思想是，从非线性系统的线性逼近稳定性质得出非线性系统在一个平衡点附近的局部稳定性的结论。

在进行电力系统的小干扰稳定分析时，我们总是假设正常运行的系统(运行在平衡点xxe或x0)在tt0时刻遭受瞬时干扰，系统的状态在该时刻由0点转移至xt0。这个xt0就是干扰消失后系统自由运动的初始状态。由于干扰足够小,xt0处x0的一个足够小的邻城内，从而使得hx在x0的邻域内是x的高阶无穷小量。因此，根据李雅普诺夫线性化理论,可以用线性化系统的稳定性来研究实际非线性电力系统的稳定性。为此，将描述电力系统动态特性的微分—代数方程式（6-1）、式（6-2）在稳态运行点x0,y0线性化,得



式中:

记R表示实数集合，Rn表示n维实向量空间，Rmn为所有m行n列实数矩阵组成的向量空间。定义Rn等于Rn1，即Rn中的元素是列向量；另一方面，R1n 中的元素是行向量。显然，上式中ARnm,BRnm,CRnm,DRnm,。

在式(7—3)中消去运行向量y，得到

式中：

矩阵ARnn,通常被称为状态矩阵或系数矩阵。

由此可见，小干扰稳定性分析实际上是研究电力系统的局部特性，即干扰前平衡点的渐近稳定性。显然，应用李雅普诺夫线性化方法研究电力系统小干扰稳定性的理论基础是干扰应足够微小。因此我们说这样的干扰为小干扰，当此干扰作用于系统后,暂态过程中系统的状态变量只有很小的变化,线性化系统的渐近稳定性能够保证实际非线性系统的某种渐近稳定性。

至此，我们知道，稳态运行情况下电力系统遭受到足够小的干扰后，可能出现两种不同的结局：一种结局是，随着时间的推移干扰逐渐趋近于零（即有扰运动趋近于无扰运动，对应于矩阵A的所有持征值都具有负实部）,我们称系统在

此稳态运行情况下是渐进稳定的，显然受扰后的系统最终将回到受扰的的稳态运行情况；另一种结局是，无论初始干扰如何小，干扰x都将随着时间的推移无限增大(对应于矩阵A至少有’一个实部为正的特征值），显然系统在此稳态运行情况下是不稳定的。对于实际运行的电力系统来说，分析临界情况下的系统稳定性并无多大意义，可以视它为系统小干扰稳定极限的情况。

最后需要说明的是，前面在研究系统的稳定性时，假设干扰是瞬时性的,即系统的状态在瞬时由x0转移至此xt0，并且引起变化的干扰消失。这同样适用于研究永久性干扰下系统的稳定性，即此时我们可以把它考虑成研究系统在新的平衡点遭受瞬时性干扰的稳定性。

另外，对一些给定的小干扰不稳定或阻尼不足的运行方式，可以通过特征分析方法得到一些控制参数和反映系统稳定性的特征值之间的关系，进而得出提高系统小干扰稳定性的最佳方案.因而进行电力系统的小干扰稳定分析显得尤为重要。

这样，电力系统在某种稳态运行情况下受到小的干扰后，系统的稳定性分析可归结为

（1)计算给定稳态运行情况下各变量的稳态值。

（2）将描述系统动态行为的非线性微分—代数方程在稳态值附近线性化，得到线性微分-代数方程。

（3）求出线性微分—代数方程的状态矩阵A，根据其特征值的性质判别系统的稳定性。

以上讨论的小干扰稳定问题主要涉及发电机组之间的机电振荡，这时我们将发电机组看成是集中的刚体质量块。然而,实际的大型汽轮发电机组的转子具有很复杂的机械结构，它是由几个主要的质量块,如各个汽缸的转子、发电机转子、励磁机转子等,通过有限刚性的轴系联接而成.当发电机受到干扰后，考虑到各质量块之间的弹性，它们在暂态过程中的转速将各不相同,从而导致各质量块之间发生扭(转)振（荡）(Torsional Oscillation)。由于各质量块的转动惯量小于发电机组总的转动惯量,因此各质量块之间扭振的频率要高于发电机组之间机电振荡的频率，这个频率一般在十几到四十几赫兹之间,因此也常将这种振荡称为次同步振荡(Subsynchronous Oscillation,SSo)。

次同步振荡发生后，在发电机组轴系中各质量块之间将产生扭力矩．轴系反复承受扭力矩会造成疲劳积累，从而降低轴系的使用寿命;当扭力矩超过一定限度后会造成大轴出现裂纹甚至断裂。系统出现的次同步振荡主要与励磁控制、调速器、HVDC控制及串联电容器补偿的输电线路的相互作用有关。进行电力系统的次同步振荡分析时,首先应建立汽轮发电机组的轴系模型；另外，由于扭振的频率较高，故系统中各元件不能再采用准稳态模型，而应计及系统的电磁暂态过程。对次同步振荡的详细分析已超出了本书的既定范围，有关电力系统次同步振荡分析的模型及方法,有兴趣的读者可参阅文献［5，6］.

本章首先推导出电力系统各动态元件的线性化方程，并给出了全系统线性化方程的形成方法和小干扰稳定计算的基本步骤，接着讨论了小干扰稳定分析中的特征值问题和电力系统振荡分析方法，最后介绍了大规模电力系统小干扰稳定分析的几种持殊方法.

7。2 电力系统动态元件的线性化方程

在进行电力系统小干扰稳定分析时，需要将各动态元件的方程线性化，下面我们推导各动态元件的线性化方程。在进行线性化时，通常不考虑所有控制装置中环节的作用。其原因是，在正常的稳态运行情况下，控制装置中状态变量的稳态值一般在其环节的之内.当干扰足够小时,各状态变量的变化也足够小，使得其变化范围不会超出其环节的。至于一些控制装置中的失灵区，一般认为失灵区很小，可以忽赂不计；而当失灵区很大时，可以认为整个控制系统不起作用。

7。2.1 同步发电机组的线性化方程 1. 同步电机.

对式(6—114）一(6-116)描述的同步电机方程，在给定的稳态运行情况下，系

统

各

变

量

的

稳

态

值

0,Eq0,Ed0,Ed0,Id0,Iq0,Vd0,Vq0,P0,0,Eq,P,Efq0可按式(6-74)一(6m0e0—78)和式(6-118）一(6-122)算出。将各方程在稳态值附近线性化,可得到同步电机的线性化方程

（2)励磁系统.

以图5—16所示的采用可控硅调节器的直流励磁机励磁系统为例,根据式(6—136）一（6—140），可以推导出其线性化方程。

对测量滤波环节，由于VCVjXCI。根据坐标变换式（5-63），发电机端电压和电

流用它们的d,q分量可表示为

这时显然有

将上式在稳态值附近线性化可得到

式中：

对式（6-136）线性化，并格式(7-10）代入其中，从而消去VC，即得到测量滤波环节的线性化方程

用式(6-140）模拟励磁机的饱和特性,将式(6—139)在稳态运行点线性化，可得到励磁机的线性化方程

最后,将式（6—137）、式（6—138)的线性化方程和式（7—12）、式（7—3）一起，并经整理后得到整个直流励磁机励磁系统的线性化方程

（3)PSS。

对于图5—l 4所示的电力系统稳定器，根据式（6—142)、式(6—143)，当输入为转速偏差，即VISs时,可依次列出如下线性化方程：

上式经适当整理后，可得到PSS线性化方程的状态表达式

(4）原动机及调速系统。

对如图5—24所示的水轮机及其调速系统,可以根据式（6-171)一(6-177）得到其线性化方程

2．同步发电机组线性化方程的矩阵描述及坐标变换 1)发电机组方程的矩阵描述。

当发电机组采用式（7—6)、式（7—7)、式(7-9)、式（7-15)、式（7-17）描述时，将其中的状态变量按如下顺序组成向量：

并定义

这时各发电机微分方程式的线性化方程写成如下矩阵形式:

而定子电压方程式的线性化方程表示为

以上两式中系数矩阵Ag,BIg,BVg,Pg,Zg的元素可以很容易地通过比较式(7—20）和式(7-6)、式（7—9）、式(7-15）、式（7-17）及比较式(7-12)和式（7—7)而得到，即

在同步电机、励磁系统、原动机及其调速系统等采用其他模型时，同上原理，总可以先写出各自的线性化方程,然后表示成式（7-20)、式(7-21）的形式.另外还需注意，式（7-18）中各状态变量的排序并不是一成不变的，不同的排序下有相应的矩阵。 （2)坐标变换。

式(7—20)和式（7-21)中的Vdqg和Idqg为各发电机本身d,q轴电压和电流分量的偏差，因此必须把它们转换成统一的同步旋转坐标参考轴xy下的相应分量，以便将它们和电力网络联系起来。

对于发电机端电压,由坐标变换式（5-62）可知

稳态值Vd0,Vq0,Vx0,Vy0和0也应满足式(7-22)，即

将式(7—22）在稳态值附近线性化，得

利用式（7-23)，式（7-24）可另写为

简写成

式中:

很明显，Tg0为正交矩阵，即满足

同理,对发电机电流也可得到以下关系：

式中

将式（7-26）和式（7—28）代入式（7—21）消去Vdqg和Idqg，可以得到

式中:

将式(7-26）和式（7—28）代入式(7-20）消去Vdqg和Idqg，并利用式（7-29)、式(7—30）消去Ig，可以得到

式中：

式（7-31)和式（7—29）便组成每个发电机组的线性化方程，它类似于一般线性定常系统的状态方程和输出方程. 7.2.2 负荷的线性化方程

在小干扰稳定性分析中，负荷大都采用电压静态特性模型.如果要考虑一些感应电动机负荷，可以用类似于推导同步电机线性化方程的方法得到感应电动机的线性化方程。

无论采用什么形式模拟负荷的电压静特性,负荷节点注入电流与节点电压的偏差关系总可以写成如下形式：

式中:

其中的系数可由负荷节点注入电流与节点电压的关系式求得，即

当采用二次多项式模拟负荷的电压静特性时，可以利用如式(6—48)所示的负荷节点注入电流与节点电压的关系和式（7—35）直接求出式(7-34)中的有关系数：

当采用指数形式模拟负荷的电压静特性时，可以利用如式(6-49）所示的负荷节点注入电流与节点电压的关系和式（7-35)直接求出式(7—34）中的有关系数：

特别地，当对负荷的电压静特性缺少足够的信息时,通常可以接受的负荷模型是：负荷的有功功率用恒定电流(即取m1）、无功功率用恒定阻抗（即取m2）模拟。

7。2。3 FACTS元件的线性化方程

(1）SVC。

由于V2Vx2Vy2，将它线性化,得

将上式代入式（7-38），经整理后得

式中：

另外，根据式（6-50）可直接得到SVC注入电流和节点电压间的偏差关系

式中：

这样式（7-40）、式(7-42）便组成了SVC的全部线性化方程式。 （2)TCSC。

从式（6-208）、式（6-209)可以直接得到如下线性化方程:

根据式（6—211）可以得到

将上式代入式（7—44），并经整理后得

式中：

另外，根据式（6-51）可直接得到TCSC注入电流和节点电压间的偏差关系

式中：

这样式(7-46）、式（7-48)便组成了TCSC的全部线性化方程式。 7.2。4 直流输电系统的线性化方程

当考虑直流线路的暂态过程时，直流线路以及整流器和逆变器的控制方程如式(6—222）、式(6—224)一（6—227)所示，利用式（6—53)中的第一式消去式(6—226)中的VdI,在忽略对,的情况下，可得到它们在稳态值附近的线性

化方程

整流器和逆变器交流母线电压的幅值与其x,y分量间的关系为

将上式在稳态值附近线性化，得

将式（7-51)代入式（7-50)消去VR和VI，并经整理后可得

式中:

式中的系数矩阵Ad,Bd通过对照式（7-52)和原方程容易得到。

两端直流输电系统的代数方程可以由换流器交直流两侧的功率关系及电流关系推得。对于整流器,将有功功率关系式

在稳定值附近线性化,得

另外，将式（6-52）中第三式两端平方,得

上式的线性化方程为

将式(7—51)代入式（7—55)中消去VR,并注意到整流器注入交流系统的无功功率QR0VyR0IxR0VxR0IyR0总不为零,于是可以从式（7-55)、式（7-57）中解出节点注入电流的偏差，并写成如下矩阵形式：

式中：

对于逆变器，将有功功率关系式

在稳态值附近线性化,得

同样，将式（6—53）中第三式的两端平方，得到的线性化方程为

同理，将式（7—51）代入式(7-61）中消去VI，式(7-61）、式（7-62）表示的电流、电压偏差关系式可以写成如下矩阵形式:

式中：

式（7-58）、式(7-63）组成了直流系统的代数方程

式中:

当直流系统采用其他数学模型时，用同样的方法可导出形如式（7—25）、式（7-65）所示的线性化方程。

7.3 小干扰稳定分析的步骤

7。3.1 网络方程

为了叙述方便，将网络方程式(6-36)写成分块矩阵形式,并注意到网络方程本身是线性的，因而可以直接写出在xy坐标下节点注入电流偏差与节点电压偏差之间的线性化方程

式中：

对各负荷节点，把式(7-33）给出的注入电流偏差与节点电压偏差关系代入上式，即可消去负荷节点的电流偏差。设负荷接在节点i，则消去该负荷后的网络方程仅是对原网络方程（7-67)的简单修正：节点i的电流偏差变为零，导纳矩阵中的第i个对角块变为YiiYli，而其他内容不变。

不失一般性，假定网络中节点编号的次序为：先是各发电机所在节点，然后是各SVC

所在节点,接下来是各TCSC的两端节点，再是各直流输电系统交流母线节点（先编整流侧节点，后编逆变侧节点），最后是其他节点.消去所有负荷节点的电流偏差后，网络方程可写成如下分块矩阵形式:

式中：IG和VG分别为由全部发电机节点注入电流和节点电压偏差组成的向量;IS和

VS分别为由全部SVC节点注入电流和节点电压偏差组成的向量;IT和VT分别为由全部TCSC节点注入电流和节点电压偏差组成的向量；ID和VD分别为由全部换流器交流母线节点注入电流和节点电压偏差组成的向量;VL为其他节点电压偏差组成的向量。这些向量可表示为

7.3。2 全系统线性化微分方程的形成

由各发电机组的方程式（7-31）、式(7—29)可以组成全部发电机组的方程式

式中:

由各SVC的方程式（7—40）、（7-42）可以组成全部SVC的方程式

式中：

由各TCSC的方程式（7-46）、式(7-48）可以组成全部TCSC的方程式

式中:

由各两端直流输电系统的方程式(7-52）、式(7-65）可以组成全部两端直流输电系统的方程式

式中：

将式（7-72）、式（7-75）、式（7—78）、式(7—81）代入式(7—69)消去IG、IS、IT、ID，所得结果与式（7-71）、式(7—74)、式(7—77)、式(7—80）一起组成如式(7—3）所示的矩阵关系式,其中：

显然，A,B,C分别为分块稀疏矩阵，而及和导纳矩阵具有同样的稀疏结构. 用式(7-83)中的矩阵A,B,C,D，根据式(7—5）即可得到状态矩阵A。至此已经得到电力系统在稳态运行点的线性化方程式。 最后，有必要说明以下几个问题：

（1)如果线性化后的系统渐近稳定，即如果A的所有特征值的实部均为负，那么实际的非线性系统在平衡点是渐近稳定的.

(2）形成矩阵A的方法，在已有的各种商业化程序中可能各不相同,以上仅给出其中的一种形成方法，皆在介绍形成A阵的原理和技巧[4—7,9,10]。式（7-83)中矩阵A,B,C,D的形式多种多样,它与状态变量的次序安排、网络方程的形式、各动态元件的代数方程和网络方程的协调处理方法等有关.不同的方法将影响到程序实现的复杂性和灵活性，但并不影响其特征值的计算结果。

(3）以上方程的形成中考虑了发电机组、SVC、TCSC、两端直流输电系统，

对电力系统中的其他动态元件可作类似处理。例如，对并联动态元件(如感应电动机负荷等)，可以仿照以上对发电机的处理方法得到其线性化方程；对多端直流输电系统，可以仿照以上对两端直流输电系统的处理方法得到其线性化方程。然后按规定的顺序将它们安排在整个系统的方程中。

（4）按照以上方法形成的系数矩阵A将必定有一个零特征值。其存在的理由是各发电机转子的绝对角度不是惟一的，换言之，系统中存在一个冗余的转子角度.事实上，由于各发电机间的功率分配取决于各发电机转子角度的相对值，如果各发电机转子的绝对角度都加上一个固定的值，并不改变各发电机间的功率分配,因而不影响系统的稳定性.若要摒除零特征值,只需选定任意一台发电机的转子角度作为参考，用其余机与该机转子的相对角度作为新的状态变量即可，这时矩阵A和相应的状态变量都将降低一阶。

(5)另外还需注意，当系统中所有发电机的转矩都与转速的变化有关,即摇摆方程的右端无阻尼项且不考虑调速器的作用时,矩阵A还将存在一个零特征值。同样，要摒除这个零特征值，只需选定任意一台发电机的转速作为参考，用其余机与该机转速的相对值作为新的状态变量即可,这时矩阵A和相应的状态变量也都降低一阶。

在明白了零特征值的来历后，可以在后面的计算步骤（5）和（6）中不作任何处理,仅需在计算结果中去除零特征值即可。然而应当注意，由于潮流和特征计算的误差，理论上的零特征值在实际上是很小的特征值. 7.3.3 小干扰稳定分析程序的组成

按照前面介绍的内容和方法，可以构成含有FACTS（例如SVC、TCSC）的交直流系统的小干扰稳定分析程序。其基本计算过程如下:

（1)对给定的系统稳定运行情况进行潮流计算，求出系统各节点电压、电流和功率。

（2）形成式（7—67)中的导纳矩阵。

（3）已知各负荷的功率及负荷节点电压的稳态值为P(0),Q(0),Vx(0),Vy(0)。根据负荷电压静特性参数,应用式(7—36）或式（7-37）求出式(7-34）中的矩阵元素

Gxx,Bxy,Byx,Gyy用它们修改导纳矩阵中对应于各负荷节点的对角子块.

(4）首先由式（6—74)—（6-78）和式(6—118）—（6—122)计算出各发电机组中所有变量的初值,然后分别形成式（7—20）和（7-21）中的矩阵

Ag,BIg,BVg,Pg,Zg及式（7-28）中的矩阵Tg0,RVg,RIg。然后应用式(7-30)和式（7—32）求出Cg,Dg,Ag,Bg,从而得到各发电机组的线性化方程。对其他动态元件，同同样的方法可得到其线性化方程中的系数矩阵.其至得出系统中所有动态元件的线性化方程。

（5）按照式（7—71）-(7-83)形成矩阵A,B,C,D，再应用式（7-5)计算出系统的状态矩阵A。

（6）应用QR法计算矩阵A的全部特征值[2—4],从而判断系统在所给定的稳态运行情况下的小干扰稳定性。计算矩阵A全部特征值的QR法将在下节介绍。

【例7-1】 9节点电力系统的单线图、支路数据、发电机参数、正常运行情况下的系统潮流分别如图6—12、表6-5、表6-6、表6-7所示。系统频率为60 Hz。

各负荷均用恒定阻抗模拟。发电机1采用经典模型,发电机2和3采用双轴模型。发电机2和3均装有自并励静止励磁系统，其参数如下：

XC0,KA200,TR0.03s,TA0.02s,TB10.0s,TC1.0s 另外，各发电机的阻尼系数Di均取为1。0。

下面研究在正常运行情况下系统的小干扰稳定性。为简单起见，在以下的矩阵中“空白”表示数0或适当维数的0矩阵。

【解】 1)利用潮流结果，根据式(6—74)一（6-78）和式(6—118)一（6-122）计算出各发电机变量的初值，如表7-1所示。负荷的等值导纳见例6—1，直接并入电力网络。

2)根据7.2。1节的方法得到各发电机组的线性化方程.

发电机1：不难算出式(7—20）、式（7—21）、式(7-26）、式(7—28)中的系数矩

阵为

最后,根据式(7-32)、式（7-30）和以上矩阵计算出发电机组线性化方程（—31）、（7—29）中的矩阵：

发电机2同上原理得到发电机2线性化方程的系数矩阵：

发电机3：同上原理得到发电机3线性化方程的系数矩阵:

3)系统的线性化方程。

显然，式(7—3)中的矩阵[见式（7—83)]为

根据式(7-5）可得到状态矩阵

4）状态矩阵A的特征值和相应的特征向量. 应用QR法求得A的全部特征值为

显然，除了我们已知的零特征值外，系统的其他所有特征值都具有负实部，因此系统在给定的运行方式下是小干扰稳定的。

7.4 小干扰稳定分析的特征值问题

既然遭受小干扰后非线性系统的稳定性可由其线性化系统的稳定性决定，而线性系统的稳定性又由状态矩阵A的特征值决定，因此下面我们简单介绍状态矩阵A的特征分析方法［2，3，6］，从而为进行电力系统的小干扰稳定性分析打下基础。

由上节可见,状态矩阵A是一个实不对称矩阵，因此后面的讨论一般仅限于

ARnn。另外，在下面的论述中要涉及到复数及复矩阵的计算，记C表示复数

集合，Cn表示n维复向量空间(列向量)，Cmn为所有m行n列复数矩阵组成的向量空间。复矩阵的标乘、相加、相乘是与实矩阵完全相对应的.但是,转置在复情形下是转置共扼(用上标H表示），即CAcijaji。n维复向量x和y的点积是sxy指满足于xHHxyii1i。另外，在p范数意义下的单位向量(或规范化向量）是

p1的向量x.例如在1、2和无穷范数意义下的单位向量x分别为

将任一向量变为单位向量的过程称为向量的归一化。 7.4。1 状态矩阵的特征特性

1. 特征值

对于标量参数C和向量vCn,如果方程

有非退化解（即v0)，则称为矩阵A的特征值。

要计算特征值,方程（7—85）可写成如下形式:

它具有非退化解的充分必要条件是

展开上式左端的行列式，得到显式的多项式方程

该方程称为矩阵A的特征方程，方程左端的多项式称为特征多项式。因为n的系数不为零,所以此方程共有n个根,这里根的集合称为谱，记为(A)。如果

(A){1,,n},则有

而且，如果我们定义A的追迹为

可以证明tr(A)12n.

实不对称矩阵A的特征值既可能是实数，也可能是复数，并且复特征值总是以共扼对的形式出现。另外,相似矩阵的特征值相同,转置矩阵的特征值不变。

2。 特征向量

对任一特征值i，满足方程

的非零向量viCn称为矩阵A关于特征值i的右特征向量。由于该方程为齐次方程,因而kvi （k为标量)也是以上方程的解，即同样是矩阵A关于特征值i的右特征向量。除非特别声明,以后提到的“特征向量”均指“右特征向量”。一个特征向量定义了一个一维子空间,这个子空间用矩阵A左乘保持不变性。

同样，满足方程

的非零向量uiCn，称为矩阵AT关于特征值i的右特征向量.方程（7—90)两边转置,得

称行向量uiT为矩阵A关于特征值i的左特征向量。

为了简明地表达矩阵A的特征持性，将A的所有特征值组成对角矩阵,相应的右特征向量按列组成矩阵XR，相应的左特征向量按行组成矩阵XL，即

以上三个n阶方阵称为模态矩阵,

利用式（7—92)，方程（7—)和(7-91）可表示成如下矩阵形式:

在上式中，前一式两边左乘XL，后一式两边右乘XR，即可得到以下关系式：

或写成

显然,相应于不同特征值的左、右特征向量是正交的；相应于同一特征值的左、右特征向量的乘积为一非零常数，通过对左、右特征向量的归一化处理总可以使这个常数为1。即有

注意,uiTvj并不是通常理解的内积。上式的矩阵形式为

根据式(7—93）和式（7—96)可得

3. 动态系统的自由运动

由状态方程(7-4）可看出，每个状态变量的变化率都是所有状态变量的线性和。由于状态之间的耦合，很难对系统的运动有一个明晰的概念.

为了消去状态变量间的耦合，引入一个新的状态向量z，它和原始状态向量x间的关系定义为

把上式代入方程（7—4）,并考虑式（7—97），状态方程即可改写为

它和原方程的差别在于：为对角阵，而A一般不是对角阵。方程（7—99）表示n个解耦的一阶方程

它的时域解为

式中：zi的初值zi(0)可根据式（7-98）用uiT和x(0)表示，即

将式（7—101)、式（7-102)代入变换式（7-98),可得到原始状态向量的时域解

其中第i个状态变量的时域解为

式中:vik表示向量k的第i个元素。上式给出了用特征值、左特征向量和右特征向量表示系统自由运动时间响应的表达式。特征值i对应于系统的第i个模态（Mode)，与之相应的时间特性为eit,这样系统自由运动的时间响应就可以说成是n个摸态的线性和。

因此，系统的稳定性可以由特征值决定:

（1)一个实特征值相应于一个非振荡模态。负实特征值表示衰减模态，其绝对值越大,则衰减越快；正实特征值表示非周期性不稳定。与实特征值有关的持征向量和z(0)都具有实数值。

(2)复特征值总是以共扼对的形式出现，即

每对复特征值相应于一个振荡模态.与复特征值有关的特征向量和z(0)都具有复数值。 因此

具有这样的形式

显然，特征值的实部刻画了系统对振荡的阻尼，而虚部则指出了振荡的频率。负实部表示衰减振荡;正实部表示增幅振荡。振荡的频率(Hz)为

定义阻尼比为

它决定了振荡幅值的衰减率和衰减特性。 7。4.2 线性系统的模态分析

1. 模态与特征向量

前面已经讨论了系统的时间响应，向量x和z的关系为

变量x1,x2,,xn式表示系统动态性能的原始状态变量。变量z1,z2,,zn 是变换后的状态变量，每一个变量仅对应于系统的一个模态。

从方程（7-107)的第一式可以看出，右特征向量呈现出模态的表现形式，即当特定的模态被激活时，各状态变量的相对活动情况。例如,右特征向量vi中的第k个元素vki给出了状态变量xk在第i个模态中的活动程度。

而各元vi中各元素的模值表征了n个状态变量在第i个模态中的活动程度，素的角度则表征了各状态变量关于该模态的相位移.

从方程(7-107)的第二式可以看出,左特征向量uiT确定呈现第i个模态时原始状态变量的组合方式。这样，右特征向量vi中的第k个元素度量变量xk在第i个模态中的活动，而左特征向量uiT中的第k个元素加权这个活动对策i个模态的贡献。

2．特征值灵敏度

我们首先考查特征值对状态矩阵A各元素akj (A的k行、j列元素)的灵敏度.将(7-)两边对akj求偏导数，得

上式两边左乘行向量uiT,并考虑式（7-91）、式（7-95），可得到

显然，A/akj中除第k行、j列的元素为1外,其余元素都为零，因此

式中：vji表示向量vi的第j个元素，uki表示向量ui的第k个元素。

设是标量，A()是由元素akj()组成的n阶方阵，如果对一切的k和j，

akj()都是可微函数，则

这样，仿照以上推导可得到特征值对标量的灵敏度

3。 参与因子(Participation Factor）

为了确定状态变量和模态之间的关系，把右特征向量和左特征向量结合起来,形成如下的参与矩阵（Participation Matrix）P，用它来度量状态变量与模态之间的关联程度：

［12]

称参与矩阵P的元素pkiukivki为参与因子,它度量了第i个模态于第k个状态

变量xk的相互参与程度。称矩阵P的第i列为pi第i个模态的参与向量。由于vki度量xk在第i个模态中的活动状况，而uki加权这个活动对模态的贡献，因此它们的乘积pki即可度量净参与程度。左、右特征向量相应元素的乘积导致pki是无量纲的，即于特征向量单位的选择。

设x(0)e，即xk(0)1且xjk(0)0,则由式(7—102）得zi(0)uki，根据式(7—103)可得

这个方程表明，被初值xk(0)1激活的第i个模态,以系数pki参与在响应xk(t)中，参与因子由此而得名。

对于所有的模态或所有的状态变量，有

上式不难得到证明。在式（7-114）中令t0，容易得到矩阵P的第k行元素之和为1；而矩阵P的第i列元素之和等于uiTvi，根据式(7—95)可知其值为1.

另外，参看式(7—110）可知，参与因子pki实际上等于特征值i对状态矩阵A的对角元素akk的灵敏度：

7.4。3 特征值的计算 1。 QR法

计算一般矩阵全部特征值的数值方法中,当首推双位移QR法，它是J。G.F。Francis在1962年提出来的.这种方法具有鲁棒性强、收敛速度快等特点，是迄今为止最有效的特征求解方法.

对于给定的ARnn和正交阵Q0Rnn，考虑如下达代：

其中，每个QkRnn都是正交阵，每个RkRnn都是上三角阵。由归纳法可得

这样，每个Ak都与A相似。由于A有复特征值，因此Ak不会收敛到严格的、“特征值暴露”的三角阵，而只能满足于计算称为实Schur分解的另一种分解。

对角均为1×1块或2×2块的分块上三角阵称为拟上三角阵（Upper Quasitriangular）。实Schur分解相当于将矩阵实归约为一个拟上三角阵。若

ARnn，则存在一个正交阵QRnn使得

其中每个Rij是1×1矩阵或2×2矩阵。若是1×1的，其元素就是A的特征值；若是2×2的，Rn的特征值是A的一对共轭的复特征值。

为了有效实现实Schur分解，选取式(7-117）中的初始正交相似变换阵Q0，使得A0为上海森伯矩阵，这样一次迭代的计算量将从O(n3)减小到O(n2)。

一个上海森伯矩阵是这样的矩阵，在其对角线下，除去第一个次对角线上元素之外的其他所有元素为零.例如,对于6×6的情形，非零元素是

至此，能够一眼看出，这样一种形式可以通过一系列豪斯霍尔德（Householder)变换得到，每一次变换将矩阵一列中的相应元素化为零.由于豪斯霍尔德变换为对称正交相似变换，所以得到的上海森伯矩阵与原矩阵具有相同的特征值。 要详细了解QR法的原理和算法实现，可参阅各种有关数值分析的教科书和专著。目前,用双位移QR法计算一般矩阵全部特征值，在一般的大、中型计算机中也都有标准的库程序供用户调用。

最后,我们应注意到，如果A中的元素数值差别很大，当实施某种迭代算法时，可能导致计算的特征值误差过大。特征值对舍入误差的敏感程度可以通过平衡来减小。由于数值过程中导致的特征系统的误差一般是与矩阵的欧几里得范数成正比，而平衡的思想就是：用相似变换将矩阵对应的行和列的范数变得相接近，从而在不改变特征值的前提下使矩阵的总范数减小.

平衡的实现是通过O(n2)运算确定对角阵D，使得若

则riciini1i2，。对角阵D选成具有形式i1,2,,nDdiag{b,b,,b}，其中b是浮点基数，这样计算A就可以没有舍入误差.当A被平衡后，计算的特征值常常会更精确。

2。 幂法（The Power Method）

一般矩阵特征值问题的库程序都能够把矩阵的全部特征值和特征向量一并求出。但在实际应用中，往往不需要计算矩阵A的全部特征值，而只需求出模数最大的特征值（通常称为主特征值，Dominant Eigenvalue）。幂法是计算矩阵主特征值和相应特征向量的一种非常有效的迭代方法。

设ACnn可对角化且X1AXdiag(1,2,,n),其中X[x1,x2,,xn]，幂法产生如下向量序12n给定2范数下的单位初始向量v(0)Cn，列v(k)：

~

显然，以上迭代中得到的向量序列v(k)都是2范数下的单位向量. 由于

并且

显然，只要2/11,当k时，

我们说幂法是线性收敛的方法，其有用性取决于比值2/1，因为它反映收敛速率。

在用幂法求得A的主特征值后，我们可以利用“收缩技术”(Deflation Technique)得到A的剩余特征值.收缩方法有多种，但数值稳定的方法并不是很多，下面我们介绍一种以相似变换为基础的收缩方法。

设1、v1已知，可以找到一个豪斯霍尔德矩阵H1,使H1v1k1e1，且k10。由A1v11v1,给出H1A1(H11H1)v11H1v1，显然H1A1H11e11e1，即H1A1H11的第一列为1e1，记

式中：B2为n1阶方阵，它显然有特征值2,,n.在23条件下，可用幂法求出B2的主特征值2和相应的特征向量y2，其中B2y22y2.设A2z22z2，为求出路z2，设为待定常数，y为n1维向量,则

b1Ty因为12，可取yy2，这样就求出z2.而v2H11z2就是A关于221的特征向量。

按以上方法和豪斯霍尔德矩阵的应用，应有

求出2、v2后，可对B2继续收缩，计算其他特征值和特征向量。理论上只

要A按模排列的前若干个特征值按模分离便可求出它们。这种收缩方法的缺点是改变了原矩阵的元素，使得当A是稀疏矩阵时，收缩过程将使矩阵的稀疏性不再保持。

最后应该指出的是，用幂法计算矩阵A的主特征值和相应特征向量的程序实现并非以上叙述的那么简单。在前面，我们仅时论了A的主特征值1是实数且是单重的情形。实际上,1还可能出现其他几种情形；1为r重实特征值;1和

2为模数相同但符号相反的实特征值;1和2为一对共轭复特征值。A的特征值

出现这些不同的情形时，幂法的实现将稍有不同，其详细讨论可参阅文献[2]。

3。 反幂法（The Inverse Power Method)

我们知道，非奇异矩阵A的逆阵A1的特征值是A的持征慎的倒数。因此，

A1的主特征值的倒数便是A的模数最小的特征值.把幂法用于矩阵A1,得到的

方法称为反幂法(或逆迭代法），它可以用来计算非奇异矩阵A的模数最小的特征值及相应的特征向量。

给定2范数下的单位初始向量v(0)Cn，则反幂法产生如下达代序列：

当k时，

反幂法的一种更有用的形式是把幂法用于矩阵(AI)1,其中为实常数或复常数.给定2范数下的单位初始向量v(0)Cn,其迭代过程如下：

当k时，

式中：p为矩阵A的所有特征值中最靠近数的特征值；而xp为相应的特征向量。对式(7—128)需要作如下解释：

由于非奇异矩阵AI的特征值为J(j1,2,,n),因此与矩阵(AI)1相对应的特征值为

1(j1,2,,n)。对矩阵(AI)1应用幂法即可得到模数J最大的特征值

1，这意味着p的模数最小,或p最靠近数。 p 这样，如果我们想要得到矩阵A的离最近的特征值和相应的特征向量，可应用式（7-127）所示的反幂法。反幂法的另外一个用途是，如果知道矩阵A的某特征值的近似值，应用反幂法可求出相应的特征向量，并能改善特征值的精度。

在实际求解式(7-127）时，可先对AI进行三角分解：

式中:L为单位下三角矩阵；U为上三角矩阵。此时求解方程成为

从而,每次迭代只需做简单的前代和回代即可. 7.4。4 稀疏特征求解方法

在小干扰稳定分析中，电力系统的动态特性由式(7-3)所示的线性微分—代数方程描述，由式（7-83)可知,其中的矩阵A、B、C、D都为稀疏矩阵。当用式(7—5）得到矩阵A,进而计算其特征值时，我们发现矩阵A几乎完全失去稀疏性.由于著名的QR法不能稀疏实现，因此在用QR法计算A的特征值时,A是否稀疏无关紧要.但是,在用其他迭代方法（例如，幂法、反幂法或后面介绍的子空间

~~~~法）计算矩阵A的部分特征值时，如果能够充分利用原始矩阵的稀疏性，直接从式（7-3)中求得A的特征值，则可提高特征值计算的效率。

对于A的特征值，满足方程

的非零向量vCn就是矩阵A关于特征值的右特征向量。上式最左侧的矩阵称为增广状态短阵.

上述结论不难得到证明。事实上，在式（7-129)中可得wDCv,然后消去w，并利用式（7-5），容易得到

~~1

这样，我们就可通过对式（7—129）所示的增广状态矩阵方程直接求解，在不破坏系统稀疏性的前提下,得到矩阵A的特征值和特征向量.下面描述幂法和反幂法的稀疏实现，并给出特征值对标量灵敏度的稀疏表达式。

1. 幂法迭代式(7-121)的稀疏实现 由于方程

所表达的z(k)与v(k1)之间的关系等价于z(k)Av(k1)，因此，可将式(7—121）中第一式的计算用下式替换：

在实施式(7-121）的迭代前，仅对D作一次稀疏三角分解，即DLU，则每步迭代中对式(7-131）的计算仅是一些稀疏矩阵与向量相乘及两个三角方程的求解.

2. 反幂法迭代式(7-127)的稀疏实现 由于方程

~~

（AI）z(k)v(k1)，所表达的z(k)与v(k1)之间的关系等价于因此，可将式(7—127)

中第一式的求解用求解式(7—132)替换，从而得到所需要的向量z(k).

对于给定的数,首先计算矩阵DDC(AI)B,并作稀疏三角分解

~~*~~~1~DLU。注意到(AI)为分块对角距阵(-个动态元件对应于一个对角块),因此

*~(AI)可通过对各对角块分别直接求逆得到.另外不难看出，D与D有同样的

~1~*~稀疏结构(2×2的分块稀疏阵)。这样，方程（7-132）的求解步骤可以归纳如下：

(1) 求解方程得到w(k)即

(2) 计算

3。 特征值对标量灵敏度的稀疏表达 类似于式（7—129），对左特征向量，有

这样，仿照前面的推导容易得到

7.4.5 特征值灵敏度分析的应用

在电力系统运行情况分析和控制器的设计中，往往需要分析系统中的某些参数（诸如控制器的放大倍数和时间常数等)对系统稳定性的影响,以便适当选择或调整这些参数，使系统由不稳定变为稳定，或进一步提高系统的稳定程度。 由于系统的状态矩阵A是系统中某一参数的函数,即A()，因此，A的任

一特征值i也是参数的函数,即i()，i1,2,,n。当改变参数时，i()将发生相应的变化，i()的变化即反映了参数的变化对系统稳定性的影响。

设参数从)变为（0）,则对应的系统特征值从i(（0）（0）变为泰勒展开，即 i(（0）)。将i(（0）)在（0）

在很小时，i的变化可近似为

式中:偏导数i/就是特征值i对参数的一阶灵敏度，简称特征值灵敏度.这样,如果能求出i/，则可以根据所希望得到的特征值变化i来近似地决定。

特征值i对参数的一阶灵敏度计算可归纳如下： （1)置，形成状态矩阵A()。 （0）（0）（2)计算A()的特征值i和相应的左、右特征向量uiT和vi，使uiTvi1。 （0）（3)计算

A(）。

（0）（4）

i(）TA(）uivi （0）A(）的计算，下面以为PSS的放大倍数KS、相位环节的时间其中对于

常数T1、T2、T3、T4为例加以说明。在发电机组g的方程式(7-20）、式（7-21）中，除了Ag与有关外，BIg、BVg、Pg、Zg都与无关。另外,RIg、RVg、Tg(0)也都与无关.因此,根据式（7—30)、式（7—32）可得



其中的矩阵

Ag不难由式(7-20）所示的矩阵Ag求出。

显然，在全系统的方程式中,由式（7-83)可求出

其中的AG/可由式（7—73)求出。另外，根据式（7-5）和上式可得

矩阵A对于其他参数的偏倒数也可以类似地求出。值得注意的是，在推导过程中必须仔细观察参数与该元件方程有关矩阵Ai、Bi、Ci、Di之间的关系，从而防止遗漏和错误。

在特征值灵敏度分析方面，除了以上介绍的特征值对参数的灵敏度外,目前还提出了特征值对运行方式的灵敏度；为了提高特征值对参数灵敏度的计算精度，还提出了特征值对参数的二阶灵敏度,并提出了一些有效的计算方法.有关这些方面的成果可参阅文献［34~371。 7.5 电力系统的振荡分析

没有控制的电力系统是不能运行的。运行调度人员可以通过对发电机功率的控制(包括一些其他控制装置的投切)来满足各种预定的负荷需求。而系统中的一些自动控制装置(例如,发电机的调理器、励磁调节器，HVDC的控制,FACTS元件的控制等）则是在系统遭受到干扰后承担着快速调节的任务．从而使得系统的频率和电压保持在设计的限值之内.

20世纪中叶，电力工业界发现将各区域电力系统互联可以获得更高的可靠性和经济性，至今50年来，电力系统规模发展得越来越大。在20世纪60年代，北美刚刚建立起来的大型互联系统便遭受到增幅振荡，从而破坏了大型系统间的并联运行.研究发现，各区域电力系统之间大多通过长距离输电线路互联，这种大型的弱耦合系统本身固有的对振荡的薄弱阻尼是产生这种现象的主要原因，而

高放大倍数的快速励磁系统则进一步加重了负阻尼的状况。工程师们认识到，通过给励磁系统引入适当的控制信号(PSS）可以使系统的阻尼得到加强。北美电力系统的经验表明，这个方案对于克服增幅振荡是非常成功的。增幅振荡妨碍了各区域电力系统的互联.一些互联的电力系统，为了避免增幅振荡的发生，不得不输送很少的功率，从而使得互联变得没有多少实际意义;一些互联的电力系统，为了避免增幅振荡的发生，不得不采用低放大倍数的励磁调节器.因此，直到异步HVDC互联方案产生前，一些系统干脆放弃互联.

在20世纪40年代，人们已经认识到励磁控制可以增加同步发电机的稳定极限。从另一个角度去看这个问题，就是允许发心机在较高的电抗下运行.励磁系统在某些情况下改善电力系统动态性能所取得的成功,加之其控制速度快、效率高,使得电力系统工程师们对它的控制能力寄予了更高的期望。但是应该明白,励磁系统的控制效果是有限的，Concordia早就警告说：“我们不能指望用越来越强大的励磁系统来无地继续补偿电抗的增加”[7］。快速励磁系统的确能够改善同步转矩，从而提高系统在第一摇摆周期的暂态稳定性.然而，快速励磁系统一般是高放大倍数的负反馈系统，它对第—摇摆周期以后系统振荡的阻尼影响很小,有时甚至减小系统对振荡的阻尼。在系统呈现负阻尼特性时，快速励磁系统（特别是高放大倍数时）通常是增大负阻尼,从而恶化系统的运行情况。

对于由m台发电机组成的互联电力系统来说,一般认为系统中机电振荡模态的总数为m1.根据对实际系统振荡的现场记录［48］和大量的仿真结果，将电力系统出现的振荡按振荡所涉及的范围及振荡频率大小大致分为两种类型［6，8]：局部模态(Local Modes)和区域之间模态（Interarea Modes）：

(1)局部模态涉及一个发电厂内的发电机组与电力系统其他部分之间的摇摆。由于发电机转子的惯性常数较大,因此这种模态振荡的频率大致在1~2Hz范围内.

(2)区域之间模态涉及系统中一个区域内的多台发电机与另一个区域内的多台发电机之间的摇摆.联系薄弱的互联系统中接近耦合的两台或多台发电机之间常发生这种振荡。由于各区域的等值发电机具有更大的惯性常数，因此这种模态要比局部模态振荡的频率还要低，大致在0。1～0。7Hz范围内。当系统表现为两群发电机之间振荡时,振荡的频率大致在0.1～0。3Hz范围内；档系统表现为

多群发电机之间的振荡时，振荡的频率大致在0。4~0.7Hz范围内。

这两种类型的机电振荡，由于振荡频率较低,因此，也常称为电力系统的低频振荡。遭受小干扰的电力系统，除了机电振荡模态外,系统中的振荡还涉及控制模态和扭振模态。扭振模态在前面已经提及；控制模态涉及系统中的各种控制装置。由于控制装置的调节速度较快,时间常数非常小．因此控制模态的振荡频率一般很高。这里我们仅关注电力系统的机电振荡模态，有关控制模态和扭振模态的分析和控制已超出本书的研究范围。

小的干扰能够引发电力系统发生低频机电振荡，只要所有的振荡模态是衰减的,系统就是小干扰稳定的。但是，在实际电力系统中,一般认为机电振荡模态的阻尼比大于0.05我们才可以接受系统的这种运行状态。然而这并不是硬性的、固定不变的原则。如果随着系统运行方式变化，模态变化很小，较低的阻尼比（例如,0.03）也是可以接受的［8］。

当今实际电力系统的小干扰不稳定大多数是由于阻尼不足而引发的振荡.1969年,Demello和Concordia[37]针对单机无穷大系统研究了可控硅励磁系统的作用，得到了系统小干扰稳定的条件，即电压调节器增益和引入发电机转速变化信号的传递函数特性应满足一定的要求.它清楚地揭示了单机无穷大系统增幅振荡发生的原因，为电力系统稳定器(PSS）的设计提供了坚实的理论基础。基于这种分析的原理和思想,研究人员将其推广到多机系统的局部振荡模态分析，并进—步推广到互联大系统区域之间振荡模态的分析。然而，应该特别指出的是,这种简单的推广经常是不合适的。

事实上，由于干扰所引发的增幅振荡受着许多因素的影响，大规模多机电力系统是—个典型的非线性动力学系统.网络的拓扑和参数、动态元件的特性、运行方式、各种控制器的控制策略从其参数等都对振荡过程起着重要的作用。要分析清楚电力系统发生增幅机电振荡的原因，并得到克服增幅振荡的有效措施是相当困难的。随着现代社会人们对经济性提出更高的要求，特别是电力系统正朝着市场化方向发展，要求电力系统在现有网络上承担着越来越多的负荷。然而，系统运行的经济性和安全性构成了一对矛盾。当受扰前的系统在较轻的负载下运行时,发电机的阻尼绕组能够产生与转子速度变化成正比的转矩，一般情况下它们能够完全吸收系统振荡所产生的能量,使得振荡幅值不断衰减,这时电力系统是小

干扰稳定的.而当干扰前的系统在较重的负载下运行时，发电机的阻尼绕组不能完全吸收系统振荡所产生的能量，这时振荡幅值不断增长，电力系统是小干扰不稳定的。另外，为了提高功率传输能力并改善暂态稳定性和一些其他性能，电力系统中加装了大量的、各种类型的控制器，其中的一些控制器由于采用了不适当的控制策略和参数，各控制器之间的配合不协调使得它们对振荡的作用发生冲突等等，这些都将导致系统发生增幅振荡。电力系统的振荡分析的目的皆在得到影响振荡模态的关键因素，从而为有效地抑制振荡打下基础.

【例7-2】例7—1已计算出状态矩阵的所有特征值，下面我们将对系统的振荡模态进行分析。几个振荡模态的振荡频率和阻尼比如表7-2所示，相应的左、右特征向量和参与向量如表7-3和表7-4所示.

【解】 下面我们根据表7—3和表7-4（其中的所有向量都规范化为无穷范数下的单位向量）进行摸态分析。首先利用指定模态的参与向量辨识机电振荡摸态:如果参与向量中模值最大的分量与发电机的转速有关，我们就说该模态为机电振荡模态.然后，我们依赖于右特征向量观察模态的表现形式：考察右特征向量与各发电机转速有关的分量，模值相差不大、方向基本相同的一组分量说明相应的发电机同群，不同群之间相应的分量方向基本相反。局部模态的右特征向量由与一台机或连接很近的一群发电机有关的几个变量所支配，区域之间模态的右特征向量分量更平均地分布于系统的所有区域。

对于振荡模态5，因此属于机电6,其参与向量中模值最大的元素对应于3，振荡模态。另外，在右特征向量中,对应于1、2的分量模值较小（分别为0。000 18和0。001 21），方向基本相同(辐角分别为-170。62°和—166．98°）;对应于3的分量，模值较大(为0。004 11)，方向（辐角为10。52°）和上述两个分量的方向基本相反.因此，该模态表现为发电机1和2与发电机3之间的机电振荡，振荡的频率为2.047 32Hz，属于局部振荡模态.

同样,对于振荡模态7，8，其参与向量中模值最大的元素对应于2，因此属

于机电振荡模态。另外，在右特征向量中，对应于2、3的分量模值较大（分别为0。002 61和0.001 41),方向基本相同(辐角分别为.08°和67。20°);对应于1的分量，模值较小（为0.00093）,方向（辐角为—113 60°)和上述两个分量的方向基本相反.因此，该模态表现为发电机1与发电机2和发电机3之间的机电振荡，振荡的频率为1。380 07Hz，也属于局部振荡模态.该模态虽然是稳定的，但阻尼比（0.017 47）不足,表现为振荡衰减的动态品质较差。

对于振荡模态11，12，其参与向量中模值最大的元素对应于VR2；对于振荡

'模态13,14,其参与向量中模值最大的元素对应于Eq因此说明这些振荡模态不属

于机电振荡模态，而是控制模态.

【例7—3】我们以New Eng1and 39节点、10机简化系统为例说明电力系统振荡分析的过程[40].其中10台发电机分别在30～39号节点，39节点的发电机为一台等值机,30~38号节点的发电机分别装有快速静止励磁系统。

【解】通过前面介绍的方法得到系统在稳态运行点的线性化方程，进而可计算出状态矩阵的所有持征值。9个与机电振荡有关的特征值阻尼不足，其中的—些特征值具有正实部。对于两个特征值0.1022j7.215（模态1）和

0.037j4.301(模态2）,其相应的右特征向量中对应于各发电机转速的分量如表

7-5所示.

由表7—5可以看出，第一个模态的特征向量中，有三个分量的模值较大（用黑体标出），并且前一个分量的方向（辐角为0°)和后两个分丝的方向(辐角约为170°）基本相反，从而说明该模态主要表现为发电机30与发电机35和36之间的机电振荡,振荡的频率为7.215/21.148Hz，因此属于局部振荡模态。第二个模态的特征向量中，除了发电机39对应分量的模值较小外，其他各发电机对应分量的模值差别不大，并且前9个分量的方向（辐角约为0°)和后一个分量的方向（辐角约为-180°)基本相反，从而说明该模态主要表现为发电机30~38与发电机39（等值机,可以看作一个区域)之间的机电振荡,振荡的频率为

4.301/20.685Hz，因此属于区域之间振荡模态。

参与向量除了用于辨识机电振荡外，还表示发电机控制机电振荡的相对价值。例如,参与向量的转子速度分量表明了特征值对作用于发电机轴系机械阻尼变化的灵敏度.如果参与向量中某台发电机速度分量为零,那么安装在该发电机上的PSS对该模态的阻尼无任何影响；如果参与向量中一些发电机的速度分量为正且数值较大时，说明这些发电机是安装PSS的很好候选位置，这样能够非常有效地增加对该模态的阻尼。表7—6是参与向量中对应于各发电机转速的分量。从表7-6中可以看出,对于局部模态，发电机30的分量最大，但发电机35、36的分量之相与发电机30的分量几乎相等；对于区域之间模态，发电机39的分量最大。在发电机30施加控制和在发电机35、36同时施加控制对于模态1阻尼的影响几乎是等价的；对于模态2，在发电机39施加控制似乎效果最好，但出于该机为等值机,因而无法施加控制。发电机30~38的分量之和大致与发电机39的分量相等，这意味着在发电机30~38同时施加控制可以提供与在发电机39施加控制同样的阻尼.另外，值得注意的是，虽然有些发电机的参与因子较大，但当该机容量很小时，在该机上施加控制对机电振荡阻尼的影响甚微。因此，在容量较大的发电机上施加控制比在容量较小的发电机上施加控制，对于增加对振荡的阻尼来说，效果要好得多。

7。6 大规模电力系统小干扰稳定分析的特殊方法

在电力系统小干扰稳定分析中，状态矩阵A一般是实不对称的稠密矩阵。现今大多数小干扰稳定分析都是使用著名的QR法计算A的全部特征值，进而判断系统的稳定性。

基于计算A的全部特征值的传统小干扰稳定分析方法有以下优势： （1）根据全部特征值能够清楚地分离并确定所有的模态.

(2）用特征向量能够容易地确定各模态的表现以及模态和状态变量之间的关系。

虽然著名的QR法鲁棒性好、收敛速度快，但在用于电力系统小干扰稳定分析时有许多固有的缺点。首先，QR法是一种基于稠密矩阵实现的特征求解方法,占用计算机内存大；另外,当矩阵A的维数特别大时，受计算机字长而产生的舍入误差的影响,QR迭代可能不收敛,有时即使收敛，误差也可能淹没真解，这些都将导致求解失败。计算经验表明，用QR法有效地进行特征求解的矩阵一般为几百至一千阶，因而其应用只能局限于小型电力系统.

现代电力系统不断朝着大规模的方向发展，系统规模可达上万个母线和几千个动态元件.以每个动态元件平均用10个状态变量描述来考虑，系统中总的状态变量可达几万个，这显然已超出传统的QR法有效地进行特征求解的范围。因而，需要一些针对超大规模系统特征求解的特殊方法，即能够求得系统全部特征值中我们所选择的特征子集，从而判断系统的稳定性或对系统的一些模态进行详细的分析，找出不稳定的机电振荡产生的原因,为抑制振荡的控制器设计打下基础. 目前，选择特征子集的求解方法大致分为两类：一类是降阶方法［11~16]。它首先将原系统降阶，使得降阶系统的特征值属于原系统的特征值,从而以较小的计算量对降阶系统进行特征求解，得到我们需要的特征子集。另-类则是直接在原系统中求得选择特征子集的特殊方法.如式（7-3)和式(7-83）所示,大规模电力系统小干扰稳定分析中所形成的增广状态矩阵是一个非常稀疏的矩阵,因此在计算选择特征子集时，我们应该充分利用稀疏矩阵技术,从而减少所用的计算机内存和计算费用。几年来,计算数学界对大型稀疏实不对称矩阵选择特征子集计算问题进行了广泛深入的研究，并提出了许多有用的方法［3,18~20］。近十几年来，电力工业界已将这些方法用于超大规模系统的小干扰稳定分析特别是振荡分析[5,6，

9，17,21~32]

, 并借助于线性系统模态分析的成果研究振荡发生的原因和影响振荡的

主要因素以及有效地抑制电力系统振荡的理论和法［33~47］。

本节首先介绍一种降阶选择模态分析法，然后介绍几种大型稀疏实不对称矩阵选择特征子集直接计算的数学方法，并将它们和电力系统的小干扰稳定分析，特别是振荡分析联系起来。 7。6.1 降阶选择模态分析法[5,6,11~13,16］

对于式(7-3）所描述的线性微分—代数方程组,将状态变量x按预定的原则分成x1，（维数为n1)和x2(维数为n2)，相应的方程改写为如下分块矩阵形式:

对上式进行拉普拉斯变换，然后消去x2和y得

式中：

如果将D1(s)当成一个常数矩阵，降阶矩阵（7-137)必须满足

其特征方程为

下面我们分析式（7-140）的根与原系统特征根之间的关系。为此，首先引出分

块矩阵行列式的Schuu’s公式：

AB(1)如果方阵J被分割成J,其中A、D都为方阵,并且如果detD0，CD则有detJdetDdet(ABD1C)。

(2)如果A、B为同阶方阵，则detABdetAdetB。

据此可得

这样,只要

式（7—140)的根必定就是原系统的特征根。这意味着求解降阶系统式(7—137）~（7-140）可以得到原系统的部分特征值。由于Ar(s)是随s的变化而变化的,因此式（7—140)的求解将不得不采用迭代方法。 2。 AESOPS算法

AESOPS是英文Analysis of Essentially Spontaneous Oscillations in Power Sysyems的缩写，它属于一类基于合理的工程判断而非基于严格数学理论的进行机电振荡分析的巧妙方法。AESOPS的基本思想是，在一台指定发电机（该发电机有时也被称为“受扰发电机”)转轴上施以具有复正弦形式的附加转矩，迫使线性动态系统产生振荡。类似于频率响应法的原理，求出关于机电振荡模态的一对复共扼特征值，附加转矩表示特征值的初始估计。在稳态情况下，线性化系统中所有变量与附加转矩具有同样的复正弦形式，通过求解相应的复代数方程组就可计算出系统的复频率响应.一般来讲,对指定发电机施以附加转矩所得到的特征值，表示指定发电机有较大参与的一个系统振荡模态.如果指定发电机有几个主导振荡模态，各模态的获得与附加转矩的初值有关。

AESOPS最初提出于1982年［13]，后经对算法的不断改进[16］，使其逐步趋于完善，现已成为EPRI开发的小干扰稳定程序(SSSP）包的一个重要组成部分［9］。

设式(7-136)中的x1仅包含指定发电机的状态变量和.，指定发电机施以附加转矩Tm后，与式（7-136)对应的方程变为

式中：

对式(7—143）作拉普拉斯变换，然后消去x2和y得

式(7—6）的第一式经拉普拉斯变换后得

它应该就是式(7-145）的第一式.因此,式(7—145）可表示为

消去上式中的变量，得

容易证明，当0时,方程

的全部零点就是降阶系统的特征值.因此，降阶系统的特征值就可以看作是当

0时，迫使Tm为零的s的取值.AESOPS算法就是要确定式（7-149)的零点。

取1.0j0，通过迭代求解,可以得到Tm趋近于零时s的值，即相应的特征值。

当使用牛顿法时，迭代公式为

其中的导数可由式（7—148）导出:

从式(7—149)可得Ks(s)/s2TJKD(s)/s，代入上式，并认为KD(s)、

dKD(s)/ds、dKS(s)/ds足够小而将它们忽略，这样式（7—151)所表示的导数可近似表示为

对于涉及多台发电机的振荡模态，dKD(s)/ds、dKS(s)/ds并不是很小，这就导致在每次迭代中对特征值的过修正。因此,常使用一个等值惯性常数TJeq，使得系统中所有发电机因速度变化产生的动能之和等于TJeq与指定发电机速度变化平方的乘积，即

这样，迭代式(7—150）变为

下面我们讨论AESOPS算法的具体实现。按式（7—1）进行迭代时,需要计算

(k)出指定发电机的附加转矩T（)和由该附加转矩引起的所有其他发电机的速ms(k)度变化(从而得到TJeq).这些计算取决于系统中所有动态元件的状态方程和电力

网络方程。由7。2节可知，任意一个动态元件的线性化微分-代数方程都可统一表示成

对上式中的第一式作拉普拉斯变换，得

将式（7-156）代入式(7—155)的第二式，得

式中:

式中:i为矩阵Ai的所有特征值组成的对角模态矩阵；XRi和XLi分别为相应的右特征向量按列及相应的左特征向量按行组成的模态矩阵。只要得到矩阵Ai的所有特征值和相应的左、右特征向量，式（7—158）的计算将变得非常简单。

对指定的发电机k施以附加转矩Tm，将状态变量划分为、及xr三部分，则系统的微分方程可表示为

发电机k的代数方程可表示为

对式(7-159）进行拉普拉斯变换,并注意到1.0j0.然后消去和xr,得到

Tm(s)和Ik(s)表达式：

式中:

其中，式(7—161)、式(7—163)、式(7-1）中(sIArr)1的计算可仿照式(7-158）中(sIAi)1的计算。

根据式（7—157）和式（7—162)，将系统中全部动态元件的代数方程统一表示为

式中:

电力网络方程可表示为

式中：下表D和小标L分别表示有动态元件和无动态元件的系统节点。将式(7-165）代入网络方程式(7—167)，从而消去ID,得

现将AESOPS算法的计算步骤归纳如下：

（1）形成各动态元件的模型Ai、Bi、Ci、Di。用QR法计算Ai的全部特

征值和相应的左、右特征向量。

(2)给定s的初始估计值s(0)。

（3)用式（7-163）、式(7-1）计算指定发电机的Ike(s)和Yke(s)，并用式（7-158)计算其他动态元件的Yie(s)。

(4)求解方程(7-168），得到各节点电压V.

（5)用式（7-161)计算指定发电机的附加转矩Tm，用式（7—156）计算所有其他发电机的i，最后用式（7-153)计算TJeq。 (6）用式(7—1）计算s的改进值.

（7)如果修正量s小于指定的误差，则迭代收敛；否则，转向步骤（3）. 由上述计算步骤可以看出，AESOPS算法每经一轮迭代，便可得到一个机电振荡模态，其中仅求解惟一的代数方程（7—168）.它不需要形成全系统的状态矩阵A，仅形成每个动态元件的状态方程，并可通过很小规模的特征求解计算出每个动态元件的Yke(s)。代数方程(7—168）的系数矩阵是高度稀疏的，因此可以利用稀疏求解技术高效地进行求解.显然，这种方法计算速度快，且不受系统规模的。

AESOPS算法的缺点是，受扰发电机的选取和s初值的确定没有系统化的方法。因此，除非关键模态的-般特性事先知道,否则要发现所有的关键模态，需要进行大量的试探，这样就不能保证一些关键模态不被遗漏。这种方法最适合于跟踪指定模态随系统运行方式的改变而变化的情况。 7.6.2 选择特征分析的序贯法和子空间法

本小节首先介绍两类计算矩阵ARnn的一个或一组模数最大的特征值的方法：序贯法(Sequential Method)和子空间法(Subspace Method);然后介绍将矩阵A变换为矩阵S的两种预处理技术,从而通过计算S的一个或一组模数最大的特征值得到A的一个或一组选择的特征值.

1。 序贯法

通过对矩阵不断的收缩处理，按模数从大到小的顺序一个一个地求出矩阵特征值和与之相应的特征向量的方法称为序贯法。序贯法主要有：幂法,Rayleigh

商迭代（RQI，Rayleigh Quotient Iteration)法及牛顿法。幂法在前面已有详细论述,这里主要介绍Rayleigh商迭代法和牛顿法。

（1）RQI［2,3，27]。

设v是一给定非零n维向量,我们可以选择使得Avv2达到最小。事实上

显然,当

时，达到最小值。纯量r(v)称为对应于v的Rayleigh商.显然,若v是A近似的特征向量，则r(v)是其对应特征值的较好估计。把这个思想与反幂法结合就得到RQI:

RQI是一种点位移方法，在每次迭代中需要矩阵的重新分解,因而其计算量将大大增加。RQI几乎总能收敛，并且可以证明，如果A是对称的，则RQI在解的邻域内立方收敛。同样，如果A非对称，设uH、v分别为A关于同一特征值近似的左、右特征向量，我们期望对于定义为

的r(uH,v)，向量Avr(uH,v)v和AHur(uH,v)u都是小量，r(uH,v)式特征值更

好的近似.一般将r(uH,v)称为对应于uH、v的“广义Rayleigh商”.这样，RQI也能立方收敛。 (2）牛顿法[21］。

我们还可以将特征分析问题看成如下方程的解：

其中对v有一些规范化约束，从而保证特征对有惟一解.可以用牛顿法求解上述

（0）（0）（0）方程组。给定和v，规定向量v中的最大元素为1（例如，第一个元素）,

并且这个元素在以后的迭代中保持不变。将式(7-172)线性化，即得到牛顿迭代式

或

它可以写成如下矩阵形式：

其中，向量v~(k)除第一个元素用(k)替代外，就是向量v(k)；矩阵B(k)除第一

列用v(k)替换外,就是矩阵A(k)I。和v的修正式为

用牛顿法求解特征值问题也可以稀疏实现，此时式（7-172）变成

其线性化方程为

将上式中最左边的矩阵简写为J(k)，同上理由，规定向量v(0)中的最大元素为1（例如，第一个元素），并且在以后的迭代中保持不变.这样，式（7-176）可写成如下矩阵形式：

其中，向量v~(k)除第一个元素用替换外，就是向量v(k)(k);矩阵J~(k)除第一列

用式（7—176)中(k)的系数向量替换外，就是矩阵J(k).

由于J表示成

式中：向量g(k)为(k)的系数向量与J(k)中的第一列向量之差；eT[1,0,,0]。利用Sherman—Morrison公式，可得式(7-177)的解

~(k)破坏了J(k)的稀疏结构，为了提高式（7-177）的求解效率，将J~(k)

只要对J(k)作稀疏三角分解,上式容易通过前代、回代计算之。

以上介绍的三种序贯法中,幂法鲁棒性好，但收敛很慢;与之相反，RQI和牛顿法收敛很快，但由于每步修正特征值,因此需要矩阵的多次三角分解,使得计算量很大.牛顿法由于其二次收敛性而颇具吸引力，但它需要良好的初值,用幂法进行几步这代的结果作为牛顿法的初值，不失为一种很好的选择. 2．子空间法

同时计算矩阵的一组模数最大的特征值和相应特征向量的方法称为子空间法。对于非对称矩阵的特征子集计算，下面介绍两种性能良好的子空间法：同时迭代法(Simultaneous Iteration Algorithm)和Arnoldi方法。 (1）同时迭代法［18，19，22]。

同时迭代法是幂法的一个直接推广，它可用来计算高维不变子空间。设r为指定的整数,满足1rn。给定n×r阶列正交矩阵Q0，如下迭代产生一系列的列正交矩阵{Qk}Cnr:

因此这种迭代也称为正交迭代.在合理的假设下，上述迭代中产生的子空间按

r1/r正比例收敛于A的按模数递减的前r个特征值所对应的不变子空间。 上述正交迭代中每次需要矩阵的正交分解，因而计算过很大。下面给出一种不需要向量正交化的不对称同时这代法（Lopsided Simultaneous Iterationn Algorithm,LOPSI)［18]。

LOPSI的每次迭代包括矩阵的左乘及重新定向，最后是归一化和误差检验。 设ARnn的特征值i满足

定义

k

式中：

A的n个右特征向量按列组成矩阵

这样，可得到

假设我们开始用m(mr)个的单位试探向量组成矩阵U，即

Uu1,u2,,umCnm,U左乘A，得

式中：Vv1,v2,,vmCnm可以被表示成特征向量的线性和:

式中：CaCmm,CbCnmm分别为系数矩阵。由式（7—181)一（7—183)可得

注意到Qb对V的贡献相对地小于对U的贡献，经过一定次数的迭代后，可假设Cb比Ca 小得多。

重新定向处理涉及矩阵BCmm的完全特征求解，B是如下方程的解：

式中：将式(7—183）、式（7—184）的U和V代入式（7-185），GUHU,HUHV。并忽略系数Cb，得

只要UHQa非奇异,从式（7-186）就可得到

上式表明,B的所有左特征向量按行组成的矩阵是Ca的近似，而B的特征值是a的近似.如果P是B的所有特征向量按列组成的矩阵，那么PCa1，因此

就给出了矩阵4的行特征向量的改进。

需要得到矩阵A的前r(rn）个模数递减的特征值和相应的特征向量，LOPSI的迭代步骤如下：

①指定参加迭代的子空间的维数m(mr)，构造m个的、单位初始试探向量，并组成矩阵U(1)[u1,u2,um]。置k1。 ②确定矩阵A连续左乘的次l。 ②计算V(k)AlU(k)。

④计算G(k)[U(k)]HU(k)和H(k)[U(k)]HV(k) ⑤求解方程G(k)B(k)H(k)得到B(k)。

⑥用QR法计算B(k)的特征值和特征向量，并将B(k)的特征值按模数递减排序,即B(k)(k)k)(k)k)(k)(k)(k)(2k)(mdiag{1,(,2,,r,r1,,m}，其中1相应的特征向量按特征值的顺序逐列排放组成矩阵P(k)。 ⑦计算U(k1)V(k)P(k)，并将U(k1)逐列归一化.

⑨收敛判断。比较B(k)和B(k1)(B(k)0)的前r个特征值，或比较U(k1)和

U(k)的前r个列向量。如果

或

则收敛.否则，置kk1,返回步骤②。

对LOPSI算法的实施作如下说明：

1）注意到参加迭代的子空间维数m要大于需要计算的模数递减的特征值的个数r，其中多出的mr个向量称为防卫向量（Guard Vector)。由于前r个特征值中的第r个特征值收敛最慢，即决定收敛速率的值m1/r最大。要使收敛加快，应使m1/r减小，相应地应使m增大。显然m越大,m1/r越小，但过大的m将使计算量大大增加，一般取m2r，从总体上来说可达到计算量和收敛速度的较好权衡。

2)一般在[-1.0,1.0]间产生n×m个均匀分布的随机数，并由此组成矩阵，各列经归一化后即得到U(1),这样基本上能够保证m个试探向量的性。 3)由于算法的收敛率主要由左乘淘汰下部特征向量的速度支配，因此大部分迭代可以省略重新定向处理，从而减少总的计算量.连续左乘的次数

1llmax,l1时为基本迭代，l1时为快速迭代，一般取lmax10.l可以取固定值，也可以按文献[18］中给出的方法决定。

4)一般总期望主导特征向量首先收敛。因此，开始总是对u1(k)进行收敛判断，一旦u1(k)收敛,就将它“锁定”(Locked），然后计算v1(k),并在后续的迭代中保持

(k)进行收敛判断。一般地说，当判定U(k)的u1u1(k)，v1v1(k)不变。下次将仅对u2前j1个向量已经收敛，它们将被锁定，仅对u(jk)进行收敛判断.这样一直迭代下

去，直到jr，并且ur(k)收敛，计算结束。使用这种被称为“锁定向量\"（Locking Vector）的技术,使得步骤③和步骤④中一些向量的内积不必重新计算，因而可节省计算时间。

(2）Arnoldi方法［3，20]。

如果矩阵ARnn非对称，则正交三对角化QTAQT一般不存在. Arnoldi方法的基本思想是,一列一列地产生正交阵Q，使得QTAQH为上海森伯矩阵。设Q[q1,q2,,qn]，比较AQQH两边的列可得

把上式右端的最后一项分离出来可得

式中：himqiTAqm，i1,2,,m。由此可知，若rm0，则qm1可通过下式来确定:

式中：hm1,mrm2。这些方程定义了Arnoldi算法。假定q1为给定的单位•2范数初始向量，称q1,q2,,qm为Arnoldi向量，它们构成了Krylov子空间K(A,q,m)的一组标准正交基：

迭代m步以后，算法的运行情况可用第m步Arnoldi分解来概括：

式中：

如果rm0，则由Qm的列组成了一个不变子空间，且(Hm）（A).问题是如何从海森伯阵Hm和Arnoldi向量组成的矩阵Qm中提取有关A的特征信息。

若y为Hm的单位特征向量,有Hmyy，则从式（7—193)可得

Tyrm2的大小可式中:xQmy。我们称为里茨值，x为相应的里茨向量。数em用来衡量误差界.

需要得到矩阵A的前r(rn)个模数递减的特征值和相应的特征向量，Arnoldi方法如下：

①初始化。选取步数m(mr),给定2范数单位初始向量q1.

②A rnoldi步骤：

③特征值计算.用QR法计算上海森伯矩阵Hm{hij}的特征值和特征向量。将Hm的特征值按模数递减排序，即：Hmdiag{1,2,,r,r1,,m}，其中

12m,相应的特征向量按特征值的顺序逐列排放组成矩阵P。

④置QmQmP,并将Qm逐列归一化为2范数单位向量。

⑤收敛判断。如果max{Aqiiqi,i1,2,,r},则收敛，退出；否则继续。 ⑧重启动。新的初始向量取为q的归一化系数。返回步骤②.

对Arnoldi方法的实施作如下说明:

1)一般在[—1.0，1。0]间产生n个均匀分布的随机数并组成向量，经归一化后即可作为初始向量q1。m大于r的理由与LOPSI算法中的理由类似.

2)“锁定向量”技术在Arnoldi方法中也可使用，从而节省计算时间。 3）迭代过程中，矩阵Qm容易失去正交性，这将使方法的收敛性很差。因此需要对Qm重新正文化,按步骤⑥的公式计算重新正交化的初始向量,可使Arnoldi方法的收敛性大大提高.

3。 预处理(Preconditioning）

预处理是通过某种线性变换将矩阵A变换为矩阵S，将A的一些关键特征值映射为S的一些主特征值,并且特征向量保持不变。这样,我们就可以利用前面介绍的方法计算S的一些主特征值，并通过反变换得到A的一些关键特征值.下面介绍两种线性变换、即位移求逆变换（Shift—invert Transformation)和凯莱（Cayley)变换［27]。 （1)位移求逆变换。

这是一种最常用线性变换，即

式中：s称为点位移量，它可以是实数或复数。

我们首先分析矩阵A和矩阵S的特征值和相应特征向量之间的关系.设矩阵

new1iqi,其中iAqiiqi,为向量

i1rA的特征值和与之相应的特征向量分别为A、vA，矩阵S的特征值和与之相应的特征向量分别为S、vS，显然

式（7—196）两端同时减去向量svA，并整理后可得

比较式(7-197)和式（7—198）不难看出

显然，这种变换将A的最靠近s的特征值映射为S的模数最大的特征值,并且特征向量不变。如果用序贯法或子空间法能够求出S的前r个模数递减的特征值,则它们是A的到s的距离由近到远的r个待征值.这显然适应于我们所需要的选择特征分析。

很明显，将S用于幂法就是前面介绍的反幂法。这种变换也可以用于RQI和牛顿法的每次迭代中，不过这时位移量s是变化的。由于采用这种变换增大了映射的特征值之间的距离，因而提高了特征求解方法的收敛性,这是它的优点。其缺点是，当A是实矩阵，而位移s为复数时，需要复数运算。另外,在用于电力系统时，为了得到所有的关键特征值,需要多次位移从面扫描整个虚轴，这就不可避免地需要做一些冗余计算,并且在理论上不能保证一些关键特征值不被遗漏.

（2)凯莱变换。

取复常数s1和s2（s1s2),凯莱变换是这样的线性变换:

为简单起见，在实际应用中一般取s1和s2为实常数（以后如无特殊说明，均设s1和s2为实常数且s2s1。

不难得到采用凯莱变换后的矩阵S和原始矩阵A的特征值和相应特征向量之间的关系.式(7—196)两端分别同时与向量s1vA相减和与向量s2vA相减，得

由上式容易得到

比较式（7—197）和式（7-202)不难看出

令Aj，将它代入式（7-203），得

显然

因此，凯莱变换将复平面上横坐标为

s1s2的直线（称为对称轴)映射到单2位圆上,对称轴右侧的复平面映射到单位圆以外，对称轴左侧的复平面映射到单位圆以内，如图7—1所示.并且特征向量保持不变.

现在,计算矩阵A的位于对称轴右侧的所有特征值可以这样进行：用序贯法一个一个地或用子空间法一次求出矩阵S的模数大于1的全部特征值，然后可通过式(7—203)计算出矩阵A的对应特征值。和位移求逆变换相比，只需—次凯莱变换就可计算A的所有关键特征值，相应地仅需一次矩阵分解。

当取s2s1h(0)，即得到所谓的S矩阵法[17］。显然这时的对称轴就是虚轴，变换将虚轴映射到单位圆上，右半复平面映射到单位圆以外，左半复平面映射到单位圆以内.设Smax为矩阵S模数最大的特征值，由式(7-205)可知，当

Smax1时，说明短阵A的全部特征值都具有负实部；当Smax1时，说明矩

阵A的全部特征值中至少有一个特征值具有正实部；当Smax1时，说明矩阵A无具有正实部的特征值，但至少有一个实部为零的特征值.以上三种情况分别对应于系统稳定、不稳定和临界。因此，只要求出短阵S模数最大的特征值Smax，便可判断系统的小干扰稳定性。由此可见，S矩阵法对于单纯地进行小干扰稳定性判断来说是非常有效的.

凯莱变换的诸多优点并不能掩盖其缺点。首先，这种变换减小了映射后特征值之间的距离，因而降低了特征求解方法的收敛性;其次,由于A的一对复特征值模数相同，一些实部、虚部完全不同的特征值也可能模数相同,因而导致S中可能有很多模数相同的特征值，使得S的按模数递减的特征值求解困难。 在实际应用中，我们不仅关心系统是否稳定，还想要清楚地知道A的最右端(实部最大)的特征值，特别是在进行低频振荡分析时，关心那些具有负阻尼(实

部为正)或阻尼不足（实部为负，但接近虚铀)的特征值.对此，一种自然的想法是用某种方法能够得到S的模数递减的部分特征值，并用式(7-20 3)计算出A的相应特征值。但遗憾的是，A的从右到左的特征值和S的按模数递减的待征值并不存在一一对应的关系。对此,一种可行的方案是选择合适的s1和s2,使对称轴右侧的全部特征值包含我们所关心的特征值，这样只要求出S的模数大于1的所有特征值，通过反变换就可得到我们所关心的待征值。显然,这种方案可能计算出许多很大的待征值，对于进行低频振荡分析来说，这些待征值并不是我们所关心的,因而这种方案可能包含大量的冗余计算。另外，如果想要得到A的最右端（实部最大)的特征值，通过试探总可得到合适的s1和s2，使得对称轴右侧仅有一个特征值，这样就将A的实部最大的特征值映射为S的主特征值。

对凯莱变换后的矩阵S,用幂法求其主特征值和相应的特征向量时有以下迭代式：

直接计算出S,进而实施式（7-206）的迭代是不经济的。受参考文献［24]的启发，下面介绍一种简单、有效的迭代方案。根据式(7—200）可得

因此，式（7-206）可改写为

上式两端左乘(As1I)得

经整理后可写成

这样,式（7—206）的求解可归结为

另外，由于

因此，在幂法的迭代式(7-121）中：

用幂法计算矩阵S的主特征值和相应特征值向量的迭代式如下：

值得注意的是,在式（7-211）每步迭代的最后,计算出的向量（s2s1）(As2I)1v(k)就是下一步所要计算的向量y(k1)。另外，利用7。4。4节介绍的方法，容易得到迭代式（7—211）的稀疏实现。