第六章 电力系统稳定性

1、天津大学硕士学位选修课电力系统稳定性分析授课人：魏 炜2022年6月19日Introduction to Power System Stability Analysis第6章 暂态稳定性分析第2页电力系统稳定性分析引言根据IEEE建议，电力系统稳定性可以根据受扰大小进行分类： 大扰动稳定(Transient Stability)：系统在经历一些大的扰动后(线路出现三相短路、变压器突然停业、发电机和重要负荷的突然退出等)，系统是否能够过渡到一个新的运行点(平衡点)； 小扰动稳定(Small Signal/Steady State Stability)：系统在缓慢变化过程中，系统的关键参量不出现急

2、剧的变化(自发的振荡或单调的发散)；等效定义：受到微小扰动后可恢复原来运行状态，也常被称为是运行点的稳定性。静态稳定(Static Stability)：系统在运行过程中，潮流方程存在可行解，同时各种关键设备的限制(发电机出力上下限、线路或变压器的传输极限等) 能够得到满足，也称为Static Power Flow Stability，近年来这方面的很多研究多侧重研究电力系统静态电压稳定性问题。第3页电力系统稳定性分析引言此外，通常还可以根据系统最后失稳的场景和性质进行划分，将系统稳定问题分为： 功角稳定(Power-Angle Stability, Transient Stability)：

3、系统在经历扰动后，某些发电机的转速出现持续加速/减速或与其他发电机产生持续的振荡，系统无法正常运行； 电压稳定(Voltage Stability)：系统在扰动后，系统关键节点的电压出现持续振荡或单调的下降，导致系统无法正常运行；第4页电力系统稳定性分析引言a) 功角失稳和电压失稳的场景，可能是大扰动引起的，也可能是由于小扰动引起的，因此也可以进一步细分为大扰动(暂态)的功角/电压稳定性、小扰动(暂态)的功角/电压稳定性。也即：上述两种分类方法是相互补充的；b) 大扰动情况下系统的功角稳定性是人们长期研究的一个问题，早期城市规模较小，负荷分散，单纯的电压稳定性问题很少发生，因此在学术/工业界，

4、常将暂态稳定性和大扰动的功角稳定性等同，在阅读文献时，需要注意；c) 在我国实际电力调度运行部门习惯上称小扰动稳定(负荷的缓慢变动引起的单调或振荡情况)为动态稳定性；d) 系统稳定是一个统一概念，只是出于研究的方便，才人为地划分了各种稳定性问题。大负荷区暂态稳定第5页电力系统稳定性分析第六章 电力系统暂态稳定分析 提 纲6.1 简介6.2 时域仿真分析方法6.3 网络方程及其相关模型6.4 简单系统稳定性分析6.5 暂稳分析中FACTS元件的考虑6.6 小结第6页电力系统稳定性分析第6章 暂态稳定分析简介 首先强调：暂态稳定分析，特指系统在受到大扰动后能否过渡到一个新的运行点的稳定性问题。因此

5、根据前面的分类，它涵盖了大扰动的功角和电压稳定性。 大扰动的功角稳定性：一般指系统在出现大的扰动后的单摆稳定问题，研究的时间较短(几秒之内) 。 大扰动的(中长期)电压稳定性：需要同时考虑快速变化元件(发电机机电-暂态过程、励磁调节器、PSS装置、HVDC、FACTS元件等)和慢速调节设备(原动机及其调速器、恢复性动态负荷)的动态过程，时间跨度几十秒到几分钟。 现在，大扰动后的多摆功角稳定性也成为业界关注的一个热点问题，这一问题的时间跨度和大扰动后的电压稳定性基本吻合(中国国电公司修改了运行导则，仿真时间50秒)。基于WAMS的系统暂态稳定在线评估研究是一个热点。第7页电力系统稳定性分析第6章

6、 暂态稳定分析简介单摆稳定单摆不稳定多摆不稳定第8页电力系统稳定性分析第6章 暂态稳定分析简介暂态稳定的一般场景系统出现较大扰动 故障点出现短暂的暂态过程(电磁暂态) 线路潮流出现变化 发电机输出的电磁功率受到影响 因发电机原动机系统调节特性较慢(几秒到几分钟) 发电机转子转矩出现失衡 导致部分发电机加速或减速 进一步导致网络中的潮流出现振荡型变化 引起FACTS/HVDC/Exciter/Governor等的进一步变化，在此过程中，故障元件可能被切除，人为的或预先设定的一些控制环节会动作。急剧变化系统部分节点电压出现导致失步速,部分发电机持续加/减处稳定下来系统在一个新的运行点第9页电力系统

