第六章-电力系统暂态稳定分析

6.1概述

在正常的稳态运行情况下，电力系统中各发电机组输出的电磁转矩和原动机输入的机械转矩平衡，因此所有发电机转子速度保持恒定。但是电力系统经常遭受到一些大干扰的冲击,例如发生各种短路故障,大容量发电机、大的负荷、重要输电设备的投入或切除等等。在遭受大的干扰后,系统中除了经历电磁暂态过程以外，也将经历机电暂态过程.事实上，由于系统的结构或参数发生了较大的变化,使得系统的潮流及各发电机的输出功率也随之发生变化，从而破坏了原动机和发电机之间的功率平衡，在发电机转轴上产生不平衡转矩，导致转子加速或减速。一般情况下,干扰后各发电机组的功率不平衡状况并不相同，加之各发电机转子的转动惯量也有所不同、使得各机组转速变化的情况各不相同。这样，发电机转子之间将产生相对运动，使得转子之间的相对角度发生变化，而转子之间相对角度的变化又反过来影响各发电机的输出功率，从而使各个发电机的功率、转速和转子之间的相对角度继续发生变化。

与此同时,由于发电机端电压和定子电流的变化，将引起励磁调节系统的调节过程;由于机组转速的变化，将引起调速系统的调节过程；由于电力网络中母线电压的变化，将引起负荷功率的变化；网络潮流的变化也将引起一些其他控制装置(如SVC、TCSC、直流系统中的换流器）的调节过程，等等。所有这些变化都将直接或间接地影响发电机转抽上的功率平衡状况.

以上各种变化过程相互影响，形成了一个以各发电机转子机械运动和电磁功率变化为主体的机电暂态过程。

电力系统遭受大干扰后所发生的机电暂态过程可能有两种不同的结局。—种是各发电机转子之间的相对角度随时间的变化呈摇摆(或振荡)状态，且振荡幅值逐渐衰减，各发电机之间的相对运动将逐渐消失，从而系统过渡到一个新的稳态运行情况，各发电机仍然保持同步运行.这时，我们就称电力系统是暂态稳定的.另—种结局是在暂态过程中某些发电机转子之间始终存在着相对运动，使得转子间的相对角度随时间不断增大、最终导致这些发电机失去同步.这时称电力系统是暂态不稳定的。当一台发电机相对于系统中的其他机失去同步时，其转子将以高于或低于需要产生系统频率下电势的速度运行,旋转的定子磁场（相应于

系统频率)与转子磁场之间的滑动将导致发电机输出功率、电流和电压发生大幅度摇摆，使得一些发电机和负荷被迫切除，严重情况下甚至导致系统曲解列或瓦解.

电力系统正常运行的必要条件是所有发电机保持同步.因此，电力系统在大干扰下的稳定性分桥，就是分析遭受大干扰后系统中各发电机维持同步运行的能力，常称为电力系统的暂态稳定分析。

上述对电力系统的暂态稳定分析通常仅涉及系统在短期内(约10s之内)的动态行为，然而有时我们还必须分析系统的中期(10s直至几分钟)和长期(几分钟直至几十分钟)动态行为,这就涉及到电力系统的中期和长期稳定性分析.

中期和长期稳定性主要关注在遭受到严重破坏时电力系统的动态响应.当电力系统遭受到严重破坏时，将导致系统的电压、频率和潮流发生重大偏移，因此必然涉及到一些在短期暂态稳定分析时未曾考虑的慢过程、控制及保护的行为.对电压和频率发生大的偏移起作用的装置，其响应过程从几秒（如发电机控制与保护装置的响应）到几分钟（如原动机能量供应系统和负载电压调节器等装置的响应)。

进行长期稳定性分析的重点是与大范围系统破坏同时发生的较慢的、持续时间长的现象.以及由此引起的发电机与负荷的有功功率和无功功率显著的持续性失配。这些现象包括：锅炉的动态，水轮机的进水口和水管功态,自动发电控制（AGC），电厂和输电系统的控制与保护,变压器饱和，负荷和网络的非正常频率效应等。长期稳定通常关心系统对特大干扰的响应，这些干扰不属于正常系统设计准则的预想事故。在这种情况下，可能引发连锁事故及系统被分离成几个孤立的子系统。这时稳定分析要回答的问题是如何在负荷损失的情况下各孤岛能达到可以接受的平衡状态。

中期响应是指短期响应向长期响应的过渡.中期稳定研究的重点是各机之间的同步功率振荡，包括一些慢现象以及可能的大的电压和频率偏移[4]。

电力系统遭受大干扰是人们所不希望的，但事实上又是无法避免的.系统在遭受大干扰后失去稳定的后果往往非常严重,甚至是灾难性的。事实上电力系统遭受到的各种大干扰，诸如短路故障．大容量发电机、大的负荷、重要输电设备的投入或切除等都是以一定的概率随机地发生，因此系统的设计、运行方式的制

定总是需要保证系统在合理选择的预想事故下能够保持稳定，而不能要求电力系统能承受所有干扰的冲击.由于各国对系统稳定性的要求不同，因此对预想事故的选择也就有不同的标准。我国对系统稳定性的要求反映在《电力系统安全稳定导则》[3］中.

判断电力系统在预想事故下能否稳定运行,需要进行暂态稳定分析。当系统不稳定时。还需要研究提高系统稳定的有效措施;当系统发生重大稳定破坏事故时，需要进行事故分析，找出系统的薄弱环节，并提出相应的对策。

下面首先讨论电力系统暂态稳定分析所用全系统数学模型的构成[1,2,4，6,25］。 在电力系统稳定分析中，各元件所采用的数学模型，不但与稳定分析结果的正确性直接相关，而且对稳定分析的复杂性有很大的影响。因此，选用适当的数学模型描述各元件的特性，使得稳定分析的结果满足合理的精度要求并且计算简单,是电力系统稳定分析中一个至关重要的问题。对于包含众多发电机、输电线路、负荷及各种控制装置的实际电力系统．考虑到任何冲击后果的复杂性，使得各元件的建模遇到很大的困难.所幸的是，各种现象时间常数的明显差别允许我们把注意力集中在影响暂态过程的关键元件和所研究区域。

在进行电力系统稳定分析时,由于在遭受干扰后电力网络的电磁暂态过程衰减很快，因此忽略其暂态过程是合理的。采用这种简化后，电力网络的模型中就仅包含代数方程。另外,在发电机定子电压方程中，d和q反映了定子绕组本身的暂态过程，忽略这两项，意味着忽略了定子中的直流分量，因此定子中仅包合基频电气分量，定子电压方程也就变成代数方程。很明显,同时忽略发电机定子和电力网络的暂态过程，能够使得定子电压方程和网络方程保持一致，即均为代数方程,且仅包含基频电气分量，因而可以用稳态关系式描述，这样做显然还使全系统微分方程的数目大大减少，从而可提高系统稳定分析的效率。由于系统中所有的电气量在交流系统中是基波交流分量的有效值，故可用相量描述（用大写字母表示)；在直流系统中是直流分量的平均值。描述各元件电压、电流关系的方程都为代数方程（和潮流计算中的稳态方程相同）；由于系统中动态元件的存在，一些电气量表现出一定的动态特性。因此，在遭受干扰后,电力系统经历的整个暂态过程可以看成是各时刻的稳态量（正弦交流量）按一定动态特性的过渡,这时系统中的电压、电流、功率能够发生突变.这就是电力系统稳定分

析常用的准稳态模型(Quasi—steady state Model)。

图6-1给出了用于电力系统稳定分析的全系统数学模型的构架。由图6—1可以看出，全部电力系统的表达式包括描述同步发电机、与同步发电机相关的励磁系统和原动机及其调速系统、负荷、其他动态装置等动态元件的数学模型及电力网络的数学模型.很明显,系统中的所有动态元件是相互的，是电力网络将它们联系在一起。

整个系统的模型在数学上可以统一描述成如下一般形式的微分—代数方程组：

式中：x表示微分方程组中描述系统动态特性的状态变量；y表示代数方程组中系统的运行参量。

微分方程组(6—1)主要包括：

（1)描述各同步发电机暂态和次暂态电势变化规律的微分方程。 （2)描述各同步发电机转子运动的摇摆方程。

（3)描述同步发电机组中励磁调节系统动态特性的微分方程. （4)描述同步发电机组中原动机及其调速系统动态特性的微分方程. （5)描述各感应电动机和同步电动机负荷动态特性的微分方程。 （6)描述直流系统整流器和逆变器控制行为的微分方程。

(7)描述其他动态装置(如SVC、TCSC等FACTS元件）动态特性的微分方程。 而代数方程组（6—2)主要包括：

(1）电力网络方程,即描述在公共参考坐标系x—y下节点电压与节点注入电流之间的关系。

(2)各同步发电机定子电压方程(建立在各自的d—q坐标系下)及d-q坐标系与x—y坐标系间联系的坐标变换方程。

(3）各自流线路的电压方程。 （4)负荷的电压静态持性方程等.

根据对计算结果精度要求的不同，可依据所研究问题的性质，本着抓住重点、忽略次因素的原则使用相应复杂程度的元件数学模型。

目前，电力系统暂态稳定分析方法基本分为两种。第一种方法是数值积分方法，又称间接法[26～32］，其基本思想是用数值积分方法求出描述受扰运动微分方程组的时间解，然后用各发电机转子之间相对角度的变化判断系统的稳定性.数值积分法由于可以适应各种不同详细程度的元件数学模型．且分析结果准确、可靠,所以得到了广泛的实际应用，并一直作为一种标准方法来考察其他分析方法的正确性和精度。目前，利用数值积分法进行电力系统暂态稳定分析已经相当成熟,并已有许多商业性程序相继问世。如我国电力科学研究院编制的《交直流电力系统综合计算程序》,由BPA根据美国WSCC标准开发的暂态稳定分析程序。PTI开发的PSSE，美国EPRI的ETMSP，TRACTEBEL／EDF开发的EUROSTAG，巴西CEPEL的ANATEM及联邦德国的VISTA程序[30］和比利时的STAG程序[31］等。这些程序除可用于分析故障后转子的摇摆过程外，还可用于各种动态行为分析，它们已成为规划和运行人员进行离线暂态稳定分析、安全备用配置、输电功率极限估计的有力工具。

另一种方法是直接法，它不需要求解微分方程组，而是通过构造一个类似于“能量”的标量函数,即李雅普诺夫函数，并通过检查该函数的时变性来确定非线性系统的稳定性质，因此它是一种定性的方法。由于构造李雅普诺夫函数比较困难，因此目前电力系统暂态稳定分析的直接法仅限于比较简单的数学模型，或用暂态能量函数近似李雅普诺夫函数，因此其分析结果尚不能令人完全满意。

本章首先介绍暂态稳定分析中全系统数学校型的构成和微分-代数方程组的数值求解方法,然后叙述各动态元件与电力网络的连接以及网络操作及故障的处理方法.接着对简单模型和带有FACTS元件的详细模型下的电力系统暂态稳定分

析算法分别进行了详细论述.最后介绍暂态稳定分析的直接法. 6.2 暂态稳定分析数值求解方法［25］

电力系统的暂态稳定分析可以归结为微分-代数方程组的初值问题。本节我们首先介绍常微分方程的数值解法,然后讨论微分—代数方程组的数值解法，最后给出暂态稳定分析的基本流程。 6。2。1 常微分方程的数值解法[1，14～16]

1. 基本概念

考虑一阶微分方程

一般地讲,上式中f是x、t的非线性函数。在很多工程实际问题中，函数f中不显含时间变量f,因此往往表现为以下的形式：

在电力系统稳定计算中，所有微分方程都不显含时间变量t。

当式（6-4)中的f为t的线性函数时,可以很容易地得到微分方程解的解析表达式。例如,对微分方程式

可以求出它的通解为

式中：A为积分常数.式(6—6）表示了一个曲线族。

根据初始条件x(t0)x0可以确定x随t变化规律的一条曲线.例如,当x(0)1时，从式(6—6）即可确定积分常数A1，这样就得到了确定的解(或积分曲线)

工程实际问题所表现出来的微分方程比较复杂，其函数往往是多元非线性的，因此一般不能用解析的形式求出像式(6-6)那样的通解，而只能用数值解法。即从已知的初始状态（tt0,xx0）开始，利用某种数值积分公式离散地逐点

求出时间序列tnt0nh，n1，2,…（h为步长）相对应的函数的近似值：xn。对微分方程的这种数值解法称为逐步积分法。

以下我们以欧拉法为例子说明逐步积分法的基本概念。 设一阶微分方程式（6-3）在:。t00、x(t0)x0时的准确解为

这一函数曲线，即微分方程式(6-3)通过点（0，x0）的积分曲线如图6-2所示。

欧拉法又称为欧拉切线法或欧拉折线法。它的基本思想是将积分曲线用折线来代替，而每段直线的斜率都由该段的初值代入式（6-3）求得.具体推算步骤如下:

对于第一段，在点（0，x0）处曲线的斜率为

将第一段曲线用斜率为的增量为

dx的直线段来代替，则可以求出t1h（h为步长)时xdt0

因此在t1h处，x芙蓉近似值应为

对于第二段，积分曲线将用另一段直线来代替，其斜率由该段的初值［即该段的起始点(t1，x1)]代入式（6-3)而得，即

这样便可以求出在t22h处x的近似值

如图6—2所示.这样继续下去又可以推算出t3处函数近似值xd，等等。一般．对于第n+1点函数值的递推公式为

现在我们来分析利用这个递推公式由（tn，xn)点推算(tn+1，x n+1）时带来的误差.为此可把积分函数式（6—8)在该点展开为泰勒级数

'\"式中：xn、xn…分别为积分函数对自变量t的一阶导数、二阶导数……在ttn点

(r)的值。n为区间[tn，tn+1]中的某一数,x为泰勒级数的余项。当取r＝2时，式（6n

—10）变为

或者写成

'这里n'仍为区间［tn，tn+1]中的某一数，一般nn。

d2xh2显然，忽略式(6-12)中余项2以后就得到欧拉法的递推公式（6—9）。

2!dtn'因此，在由n点推算n+1点函数值时所引起的误差为

d2x设整个计算的区间[0,tm］内，2f'(x,t)的最大值为M,则误差En+1应满足

dt

式中：M值与步长h的选择无关.式(6—13）、式(6-14)中的误差En+1是由n点推算n十1点函数值时引起的误差，称为局部截断误差.欧拉法的局部截断误差与h2成比例，通常说它的局部截断误差是o(h2)阶的。

应该指出，在计算xn1以前，xn也是用同一递推公式求得的,所以xn本身就有误差，因此在用式(6—9)计算xn1时,除了忽略余项而引起的局部截断误差以外，还应加上xn误差的影响。这个误差叫做全局截断误差或简称截断误差，因此，由于欧拉法递推公式不精确而引起的误差要比式(6-13）、式(6-14)所表示的局部截断误差大。可以证明，欧拉法的全局截断误差是和步长h成比例的，或者说它是o(h)阶的。

由以上讨论可以看出，为了减小欧拉法的计算误差，应该选择较小的步长h.但绝不能由此得到步长愈小则计算误差愈小的结论，因为在以上的讨论中，我们完全没有考虑计算机本身由于有效位数的而引起的舍入误差。当取较小步长h时,将使运算量成反比地增加，从而使舍入误差的影响加大。如图6—3所示,图中hmin为最小误差所对应的步长，因此，我们不能单单用缩小步长的方法来减小误差。当计算精度要求较高时，必须选择更完善的计算方法。

在以上欧拉法的计算过程中，当计算tn+1点的函数值时，仅需利用它的前一点tn处的函数值,这种方法称为单步法。本节介绍的方法都属于这一类。与此对应的是多步法（或多值法),这类方法的精确度较高，它在推算tn+1点的函数值xn+1时需要利用前而几点的数据：tn,xn,tn1,xn1,…,tnk1,xnk1。

2。 改进欧拉法

在应用欧拉法时，由各时段始点计算出的导数值

dxf(xn,tn)被用于[tn，dtntn+1]的整个时段，即代替积分曲线的各折线段的斜率仅由相应时段的始点决定，因而给计算造成较大的误差。如果各折线段斜率取该时段始点导数值与终点导数值的平均值，我们就可以期望得到比较精确的计算结果。改进欧拉法就是根据这个原则提出来的计算方法。

对于一阶微分方程式(6-3),设给定初值为t00时x(t0)x0,以下介绍改进欧拉法的具体步骤。

为了求t1h时的函数值x1，首先用欧拉法求x1的近似值

式中：

当x1(0)由式（6-15）求得以后，即可将t1，x1(0)代入式（6-3）求出该时段末导数的近似值

dxdx然后就可以用和的平均值来求x1改进值

dt0dt1(0)

这样求得的x1(1)比单纯用欧拉法求得的x1(0)更接近微分方程的正确解x1，其几何解释如图6—4所示。

当由tn,xn点推算tn1,xn1点时，递推公式的一般形式为

由式(6-17）中第二式及第四式消去xn，可将第四式改写为

式中：

这样，也可以把改进欧拉法的递推公式归结为以下形式：

(0)当应用式（6—19）计算xn1时，其形式与xn1的公式具有相向的形式,因此可以(0)简化程序，并且在求得xn1以后不必再记忆xn,因此也节省了内存单元.

以下讨论改进欧拉法递推公式的局部截断误差。 为此，仍需要利用式（6—l0）的泰勒级数展开式

h3式中：x\"为泰勒级的余项.

n3!\"' 改进欧拉法递推公式（6-17)中第四式可以改写为

将式(6-17）中第一式代入上式,可得

把上式中右端第三项按泰勒级数展开:

因为

所以

将上式代入式（6-21）中，则得

再把上式与式（6—20)相减，可知

因此，改进欧拉法的局部截断误差o(h3)阶。同样可以证明改进欧拉法的全局截断误差是o(h2)阶的。

【例6—1】 用改进欧拉法求解微分方程

其初值为t00，x01。

【解】 步长取0。2.计算结果见下表:

这一微分方程的准确解为

当t=1时,x=1.732 05，故误差为

改进欧拉法也可以用来求解一阶微分方程组。例如,对于微分方程组

其初始值为t0，x0，y0。当选定步长h以后，对于第一时段可以求出变量的近似值为

式中：

再由t1h,x1(0)，y1(0)求出

这样，函数在t点的值应为

式中：

以此类推。

由递推公式(6—17）可以看出,改进欧拉法计算一个时段所需要的运算过比欧拉法大一倍，但是如果步长一样，改进欧拉法的计算精确度却比欧拉法高.如上所述,改进欧拉法的截断误差是o(h2)阶的，而欧拉法是o(h)阶的。如图6—5所

示，当容许误差为1时，改进欧拉法容许步长h1'和欧拉法步长h1相差不大，在这种情况下，用改进欧拉法的运算量比欧拉法要大。当容许误差为2时，改进欧

''拉法的容许步长h2比欧拉法步长h2相对大得多，显然当h22h2时,改进欧拉法的

总运算量比欧拉法要小.

3. 龙格-库塔法

改进欧拉法用［tn，tn+1]区间两点的导数（或斜率）推算xn1，拟合了积分函数泰勒级数的前三项，从而使局部截断误差达到了o(h3)阶。这就启发人们去考虑：是合可利用[tn，tn+1］区间上更多点的导数去推算xn1，以便拟台泰勒级数更多的项数?结论是肯定的。龙格—库塔法就是基于这种原理建立起来的微分方程数值解法。最常用的是四阶龙格-库塔法，这种方法用［tn，tn+1］区间四个点的导数去推算xn1，从而拟合了泰勒级数的前五项：

因此，它的局部截断误差是o(h5)阶的．其全局截断误差是o(h4)阶的。

对于一阶微分方程式（6-3)，当利用四阶龙格-库塔法求解时，可以利用递推公式

求出x1，x2，x3，…。

【例6-2】 用四阶龙格-库塔法求解例6—1中的一阶微分方程 【解】 步长取h0.2。计算过程及结果如下表所示:

由以上计算结果可知，当采用龙格-库塔法时，函数值在t1时x1.732141,和准确解相比，其误差为

和例6-1相比精确度提高很显著。

应用龙格-库塔法也可以求解一阶微分方程组。例如,对于式（6-23）所示的微分方程,可按以下递推公式进行计算：

式中：

龙格-库塔法的精度较高，但运算量较大，为欧拉法的4倍。目前，当精度要求较高时，已逐步趋向于采用运算量较小的多步法来代替龙格-库塔法.龙格—库塔法往往只作为多步法起步时的一种辅助计算方法。

4. 隐式积分法

微分方程数值解法可以分为显式解法与隐式解法两大类。目前所介绍的方法都属于显式解法.分析它们的计算公式（6—9)、式(6-17）、式(6-24）即可看出,这些公式等号右端都是已知量，因此利用这些递推公式可以直接计算出相应时段终点的函数值xn1。与此不同,微分方程的隐式解法不是给出递推公式，而是首先把微分方程化为差分方程，然后利用求解差分方程的方法确定函数值xn1。

现在我们来介绍隐式梯形积分法。

对于微分方程式（6-3)，当tn处函数值xn已知时，可以按下式求出tn1tnh处的函数值xn1：

上式中的定积分相当于求图6—6中阴影部分的面积.当步长h足够小时，函数

f(x,t)在tn到tn1之间的曲线可以近似地用直线来代替，如图中虚线所示。这样，

阴影部分的面积就近似为梯形ABCD的面积，因此式（6-25）可以改写为

这就是隐式梯形积分法的差分方程。

显然，在这种情况下已不能简单地利用递推运算求出xn1，因为式(6-26）等号的右端也含有待求量xn1。这时必须对式(6-26)采用求解代数方程式的方法去计算xn1。

一般地说，微分方程隐式解法的特点就是把微分方程的求解问题转换成一系列代数方程的求解过程。例如,当初始值t0、x0给定时,根据式(6-25)可以得到第一时段的差分方程式

上式中只有x1为未知数，因而利用求解代数方程式的方法即可求得x1.当t1、x1已知后由式(6-26）又可得到第二时段的差分方程式

由上式又可解出x2，以此类推。

如果我们把f(xn,tn)和f(xn1,tn1)理解为积分曲线在［tn,tn+1］区间始点和终点的斜率，那么就有理由把隐式梯形积分法称为隐式改进欧拉法，也就是说,差分方程式(6—26）可以理解为改进欧拉法的隐式解法.实际上，隐式解法不限于改进欧拉法，前面介绍的欧拉法、龙格—库塔法以至多步法都可以采取隐式解法。例如,把欧拉法的递推公式(6-9)改为

''即把[tn,tn+1］区间始点的导数值xn改为终点的导数值xn1，我们就得到了隐式欧

拉法。式(6-27）就是隐式欧拉法的差分方程。

差分方程式(6—26)、式(6-27)可能是非线性的，因为微分方程式（6-3)中给山的函数f(x,t）可能是非线性的，因此，隐式解法比显式解法的求解过程要复杂一些.

顺便指出,隐式梯形积分法的截断误差可以解释为是由于以梯形面积代替阴影部分面积引起的（见图6—6）。利用前面介绍的方法同样可以证明差分方程（6—26)的局部截断误差也是o(h3)阶的。

隐式解法相对于显式解法来说的优点是可以来取较大的步长。这个问题涉及到微分方程数值解的稳定性问题，读者可参看有关文献。我们在这里只用一个简单的例子来直观地说明这个问题。

设有一阶微分方程

初值为t0时x01。

对于这个微分方程，不难求出她的解析解为

这是一个按指数曲线衰减很快的函数，如图6—7所示。

当步长取h0.025时，用欧拉法计算结果如下表所示：

可以看出，上表所列函数值随时间在作振荡的变化，而且振荡的幅值愈来愈大，如图6—7所示。从数学上来说,这种情况表示欧拉法数值解本身已经不稳定。

当采用隐式欧拉法计算时，就不会出现这种情况.首先将式（6-28)化为差分方程

因此可以得到

当h0.025时上式变为

因此可以得到下表计算的结果：

上表所列函数值随时间单调衰减,如图6-7所示。

为了进一步说明以上现象与步长的关系，我们把微分方程式（6—28)写成更一般的形式

式中：常数T具有时间的量纲，在工程上叫做时间常数.

当采用欧拉法时，将式(6—29）代入递推公式（6—9）,得

因此

显然要使函数x成为单调衰减的函数,式（6-30)右端括弧中的值必须满足

因此，步长的选择应满足

对式（6-28）来说，在采用显式解法时，要得到合理的解步长的选择则必须满足

h0.01

移项整理后可得

因此我们有

由上式可以看出，当采用隐式欧拉法时，选择任何步长都可以满足使函数x成为单调衰减函数的要求。

一般地说,在采用显式积分法时，步长的选择要受到微分方程中最小时间常数的，否则就会导致错误的计算结果.隐式积分法的步长则没有这个,容许选择较大的步长。

5. 常微分方程求解方法的选取

如前所述，常微分方程初值问题的数值解法可分为显式法或隐式法及单步法或多步法。在显式法中，积分公式可直接用于对每个微分方程进行求解,因此计算虽小，但数值稳定性差；在隐式法中，各微分方程先被差分化,然后对导出的代数方程组联立求解,显然它的计算更复杂，然而可获得更高的数值稳定性.单步法仅用前一步的信息得到本步的解，因而它是自启动的，便于处理不连续的情况；多步法用前几步的信息得到本步的解，因而原理上效率更高，然而当发生不连续时需要重启动。至于采用何种数值积分公式求解微分方程，至少应从以下三个方面考虑:

（1）方法的精度。

在用数值方法求解微分方程时至少存在两种误差，即舍入误差和截断误差。舍入误差的存在是由于计算机不能精确地表示浮点数所致，因而要减小舍入误差,

惟一的办法是使用大字长的计算机。

解的真值与计算值之差主要由截断误差决定，它体现了方法的精度。截断误差的大小与所用的数值积分公式有关，显然方法的阶数越高，在同样的步长下其计算精度越高.另外,随着积分的逐步进行，局部截断误差也在不断积累，因此把某时刻解的真值与计算值之差称为累积截断误差。 （2）方法的数值稳定性。

粗略地讲，积分方法的数值稳定性是关于在逐步积分过程中误差的传播问题。一个不稳定的方法是指误差趋向丁不断地积累甚至放大，以至于可能淹没真解。一个稳定的方法对累积的误差不但不放大，有时甚至缩小。

有关积分方法数值稳定性的定义有很多，这里不作详细讨论。为比较起见，我们只粗略地将方法分为“很稳定”或“不很稳定”.一般情况下，隐式法相对于显式法有更好的数值稳定性。 (3）对刚性方程的适应性。

微分方程组的刚性（Stiffness)是指微分方程类似于代数方程病态的一种特性。一般情况下，如果微分方程组中最大时间常数与最小时间常数之比较大，我们就称该微分方程组为刚性的。更精确地讲,刚性是用线性化系统的最大特征值与最小特征值之比度量.

对于刚性问题,为了保证截断误差在安全的低数值，相对不太稳定的积分方法将需要很小的步长，从而准确地跟踪系统响应中变化快速的分量。而更稳定的积分方法由于每步可以承受更大的误差,为得到同样精度的解,可以采用较大的步长。

除了经典模型外，电力系统稳定分析中的微分方程组一般是刚性的,而且随着同步机模型详细程度的增加，时间常数的范围增大,刚性愈发明显。代数方程中也隐含刚性，特别是负荷为非阻抗形式时。另外，存在微分-代数方程的不连续性、发电机组中调节器的限幅环节等.

常微分方程初值问题的解法很多,但适合于电力系统应用的却较少，这主要由稳定分析中微分方程的特点决定。数值积分方法的选取主要从方法的计算速度、精度、数值稳定性、对刚性微分方程的适应性及计算的灵活性（处理不连续和限幅比较容易）等几个方面考虑。这方面的研究工作很多，而且已有许多方法

得到应用，目前认为较满意的方法包括欧拉法、改进欧拉法、龙格—库塔法等显式方法，及隐式梯形法。这些方法大多出现在当前流行的生产级商业软件包中. 6.2.2 微分—代数方程组的数值解法

在进行电力系统暂态稳定分析时，需要寻求的是微分—代数方程组的联立解，这里的关键问题是微分方程组（6—1)和代数方程组（6-2)的交接处理。为此,我们可以来用交替求解法或联立求解法，现分别介绍如下。

1. 交替求解法

在这种方法中，数值积分方法用于微分方程组,可地求出x，单独求解代数方程组得到y。显然，积分方法和代数方程的求解方法可以相互。一般情况下，x和y的求解按某种指定方式交替进行。在交替求解法中，微分方程组用显式法和隐式法求解也有所不同。下面给出在已知t时刻的量x(t)和y(t)，求

tt时刻的量x(tt)和y(tt)的两个例子［按电力系统计算中的惯例,用下标（t）

表示t时刻的计算值，而积分步长用t表示，下同］。

当用式（6—24)所示的显式四阶龙格—库塔法求解微分方程组时，计算步骤如下:

（1）计算向量k1tf(x(t),y(t)). （2）计算向量x1x(t)1k1，然后求解代数方程组0g（x1,y1）得到y1，2最后计算向量k2tf(x1,y1)。

（3）计算向量x2x(t)1k2，然后求解代数方程组0g（x2,y2）得到y2,2最后计算向量k3tf(x3,y3)。

（4)计算向量x3x(t)1k3，然后求解代数方程组0g（x3,y3）得到y3：，2最后计算向量k4tf(x4,y4)。

1(5）最后得到(xtt)x(t)(k12k22k3k4)，相应地求解代数方程组

60g(x(tt),y(tt))得到y(tt)。

当用式（6-26）所示的隐式梯形法求解微分方程组时，整个计算工作为求如下方程的联立解：

对此,非线性方程组的交替迭代求解步骤为：

]]（1)给定y(tt)的初始估计值y([t0应用式（6—33）得到x(tt)估计值x([t0t)，t)，

即求解方程

][1]（2)用x([t0t)和式(6-34)得到y(tt)估计值的修正值y(tt)，即求解代数方程

][0] （3）用y([1tt)代替y(tt),返回步骤(1),继续迭代,直至收敛.

]为了给出良好的初值，从而使迭代次数减少y([t0t),可以取前一步的值，也

可以用前几步的值通过外推得到。从以上的迭代序列不难看出，除非迭代次数无限大，否则最终得到的解x(tt)和y(tt)不会相一致，即它们不会以同样的精度同时满足式（6—33）和式(6—34)，由此造成的误差称为交接误差.显然，减少交接误差的惟一方法是增加迭代次数，但相应地增加了计算量。

2.联立求解法

联立求解法一般针对微分方程用隐式积分法求解的情况。其基本过程为，先用隐式积分公式将微分方程组代数化，它和代数方程组—起形成联立非线性方程组，然后求解此非线性方程组，即可得到所要的解。显然,这种求解方法不存在交接误差。当应用隐式梯形积分公式时，联立求解法就是对式（6—33）和式(6—34)联立求解。联立求解的方法一般采用牛顿法，在求解中，为提高计算效率，应充分考虑方程的稀疏性。 6.2.3 暂态稳定分析的基本流程

分析电力系统暂态稳定的主要途径是通过对遭受大干扰后系统动态响应的计算得出系统是否稳定的结论。

事实上，在系统遭受干扰后的整个暂态过程中，描述系统动态特性的微分—代数方程组[式（6—1)和式（6—2）］实际上是非自治、不连续的。微分方程和

代数方程的组成或/和内容在暂态过程中可能发生变化,即它们是“故障或操作”的内容及其发生时刻t的函数。系统可能发生的“故障或操作”有很多，例如，发生短路故障、切除输电设备、输电线路继电保护及自动重合闸的动作、串联电容的强行补偿以及制动电阻的投入或退出等，这些情况下电力网络的结构或／和参数将发生变化,因此需要在计算过程中相应地改变代数方程。又如，切除发电机、投入强行励碰、进行快速汽门控制等，特使发电机组有关元件的结构或／和参数发生变化，因此需要改变相应的微分方程.除了“故障或操作”外,一些调节系统限幅环节的存在也导致在暂态过程中微分方程和代数方程的不连续。 由于在不同时刻发生的各种“故障或操作”将导致微分-代数方程组不连续，这就使得运行参量y（t）在“故障或操作”时刻发生突变，但根据微分方程解对初值的连续依赖性[13]可知，状态变量x（t）在整个暂态过程中总是连续变化的。因此，在进行暂态稳定分析时，可以根据“故障或操作”发生的时刻把整个暂态过程自然地划分为多个时段。在一个时段内,函数f和g的结构和形式是不随时间变化的，因而微分-代数方程组是自治的。显然，在一个时段的计算结束后(t0时刻)和下一时段的计算开始前（t0时刻)，应根据发生的“故障或操作\"修改式(6—1)和式(6—2）的形式和内容，由于xt0xt0，这样就可以根据xt0重新求解修改庸的网络方程，从而得到yt0。在得到新时段的微分—代数方程组及其初值xt0、yt0后，就可以用6。2。2节中介绍的方法求解微分-代数方程组。