7、稳定性分析第6章 暂态稳定分析简介暂态稳定涉及的系统参数变化 网络和相关设备的电磁暂态过程 网络潮流参数的动态变化(包括长线的波过程和谐波) 继电保护装置动作引起的电磁暂态过程 发电机的机电暂态过程(转子在电磁功率变动过程中动态) Exciter/FACTS元件的控制过程 动态负荷的变化 慢调节特性元件的控制和调整(Governor/OLTC/可投切补偿装置)等 时间跨度为毫秒级，一般简化处理为代数变量，即 利用其 稳态结果参与稳定计算(波过程、谐波衰减、定子及线路电 磁暂态过程) 时间跨度在几十秒到几分钟之间，1) 近似处理为代数量； 2) 考虑为动态参数 在传统的暂态稳定研究中，重点关注的

8、动态变量，此方法 称为Quasi Steady State Model/Approach第10页电力系统稳定性分析第6章 暂态稳定分析简介暂态稳定分析模型),(0),(yxgyxfx 暂态稳定分析的电力系统模型，经过部分简化后，可以处理为如下的微分-代数(Differential Algebraic Equation)方程 发电机暂态/次暂态电势； 发电转子角速度/运行角； 励磁系统和PSS动态参数； FACTS元件控制参数 动态负荷参数； Two-Term/Multi-Term HVDC动态参数 网络潮流方程； 发电机电压/电势方程； 不计及Governor时原动机输出功率； HVDC的电压

9、方程； 静态负荷变量； 第11页电力系统稳定性分析第6章 暂态稳定分析简介暂态稳定分析方法总体上分为两类 时域仿真法：直接利用积分运算，对描述系统动态的DAE方程直接进行求解，并监视系统关键参量(功角/电压)，以判断系统稳定性状况。 直接法：利用Lyapunov稳定性定理，列解电力系统适合的能量函数，不通过积分或通过短暂的积分运算，直接判断电力系统的稳定性状况。 基于WAMS量测信息的在线稳定性分析，也是通过判断系统关键量的变化规律来判断系统稳定性第12页电力系统稳定性分析第6章 暂态稳定分析简介(一) 时域仿真法 电磁暂态：需要考虑系统的高次谐波、波过程等，因此仿真步长较短，模型要求精确，对

10、积分算法的数值稳定性要求较高。主要用于电气设备(保护/断路器/HVDC)的控制环节分析 EMTP (Electomagnetic Transients Program)：最著名的电磁暂态仿真程序(Fortran)，最早由德国的Hermann W. Dommel设计，后由BPA (Bonneville Power Administration)组织系统开发，开始主要在大/小型机上运行，后形成微机版EMTP-ATP (Alternative Transients Program)，加入ATP-EMTP协会后，可以免费使用，网址：http:/、 。第13页电力系

11、统稳定性分析第6章 暂态稳定分析简介 EMTDC/PSCAD (ElectroMagnetic Transient in DC System)/(Power System Computer Aided Design)：由加拿大Manitoba HVDC Research Center开发完成，最初主要是在仿真中精确地考虑直流线路。现在功能强大，可以进行用户自定义模型的设定(For/C /C+)，可以进行AC/HVDC组合仿真等，PSCAD是其前台界面，最新版本V4.2，可从http:/下载免费版。(一) 时域仿真法 电磁暂态：需要考虑系统的高次谐波、波过程等，因此仿真步长较短，模型要求精确，对

12、积分算法的数值稳定性要求较高。主要用于电气设备(保护/断路器/HVDC)的控制环节分析第14页电力系统稳定性分析第6章 暂态稳定分析简介 NETOMAC (Network Torsion Machine Control)：由德国西门子公司研制的电磁暂态仿真软件，同时具备机电暂态和稳态分析的相应功能，在欧洲电力系统中应用较为广泛，可以从其程序网站：http:/上关注其最新版本情况。最值得注意的是其二次开发功能。清华、南方电网公司等已引进。(一) 时域仿真法 电磁暂态：需要考虑系统的高次谐波、波过程等，因此仿真步长较短，模型要求精确，对积分算法的数值稳定性要求较高。主要用于电气设备(保护/断路器/

13、HVDC)的控制环节分析第15页电力系统稳定性分析第6章 暂态稳定分析简介(一) 时域仿真法 机电暂态：忽略系统的高次谐波、波过程等，主要关注系统的基频分量的变化规律、发电机转子运动方程以及相关的FACTS元件的控制过程，仿真时间一般较长 BPA (Bonneville Power Administration)：由BPA公司研制开发，曾是北美电力公司主要的机电暂态仿真程序，我国于1983.9引进，汉化并添加前台管理界面后形成中国版BPA程序，现在在华北、国网公司、南方电网等电力企业中广泛使用。但在国外，已是落日黄花！第16页电力系统稳定性分析第6章 暂态稳定分析简介(一) 时域仿真法 机电暂

14、态：忽略系统的高次谐波、波过程等，主要关注系统的基频分量的变化规律、发电机转子运动方程以及相关的FACTS元件的控制过程，方针时间一般较长 PSS/E (Power System Simulator for Engineering)：由PTI-US公司一家由Siemens控股的美国公司研制，现在已是北美电力企业主要的暂态稳定仿真程序，最新版本30.1.2，是集暂态仿真、潮流计算、系统优化、短路分析为一体的一个综合产品，已在世界上123个国家成功应用。我国中华电力、华东电网、南方电网、国电公司等已先后引进。http:/www.pti-第17页电力系统稳定性分析第6章 暂态稳定分析简介(一) 时域

15、仿真法 机电暂态：忽略系统的高次谐波、波过程等，主要关注系统的基频分量的变化规律、发电机转子运动方程以及相关的FACTS元件的控制过程，方针时间一般较长 EUROSTAG：由法国EDF开发，是一款欧洲电力公司普遍使用的电力系统稳定性分析软件，除了可以进行机电暂态仿真外，该程序的最大优点是可以进行中长期稳定性分析(由于其小扰动稳定分析见长)，具有与PSS/E、MATLAB Simulink等程序的接口，可以从其网站上了解其最新进展 http:/www.eurostag.be。国内由清华负责其汉化工作，国网公司和华北电力公司已经正式使用。第18页电力系统稳定性分析第6章 暂态稳定分析简介(一) 时

16、域仿真法 机电暂态：忽略系统的高次谐波、波过程等，主要关注系统的基频分量的变化规律、发电机转子运动方程以及相关的FACTS元件的控制过程，方针时间一般较长 PSAPAC (Power System Analysis Package)：由美国电科院EPRI研制开发，是一款集潮流计算、暂态和中长期稳定仿真、小扰动稳定分析、电压稳定性分析为一体的综合程序，在各大高校使用广泛。后被Powertechlabs继承相应程序并进一步完善，形成DSATools的全套产品，也有相当的电力用户。可在如下网址，察看相关信息http:/和或第19页电力系统稳定性分析第6章 暂态稳定分析简介(一) 时域仿真法 机电暂态

17、：忽略系统的高次谐波、波过程等，主要关注系统的基频分量的变化规律、发电机转子运动方程以及相关的FACTS元件的控制过程，方针时间一般较长 国产BPA/PSASP：由国家电力科学研究院研制，前者属引进消化吸收，后者为自主开发(周孝信/吴中习)，并获首届国家科技进步一等奖，具有潮流计算、短路分析、暂态稳定分析、小扰动稳定分析、电压稳定性计算等相应功能，最新版本情况：6.227.0。在各大电网公司广泛采用，现已有基于大型机的并行计算程序，可以实现超实时仿真。第20页电力系统稳定性分析第6章 暂态稳定分析简介(二) 暂态稳定直接法)()()(),(xxxpxfx EkVVV负定的能量函数 利用直接法进

18、行电力系统暂态稳定性分析，主要工作是寻求更为科学的能量函数和判稳准则； 典型的直接法：PEBS (Potential Energy Boundary Surface ), BCU (Boundary Controlling Unstable equilibrium point approach) Method，Hybrid Method，EEAC Method 难点问题：详细电气模型难以计及；Controlling UEP(事故后不稳定平衡点)的准确求解第21页电力系统稳定性分析第六章 电力系统暂态稳定分析 提 纲6.1 简介6.2 时域仿真分析方法6.3 网络方程及其相关模型6.4 简单系统

19、稳定性分析6.5 暂稳分析中FACTS元件的考虑6.6 小结第22页电力系统稳定性分析第6章 暂态稳定分析时域仿真法时域仿真方法 6.2.1 仿真方法概述6.2.2 常微分方程(ODE)仿真方法简介6.2.3 微分代数方程(DAE)仿真过程6.2.4 暂态稳定仿真计算的一般流程6.2.5 暂稳计算程序优劣比较第23页电力系统稳定性分析第6章 暂态稳定分析时域仿真法Simulation and Integration 微积分是人类数学史上的一个奇迹，它将人类对世界的认识从呆板的静态(Static)世界引导到多姿多彩的动态世界(Dynamic System)中，使人类科学技术上升了一大台阶，现代科

20、学的数学模型几乎都与微积分有关。 谁发明了微积分是人类科学史上一段著名公案牛顿1665.5. 20在一手稿中提到“流数术”(微积分初步思想，很多文献认为此日是微积分创生日)，但因一次离奇火灾，其完整思想在1736出版的流数法和无穷级数才体现的(Newton 1671)；而德国数学家莱布尼茨在1684从几何学角度提出微积分思想，并在1686年给出现在沿用的积分算法及微积分表示及运算符号。第24页电力系统稳定性分析第6章 暂态稳定分析时域仿真法Simulation and Integration 牛顿-胡克(万有引力) 牛顿-莱布尼兹(微积分) Newton-Leibniz公式： 动力系统(Dyn

21、amic System)可以用微分方程描述的动态系统： Time Varying System Time Invariant Systemttdtt)()()(zzz),(),()(ttttpxfx),(pxfx 第25页电力系统稳定性分析第6章 暂态稳定分析时域仿真法Simulation and Integration对于定常动力系统(Ordinary Differential Equation)：(*) )(),(00 xxpxfxt当系统的模型和初值确定后，系统的动态行为将完全确定，由高等数学的知识可知，理论上解析式子：ttdttt0),(),()(0pxfpxFxtx轨迹(Flow,

22、Trajectory)初值1初值2实际上，很多实际的物理系统，其动态方程往往是复杂的非线性(nonlinear)方程，难以直接得到对应的解析式，例如：xxxsin只能通过数值(numerical)方法加以研究(time domain)第26页电力系统稳定性分析第6章 暂态稳定分析时域仿真法Simulation and Integration影响常微分方程性态(包括稳定性)的两个关键因素：1) 方程(模型)本身；2) 方程的初值电力系统暂态仿真成败的关键，与上述两个因素息息相关：1) 方程(模型)本身性态的保持：选用具有足够精度、数值稳定性好且求解速度较快的数值积分算法在模型准确的前提下，仿真轨