通常将系统遭受大干扰的时刻定为初始时刻（即t=0s），在对微分-代数方程组用某种数值方法的求解过程中,可以根据系统的运行状态利用适当的判据判断系统的稳定性.暂态稳定分析的其本流程如图6—8所示.

在进行暂态稳定分析前，首先应利用潮流计算程序算出干扰前系统的运行状态，即由潮流计算得到各节点的电压及注入功率，然后算出系统的运行参量

，并由此计算出状态变量的初始值x(0)。见图6—8中的①、②两框。 y(0)，

框③是根据各元件所采用的数学模型形成相应的微分方程,并根据所用的求解方法形成相应的电力网络方程。应当注意的是,在暂态稳定计算中的网络模型和潮流计算中有所区别，前者应考虑发电机和负荷的影响。关于这一点在本章稍后有详细论述。

从框④开始，进入暂态过程计算。目前的大多数程序中,积分步长t取为固定不变的常数。假定暂态过程的计算已进行到t时刻,这时的x(t)和y(t)为已知量.在计算x(tt)和y(tt)时,应首先检查在t时刻系统有无故障或操作，如果有故障或操作，则需对微分或/和代数方程式进行修改，见图6—8中的⑤、⑥两框。而且当故障或操作发生在电力网络内时，系统的运行参量y(t)可能发生突变,因此必须重新求解网络方程,以得到故障或操作后的运行参量y(t0).见图6—8中⑦、⑧两框，由于状态变量不会发生突变，因此故障或操作前后的x(t)和x(t0)相同。 框⑨是微分-代数方程组一步的计算，根据x(t)和y(t)采用交替求解法或联立求解法得到x(tt)和y(tt)的值。然后在框⑩中利用适当的判据(例如,可以采用任意两台发电机转子间的相对摇摆角超过180°作为系统失稳的判据）进行系统稳定性的判断.如果系统失去稳定，则打印计算结果，并停止计算（框⒀)；否则，经框⑾将时间向前推进t，进行下一步的计算，直至到达预定的时刻tmax(框⑿)。

tmax大小与所研究问题的性质有关.当仅关心第一摇摆周期系统的稳定性时,通常取tmax1~1.5s.这时的暂态稳定计其容许采用较多的简化。例如，可以忽略调速器的作用而假定原动机的机械功率保持不变；可以把励磁调节系统的作用近似考虑为在暂态过程中保持发电机暂态电势不变.这种简化模型下的暂态稳定计算将在6.4节中详细讨论.对于大规模互联电力系统，系统失去稳定的过程发展较慢，往往需要计算到几秒甚至十几秒才可能判断出系统是否稳定。这种情况下，必须用更复杂的模型来模拟系统的暂态过程，例如计及发电机组励磁调节系统和原动机调速系统的作用，考虑直流输电系统，考虑其他控制装置的作用等，这些将在6.5节中详细讨论。

最后需要指出的是,一个商业化的暂态稳定分析程序至少应满足以下基本要求：

（1）有足够的准确度。整个暂态过程中发电机转子角度的最大相对误差应小十几个百分点。

（2)算法可靠。数值积分方法的数值稳定性和任何迭代过程的收敛性要好. （3)占用内存少。使得一定容量的计算机可以进行大系统的计算。 (4)使用灵活且容易维护。可以根据不同的需要组织相应的模型进行计算，模型修改容易.

这样，在程序的构成上需要在计算速度、精度、可靠性、内存占用、灵活性等方面之间进行综合权衡.

6。3 暂态稳定分析的网络数学模型及其求解方法

和潮流、短路计算中一样，电力网络的节点电压方程可用相量表示成

式中：I、V分别为电力网络节点注入电流和节点电压组成的列向量;Y为节点导纳矩阵。式(6—35）所描述的网络方程在形式上为线性方程组,其中的导纳矩阵y仅由电力网络的结构和参数所决定.

在以后的计算中，经常把电力网络方程写成实数形式；

式中：n为电力网络的节点数；Gij段Bij分别表示网络导纳矩阵元素Yij的实部和虚部；Ixi 、Iyi和Vxi、Vyi分别表示节点注入电流和节点电压的实部相虚部。 在电力系统中,电力网络将系统中看起来相互的所有动态元件联系在一起。在暂态过程中的任一时刻,各动态元件注入网络的电流不但由其自身的持性决定,而且整个电力网络必须满足基尔霍夫定律。其中前者由各动态元件自身的

代数方程描述,后者反映在电力网络方程中.因此，为了求解网络方程,需要列出各动态元件自身的代数方程，并对其进行处理，从而可以和网络方程联立求解。一般来说，各动态元件注入网络的电流是描述其动态行为的状态变量和相应节点电压的函数，推导出该函数的表达式，即是本节的主要工作。

此外，在暂态过程中系统发生的故障或操作会引起电力网络结构和参数的变化。特别地，当故障或操作为三相不对称时，网络方程除反映电力系统的正序网络以外还可能与负序网络及零序网络有关.这样在故障或操作情况下如何处理网络方程也成为暂态稳定计算中需要解决的问题. 6。3。1 各动态元件与网络的连接 1。发电机于网络的连接

对于5。2。2节和5。2.3节中描述的各种同步电机模型，d—q坐标系下的定子电压方程都可统一表示为

式中：Ed、Eq、Xd、Xq分别表示同步机d轴和q轴的电势和电抗,随着同步电机所用模型的不同而不同，他们的取值可通过对比式(6—37）与原始定子电压方程而得到，如表6-1所以。



把坐标变换式（5—62）用于式（6-37)，可得到x—y坐标下的定子电压方程

sincos(6-38)

cosVxEdRaVsinyEqXdXqsincosRacosIxIsiny

将式（6-38)加以整理,即为发电机节点注入电流的表达式

IxgxIbyy式中：

bxEdGxgyEqByBxVx （6-39) GyVygxbyGxBxRasinXdcosRXdXq2aRacosXdsinRXdXq2aRa(XdXq)sincosRXdXqXdsinXqcos2RXdXq2a22abx2RaXdXqRasinXqcosgy2RaXdXq （6-40） Xdcos2Xqsin2Bx2RaXdXqR(XdXq)sincosGya2RaXdXqRacosXqsin将式（6—39）得到的注入电流表达式代入网络方程式（6-35），经整理后可以看出，一台发电机接入系统相当于在网络中的相应母线上接入一个电流源

'Ixgx'IybybxEd gyEq这个电流也称为发电机的虚拟电流，并且网络导纳矩阵中的相应对角块上应加上矩阵

GxByBx Gy由此可见，发电机接入系统后，在暂态过程中的任何时刻网络方程仍为线性方程，但其中的发电机虚拟注入电流及相应的导纳矩阵式发电机本身的状态量Ed、

Eq、的函数,因此这个线性方程是时变的。

在同步发电机采用简化模型时，相应的网络方程能够得到简化，方程可保持为n阶的复型方程。除非网络发生故障或操作，否则网络方程的系数矩阵为定常矩阵。这样在整个暂态过程中，为了求解网络方程,只需在网络发生故障或操作时对方程的系数矩阵重新进行三角分解。下面讨论两种同步电机采用简化模型时的网络方程.

''当不计阻尼绕组影响时,同步电机模型对应于表6—1中Eq变化或EqC的

情况，此时式（6—39）可改写为

'XdXqRasin222'IxRaXdXqIX'X'XXqydqdcos2222'RXadXq'XdXqcos222'2'RaXdXqXqcosVx'X'sinV XdXqyqRasin222'RaXdXq'XdXq（6—41）

由此可以得到接在节点i的发电机注入电流的复数表达式

IiIiYiVi （6—42）

式中:

..''.

式(6-42)可以用图6—9所示等值电路表示。其中Yi'称为发电机的虚拟导纳，其值仅由发电机本身的参数决定，可以把它并入电力网络中；Ii'称为发电机的虚拟注入电流，它与发电机的端电压Vi有关.这就使得网络方程成为非线性方程，其求解必然要用迭代法，即首先给出电压的初值Vi，用式（6-44)算出Ii'，然后以Ii'为注入电流求解网络方程，得到Vi的修正值。继续迭代直到收敛为止。在正常时段计算时,一般迭代2～3次即可收敛，而在故障或操作的瞬间迭代次数要多一些

［26］

......。

当同步电机采用经典模型描述时，对应于不计阻尼绕组的影响并忽略凸极效应，同时认为暂态电抗Xd后的虚构电势E'在暂态过程中保持不变。这对应于表6—1

''中E'Eq的情况。相应地，根据式（6—42）～(6—44）可以得C，XqXd到

显然这时发电机的虚拟注入电流Ii'与发电机的端电压Vi无关。因此只要把虚拟导纳Yi并入电力网络中，在任何积分步,Ii'为已知量,可通过直接法求解网络方程得到节点电压.

2。负荷于网络的连接

如果负荷接入网络，负荷性质的不同将会使对负荷与网络连接的处理有所区别。

(1)当负荷用恒定阻抗模拟时,可将相应的等值导纳直接并入电力网络中的节点。

（2）当按动态持性模拟负荷时，如果只考虑综合负荷中感应电动机的机械暂态过程，负荷仍可用阻抗来模拟。但这个阻抗不是恒定的,而是随感应电动机滑差s的变化在变化。因此，在暂态稳定计算过程中的每一时刻．都必须根据当

'...时感应电动机的滑差重新计算综合负荷的等值阻抗.这就使得网络导纳矩阵的对角元素在计算过程中不断变化。所以在求解网络方程时，每个时段都要对导纳矩阵重新进行三角分解。

（3)当按动态特性模拟负荷时、如果只考虑综合负荷中感应电动机的机电暂态过程，可以按5。5.2节中所叙述的负荷动态模型将动态负荷用如图6—10所示的诺顿型等值电路描述.即把与负荷有关的阻抗RjX和KM(r1jx')并入网络，从而使负荷变为一个简单的电流源.这和前面对发电机与网络连接的处理类似。

以上三种负荷模型下的电力网络方程都是线性的。

(4）当按电压静态特性模拟负荷时，相应的节点注入电流是节点电压的非线性函数，导致电力网络方程成为非线性方程。根据式（5-148)和式（5-153）可知，负荷电压静特性的模拟有二次多项式和指数两种形式：

注意，式中的功率为负荷从电网吸收的功率。

根据节点电压、注入电流和注入功率的关系式

容易得到负荷注入网络电流和相应节点电压之间的关系式.

当用二次多项式模拟负荷的电压静特性时,负荷注入网络电流的表达式为

其中,也可以将负荷中与电压平方成正比的项以恒定导纳的形式并入网络中的节点i，这时负荷注入网络电流的表达式就只剩下式（6-48)中的最后两项.

当用指数模拟负荷时，负荷注入网络电流的表达式为

3。 FACTS元件与网络的连接

这里我们仅对SVC和TCSC与网络的连接加以叙述，其他FACTS元件与网络的连接可按照同样的方法加以处理. (1）SVC.

一般情况下,SVC经升压变压器接入所控制的高压母线上（设其节点号为i)，

.BSVC这样，该接地支路的导纳为j.根据高压母线电压Vi和注入其节点电流

1XTBSVC.Ii的关系式不难得到注入电流的实部和虚部

式中：XT为升压变压器的电抗；BSVC为SVC的等值电纳；Vxi和Vyi分别为高压母线电压的实部和虚部。 （2)TCSC。

无论TCSC串联在线路中间的什么位置，总可以在TCSC两端设置两个虚拟节点i和j这样，TCSC的作用相当于在节点i和节点j分别接入一个大小相等、方向相反的电流源，容易得到两端注入电流的表达式

式中：BTCSC为TCSC的等值电抗;Vxi、Vyi、和Vxj、Vyj分别问两端母线电压的实部和虚部。

4。 两端直流输电系统与网络的连接

用下标“d\"代表直流侧,下标“R”和下标“I'’分别代表整流侧和逆变侧（下同），借助于式(4-52）～(4—54）和式(4-57)（其中取k1），可以分别写出整流器的稳态方程

和逆变器的稳态方程

这样利用式（6-52）、式（6-53）可以把直流系统注入交流系统的功率表示为变量Id、、、VxR、VyR、Vx1、Vy1的函数。 注入整流侧交流母线的功率可表示为

注入逆变侧交流母线的功率可表示为

求得注入整流侧交流母线和注入逆变侧交流母线的电流为

将式（6—54)、式(6—55）代入式(6—56)消去PR、QR、PI、QI后可知，注入电流IxR、IyR是变量IdR、、VxR、VyR的函数，而IxI、IyI是变量IdI、、

Vx1、Vy1的函数。

6.3。2 网络操作及故障的处理[1］

电力网络中发生故障或操作时,需要通过修改网络导纳矩阵来反映相应的故障或操作。当故障或操作为对称时，例如，三相短路、元件的三相开断、串联电容的强行补偿和电气制动的投入或退出等，对应于网络中某些接地或不接地支路的参数发生变化，从而可以很容易地修改网络导纳矩阵.

电力网络中发生的短路和开断大部分是不对称的,因此需要用到对称分量法加以分析，这样除了涉及正序网络外，还涉及到负序网络和零序网络的问题。但是,在电力系统稳定分桥时，我们仅关心网络中节点电压和电流的正序分量,而对负序网络和零序网络中的电压和电流不感兴趣，它们的影响可以用在故障端口看进去的等值阻抗来模拟。

在用对称分量法分析不对称问题时，一般取A相作为基准相，这样各种短路或开断的边界条件都以A相的序分量表示。另外在发生短路或开断时，把和其他两相表现不同的相称为特殊相，即，单相接地短路时的特殊相为短路相；两相接地短路和两相短路时的特殊相为非短路相;单相开断时的特殊相为开断相;两相开断时的特殊相为非开断相。当发生短路或开断的特殊相就是A相时,根据边界条件可以将三个序网直接连成所谓的复合序网,这相当于在正序网络的故障端口处接入附加阻抗。附加阻抗的取值与故障类型有关．如表6-2和表6-3所示。这里所说的故障端口，在短路时指短路点和大地之间形成的端口，在开断时指开断形成的断口本身。这时同样可以很容易地修改正序网络导纳短阵.

[2][0]注：Z为负序网中短路点的的自阻抗，Z为零序网中短路点的自阻抗。

注:Z(2)负序网中断口处的等值阻抗，Z(0)为零序网中断口处的等值阻抗。 当发生短路或开断的特殊相不是A相时，边界条件中出现了复数算子

aej120，因此不能将三个序网直接连成复合序网，但我们可以通过在零序网络、

正序网络和负序网络各设置1:n'(0)、1:n'(1)及1；1:n'(2)的理想变压器后将它们相连。这里的理想变压器，其两侧的电压和电流具有同样的变比，从而使得变压器没有损耗。对不同的特殊相来说，理想变压器在不同的序网中有不同的变比,如表6-4所示。

在引入理想变压器以后，可以将电力系统中发生的各种不对称短路和开断依三序网的连接方式归纳为串联型故障和并联型故障。属于串联型故障的有单相接地短路、两相断线和串联电容单相击穿，这类故障的边界条件为：在理想变压器的非标准侧，各序网电压之和为零．各序网电流相等。属于并联型故障的有两相接地短路、单相断线和串联电容两相击穿，这类故障的边界条件为：在理想变压器的非标准侧，各序网电流之和为零,备序网电压相等。

在发生多重短路或／和开断，并且它们发生在非同名相时，可用类似于单重故障的处理方法修改正序网络导纳矩阵，但这时的附加阻抗要推广为综合阻抗矩阵。下面以单相接地短路和单相断线同时发生的情况为例说明综合阻抗矩阵的基本概念和处理方法。

设电力网络中k点发生单相接地短路(编号为1），i、j点之间发生单相断线(编号为2)．并且这两个故障发生在非同名相上。根据故障处三序分量的边界条

'(1)'(2)件，可得到复合序网如图6—11（a)所示。图中n1'(1)、n2、n1'(2)、n2为理想变

压器的变比，与故障的特殊相有关.为了数学上处理方便，将图6-11(a）变为图6—11(b）所示的复合序网。不难看出，两个复合序网中的变比有以下关系：

下面我们根据复合序网图推导从正序网络故障端口向负序和零序网络看进去的阻抗矩阵Zf，称此矩阵为多重故障的综合阻抗矩阵。

在图6—11（b）中．左边单相接地短路部分形成一个回路，设其回路电流为I1，右边单相断线部分形成两个回路，设其回路电流分别为I2和I2。这样负序、零序网络中故障端口电流I1、I2、I1、I2可以用这些回路电流来表示：

式中:C为与故障条件有关的关联矩阵；各向量和矩阵为:

.(2).(2).(0).(0).(1).(1).(0)

另外，根据回路电压方程可以得到各序网故障断口电压关系

式中:CT为C的转置，各向量和矩阵为

从变压器的非标准侧看进去，负序和零序的端口电压、电流间的关系可以表示成

由于负序网和零序网中包含了理想变压器,所以式(6-61）．式(6-62)中的阻抗矩阵一般是不对称的，阻抗矩阵元素的求法将在后面介绍。将式（6—61)、式(6-62)合并后得到

或简写为

利用上式中的Z及关联矩阵C就可以求出正序网电压、电流间的关系式.为此，将式(6-)和式(6—57）代入式(6-59）中得到故障部分回路电压与电流之间的关系

式中:ZL称为回路阻抗矩阵：

在本例中

式（6—65）中消去电流I2后得到

.(0)

式中:阻抗矩阵各元素Zmn（m、n可取1或2)可用下式求出：

式（6-68）可简化为

至此就得到了从正序网络故障端口向负序和零序网络看进去的阻抗矩阵Zf。

式（6-70）也可以写成综合导纳矩阵的形式

1式中:YfZf,这样得到Yf后,就可以将其中的元素追加到正序网导纳矩阵的相

应位置。

在本例中,利用关系

和式(6-71）就可得到正序网络中节点k、i、j电压和注入电流间的关系

总结上述处理，综合阻抗矩阵的计算包括以下二个过程:

（1）形成负序、零序网络的故障端口阻抗矩阵[见式(6-63）］。

(2)利用表示多重故障边界条件的关联矩阵，形成回路阻抗矩阵ZL［见式（6—66）、(6-67)]。

（3）消去闭合回路，形成综合阻抗矩阵Zf[见式(6-68）、式（6-69）］。 现将这三个过程分别叙述如下：

(1）负序、零序网络故障端口阻抗矩阵的形成。

在进行电力系统暂态稳定分析时,应首先形成各序网的导纳矩阵，并对其进行三角分解得到相应的因子表。这样就可很容易地根据故障信息得到各序网故障端口的阻抗矩阵.

对负序网，从图6-11（b)可以看出，只要在理想变压器的非标准侧向节点k注入单位电流，其他节点注入电流为零，即Ik1，Im0(m属于除节点k外的所有节点)，求解包括理想变压器的负序网络方程，得到的电压Vk和

.(2).(2).(2)VijViV(2)。 Z21.(2).(2).(2)j(2)在数值上即分别等于式（6-61)中阻抗矩阵第一列元素Z11和

但具体求解是这样的：在理想变压器的非标准侧向节点k注入单位电流等价于直接在负序网向节点k注入的Ikn1.'(2).(2)(2)的电流，利用负序网的因子表进行稀

.'(2)j疏前代和稀疏回代运算,求得Vk和VijViV的变比,得到VknVk，VijnVij。

.(2)(2)1.'(2).(2)(2)2.'(2).(2).'(2)，然后乘以相应理想变压器

同理，直接在负序阿向节点i和节点j分别注入n2(2)和n2(2)的电流，利

.'(2).'(2)j用负序网的因子表进行稀疏前代和稀疏回代运算，求得Vk和Vi然后乘以相应理想变压器的变比，得到VknVk，Vijn(2)(2)别在数值上等于式(6—61）中阻抗矩阵第二列元素Z12和Z22。

.(2)(2)1.'(2).‘(2).'(2).(2)ViV.'(2)，

(2)2Vij，它们分

对零序网，可用仿照负序网计算的同样方法得到式（6-62)中阻抗矩阵各元素。

(2)由关联矩阵形成回路阻抗矩阵.

如前所述,串联型故障的复合序网由三个序网串联而成，因此只有一个回路.并联型故障的复合序网由三个序网并联而成，因此形成了两个回路.另外，我们把两相短路也作为一种特殊的并联型故障来处理。

由式(6-57)、式（6—59)可知，关联矩阵C表示复合序网边界回路电流与负序、零序网故障端口电流之间的关系，因此关联矩阵的行数应等于IS的维数，即两倍故障重数（当发生两相短路时,在零序网设置一个空着的故障端口）。关联矩阵的列数应等于IL的维数。一个串联型故障在关联矩阵C中占一列，其元素排列为

其中非零元素的列号为该故障在所有故障中的排列顺序号。并联型故障在关联矩阵C中占两列,其元素配置为

其中第一列描述负序网络与正序网络的连接情况，第二列描述零序网络与负序网络的连接情况。其中非零元素的列号为该故障在所有故障中的排列顺序号。对于两相短路来说，由于负序网络与零序网络之间没有回路相联系，因此关联矩阵只有上述第一列元素。

根据以上原则，我们可以很容易地根据故障的类型形成任意复杂故障边界条件的关联矩阵.例如，当系统发生三重故障,其编排顺序为单相接地短路、单相断线、两相短路,这时的关联矩阵为

有了描述复杂故障边界条件的关联矩阵以后，就可以根据式（6-66)、式（6—67)求出复合序网的回路阻抗矩阵,这种矩阵乘积可经过一些简单的加减运算实现.

（3）消去并联型故障的闭合回路，形成综合阻抗矩阵。

回路阻抗矩阵的阶数等于复合序网中的回路电流数。为最终形成综合阻抗矩阵，还必须把非正序电流消去［见式（6-68)、式（6-69）］。 6.4 简单模型下的暂态稳定分析

对于地区性的电力系统来说，一般失去暂态稳定的过程发展很快，通常分析系统遭受干扰后第一摇摆周期（1～1．5s)的机电暂态过程就可以判断系统是否能够维持稳定运行。这种情况下的暂态稳定分析中，由于调速系统的惯性,使得在短时间内原动机的功率不会发生很大变化，因此可以忽略调速系统的作用而假定原动机的功率保持不变；此外由于发电机励磁绕组的时间常数较大，这样在短时间内其磁链也不会发生显著变化，而对励磁调节系统的作用，可以用发电机暂

'态电势Eq或E'保持恒定来近似模拟，即认为在第一摇摆周期内，励磁绕组中自

由电流分量的变化由励磁调节系统的作用所补偿,从而使励磁绕组的磁链f在这段时间内保持不变.相应地，阻尼绕组的影响将略去不计.

简单模型下的暂态稳定分析程序在电力系统运行与规划中获得了广泛的应用，用它可以验证电力系统接线方式和运行方式的合理性，计算输电线路的最大输送功率，确定系统故障切除的临界时间，以及研究某些提高电力系统稳定措施的效果,等等。

对发电机、负荷及电力网络采用不同的数学模型可以构成各种简化的暂态稳

定分析程序，采用何种组合要根据所研究问题的性质而定。为说明简化暂态稳定分析程序的原理和方法，下面介绍的简化暂态稳定分析程序，采用如下的数学模型和计算方法:

' 发电机：用发电机暂态电势Eq，保持恒定来模拟。

负荷：较小的负荷用恒定阻抗模拟,较大的负荷考虑综合负荷电动机转子机械暂态过。

电力网络：用导纳矩阵描述。 微分方程：用改进欧拉法求解。 网络方程:用直接法求解。

整个暂态稳定计算的基本过程仍如图6—8所示.下面给出机电暂态过程的计算的实现。 6.4。1 初值计算

在进行暂态稳定分析前，首先应根据潮流程序算出的干扰前系统运行状态确定微分方程求解所需的初值.在简单模型下的暂态稳定分析中，初值的计算包括干扰前瞬间发电机的暂态电势、转子角度、原动机的机械功率以及综合负荷电动机的滑差、等值电纳等.这些参数在系统受到干扰后的瞬间是不会发生突变的。以下各变量的下标(0)表示初值。

首先，我们介绍发电机的初值计算。

由潮流计算可得到干扰前各发电机的端电压V(0)Vx(0)jVy(0)和各发电机注入网络的功率S(0)P(0)jQ(0)，进而可以计算出发电机注入网络的电流为

.这样，根据式（5—61)就可求出虚构电势EQ(0),即

依此就可以确定发电机转子角度的初值

在系统稳态运行时，发电机转子以同步转速旋转，于是有

.

利用坐标变换公式(5—62)，可以求出发电机定子电压和电流的d、q分量

根据式(5—）可以算出暂态电势的值

另外，稳态运行时发电机的电磁功率Pe（0）等于原动机的机械功率Pm（0），即有

负荷初值的计算比较简单。

由潮流计算结果可知干扰前各负荷节点电压V(0)和负荷所吸收的功率S(0)，据此容易得到负荷的等值导纳

.

当按恒定阻抗模拟负荷时，其等值导纳在整个暂态过程中保持不变，如前所述,可以将它包括在网络导纳矩阵里.对于考虑综合负荷电动机机械暂态过程的负荷，由于在扰动瞬间电动机的滑差不能突变，因此负荷的等值导纳也不应突变，即扰动后瞬间负荷的等值导纳与正常运行情况下的等值导纳相同。 6。4.2 用直接法求解网绕方程

在这种求解法中，网络方程用实数形式描述，如式(6—36）所示。进行机电暂态过程计算前，对于用恒定阻抗模拟的负荷，应首光将其等值导纳并入电力网络，从而得到考虑了恒定阻抗负荷后的电力网络方程,它在整个暂态过程中是保持不变的。

设电动机负荷接在网络中的节点j。在暂态过程中各电动机的转差sj是随时间变化的．可以根据某时刻的sj由式(5-160）求出相应时刻电动机约实际阻抗

式中:ZMj(0)和ZM(0)别为正常运行状态下全部感应电动机的等值阻抗和典型电动机的等值阻抗。与实际阻抗相应的导纳可写为

'设发电机接在网络中的节点i。当发电机采用Eq变化的模型时，可参照表

6-1在式(6—39）中取Edi0,EqiE，XdiX，XqiXqi得到发电机节点注入电流的表达式

'qi'di

其中的元素根据式（6—40）写为

将发电机节点注入电流的表达式（6-84)代入考虑了恒定阻抗负荷后的电力网络方程．并把电动机的等值导纳[式(6-83）]并入网络，即得到新的网络方程.很显然,新的网络方程只是对原网络方程的简单修改；导纳矩阵的相应对角块发生变化，电流向量中仅发电机节点有虚拟注入电流，其余节点的注入电流为零，即：

导纳矩阵的第i个对角块变为

第j个对角块变为

发电机节点的虚拟注入电流为

这样在每个积分步得到的线性方程组可以用高斯消去法或三角分解法直接求解，从而解出此时刻网络各节点电压的实部和虚部Vx、Vy。在得到发电机节点的电压后,即可按式（6-84）算出发电机节点的注入电流Ix、Iy. 6.4。3 用改进欧拉法求解微分方程

简单模型下的暂态稳定分析中，全系统的微分方程式包括各发电机的转子运动方程式(5—76）和各典型综合负荷电动机的转子运动方程式（5-155）：

假设电力系统的机电暂态过程已经计算到t时刻，现在我们围绕微分方程的求解方法来讨论tt时刻系统运行状态的计算过程。新的时段计算前总是先判断系统在t时刻有无故障或操作发生.若无故障或操作发生，则直接以t时刻系统状态作为初值求解微分方程；否则需要计算出故障或操作后电力网络的运行参数，用它和t时刻的状态变量作为初值求解微分方程。用改进欧拉法求解微分方程的步骤如下：

（1）由t时刻各发电机的x(t）和各电动机的sj(t）按6。4.2节的方法计算出系统所有节点的电压Vx(t)、Vy(t)和发电机节点的注入电流Ixi(t)、Iyi(t)。

(2）根据式（6-17)，首先应求出状态量在t时刻的导数值，即

其中，各发电机的电磁功率Pei(t)按下式计算：

各电动机的机械转矩TmMj(t)和电磁转矩TeMj(t)可根据式(5—157）和式(5—156）求出

式中:Vxj(0）、Vyj(0）分别表示节点j在干扰前正常运行状态下电压的实部和虚部.

（3)然后求出tt时刻状态量的初值估计

](4）类似于第（1)步,由各发电机的i[(0t]t)和各电动机的[j0(tt),按6。4。2节

的方法计算出系统所有节点的电压Vx[(0t]t)、Vy[(0t]t)和发电机节点的注入电流

[0][0]Ix(tt)、Iy(tt)。

di[0](5)类似于第（2)步，应求出在tt时刻状态量初值估计的导数值

dt、

ttdi[0]dtdsi[0]、

dttt。对此，只需在式（6-90)~(6-92）中将i(t)、Pei(t)、MmMj(t)、

tt[0][0][0]MeMj(t)分别换成i[(0t])、Pei(tt)、MmMj(tt)、MeMj(tt)即可。而为了得到它们,还应[0][0][0]把Vxi(t)、Vyi(t)、Ixi(t)、Iyi(t)、sj(t)、Vxj(t)、Vyj(t)分别换成Vxi(tt)、Vyi(tt)、Ixi(tt)、0][0][0]、、、sI[yiVVj(t)(tt)xj(tt)yj(tt)。

（6）最后，求出各状态变量在tt时刻的数值，即

【例6-3】 考虑图6—12所示的9节点电力系统［5]。该系统有3台发电机、3个负荷以及9条支路。支路数据和发电机参数分别列于表6—5和表6—6，正常运行情况下的系统潮流如表6-7所示，系统频率为60Hz.

注：表中所有时间常数的单位为“s”，阻尼系数D及所有电阻、电抗均为“标幺值”。

【解】下面针对该系统进行简单模型下的暂态稳定分析。干扰是在零秒线路5-7靠近母线7处发生三相接地短路，故障在5个周波(约0。083 33s）由断开线路5—7而被消除。

' 各发电机用暂态电势Eq保持恒定来模拟，各负荷用恒定阻抗模拟，电力网

络用导纳矩阵描述，微分方程用改进欧拉法求解,网络方程用直接法求解. 根据图6-8中暂态稳定分析的基本过程和上节所介绍的方法及计算公式,暂态稳定分析可归纳如下: 1）初值计算.

根据式（6-81）计算各负荷的等值并联导纳,结果如下： 负荷（节点5)：1.260 99-j0。504 40 负荷（节点6):0。877 65-j0.292 55 负荷（节点8）: 0.968 98—j0。33914

'根据式（6—74)~(6-80）计算各发电机的暂态电势Eq、初始功角（0）及其输

入机械功率Pm(0)，结果如表6—8所示.另外各发电机角速度的初值为

'在以下计算中,不计发电机凸极效应相当于在Eq1(0)2(0)3(0)1。C模型'中令XqXd,即为发电机的经典模型。

2)故障系统与故障后系统描述。

故障期间的电力网络相当于在7号母线处并联一条阻抗为零的接地支路，这时只要将正常情况下导纳矩阵Y中的对角元素Y77改为无穷大(可以用大数模拟，例如在实际计算中取为1020。），即可得到故障期间系统的节点导纳矩阵YF。

故障后的电力网络是切除线路5—7后的情况。出于线路5-7对导纳矩阵的贡献为

式中：r0.032，x0.161，b0.153，因此故障后系统的节点导纳矩阵

YPYYl(57)。

3)微分-代数方程的数值积分。

我们仅计算从短路故障开始系统在2s内的暂态过程.这样,0~2s的暂态过程可划分为两个自治系统，即0～0。083 33s的故障系统和0.083 33～2s的故障后系统。数值积分采用0.001s的步长。表6—9列出了不计凸极效应和计及凸极效应时各台发电机的及最大相对摇摆角，后者也如图6-13所示. （t）

由图6-13可以看出,无论计及凸极效应还是不计凸极效应，系统都是暂态稳定的.在计及凸极效应的情况下，最大相对摇摆角为21151.48396 (t0.80133s）。在不计凸极效应的情况下，最大相对摇摆角为2185.65788。（t0.44633s)，而第二摆的角度2185.43378 (t1.53433s)比第一摆的角度小。

最后，对不计凸极效应和计及凸极效应情况下的临界切除时间进行了计算。得到前者所对应的临界切除时间在0。162~0。163s之间,后者所对应的临界切除时间在0。085～0.086s之间。它们对应的摇摆曲线分别如图6-14和图6-15所示。

6.4.4 经典模型下暂态稳定分析的数值积分方法［49,50］

在当今的EMS(Energy Management System）中,为了评价系统运行的安全性,需要在限定的时间内对各种预想事故下系统的暂态稳定性进行在线评价.由于预想事故较多，为了满足在线动态安全评价的时间要求，需要每次暂态稳定分析有很快的速度.显然暂态稳定分析中传统的数值积分方法由于速度的已不能直接应用，而需要开发一些特殊的快速分析方法。

动态安全评价方法是近年来电力系统稳定分析领域研究的热点。早在1983年IEEE就成立了快速暂态稳定分析工作组,一直注视和引导着这一领域的研究工作。动态安全评价方法对计算速度要求很高,特别是在在线环境之下，虽然对结果的精度要求可适当放宽，但仍要求结果具有足够的可靠性。目前，提高在线动态安全评价速度的途径无非来自于以下几个方面：首先是对描述系统的数学模型加以简化，使得暂态稳定分析问题变得相对简单;其次是提出快速的暂态稳定分析方法。下面介绍一种在经典模型下暂态稳定分析的快速数值积分方法。

1．电力系统的经典数学模型

对描述系统的数学模型作如下假设后即得到所谓的“经典模型”： (1)假定发电机输入机械功率在暂态过程中保持恒定,并忽略阻尼的作用。 （2）认为在暂态过程中发电机暂态电势E'保持不变，并假定该电势的相角与转子角度相一致。

（3)负荷用恒定阻抗模拟。

在此情况下，第i台发电机的转子运动方程为

由潮流结果可计算出各发电机的暂态电势

而由稳态运行情况可知

根据6.3。1节中式(6-42)、式(6-45）、式（6-46)对发电机接入网络的处理，当把各发电机的虚拟导纳[如式(6—45)所示]并入电力网络，并将各负荷的等值导纳[如式（6—81）所示]并入电力网络后,网络方程式（6-35)将有所变化：导纳矩阵的对角元素应加上相应节点的发电机虚拟导纳或负荷等值导纳,右端电流向量中仅发电机节点有虚拟注入电流，其表达式如式（6-46)所示，而其余节点的注入电流为零。

容易得到发电机电磁功率的表达式为

2. 电力网络方程的求解

首先对导纳矩阵Y(对称矩阵)进行三角分解：

式中：U为单位上三角矩阵，D为对角矩阵。在做完以下前代和回代后,即可得到电压向量：

由于I是稀疏向量，并且为了得到各发电机的电磁功率,只需要求出各发电机节点的电压，这样所要求的V也为稀疏向量，因而网络方程的求解可以应用稀疏向量法进行快速前代和快速回代运算.一些中等规模系统的计算表明，采用稀疏向量法求解网络方程相对于执行满前代和满回代来说，可以节省大约1／3的计算量.

应用稀疏技术求解网络方程时，相当大的计算量在于对导纳矩阵的三角分解。在电力系统动态安全评价中，对于各种预想事故，故障期间和故障后的网络导纳矩阵各不相同。如果对这些导纳矩阵部分别进行三角分解，则需要花费大量的计算时间。然而，一般情况下,故障期间和故障后的网络导纳矩阵相对于故障前而言，只有少数矩阵元素有所变化,这样就可以采用补偿法求解修正后的网络方程，避免导纳矩阵的重复三角分解，从而可以大大减少计算量。

考虑网络方程

式中:Y为故障前的网络导纳矩阵;Y是由于故障或操作引起的对Y的修正量，它可以进一步表示为

式中；y为qq阶矩阵，它包含对Y的修正信息,q通常为1或2；M为nq阶与故障或操作有关的关联矩阵。

根据逆矩阵修正引理,对式(6—102)、（6—103），有

式中:qq阶矩阵C为

而qq阶矩阵Z为

这样，根据式(6-104），并考虑式式（6-99），用中补法求解网络方程式（6—102）的步骤如下：

前置计算：

求解网络方程：

式(6—107)、式(6-108）中的前代和回代运算都采用稀疏向量法.

3．二阶保守微分方程组的数值积分方法

式(6—95)所描述的微分方程组可以写成如下向量形式：

式中：

式(6-109）所示的二阶微分方程的右端函数中不显含一阶导数项,因此被称为二阶保守方程.与对两个一阶方程的求解相比，通过对方程的直接差分可使求解效率提高一倍.考虑Stormer和Numerow数值积分公式[16]：

式（6—111）为显式二阶方法，而式(6—112）为隐式四阶方法。用式（6—111）求解微分方程组（6—109)时，由于方法阶数较低且数值稳定差，因而积分时需要较小的步长。用式（6—112)求解微分方程组（6-109）时,虽然方法有较高的阶数和较大的绝对稳定域，可以采用较大的积分步长,但每步求解一个非线性方程组，因而需要过多的运算量。然而,在对式(6-112)进行迭代求解时,如果能够提供

[0]良好的初值k，2,则可使收敛加快。因此，可以采用预估-校正法求解式(6-109）

显示法[式（6—111)］用作预估器，而隐式法[式(6—112）］用作校正器。 设P、C分别表示预估器和校正器的一次使用，E表示计算函数f()，则可

[0][0][0]以构造预估校正对PECE.即从预估式中求出k，计算代入校正ff(2k2k2)，[1][1][1]式得到k2，最后计算出fk2f(k2)。

上述方法属于多步法，其起步计算可以采用如下特殊形式的四阶龙格-库塔法[15］:

电力系统的经典模型仅在机电暂态的“第一摇摆周期”是合适的(大约干扰后的1。5s之内)，并且几乎不存在刚性问题，因此数值积分可以采用较大的步长(0。1~0。2s)。

6.5含有FACTS的复杂模型暂态稳定分析

为了详细分析在遭受各种大干扰下大规模互联电力系统的暂态稳定性，并仔细分析各种控制装置对系统稳定性的影响，进而寻求改善电力系统稳定性的措施，在进行电力系统暂态稳定分析时需要采用详细的元件数学模型。

随着直流输电技术的发展，在—些远距离输电、海底电缆输电及系统互联中广泛采用直流输电系统；而近年来发展起来的柔性交流输电系统(F1exible AC Transmission System，FACTS)装置也在电力系统中得到广泛应用，它们不但能够改善系统的稳态运行特性,也在不同程度上明显改善电力系统的稳定性，使得输电线路的传输能力大大提高。大型发电机组原动机及其调速系统、励磁系统、PSS和其他的系统控制装置等对系统的稳定性也有明显的影响。这种规模不断变大、动态装置的类型和数量不断增多的电力系统，在遭受干扰后的动态行为将更为复杂,系统经历的机电暂态过程持续时间更长,系统失去稳定以前所经历的“振荡”过程可能长达几秒甚至十几秒.

在本节中，给出了当系统中包含多种动态元件，而各元件采用较为详细的数学模型时。大规模互联电力系统暂态稳定分析的基本方法.应当注意的是，它不是一个商用程序的实现，而是一个计算基本思想的介绍。

'\"'\" 本节中考虑的动态元件的数学模型如下:同步发电机采用Eq、Eq、Ed、Ed变化，加上转子运动方程的六阶模型；水轮机及其调速系统;采用可控硅调节器

的直流励磁机励磁系统；用发电机转速偏差作为输入的PSS；两端直流输电系统；FACTS中的SVC、TCSC；恒定阻抗或具有二次电压静持性的负荷等。当系统中的元件采用其他数学模型时，可用类似的方法加以考虑.

这种动态电力大系统中．由于各动态元件时间常数的差别很大.使得它成为一个典型的刚性系统。在用数值积分方法求解这类微分方程时，如果用显式方法，由于其稳定域较小，因而需要很小的积分步长;隐式梯形法是-个二阶方法，其稳定域在整个左半平面，因而采用隐式梯形法时可以取较大的积分步长。较早的暂态稳定分析程序中大多采用显式方法，例如四阶龙格—库塔法等；由于二阶隐式梯形法良好的数值稳定性和对刚性向题的适应性，加之控制装置对一些变量的使得数值积分不宜采用太大的积分步长，因而自从20世纪70年代以来,隐式梯形法得到了广泛的应用[26］，并成为迄今为止公认的最好方法。许多商用程序几乎都采用隐式梯形法,例如美国邦捏维尔电力公司（Bonneville Power Administration,BPA)的暂态稳定分析程序，我国电力科学研究院的《电力系统分析综合程序》等。在进行暂态稳定分析时，常采用定步长的隐式梯形法，步长一般可达0.01～0。02s甚至更长,而差分方程和代数方程既可采用联立求解法[1]，也可采用交替求解法[2]。

以下将介绍一个由隐式梯形法求解微分方程,而差分方程和代数方程采用牛顿法联立求解的复杂模型下大规模互联电力系统暂态稳定分析的数值方法。 6.5.1 发电机组的初值及差分方程 1．发电机

同步发电机的数学模型包括转子运动方程和转子电磁暂态方程等微分方程，以及定子电压方程相电磁功率的表达式。根据式（5-1)～(5—4)，将这些方程重写如下。

转子运动方程：

转子电磁暂态方程：

式中：

定子电压方程:

电磁功率等于输出功率加上定子铜耗：

（1）初值计算。

可以根据潮流结果用式（6-74）～（6-78）计算出发电机的部分初值.并注意到在系统稳态运行时各阻尼绕组的电流都为零,于是可根据式（5—60）、式（5-）、式(5—65)分别求出发电机空载电势、暂态电势和次暂态电势的初值：

另外，稳态运行时发电机的电磁功率Pe(0)可以由式（6—117)直接得到：

原动机的机械功率Pm(0)，可在式（6—114）中令

d0得到： dt

（2)差分方程。

首先对发电机转子运动方程式（6—114）应用梯形积分公式，得

从式(6—124)求出tt的表达式，然后将其代入式(6—123)中后可以得到

式中：

式(6—126)的J与步长t和一些常数有关,如果取固定的步长，则在暂态工程中J为常数。至于式（6—127）中的0,它仅在差分方程（6-125)中为常数，而在不同时刻的值是不相同的.

在得到tt的值后,可依式（6—123)得到tt的值，即

然后对电磁暂态过程方程式(6—115）应用梯形积分公式，得

''在式（6—129）、式（6-130）中分别消去变量Eq、Ettdtt后可以得到

式中:

\"\"式(6—134）中的系数d1、d2、d、q1、q2、q在步长t固定时都为常数.\"\"而式（6-133)中的Eq0、Ed0都为已知的t时刻量，但它们在不同时刻的值是不相

同的。

\"\"在求得Eqtt、Edtt的值以后，可分别依式（6-129）和式（6-130）求出\"'Eqtt和Edtt的值，即

2。 励磁系统及PSS

下面以图5-16所示的采用可控硅调节器的直流励磁机励磁系统为例,根据其传递函数框图列出相应的微分和代数方程。其中，忽略Rc的作用，忽略模拟电压调节器固有等值时间常数的量TB、TC。另外,在“单位励磁电压／单位定子电压”基准值系统下，根据式(5—51）可知VfEfq.

测量滤波环节：

软负反馈环节:

综合放大环节：

励磁机：

其中励磁机的饱和系数SE，既可以按式(5—101）所示的指数形式模拟，当采用“单位励磁电压／单位定于电压”基准系统时，式(5-101）可简写为

也可以将其饱和特性逐段线性化,即用下列线性关系表示饱和特性：

根据图5—14可列出PSS的微分方程

PSS输出的

（1)初值计算.

励磁系统中各状态变量的初值,可以在其传递函数框图中令s0或可通过令相应的微分方程左端等于零而得到。由于在正常运行情况下调节系统或控制装置中的各受量一般不会超出其上、下限，因此在计算初值时可以不考虑各种的作用。下面我们给出上述励磁系统初值的计算方法,其他励磁系统的初值计算类似.

令式（6—139)左端等于零从而得到放大器输出的初值

NE1式中的饱和系数用式（6-140）求得,即SE(0)CEEfq(0)。

令式（6-136）～（6—138）中微分方程左端等于零,可得到

令式（6-142）中各微分方程左端等于零，并考虑式（6-143），即得到PSS的初值

式中,由于VIS一般取速度或功率的变化量，故初值为零。

（2）差分方程。

对式（6-136）应用梯形积分公式可得测量滤波环节的差分方程

式中:

对式（6—137）应用梯形积分公式可得软负反馈环节的差分方程

式中：

不考虑的作用时,对式（6—138）应用梯形积分公式可得综合放大环节的差分方程

式中：

将式（6-141）代入式（6-139），应用梯形积分公式可得励磁机的差分方程

式中:

对式（6—142）应用梯形积分公式，得

式中：

在式（6-160）中消去中间变量V1tt、V2tt和V3tt，得

当PSS的输入为VISs时,显然VIS(t)(t)s。将VIS(tt)(tt)s代入式（6—163）,再利用式（6—128)消去变量(tt),得

式中：

当不考虑PSS的输出时，显然有

当PSS有其他输入信号时，类似以上推导,可得到相应的表达式。

在式（6-1）、式(6-167）、式(6-147)、式（6-151）、式（6-154）、式(6—157)中消去中间变量V4tt、VStt、VMtt、VFtt、VRtt后，便可得到各

环节未起作用时励磁系统的差分方程

式中：

3.原动机及其调速系统

下面以图5-24所示的水轮机及其调速系统为例，根据其传递函数框图列出相应的微分和代数方程。

离心飞摆机构：

配压阀:

失灵区的表达式为

配压阀行程的为

伺服机构：

阀门开度的为

反馈环节：

水轮机:

式中：参数KmH为发电机额定功率与系统基准容量之比，即

一般情况下，原动机及其调速系统中的量纲均是以自身容量为基准的标幺值，引入参数KmH可以使Pm和Pe一样都为系统统一基准SB下的标幺值。 （1)初值计算。

与励磁系统的初值计算类似，原动机及其调速系统中各状态变量的初值，可以在其传递函数框图中令s0或可通过令相应的微分方程左端等于零而得到。不考虑测量失灵区以及各种的作用。令式（6—177）、式（6—176)、式(6-174）左端等于零，并利用式(6-171）、式（6-172）、式（6—173）、式（6—175）中的线性关系和式(6—77)，可得到各状态变量的初值

(2）差分方程

根据式（6—171）可直接写出离心飞摆机构在tt时刻的方程

不考虑测量失灵区，根据式（6—172）可得到

不考虑配压阀行程，根据式(6-173），显然有

对式（6-174）应用梯形积分公式可得如下差分方程：

式中：

不考虑阀门开度的，根据式（6-175），显然有

对式（6—176)应用梯形积分公式可得反馈环节的差分方程

式中：

对式（6-177）应用梯形积分公式可得水轮机的差分方程

式中：

在式(6—180)~(6—183)、式（6—186）、式（6-187)、式(6—190）中消去中间变量tt、(tt)、tt、tt、tttt并考虑式（6—128）．消去变量tt便可得到tt时刻各环节未起作用时水轮机及其调速系统的差分

方程

式中：

最后,将Pe(tt)的表达式(6—117)和Pm(tt)的差分方程式（6—193）代入式（6—l 25)，将Efq(tt)的差分方程式(6—168）代入式(6—131），它们和式(6-132）一起组成tt时刻发电机组的差分方程式。将其中的d-q坐标系下的定子电流转换到x-y坐标系下，为简单起见，并省略其中关于时间的下标tt，得

式（6-196）中包含三个方程式,第一个方程式反映了发电机的机械运动,后两个方程式反映了发电机转子绕组的电磁暂态过程.根据式(6—39)，发电机注入电

\"\"流Ix、Iy是Vx、Vy、、Eq、Ed的函数［具体可参见式（6-258）]，因此可消\"\"去Ix、Iy，则这三个方程实际包含三个状态变量、Eq、Ed和两个运行参量Vx、

Vy。

6.5。2 FACTS及直流输电系统的初值及差分方程 1。 SVC

这里我们仅给出由固定电容器（Fixed Capacitor，FC)和晶闸管控制的电抗器（Thyristor-Controlled Reactor，TCR)并联组成的SVC。为简单起见，我们以比例调节器型SVC为例,其传递函数框图如图6—16所示。

SVC一般经升压变压器接入高压系统，通过对晶闸管触发角的控制来改变TCR的等效电纳，从而改变SVC的等值电纳BSVC,使高压母线的电压V达到指定值Vref，其数学模型可按框图6-16直接写出：

SVC输出的为

式中：BCC为固定电容器的电纳;BL1/L为电抗器的电纳；输出BSVC为SVC的等值电纳。SVC输出环节的上限对应于晶闸管完全关断，而下限对应于晶闸管完全导通。BSVC在范围之内则是晶闹管处于部分导通的情况。 （1）初值计算。

虽然SVC接在升压变压器的低压母线上，但它仍可看成是接在高压母线上的无功电源，用它来功率控制高压母线的电压。因此,在潮流计算时，可以将高压母线的节点处理成PV型节点（P0，VVSP）。由潮流计算结果可知高压母线的电压V(0)VSP0以及SVC注入高压母线的功率S0jQ0。设升压变压器的电抗为XT,则SVC注入高压母线的功率可表示为

. 令式（6-197）中的两式左端等于零,再考虑式（6—198）、和式(6—199)，即可得

到SVC的初值

(2）差分方程。

对式（6-197）的第一式应用梯形积分公式，得

式中：

对式（6-197）的第二式应用梯形积分公式,并利用式（6-201)消去BS1tt得

式中：

当不考虑SVC的输出时,显然有BSttBS2tt,因此得到

2。 TCSC

晶闸管控制的串联补偿器（Thyristor—Controlled Series Compensator，TCSC）串联在输电线路上，通过控制其等值电纳BTCSC改变线路的等值电抗。这里我们仅列出由FC和TCR并联实现的TCSC的数学模型（类似于SVC）：

式中：输入信号PT为TCSC安装处流过的有功功率。 TCSC输出的为

maxmin式中：输出BTCSC为TCSC的等值电纳，BTCSC和BTCSC的取值与L与C的大小有

关，具体数值可按式(4-153）～（4—155)计算。 （1）初值计算。

由潮流计算可得到BTCSC0，PT0PSP,类似于SVC初值的计算方法可得到

(2）差分方程

如果TCSC的量测量PT从节点i流入节点j的有功功率，容易得到PT的表达式

对式（6-208）的第一式应用梯形积分公式，得

式中：

对式(6—208)的第二式应用梯形积分公式，并利用式（6—212)和式（6-211)消去BT1tt和PTtt,可得到

式中：

当不考虑TCSC的输出时，显然有BTCSttBT2tt，因此得到

3. 两端直流输电系统

在稳定分析中，交流系统的网络方程用正序分量表示，这隐含了对直流系统模型的基本，特别是不能正确预测换相失败。换相失败可能来自于逆变器附近严重的三相故障，逆变器交流侧的不对称故障，或在动态过压情况下换流变压器的饱和。

早期的直流系统模型一般考虑线路和换流器控制的动态特性。最近几年，有使用简单模型的趋势.常用的直流系统模型方以下两种，即简单模型和准稳态模型。

（1）简单模型.

遥远的直流系统对稳定分析的结果无明显的影响,故可采用很简单的模型：用在换流器的交流母线注入恒定的有功和无功功率来模拟直流系统.

更实际的模型是所谓的稳态模型（Steady—State Model）。这时根据式(4-2）,直流线路用电阻电路的代数方程表示：

式中：Rdc表示直流线路的电阻。

注意到IdRIdIId，并利用式(6-52)、式(6—53）消去式(6—219)中的VdR、

VdI得

式中：

认为极性控制行为是即时的，即对各种控制功能的模拟，只表示它们的作用，而不体现它们硬件的实际特性.这个模型以代数方程的形式出现，交、直流系统的交接与潮流计算中的交接相同.

(2)准稳态模型（Quasi Steady-State Model）。

如果直流输电系统两端中任意一端的交流短路水平较低，则直流系统中元件的动态性能将影响交流系统的稳定性,从而要求暂态稳定分析中使用较为精确的直流输电模型。

在准稳态模型中，换流器特性仍用直流分量的平均值与基频交流分量的有效值之间的关系方程表示。这时，直流线路的模拟可根据不同的精度要求采用不同的模型.最简单的直流线路模型和稳态模型中一样，如式(6—220)所示.进一步更详细的直流线路模型可用R—L电路模拟：

式中：R如式（6—221)所示，另外

其中，Ldc、LR、LI分别为宜流线路、两端平波电抗器的电感。

对于控制系统，以整流器采用定电流控制方式、逆变器采用定电压控制方式为例,由图4-18所示的传递函数框图可以列出相应的微分方程

滞后触发角的为

超前触发角的为

（1)初值计算。

SP当整流器采用定电流控制方式、逆变器采用定电压控制方式时,有Id0IdSP和VdI0VdI,由潮流计算结果可得VR0、VI0。根据式（6—224)~(6-227)，并考

虑式(6—219)或式（6—222)和式（6-52)、式（6—53)可以得到

(2）差分方程。

对式（6-224）的第一式应用梯形公式可得

式中：

对式（6-224)的第二式应用梯形公式，并利用式（6-229）消去x1tt，得

式中:

不考虑对触发角的时，显然有

对式（6—226)的第一式应用梯形公式可得

式中：

对式（6-226）的第二式应用梯形公式，利用式（6-236）消去x4tt，并考虑式（6-53)的第一式从而消去VdItt，整理后得到

式中：

不考虑对触发角的时,显然有

在准稳态模型下，不考虑和考虑直流线路的暂态过程将得到不同的差分方程。

当不考虑直流线路的暂态过程时，直流线路采用式（6—220)所示模型,可以

将Id表达成变量、、VxR、VyR、VxI、VyI的函数：

在式（6-232）、式（6-235）、式(6—243）中消去tt和Idtt，即得到不考虑对的时整流器采用恒定电流控制的差分方程



式中：

同样,在式（6-239)、式(6-242）、式（6—243）中消去tt和Idtt，即得到不考虑对的时逆变器采用恒定电压控制的差分方程



式中:

当考虑直流线路暂态过程时，直流线路采用式（6-222）所示模型，对其应用梯形公式，得

式中：

在式（6-232）、式（6-235）、式（6—248）中消去tt和Idtt,即得到不考虑对的时整流器采用恒定电流控制的差分方程



式中：

同样,在式（6—239）、式(6—242）、式（6-248）中消去tt和Idtt,即得到不考虑对的时逆变器采用恒定电压控制的差分方程



式中:

6。5。3 电力网络方程的形成

电力网络方程的实数形式如式（6—36）所示.在暂态稳定计算时，我们把电力网络中的节点分为三种类型，即有动态元件并联的节点（包括发电机节点、SVC节点、负荷节点等）、有动态元件串联的节点（包括两端直流系统交流母线的节点、TCSC的两端节点等）及联络节点或故障节点。把6。3。1节中得到的各动态元件注入电力网络的电流表达式代入网络方程，并按6。3.2节完成对故障或操作的处理，即可得到用于求解的网络方程。

1. 有动态元件并联的节点

如果动态元件并联在节点i，则节点i网络方程为

节点i注入电流Ixi、Iyi的表达式取决于所接入的动态元件。 (1)接入发电机。

'\"'\"注意到发电机采用Eq、、、变化的模型，因此对照表6-1对式（6-40)EqEdEd中元素取相应的值，这时发电机节点的注入电流表达式(6—39）重写为

（2）接入负荷.

如6。3。1节所述，如果是恒定阻抗负荷,可直接并入电力网络。当按电压的二次多项式模拟负荷静特性时,可以直接将式(6-48)所示的电流注入网络，也可以首先将非线性负荷的恒定阻抗部分并入电力网络，负荷的其余部分取式(6-48)的后两项电流注入网络；当按电压的指数形式模拟负荷静特性时，可以直接将式（6-49)所示的电流注入网络.

(3）接入SVC.

SVC的注入电流表达式如式(6-50)所示。 2. 有动态元件串联的节点

如果动态元件串联在节点i和节点j之间,则节点i、j的网络方程为

节点i的注入电流Ixi、Iyi及节点j的注入电流Ixj、Iyj，随串联的动态元件的不同有不同的表达式.

（1）串联TCSC。

这时Ixi、Iyi及Ixj、Iyj的表达式如式(6—51)所示。

(2)串联直流系统.

在采用交直流系统联立求解时，直流系统注入两端交流系统母线电流Ixi、Iyi及Ixj、Iyj的表达式如式（6-56）所示。其中的功率表达式（6-54）、式(6—55）中，直流电流Id可用式（6—243）或式(6-248）代替,从而使得注入网络的电流仅为变量、、VxR、VyR、VxI、VyI的函数.

3. 联络节点或故障节点

联络节点是指注入电流恒为零的节点.如前所述,当采用综合阻抗矩阵的概念处理系统故障时，任何形式的故障都可以用修改正序网络导纳短阵的方法模拟。这样，在扩展了的正序网络的故障点上就不会再有其他序网来的注入电流,因此我们也可以把故障节点作为联络节点处理。故障节点或联络节点(节点号为f）的网络方程为

6。5.4 差分方程与网络方程的联立求解

上面已列出了tt时刻电力系统的所有方程，包括电力网络方程和各动态元件的差分方程。其中,待求的变量包括:电力系统的运行参量，即电力网络中所

\"有节点的电压Vx、Vy；所有动态元件的状态变量，例如对每台发电机为、Eq、\"，对每个SVC为BSVC，对每个TCSC为BTCSC，对每个两端直流输电系统为、Ed.设电力网络共有n个节点，其中装有nG台发电机、nS个SVC、nT个TCSC、nD个两端直流输电系统，显然总的待求变量数为2n3nGnsnT2nD，而相应的方程数与变量数恰好相同,因此可以进行求解。

总体来讲，网络方程式(6-257)、式(6-259）、式(6—260），发电机差分方程式(6-196），SVC的方程式（6—207)，TCSC的方程式（6—218），两端直流系统的差分方程式（6-244）、式(6—246）或式（6-251）、式（6—254)一起构成了一个非线性方程组。网络方程式中的注入电流、各动态元件的差分方程式是随时间变

化的，而网络方程中的其余部分除了在干扰(对“故障或操作”的泛称）发生时刻有变化外，在两次干扰之间具有同样的形式.当干扰发生而需要求出干扰后的运行状态时，仅需要求解网络方程，此时方程中注入电流表达式中动态元件的状

\"\"态量、Eq、Ed、BSVC、BTCSC、、应取干扰前瞬间的值.这样,对由网络方

程式和各动态元件的差分方程式组成非线性方程组不断进行求解，就递推出电力系统在不同时刻的运行状态.

值得注意的是，以上在推导tt时刻各动态元件的差分方程时,没有考虑各控制器中各种环节的作用。考虑各种环节作用后各动态元件的差分方程将在6.5。6节中讨论。

上述非线性方程组一般采用牛顿法求解。由于牛顿法本身已为大家所熟知，故这里仅简要介绍求解的基本步骤：

\"\" （1)给出tt时刻各发电机状态变量、Eq、Ed的初值，各SVC状态变

量BSVC的初值，各TCSC状态变量BTCSC的初值，各两端直流系统状态变量、的初值，以及电力网络各节点电压Vx、Vy的初值.这些初值可直接取为t时刻的值，也可以由前一步或前几步的值通过外推得到.

(2)对于由各发电机组的差分方程、各SVC的差分方程、各TCSC的差分方程、各两端直流输电系统的差分方程以及网络方程织成的非线性方程组，按步骤(1)中给出的初值，先求出雅可比矩阵和残差量的值,然后解线性方程组得到各变量的修正量,进而修正各变量。

(3）判断是否收敛。若收敛，结束；不收敛时，返回步骤（2)继续迭代，直到收敛为止。

(4）在得到tt时到网络各节点电压、各参与联立求解的状态变量的值后，还要按6。5.1节和6.5。2节中导出的其他差分方程或代数方程计算该时刻各动态元件的其他变量的值，以资下个时段求解使用.注意，在计算其他变量的过程中应考虑各环节的。 6。5.5 交、直流系统的交替求解

由于直流系统准稳态模型中包含的动态比与交流系统模型有关的动态快很多,因此为了详细分析直流系统的行为，在用数值积分法求解时，需要采用比交

流系统中所用的积分步长t小得多的积分步长(例如BPA的暂态稳定程序中取mt/8）求解直流系统方程，两者在时间上的关系如图6-17所示。对于交流系统从t到tt的一个步长内,交流系统和直流系统暂态过程计算的配合可以来用交替迭代法,它在原理上类似于交直流潮流计算所采用的交谷迭代求解法.

设已知t时刻和tt时刻换流器交流母线的电压为VRt、VIt、VRtt、

VItt，在计算时可以认为它们在t内线性变化，因此可得到k时刻的电压

在已知k时刻直流系统的状态，用隐式梯形法计算k时刻的状态时，只要在式（6-229）~（6-256)中把t换成k，把t换成，就可得到k时刻两端直流系统的差分方程。

注意到在计算的状态时，、、VR、VI己知，而VR、

VI可用式(6-261）求出，因此可通过求解方程式(6-244）、式（6—246）或

式(6-251)、式(6-254）得到、。此方程组为非线性方程组，因此需要用迭代法求解.这样,经过m步后可得到tt、tt。 交、直流系统的交替迭代解法可归结为：

\"\"（1)给出tt时刻各发电机状态变量、Eq、Ed,各SVC状态变量BSVC，

各TCSC状态变量BTCSC，以及电力网络各节点电压Vx、Vy的初值。

（2）按以上方法求解直流系统方程，得到tt、tt，据此用式(6-243)

或式(6-248）算出Idtt,进而用式（6—52)、式(6-53)算出各注入功率，最后用式（6-56）求出直流系统注入交流系统的电流IxRtt、IyRtt、IxItt、IyItt.

\"\" (3)用6.5。4节的方法求解交流系统，得到各发电机状态变量、Eq、Ed，

各SVC状态变量BSVC，各TCSC状态变量BTCSC及电力网络各节点电压Vx、Vy的修正值。

（4)判断是否收敛。若收敛，结束；不收敛时,返回步骤(2)继续迭代，直到收敛为止.

(5）在得到tt时刻网络各节点电压、各参与联立求解的状态变量的值后，还要按6.5.1节和6。5。2节中导出的其他差分方程或代数方程计算该时刻各动态元件的其他变量的值，以资下个时段求解使用。注意,在计算其他变量的过程中应考虑各环节的。

6。5.6 数值求解过程中一些特殊问题的处理

在电力系统机电暂态过程的数值求解中，一些特殊问题需要妥善处理。例如，对控制器中的各种如何考虑;当采用定步长积分时，积分时刻与网络故障或操作时刻不重合时如何处理等等。本小节主要讨论这两类问题的处理方法。 1．控制器中各种的处理

控制器中的一般分为终端和非终端,在机电暂态过程计算中对它们的处理也有所不同. (1）励磁调节器。

首先以PSS为例,说明对终端的处理方法。在推导tt时刻PSS的差分方程时,根据式(6—143)可知，需要知道V4tt。然而V4tt 是待求量，由于存在环节，因此无法得到V4tt和VStt之间的关系式。但是,在实际计算中，我们可以根据微分方程式（6—142）,利用欧拉法求出V4tt的估计值：

然后根据式（6—143）判断V4[0t]t是否越限，如果越限，则VStt取相应的固定

值：

另外，式（6—138）所描述的励磁系统的综合放大环节是一个非终端环节。以此为例说明对非终端的处理方法.

当VRminVRtVRmax时，首先计算VRtt的估计值

然后利用式（6—138）判断VRtt是否达到的边界上，如果达到边界则VRtt取相应的固定值:

[C]当判断出电压放大器的作用时,将VR，便可直接得tt代入式（6-157）

到取固定值的励磁系统输出

由上式代替式（6-168），显然发电机组差分方程式（6—196）中的第二式将变为

当判断出电压放大器的不起作用，而按式（6-263）判断出PSS的输出起作用时，励磁系统的差分方程将由式（6-168)变为

]注意，式中的Efq0按式(6—170)计算时,需要将其中的VS0换为固定值VS[C这样，tt。

发电机组差分方程式（6—l 96)中的第二式将变为

（2）调速器.

对于发电机的调速系统，在暂态稳定分析中一般不考虑测量失灵区的影响。

配压阀行程和阀门开度都属于终端．因此可用前述类似方法处理。

首先根据微分方程式（6-174),利用欧拉法得到tt的估计值

然后根据式（6—175）判断tt是否越限，如果越限，则tt取相应的固定值：

[0]

]当判断出阀门开度起作用时，将[tCt代入式（6-190）便可直接得到[C]取固定值的Pmtt

相应地,发电机组差分方程式（6—196）中的第一式将变为

当判断出阀门不起作用时，就要判断配压阀行程的是否起作用。这时,先根据微分方程式（6—114),利用欧拉法得到tt的估计值

然后利用式（6-180)计算tt的估计值，并利用式（6-187）计算tt的估计值，从而按式（6—181）得到tt的估计值:

这样便可利用式（6-173）判断tt是否越限,如果越限，则tt取相应的固定值：

[0]

]如果配压阀行程的起作用，将[tC，式（5-165）、t代入式（5—162）

[C]式（5-169)可直接推得取固定值的Pmtt：

相应地，发电机组差分方程式（6-196）中的第一式将变为

(3）SVC和TCSC

与上面对终端的处理方法相同。对SVC，首先根据微分方程式（6—197）,利用欧拉法得到BS2tt的估计值

[0]然后利用式(6—198)判断BS 2tt是否越限,如果越限,则BSVCtt取相应的固定值：

同样，对TCSC，首先根据微分方程式（6-208），利用欧拉法得到BT2tt的估计值

[0]然后利用式（6-209）判断BT2tt是否越限，如果越限，则BTCSCtt取相应的固

定值：

（4)两端直流系统。

首先根据微分方程式(6-224)，利用欧拉法得到tt估计值



然后利用式(6-225)判断tt，是否越限,如果越限，则tt，取相应的固定值

[0]

进而根据微分方程式（6—226）,利用欧拉法得到tt的估计值



然后利用式（6—227）判断tt是否越限，如果越限,则tt取相应的固定值:

[0]

2。 积分时刻与网络故障或操作时刻不重合时的处理

设网络故障或操作将在t2时刻发生，而由固定步长t的隐式梯形法把机电暂态过程已计算到t1时刻,并且ttt2ttt。如果按固定步长继续积分下去,将会跨越网络故障或操作。

处理这种问题有多种方法，一种很直观的处理方法是采用新的积分步长进行机电暂态过程计算，即基于t1时刻的系统状态求解t2时刻的系统状态时，采用的积分步长为t2—t1。此后可重新用固定步长积分，直至同类情况发生为止。这种处理方法原理简单，但要注意在差分方程中一些与t有关的常数需重新计算。 另外．也有从t1时刻开始改用其他显式积分方法求解微分方程到t2时刻的处理方法.

6.6 暂态稳定分析的直接法

直接法，顾名思义就是不需要借助于各状态变量的时间响应来判断系统的稳定性，通过对特定函数的数值计算结果直接判断系统稳定性的方法。因而，一般来说直接法可避免对描述系统动态特性的微分方程进行数值积分的繁重工作。为此，本节首先介绍电力系统暂态稳定问题的数学基础,然后论述暂态稳定分析的直接法。

6。6.1 暂态稳定问题直接法的数学基础[9～12］

稳定性理论是研究动态系统中的过程(包括平衡点）相对于干扰是否具有自我保持能力的理论。定性地说，如果一个系统在靠近其期望的运行点的某处开始运动,且该系统以后永远保持在此点附近运动，那么所描述的系统就是稳定的.研究非线性系统稳定性的最有用和一般的方法是李雅普诺夫稳定性理论。

(1）非线性系统和平衡点。

考虑用微分方程描述的一般非线性自治系统

式中，x[x1,x2,,xn]T；fx[f1x,f2x,,fnx]T状。状态向量的一个特定值也叫做-个点。方程（6—288)的一个解xt通常对应于状态空间内时间t从0变到无穷大时的一条曲线，这条曲线一般也叫状态轨线或系统轨线.

定义1 如果系统在t0时刻的状态是xe，并且在无任何输入或干扰情况下，对一切tt0，有xtxe，那么称xe为动态系统的一个平衡点(或平衡状态)。 很明显，平衡点xe是式（6-288)的一个未受干扰的解，这意味对

tt0,fxe0。显然非线性动态系统可能存在多个（甚至无穷多个）平衡点。

如果我们感兴趣的平衡点为xe，那么通过引入一个新的向量

并把xxxe代入式（6—288）,可得到关于状态变量x的方程

显然,式(6—288）的解xxe对应于式（6—290)的平凡解x0(原点）。因此，不失一般性，要研究式（6-288)在平衡点xe的稳定性，今后只研究式（6-290）在原点的稳定性就够了. （2）稳定性的概念。

下面首先给出系统稳定和不稳定的基本概念。

定义2 原点是稳定平衡点。如果对于任意给定的值0，存在数

,t00，使得当xt0时，对一切tt0，系统的运动xt满足

xt。

这个稳定性的定义也叫李雅普诺夫意义下的稳定性.具有这样的稳定性企本

质上意味着，若系统在足够靠近原点处开始运动，则该系统轨线就可以保持在任意地接近原点的一个邻域内。更确切地说,假如我们想让状态轨线保持在原点的邻域内，只要限定初始干扰的范数小于.注意，这时必然有。

定义3 原点是渐近稳定平衡点。如果（a）它是稳定的、并且如果（b）存在数‘t00,使得当xt0't0时，系统的运动满足limxt0

t渐近稳定性意味着，若系统在足够靠近原点处开始运动，则该系统轨线最终收敛于原点.'t0叫做平衡点的吸引域。

以上的定义是为了表征系统的局部特性而做出的,即公平衡点附近起动后状态将怎样演变.当初始状态与平衡点有一定的距离时，局部性质几乎不能提供任何有关系统运动的信息.

在电力系统遭受到一些大的干扰后,将往往伴随着一系列干扰的发生.例如,某条输电线路上发生短路故障，保护装置的动作将使得输电线路在很短的时间后退出运行；当有自动重合装置时，可能使得输电线路相继退出、投入、再退出运行情况的发生。此外，为了使得电力系统不失去稳定或提高系统的稳定性，其间还可能伴随着切除发电机、切除负荷、投入强行励磁、快关汽门等控制措施。因此,判断系统是否暂态稳定，是判断最后一个干扰发生后，系统最终能否过渡到一个可以接受的稳定平衡点。这里，我们之所以强调稳定平衡点是可以接受的原因是，事实上,实际电力系统运行要求系统的电压、频率等电气量在指定的范围内,如果它们的值超出指定范围,一些自动装置便会有相应的动作,系统又要经历新的暂态过程。

实际上,最后一个干扰发生后的系统（称为干扰后系统）是一个自治系统,如果系统在经历了一系列干扰并在最后一个干扰发生前瞬间的状态为xtk，那么在干扰发生后，系统将以xtk为初值进入在干扰后系统中的自由运动。假设干扰后系统存在稳定平衡点xS,而且系统的自由运动最终能收敛于xS，那么电力系统是暂态稳定的，否则,是暂态不稳定的。因此,电力系统的暂态稳定性本质上属于渐近稳定性的范畴。按照上述关于系统稳定性的定义2，虽然从数学意义上有限范围内的振荡是稳定的，但是我们仍然把它排除在电力系统稳定的行列之外。这显然是符合实际的，因为一个连续振荡的系统对于供电者和用电者来说都是不

合要求的.

另外,电力系统的暂态稳定性分析也可以看成对干扰后系统稳定平衡点xS吸引域的求取，当xtk在xS的吸引域之内时，系统是暂态稳定的，否则是暂态不稳定的.这里值得注意的是,一个暂态稳定的系统，首先是干扰后系统应存在稳定平衡点(否则不可能稳定)，并且受扰系统能够逐渐过渡到该平衡点。 6。6。2 暂态稳定分析的直接法[35,37~40］ 1。 概述

设正常运行的系统在tf时刻发生故障，tcl时刻故障被切除。将tf时刻以前的系统称为故障前系统（Prefault System)，时间间隔tf,tcl内的系统称为故障系统（Fault-on System），tcl时刻以后的系统称为故障后系统（Postfault System）。假设故障后系统存在一个满足运行约束的渐近稳定平衡点(Stable Equilibrium Point，SEP) xS。暂态稳定分析的直接法就是直接判断开始于故障后系统的初始状态xcl是否落入xS的吸引域（稳定域）之内。如果xcl落入xS吸引域之内，则随着时间的推移，系统状态将会趋向于xS，这时系统是渐近稳定的.反之，如果xcl在xS的吸引域之外,则系统是不稳定的。在用直接法分析系统的暂态稳定性时,常作如下假设:

Pre (1）故障前系统的SEPxS和故障后系统的SEPxS充分接近。 Pre （2)故障前系统的SEPxS在故障后系统SEPxS的吸引域内.

 直接法是基于所建立的能量函数Vx判定系统暂态稳定性的方法。通过比较故障切除时刻的系统能量与故障后系统的临界能量得出故障后系统是否稳定的结论。系统开始运行于xSpre,当故障发生后,系统的平衡状态将被打破，导致各同步发电机加速或减速。在故障期间,系统离开xSpre并获得动能和势能，故障切除时刻系统获得的总能量称为暂态能，可表示为Vxcl。当故障切除后，系统将带着这个能量进入在故障后系统中的自由运动,系统能否稳定取决于故障后系统吸收

暂态能量的能力，如果故障后系统能够将Vxcl中的动能完全转化为势能，我们就说系统是暂态稳定的,否则是暂态不稳定的。故障后系统吸收的能量有一个最大值，称为临界能量（Critical Energy)Vcr，它与系统的不稳定平衡点（Unstable Equilibrium Point，UEP）xu有关，如果能求出xcl和xu，并进而得到Vxcl和Vcr，那么就可通过比较Vxcl和Vcr来判断系统的稳定性.差值VclVxcr是系统相对稳定程度的很好度量，称为暂态能量裕度(Transient Energy Margin，TEM).这时的直接法称为暂态能量函数（Transient Energy Function，TEF）法。

用直接法分析多机电力系统的暂态稳定性始于20世纪60年代中期。借助于Moore和Anderson针对一类非线性系统提出的构造李雅普诺夫函数的一般方法，当时曾提出了一些用于暂态稳定分析的能量函数。但同时发现，只有在忽略各发电机内电势节点之间的转移电导时，这些能量函数才是李雅普诺夫函数.迄今为止，针对多机电力系统的暂态稳定分析，当考虑转移电导时，还没有找到严格的李雅普诺夫函数。由于转移电导不容忽略，因此后来所有的直接法都不是严格的李雅普诺夫意义下的直接法。

我们知道，李雅普诺夫第二定理只是判断系统稳定性的充分条件，即如果判断出系统是稳定的，则系统肯定稳定，而当判断出系统不稳定时，不能说系统肯定不稳定,因此可以肯定这种方法是保守的。而用考虑转移电导时的能量函数来判断系统的稳定性,其结果是保守还是冒进则不得而知。

关于临界能量Vcr的求取,曾提出过许多方法。早期大多使用“最近不稳定平衡点法”,其结果更趋保守。后来曾提出与故障轨线有关的“控制不稳定平衡点法”[41］，但其求解遇到数值困难.文献［43]提出的“BCU法”是目前求取临界能量的较好方法。

用于单机无穷大（OMIB)系统时，TEF法完全等价于著名的“等面积定则”。这就促成了“扩展等面积定则”法（Extended Equal Area Criterion，EEAC）的问世［45，46］。EEAC法的基本思想是，将系统中的多台发电机按照一定的规则划分为各自同步的两群,这样就可以将多机系统等值成两机系统,进而变换成一个OMIB,然后用“等面积定则”判断系统的暂态稳定性。当暂态过程中系统表现为明显的两机群摇摆模式时，这种方法是相当有效的。在用EEAC法分析多摇摆

模式问题时，发现分群困难并且群内同调不好，从而影响EEAC法分析结果的可靠性和精度，于是出现了动态EEAC，即DEEAC［48]。DEEAC法通过大步长数值积分得到多机系统的运动轨迹，在各离散点上按实际的运动轨迹动态地修正EEAC的等值两机系统的功率曲线参数。另外,在DEEAC的工程实现中还配有辨识临界机群的智能化方法.所有这些改进都使得EEAC法不断地向良好的方向发展。虽然EEAC法迄今未得到国际学术界的普遍认可，但在国外及国内几个实际系统的多种运行方式的各种故障的穷尽式考核中，对第一摆稳定性可靠的分析结果是有目共睹的.要详细了解EEAC的原理及实现，可参阅有关专著[36］.

用直接法分析电力系统的暂态稳定性大多在经典模型下进行,电力系统的经典模型在6.4.4节中已有论述。下面我们首先给出在经典模型下以系统惯性中心作为参考时各发电机的转子运动方程和相应的其他方程，然后导出多机系统的暂态能量函数并给出了一些临界能量的求取方法，最后讨论了直接法的应用和局限性。

2．多机系统的经典数学模型

'各发电机的转子运动方程如式（6-95）所示。各发电机可以表示为RajXd.后的电压源E'E'，将其接入电力网络，将各负荷的等值导纳并入网络，形成电力网络方程。在网络方程中消去除发电机内节点外的所有其他节点，得

式中：YR是收缩后的网络导纳矩阵；EG是由各发电机暂态电势E'组成的向量；

.IG是由各发电机注入电流组成的向量.

如果系统中有m台发电机，则各发电机的电磁功率可表示为

式中:

注意,Cij、Dij在故障期间和故障后为不同的常数。

在TEF法中，常取系统的惯性中心（Center of Inertia，COI)作为参考。COI的定义如下:

各发电机相对于惯性中心的运动可表示为

显然，COI和COI在暂态过程中将随时间发生变化。根据式（6—294)、式（6-295)、式(6-95)容易推得

式中：

由式（6—295）、式（6-95）、式（6—296)可导出以COI作为参考时各发电机的转子运动方程

3。 多机系统的暂态能量函数

2，,m和变量 1，2，,m组成的向量表示为(,)。令式将变量1，（6—298）左端为零,即得到求解故障后系统平衡点的非线性方程组

~~~~~~~~

式（6—299）可能有多个解,其中一个解为SEP，其他解为不稳定平衡点UEP。由于平衡点处的0，因此系统的SEP和UEP可分别表示为 （,0)和（,0)。

根据式(6—298），在第二式的左右两端分别乘以sidt和di，再对所有发电机求和得

~~~S~u

上式中右端的第二和第三项可另行表示，第四项为零:

利用式（6-301），积分式(6-300）,得

式中:C为积分常数.上式说明沿着同一条轨迹，系统的总能量保持不变。因此，

（，）可以定义故障后系统在运动轨迹上的任意一点处相对于SEP（,0)的暂

~~~S态能量函数为

式中：第一项是动能;第二项是位能；第三项是磁能，它是网络中所有支路存储的能量；第四项是耗散能量，它是网络中所有支路消耗的能量.后三项我们统一称为势能。显然，动能Vkc仅是各发电机转速的函数，而势能Vpc仅是各发电机功角的函数。

值得注意的是，式(6—303）中最后一项仅在已知系统的运动轨线后才可计算，而直接法又恰好是希望避免对运动轨线的计算。因此，在势能计算时不得不做近似处理,一种简单的处理方法是取积分路径从到的直线，即

根据以上假设，容易得到

~S~

由式（6-305），并利用式（6—304)，得

这样就得到如下近似表达式：

将上式代入式（6-303），即得到目前比较广泛使用的暂态能量表达式

这样，暂态稳定评价的直接法包含以下步骤：

（，0） (1）将代入式(6—308)计算故障后的临界能量Vcr。

~ （2）将(cl,cl)代入式（6—308)计算故障切除时刻系统的暂态能量Vcl。 (3）计算能量裕度VtemVcrVcl。如果Vtem0,则系统是暂态稳定的.

可以应用数值积分方法计算出故障切除时刻各发电机的状态(cl,cl)；根据前面的假设,xSpre和xS充分接近，因此以xSpre为初值，用牛顿法求解方程式(6-299）容易得到xS。

4。 临界能量的求取

Vcr为系统UEP相对于SEP的势能，它的计算是TEF法中最困难的工作，下面介绍几种求取Vcr的方法。

(1）最近不稳定平衡点法(The C1osest UEP Approach）。

TEF用于暂态稳定分析的早期，大多用以下方法确定系统最小的Vcr。 1）计算出所有的UEP。

2）计算与每一个UEP有关的系统势能. 3)用最小的势能作为系统的Vcr。

显然，这样得到的Vcr与系统中发生的故障类型和地点无关,即与故障轨迹无

~~~~关，造成对系统稳定性的分析结果过于保守.

(2）控制不稳定平衡点法（The Controlling UEP Approach,CUEP)[41，35,37，38］。 用这种方法计算Vcr时考虑故障的类型和发生的地点，因而对最近不稳定平衡点法的保守性有很大程度的改进.这种方法的依据是观察到所有临界稳定情况的系统故障轨线到达那些与系统分离边界有关的UEP附近.CUEP法的本质是，用通过CUEP的恒定能量界面去近似故障轨迹指向的稳定边界（控制不稳定平衡点的稳定流形）的相关部分。因此，该方法也称为相关不稳定平衡点法(The Relevant UFP Approach）。这种方法的计算过程大致分为两步进行,即干扰模态的识别和CUEP的计算。

干扰模态（Mode of Disturbance,MOD)的识别是指辨识在给定干扰下受扰严重的发电机.如果干扰严重到足以使系统失稳，那么相对于其他发电机来说，这些发电机更容易失去同步。识别干扰模态的一种简单方法是进行数值积分，找出率先失去同步的发电机。

可以通过求解如下极小化问题:

计算系统的CUEP。其中，用基于干扰模态的近似UEP作为初值.

（3）BCU法(The Boundary of Stability-region-based Controlling UEP Method)[43,44].

CUEP法在求解UEP时面临严重的收敛性问题，特别是当初值不很接近CUEP时。BCU法能够发现与故障轨线有关的、准确的CUEP，并且具有可靠的理论基础。此外,这种方法的计算速度也较快。

在BCU法中，SEP的稳定域对于暂态稳定分析是非常重要的。将所研究的电力系统描述为如式(6—288）所示。假设xS是式(6—288）系统的一个渐近SEP［即在xS处f(x)的雅可比矩阵的所有特征值都具有负实部］，那么存在一个包含

xS的域A(xS)，在这个域内出发的任意轨线随着时间的推移将收敛于xS。称域A(xS)为xS的稳定域。A(xS)的边界,表示为A(xS)，称为xS的稳定边界。 BCU法是基于原始电力系统的稳定域边界与一个简化系统(或梯度系统）的

稳定域边界之间的关系来实现的。它把稳定域边界定义为边界上所有UEP稳定流形的交集，从边界上一点出发的任何轨线,随着时间的推移将收敛于某一个UEP。利用这个性质,BCU法通过计算相关简化系统的CUEP来计算原始系统的CUEP，简化系统的CUEP计算起来相对容易且计算量小。对于经典的电力系统模型[式(6—95）]，其简化系统(故障后）被定义为

显然，简化系统式（6-310）的平衡点与原始系统式(6-95)的平衡点相同BCU法计算CUEP的基本步骤如下：

（t),（t)))寻找离开1）对于给定的系统和预想事故，沿着持续故障轨线（(~~简化系统稳定域边界的点，即逸出点（Exit Point）。其有效的计算方法是,用

（t),（t)),直到投影轨线（t)使得VP()出数值积分方法计算系统的故障轨线(~~~~*现首次局部最大。

2）把逸出点作为初始条件，用数值积分方法求解故障后的简化系统,找到使得fi到达首次局部最小的点。

i1m~*~*0 3）以为初值，用鲁棒性好的数值方法求解方程fi0得到i1~*0m~*co

4)把(，0)作为原始系统关于故障轨线的CUEP. 5。 直接法的应用和局限性

对于6.4。4节所提及的动态安全评价，应用直接法有很多优点。直接法除了可以避免对故障后系统费时的数值积分外,还可以提供度量系统稳定程度的定量指标。当需要比较不同规划方案的相对稳定程度或需要快速算出稳定极限或需要得到提高系统哲态稳定的措施时，这些附加的信息相当有用。这些用途使得直接法具有很大的吸引力。

另外，直接法也可用于在电力系统规划、制定运行方式时,对预想事故集进行筛选，过滤掉一些明显稳定或不稳定的事故，对其余的少量事故进行详细的暂态稳定分析。

~*co 尽管直接法在近年来取得了长足的进展，但对模型的和计算结果的可靠性仍然是其工业化应用的主要障碍。众所周知，两机系统的TEF分析等价于著名的“等面积定则”,而目前一些实用的TEF法大多是两机系统TEF分析的直接或间接推广，特别是暂态能量的计算与故障轨线有关,这当然影响其计算结果的可靠性.另外，CUEP求取中繁重的计算量以及差的收敛性可能使得直接法比一般的数值积分法更慢。

当前，暂态稳定分析的混合法似乎是一种较好的选择,即在数值积分过程中计算暂态能量，从而得到更多的反映系统稳定性的信息。