23、迹与系统真实轨迹足够接近；2) 精确设定系统的初始值：往往需通过仿真场景初始时刻的系统潮流方程来决定；第27页电力系统稳定性分析第6章 暂态稳定分析时域仿真法00)( . . )(xxxfxtts 考虑如下ODE方程，数值积分的目的在于给定初始点信息，求解一段时间区间上的系统的运行轨迹，为后续系统稳定性研究提供基础。Our Target is:utttt0 ?)(x第28页电力系统稳定性分析第6章 暂态稳定分析时域仿真法 ODE方程数值积分方法有：欧拉法(Euler Method)，改进欧拉法(Modified Euler Method/ Improved Euler Method)，龙格-库

24、塔法(Runge-Kutta Method)，线性多步法(如 Adamas 插值法等)，几乎所有的数值积分算法，均是在Euler Method基础上发展起来的(本科电力系统分析中提到的分段计算法)，我们简单回顾一下各种方法的基本原理，重点放在各种方法的误差分析和比较上。第29页电力系统稳定性分析第6章 暂态稳定分析Euler Method1、Euler Methodtltutx)(xfx Our target is the area under the line 将积分区间等分为m 段，区间长度为 ，并认为每一个期间n内的导数近似等于期间起始位置的导数f(xn)，即： mtthlu)(,0,1

25、,2,),(nnntmndtdxxxfx第30页电力系统稳定性分析第6章 暂态稳定分析Euler MethodSteps of Euler Method1. n=0时，初始点的导数 f(x0) 已知，如前述，根据它就可以求解(推出)n=1时的初始点信息)(001xfxxh2. n=1时，由x1 f(x1) ，再根据f(x1)就可以求解(推出)n=2时的初始点处的信息x2 )(112xfxxh3. 可得递推公式：)(1nnnhxfxx)(1nnnnhxxdtdxfx(6-9)欧拉法又称为欧拉切线法或欧拉折线法，欧拉法又称为欧拉切线法或欧拉折线法，基本思想是将积分曲线用折线来代替，基本思想是将积分

26、曲线用折线来代替，折线曲率就对应着等分点处由微分方程折线曲率就对应着等分点处由微分方程右端项给出的广义切线，其示意图见右端项给出的广义切线，其示意图见p.296的图的图6-2。 第31页电力系统稳定性分析第6章 暂态稳定分析Euler MethodErrors of the Euler Method 对于非线性微分方程组的解曲线，在第n个积分段对时间 t 进行Talyor级数展开： dt,t)()(xfxFx)(2)( )(232321hOhhhOhhnnnnnnn xxfxxxxx第32页电力系统稳定性分析第6章 暂态稳定分析Euler MethodErrors of the Euler M

27、ethod)()(21321nnnnnnnhhOhhxfxxxxxx 可以看到，实际上Euler 法只用到了Talyor级数的线性项，二阶以上项均形成此次积分运算的截断误差(Truncation Error)，同时可以看到，除了第一次n=0时，系统积分不存在截断误差外，后续运算都存在TE误差，并产生积累效应，我们称每一步TE误差为局部TE误差 (Local Truncation Error)，累积效应为累积截断误差或全局截断误差(Global Truncation Error)。 第33页电力系统稳定性分析第6章 暂态稳定分析Euler MethodTruncation Error of th

28、e Euler Method可积先决(Lipschitz )条件：nnnnLxxxfxf11)()(积分的真值：11)(nnnnThhxfxxEuler近似值：)(1nnnhxfxx这一步误差：1111)()(nnnnnnnnThhexfxfxxxx1)()( nnnnnThhxfxfxx1 nnnnThLhexxnnTT1maxmax :令max)1 ( TheLhn第34页电力系统稳定性分析第6章 暂态稳定分析Euler MethodTruncation Error of the Euler Method这一步误差： 111nnnexxmax)1 ( TheLhnmaxmax1)1 ()1

29、 ( ThTheLhLhnmax121)1 ()1 ( ThLheLhn max011)1 (.)1 ()1 ( ThLhLheLhnnmax101)1 (1)1 (1)1 ( ThLhLheLhnn第35页电力系统稳定性分析第6章 暂态稳定分析Euler MethodTruncation Error of the Euler Methodmax1011)1 (1)1 (1)1 (ThLhLheLhennn第n步误差：max1011)1 ()1 ( ThLhLheLhnnmax101)1 ()1 ( ThLhLheLhnnLTeLhnmax01)1 ( nyneye)1 (00考虑到：LTeL

30、hnmax) 1( )(maxhOT由于因此最后可得：)(1hOen系统累积误差为O(h)第36页电力系统稳定性分析第6章 暂态稳定分析Euler Method 全局截断误差(GTE)和局部截断误差(LTE)，与步长h存在O(h)和O(h2)关系，当步长减小时，两种误差都会减小； 除了截断误差外，还存在计算机计算时的舍入误差，它与计算步骤近似成正比关系，因此需要合理选择计算程序的最小步长，如p.298所示的图6-3，理论上存在一个最佳步长，注意不是越小越好； 欧拉法属于自启动的和显式数值积分方法，只需要知道方程和初值即可计算，但其计算误差受步长影响较大，同时对于刚性系统，其稳定性较差，数值方法

31、中称为非A-稳定的。 Summary of the Truncation ErrorRound-Off Error第37页电力系统稳定性分析第6章 暂态稳定分析Modified Euler Method2、Modified/Improved Euler Method基本原理：由于Euler法只利用了积分曲线积分段起始位置的斜率信息，如果将积分区间的斜率取为该区间前后两点的斜率的平均值，则计算精度可以得到一定的提高。但积分区间后一点的精确值是未知的，因此只能利用近似结果来估计它。 第38页电力系统稳定性分析第6章 暂态稳定分析Modified Euler Method2、Modified/Imp

32、roved Euler MethodEuler Method ImprovedEuler Method截断误差截断误差一次Euler计算第39页电力系统稳定性分析第6章 暂态稳定分析Modified Euler MethodSteps of M/I Euler Method1. 利用前面计算结果，计算起始点的斜率值)(nnxfx 2. 利用Euler算法，估算积分区间n, n+1后一点的结果)()0(1nnxfxxhn3. 利用前一步估计结果，估计后一点斜率近似值 )()0(1)0(1nnxfx 4. 利用前后两点斜率平均值作为起点处的斜率进行积分运算 )()(21)0(11nnnnhxfxf

33、xx梯形面积计算公式前两步实际即是一次Euler过程第40页电力系统稳定性分析第6章 暂态稳定分析Modified Euler MethodSteps of M/I Euler Method最后的计算公式 )()(21)0(11nnnnhxfxfxx计算xn+1需要用到两个前面的已知变量:xn前一次运算的结果 本次运算的中间结果)0(1nx意味着在内存中需存储3个已知结果可以从第二步中解得)()0(1nnnhxfxx带入上式：)()(21)0(1)0(11nnnnhxfxfxx作为一个内存值第41页电力系统稳定性分析第6章 暂态稳定分析Modified Euler MethodSteps of

34、 M/I Euler Method最后的计算步骤 )(nnxfx )()0(1nnxfxxhn)()(21)()0(1nnnxfxfxf )()0(11nnnhxfxx 采用上述计算步骤，中间存储的临时变量将减少，计算机的存储量也会相应减少 p.300 (6-19)第42页电力系统稳定性分析第6章 暂态稳定分析Modified Euler MethodErrors of M/I Euler Method)()(21)0(11nnnnhxfxfxxImproved Euler Method)(21)(21 )0(1nnnhhxfxfx)(21)(21 nnnnhhhxfxfxfx一次Euler计

35、算)()(21)( 32hOdtddtdhhnnnnxfxxfxfx)()()()(2hOdtdhdtdhhnnnnnxfxxfxfxfxfdtddtddtddtd)(22xfxxfxxx 考虑到：)(21)(321hOhhnnn xxfxxTalyor展开式改进Euler法局部截断误差O(h3)第43页电力系统稳定性分析第6章 暂态稳定分析Modified Euler MethodErrors of M/I Euler Method 改进欧拉法实际是利用了Taylor展开式的前两项(线性项和二次项)，而欧拉法则只是用了展开式的线性项(一次项)。 与分析欧拉法的累计误差的过程类似，通过递推过程

36、，可以证明改进欧拉法的累计截断误差为O(h2)，因此具有比欧拉法更高的计算精度。 一般算法局部截断误差为O(hk)，全局截断误差为O(hk-1)第44页电力系统稳定性分析第6章 暂态稳定分析Modified Euler MethodSummary of M/I Euler Method 局部截断误差和全局截断误差，与步长O(h3)和O(h2)的关系，当步长减小时，两种误差都会减小； 由于增加了一个预测环节，其计算量和存储量约为欧拉法的2倍； 同样属于显式积分方法，数值稳定性与欧拉法相似。 第45页电力系统稳定性分析第6章 暂态稳定分析Runge-Kutta Method3. Runge-Kut

37、ta Method 改进Euler法的原理在于，通过增加一次“内嵌式”的Euler计算过程，实现将计算的局部截断误差由原来的O(h2)降为O(h3)， “内嵌式”的Euler计算可看作是一次内部插值过程，那么能否通过增加一些插值过程，使得算法的局部截断误差进一步降低呢？ 龙格-库塔(Runge-Kutta)法的原理第46页电力系统稳定性分析第6章 暂态稳定分析Runge-Kutta MethodRunge-Kutta Method Runge-Kutta法通过在一次“外层”Euler过程中增加一次或多次的内部插值，实现局部截断误差的显著降低，从而成为一种常见的显示积分方法，最常见的RK方法有三

38、阶、四阶方法(二阶方法实际就是改进Euler法)。其中RK4阶的运算量约为Euler法的4倍，而存储量约为Euler法的2.5倍。 第47页电力系统稳定性分析第6章 暂态稳定分析Runge-Kutta MethodRunge-Kutta Method递推公式：其中：)22(6143211kkkkxxnn)(1nhxfk)2(12kxfknh)2(23kxfknh)(34kxfknhRK-4局部截断误差为O(h4)第48页电力系统稳定性分析第6章 暂态稳定分析Runge-Kutta MethodErrors of Runge-Kutta MethodRK-n2345678局部截断误差O(h3)O

39、(h3)O(h4)O(h4)O(h5)O(h6)O(hn-2) J.C.Butcher, The Numerical Analysis of Ordinary Differential EquationsRunge-Kutta and General Linear Methods, John Wiley & Sons, 1987. 不同阶次RK方法的局部误差比较 RK-2即前面所述的改进Euler法 RK-2/3、RK-4/5以及RK-7/8局部截断误差相应，适应性不同第49页电力系统稳定性分析第6章 暂态稳定分析Runge-Kutta MethodRunge-Kutta Method in

40、MatlabCmdsDescriptionODE23Solve non-stiff differential equations, low order method.ODE45Solve non-stiff differential equations, medium order method.ODE113Solve non-stiff differential equations, variable order method.ODE15SSolve stiff differential equations and DAEs, variable order method.ODE23SSolve

41、 stiff differential equations, low order method.ODE23T Solve moderately stiff ODEs and DAEs, trapezoidal rule.ODE23TB Solve stiff differential equations, low order method.第50页电力系统稳定性分析第6章 暂态稳定分析Runge-Kutta MethodRunge-Kutta Method in Matlab第51页电力系统稳定性分析第6章 暂态稳定分析Linear multistep method4. Linear Mult

42、istep Method 欧拉法(Euler)、改进欧拉法(I/M Euler)以及龙格-库塔(RK)方法，都仅仅用到了积分区间段起始位置(n)的信息，而没有考虑已经积分完成的前面区间段的相应信息(例如 n-1, n-2, n-3)，多步法则考虑了上述结果的影响，首先利用前面多步计算结果的信息，进行后续点预测，再调用校正环节，对预测结果进行校正 第52页电力系统稳定性分析第6章 暂态稳定分析Linear multistep method4. Linear Multistep MethodAdams 3-Step Method: 预测环节： 校正环节：2115162312nnnnnhxxxxx2

43、1118512nnnnnhxxxxxAdams 4-Step Method: 预测环节： 校正环节：3211937595524nnnnnnhxxxxxx32111519924nnnnnnhxxxxxx第53页电力系统稳定性分析第6章 暂态稳定分析Linear multistep method4. Linear Multistep Method 采用上述预测校正环节，可以有效提高算法的收敛精度，可以将全局截断误差降为O(h3)，但其内存存储量较大(但较同阶RK方法要小)； 算法本身需要采用单步法进行启动； 在Adams思想的前提下，根据不同的实现思路，形成了不同的改进线性多步法，如常见的Adam

44、s-Bashforth、Adams-Moulton、Variable Stepsize Adams method，可从如下网址下载相关算法源程序：http:/ 第54页电力系统稳定性分析第6章 暂态稳定分析显式积分算法小 结 Euler method, improved/modified Euler method, Runge-Kutta method, linear multistep method (Admas method) are all belonging to the explicit integration methods. Their local and global trun

45、cation errors can be reduced when the step size turns small. Since all the explicit methods are derived based on the Taylor series expansion, they require the differential equations to be smooth enough. Trunction errors of the explicit methods can be reduced through cutting down stepsize and add mor

46、e mid-points in a iteration. However, we also have another way第55页电力系统稳定性分析第6章 暂态稳定分析Implicit Method隐式梯形法对于微分方程： 在已知tn结果xn后，下一节点tn+1=tn+h积分结果xn+1，理论上可通过求解如下积分过程得到： 当h足够小时，可认为积分区间段内的f(x)不变(p.304图6-6)，从而上式可以近似表示为： nRdtdxxfx),(dtnnttnn1)(1xfxx这样构成了n个代数方程，其中xn+1为待求的n个变量，通过迭代运算即可求解上述n+1步的积分结果。 +=+11 (*)2

47、nnnnh()()xxf xf x第56页电力系统稳定性分析第6章 暂态稳定分析Implicit Method隐式梯形法 对比一下隐式梯形法的(*)式和改进欧拉法的计算式，可以看到，两种方法的原理类似，但改进欧拉法直接利用xn预测xn+1处的值，而隐式梯形法则采用代数迭代运算对上述公式进行求解，可以验证，上述计算式(*)的局部截断误差为O(h3)。此外，可以采用类似的方法，将前面所介绍的各种显式方法，“改造”为隐式求解方法，但隐式求解方法计算更为复杂。 第57页电力系统稳定性分析第6章 暂态稳定分析Implicit Method隐式积分法优缺点 可以采用较大的积分步长，并不影响算法的数值稳定性

48、； 在进行微分代数方程(Differential Algebraic Equation, DAE)的积分运算时，可以联立求解，消除交接误差，保证算法的稳定性； 算法较为复杂； 算法的扩展性较差，尤其在DAE问题的求解时，往往需要进行相应的变通，形成所谓的类隐式方法。隐式方法的数值稳定性较显示方法要好得多第58页电力系统稳定性分析第6章 暂态稳定分析ODE积分算法小结 除上述方法外，还有很多成熟的数值积分方法，如 采用自适应步长(变步长)积分方法、有针对刚性问题或病态方程的适用求解方法、有具有更好A-稳定的类隐式方法(Quasi-implicit)等，一些专门介绍数值积分方法的专著上有更为深入的

49、介绍。 第59页电力系统稳定性分析第6章 暂态稳定分析时域仿真法时域仿真方法 6.2.1 仿真方法概述6.2.2 常微分方程(ODE)仿真方法简介6.2.3 微分代数方程(DAE)仿真过程6.2.4 暂态稳定仿真计算的一般流程6.2.5 暂稳计算程序优劣比较第60页电力系统稳定性分析第6章 暂态稳定分析DAE SystemDAE Dynamic SystemWhere: is the vector of the state variables. is the vector of the algbraic variables.In the DAE dynamic system, dynamics

50、 of the state variablesis constrained on the solution space of the algbraic equation(constraint manifold). ),(0),(yxgyxfx mRynRx第61页电力系统稳定性分析第6章 暂态稳定分析DAE Integration Method 针对DAE系统的仿真方法总体可分为两种类型：交替求解法：微分方程和代数方程分别进行求解，各自可采用不同的求解算法；联立求解法：微分方程和代数方程联立进行求解，微分方程一般采用隐式方法进行差分化后与代数方程联立进行求解。第62页电力系统稳定性分析第6章

51、暂态稳定分析交替求解法 根据t时刻的已知结果(xt, yt), t=0, h, 2h, ，代入微分方程，然后采用某一积分方法(Euler, Modified Euler, Runge-Kutta Method,)，求解出t+1动态参数xt+1；再传递给代数方程，调用相应的代数求解方法(Newton-Raphson Method, Quasi Newton Method,)，求解出t+1时刻相应的代数量yt+1，如此往复，就可以得到系统的积分解。 第63页电力系统稳定性分析第6章 暂态稳定分析交替求解法以动态方程采用Euler法为例Step1: (x0, y0)采用欧拉法，计算出xh= x0+h

52、f(x0, y0)；Step2: 以(xh, y0)作为初始值，代入代数方程，采用迭代 法，求解yh； Step3: 在(xh, yh)处，再次采用欧拉法，可计算出x2h； Step4: 以(x2h, yh)作为初始值，代入代数方程，采用迭代 方法求解y2h； 教材308页给出动态方程采用RK-4法时，一个具体求解周期的相关运算式；而309页则给出了动态方程采用隐式梯形方法时，一个具体求解周期的相关运算式。 第64页电力系统稳定性分析第6章 暂态稳定分析联立求解法基本思想 在t时刻将微分方程直接利用隐式求解方法加以差分化，处理为代数方程组，然后和代数方程一起形成一个增广的代数矩阵，再调用代数方

53、程求解方法(NR Method, Quasi Newton Method, ) ，直接求解出t+1时刻系统的动态量和代数量。 第65页电力系统稳定性分析第6章 暂态稳定分析DAE求解方法比较例子：利用交替法求解，微分方程部分至少需存储一个的nn雅可比矩阵，代数方程求解时至少需要存放一个mm的矩阵，因此内存开销至少为n2+m2；而采用联立法求解时，内存至少需要存储一个(n+m)2的矩阵，注意到(n+m)2= n2+m2+2mn. 当n=500, m=1000，数据采用double(real*8)时：交替求解法的开销：8*(500*500+1000*1000)/1024/1024=9.54MB 联

54、立求解法的开销：8*(1500*1500)/1024/1024=17.17MB 交替求解法较后者，计算机实现存储量较少； 交替求解法较后者实现灵活：微分方程和代数方程部分的求解均可采用最适合的办法；交替求解法较后者易于扩展：当有新的(动态/静态)元件时，只需修改相应的求解模块(微分/代数)部分即可；而联立求解需同时兼顾，实现复杂；交替求解法(尤其是微分方程采用隐式法时)，会产生交接误差，而联立求解则可以避免。 第66页电力系统稳定性分析第6章 暂态稳定分析DAE方法比较(误差)数值算法的误差产生原因主要包括：截断误差(Truncation Error)实际的物理系统均是非线性的(nonline

55、ar)，在某一点处对其进行Taylor级数的展开，往往对应的级数项为无穷大的情况，而每一种数值积分方法(无论其精度多高)，实际上都是利用积分区间段内有限项Taylor级数来近似表征真解的，由此产生的误差就对应着截断误差。截断误差会产生积累情况，对应着积分过程的全局截断误差。 第67页电力系统稳定性分析第6章 暂态稳定分析DAE方法比较(误差)数值算法的误差产生原因主要包括：舍入误差(Rounding/Round-Off Error)由于计算机位数有限，在表示一个无理数时，只能取其有限次项 (注意，因为计算机均采用二进制运算器和存储器，即使是十进制的有理数，在用二进制表示时，也可能对应着无穷级数

56、)。在计算程序中采用高精度的数据类型，可以一定程度上减少舍入误差，但存储量会随之增大。 C/C+：Single(float, int) Double(double, long)Fortran：Integer 2,4,8； Real 2,4,8第68页电力系统稳定性分析第6章 暂态稳定分析DAE方法比较(误差)数值算法的误差产生原因主要包括：交替误差(Alternation Error)在求解微分-代数方程时，由于代数量和微分量取不同时刻的值，由此会给系统引入相应的误差，成为交接误差，可以通过将微分方程差分化后与代数方程一同求解解决。 大型商业化的仿真计算程序(BPA/ETMSP/PSASP/

57、PSS-E)均采用联立求解法第69页电力系统稳定性分析限制误差(Limit Induced Error)实际系统中很多变量(代数量：OLTC变比上下限，发电机有功/无功上下限；微分量： Efd，SVC的Bc等)，如果限制并非刚好落在积分区间段节点上，会产生相应的误差。 第6章 暂态稳定分析DAE方法比较(误差)数值算法的误差产生原因主要包括：近似/模型误差(Model Error)数学模型是对真实物理系统的一种近似，因此必然存在一定的误差。 第70页电力系统稳定性分析第6章 暂态稳定分析DAE方法比较(误差)数值算法截断误差PK舍入误差ErrStepsizeRETETE第71页电力系统稳定性分析第6章 暂态稳定分析DAE积分方法方法精度系统考虑方法的截断误差、舍入误差和交替误差等因素，并在此基础上合理的选择求解算法及其步长。例如：进行电磁暂态仿真的计算程序(ElectroMagnetic)和进行机电暂态仿真的计算程序(ElectroMach