基于显式积分方法(龙格库塔法)的电力系统暂态稳定计算.doc

第1章 绪论1.1本课题研究的依据和意义随着电力系统规模的不断扩大，暂态稳定性问题日趋严重。电力系统一旦失去稳定，往往造成大范围、较长时间停电，给国民经济和人民生活造成巨大损失和严重危害，在最严重的情况下，则可能使电力系统崩溃和瓦解。长期以来，国内外的专家、学者对如何保证和提高电力系统的暂态稳定性进行了大量的研究工作，并且至今仍将其作为电力系统方面的一个重要研究课题。特别在我国，由于目前输电系统建设滞后于电源的建设，高低压电磁环网结构较多，且电网间联系薄弱，从而更易发生暂态稳定性破坏事故。目前，由于我国经济发展较快，电力增加难以满足国民经济增长的需要，难以在短时间内通过网络建设、备用容量的增加以及新型装置的投入来解决系统的稳定运行问题，因此在电网实际运行中，如何运用各种控制调节手段提高系统的稳定运行水平就显得更为重要。电力系统遭受大干扰是人们所不希望的，但事实上又是无法避免的。系统在遭受大干扰后失去稳定的后果往往非常严重，甚至是灾难性的。事实上电力系统遭受到的各种大干扰，诸如短路故障，大容量发电机、大的负荷、重要输电设备的投入或切除等都是以一定的概率随机地发生，因此系统的设计、运行方式的制定总是需要保证系统在合理选择的预想事故下能够保持稳定，而此系统的设计、运行方式的制定总是需要保证系统在合理选择的预想事故下能够保持稳定，而不能要求电力系统能承受所有干扰的冲击。判断电力系统在预想事故下能否稳定运行，需要进行暂态稳定分析。当系统不稳定时，还需要研究提高系统稳定的有效措施；当系统发生重大稳定破坏事故时，需要进行事故分析，找出系统的薄弱环节，并提出相应的对策。暂态稳定计算，是电力系统最基本的计算之一。它的主要目的是确定电力系统在受到一定扰动(如输电系统故障、发电机切除、大负荷突变等)后，全系统维持同步运行的能力，求出系统状态变量(如发电机功角、母线电压、线路潮流等)对扰动的响应，并在此基础上分析影响电力系统动态稳定的各种因素，研究提高电力系统稳定性的措施。1.2 电力系统暂态稳定概述电力系统中各台同步发电机在同步运行状态下，输出的电功率为定值，系统中各节点电压和各支路的潮流也是定值，这就是电力系统的稳定运行状态。当电力系统在某一正常稳定运行状态下受到某种干扰后，如果能够经过一定的时间后回到原来的运行状态或者过渡到一个新的稳态运行状态，这时我们就认为系统在该正常运行状态下是稳定的；反之，若系统不能回到原来的运行状态或者不能建立一个新的稳定运行状态，以至最后使各发电机间失去同步，则说明系统是不稳定的。一般是表示发电机组对系统或系统对系统间的同步运行稳定性 。 电力系统稳定性通常是按照干扰的大小和控制系统的控制方法来分类的。电力系统按照所承受的干扰大小可分为静态稳定和暂态稳定两大类 。电力系统的静态稳定是指电力系统受到干扰后，不发生非周期性的失步，自动恢复到起始运行状态的能力，对于简单电力系统，可以用 作为静态稳定的判据；电力系统的暂态稳定是指电力系统受到大干扰后，各同步发电机保持同步运行并过渡到或恢复到原来稳定运行状态的能力，通常指第一或第二振荡周期不失步，一般用各发电机的转子相对角度随时间的变化曲线（摇摆曲线）的变化规律来判断系统的稳定性。 在正常的稳定运行情况下，电力系统中各发电机组输出的电磁转矩和原动机输入的机械转矩平衡，因此所有发电机转子速度保持恒定。但是电力系统一旦发生短路事故等大干扰时，便发生了电力振荡。这种振荡使发电机组间发生相对运动，相对角不断变化，从而影响发电机的输出功率特性、发电机端子电压的励磁系统控制特性、负荷电压频率等，并相互影响引起很复杂的振荡。控制系统对此采取了抑制措施，经过一个暂态机电过程恢复原来的稳定状态或达到新的稳定状态。如果控制措施不当，则可能扩大振荡，系统就不能稳定，从而引发发电机失步、失磁等异常现象。电力系统稳定性与系统结构、运行方式、调节装置的参数和干扰的大小、地点及延续时间等有关。在一种干扰下时稳定的系统，在另一种更大的干扰下可能是不稳定的，所以说没有绝对稳定的系统。也就是说，一个系统的暂态稳定情况和系统原来的运行方式以及干扰方式有关，同一个系统在某个运行方式和某种扰动下是暂态稳定的，但是在另一个运行方式和另一种扰动下它可能就是不稳定的。因此，在分析一个系统的暂态稳定性时，首先必须结合系统的实际情况定出系统的初始运行方式。电力系统受到大扰动，经过一段时间后，或是逐步趋向稳态运行，或是趋向失去同步。这段时间的长短与系统本身的状态有关，有的持续约一秒钟(例如联系紧密的系统)，有的则要持续几秒钟甚至几分钟。也就是说，在某些情况下只要分析扰动后一秒钟左右的暂态过程就可以判断系统能否保持稳定，而在另一些情况则必须分析更长的时间。由于在扰动后的不同时间里系统各部分的反应不同，在分析暂态稳定时往往按下面三种不同的时间阶段分类:(1)起始阶段 即故障后约一秒钟内的时间段。在这期间系统中的保护和自动装置有一系列的动作，例如切除故障线路和重合闸，切除发电机等。在这个时间段中发电机的调节系统起不到明显的作用。(2)中间阶段 在起始阶段后，持续5秒钟左右的时间段。在此期间发电机的调节系统将发挥作用。(3)后期阶段 在故障后几分钟时间内。这时热力设备(如锅炉等)中的过程将影响到电力系统的暂态过程，另外，系统中还将发生由于频率的下降而自动切除部分负荷等操作。分析电力系统暂态稳定时要采用一些假设，这里介绍几个最基本的假设(1)由于发电机组惯性较大，在所研究的短暂时间里各机组的电角速度相对于同步角速度(314ard/s)的偏离是不大的。所以，在分析系统的暂态稳定时假定在故障后的暂态过程中，网络中的频率仍为50Hz。(2)忽略突然发生故障后网络中的非周期分量电流。一方面是由于它衰减较快，另一方面，非周期分量电流产生的磁场在空间不动，它和转子绕组电流产生的磁场相互作用将产生一同步频率交变、平均值接近于零的制动转矩。此转矩对发电机的机电暂态过程影响不大，可以略去不计。(3)当故障为不对称故障时，发电机定子回路中将流过负序电流。负序电流产生的磁场和转子绕组电流的磁场形成的转矩，主要是以两倍同步频率交变的，平均值接近于零的制动转矩。它对电力系统的机电暂态过程也没有明显影响，也可略去不计。如果有零序电流流过发电机，由于零序电流在转子空间的合成磁场为零，它不产生转矩，完全可以略去。(4) 忽略负荷的动态影响。此外，在简化计算中，还忽略暂态过程中发电机的附加损耗。根据计算结果精度的不同要求，以及分析方法本身的限制，还将对各元件的数学模型采用不同程度的简化，有时甚至对一部分发电机或系统中的某些部分进行动态等值的简化处理。根据以上几个假设，网络中的电流、电压只有频率为50Hz的分量，也就是说描述网络的方程仍用代数方程。1.3电力系统暂态分析计算方法的比较与选择目前，电力系统暂态稳定分析方法基本分为两种。第一种方法是数值积分方法，又称间接法，其基本思想是用数值积分方法求出描述受扰运动微分方程组的时间解，然后用各发电机转子之间相对角度的变化判断系统的稳定性。其优点是系统模型精致，计算结果准确，能提供系统中各种变量的时间响应；缺点是计算量非常大、耗时长，而在线应用感兴趣的时间仅在几秒内，而且不易导出稳定程度的定量信息以及对系统关键参数的灵敏度分析。数值解法是目前广泛应用的分析方法，已发展得比较成熟，并基本上能满足电力系统规划、设计和运行过程中所进行的离线暂态稳定分析对计算速度和精度的要求。另一种方法是直接法，它不需要求解微分方程组，而是通过构造一个类似于“能量”的标量函数，即李雅普诺夫函数，并通过检查该函数的时变性来确定非线性系统的稳定性质，因此它是一种定性的方法。其优点是能提供系统稳定程度的定量信息；提供系统稳定与对系统关键参数或运行条件变化的灵敏度分析；对极限参数计算速度快，可快速扫描系统暂态过程。缺点是模型能力差，计算结果具有保守性。而且构造李雅普诺夫函数比较困难，因此，目前电力系统暂态稳定分析的直接法仅限于比较简单的数学模型，或用暂态能量函数近似李雅普诺夫函数，因此其分析结果尚不能令人完全满意。其中，数值解法又分为两类。一类是在电力系统稳定分析中广泛应用的隐式梯形积分法。梯形积分法具有良好的数值稳定性和对刚性微分方程组的适应性，从而可以采用较大的步长（最大可达5周波），并且在发生不连续时无需重新起步。另一类是包括欧拉法、改进欧拉法和龙格-库塔法的显示数值积分法。本毕业设计采用的即是龙格-库塔法，其运算精度高，数值稳定性较好，但是运算量大，为欧拉法的四倍。第2章 模型的建立随着电力系统的规模持续增大，结构日益复杂，元件不断更新，电力系统运行对电力系统的分析、规划和控制的方法不断地提出新的、更高的要求。与此相适应的是计算工具和计算数学以及其他技术领域也在不断地进步，为研究电力系统提供了新的手段。现代电力系统分析目前大多是以电子数字计算机为计算工具，因而，建立描述电力系统的数学模型是研究分析电力系统备种专门问题的基础。把数学与客观物理系统联系起来的过程就是通常所说的建模过程。数学模型的正确性和准确性是保证计算结果的正确性和准确性的基本前提。2.1同步发电机的数学模型电力系统中的电源是同步发电机。本毕业设计中只讨论为常数的简单模型。不计阻尼绕组且忽略励磁绕组的暂态过程，认为励磁调节器的控制作用使得。对于简单电力系统，用标幺值描写的发电机转子运动方程为 （2-1）2.2负荷数学模型负荷是电力系统的一个重要组成部分。分析电力系统在各种状态下的行为必须建立负荷的数学模型。 最简单的负荷模型是将负荷用恒定阻抗模拟，即认为在暂态过程中负荷的等值阻抗保持不变，其数值由扰动前稳态情况下负荷所吸收的功率和负荷节点的电压来决定。这种模型十分粗略，但由于其简单，在计算精度要求不太高的情况下仍有广泛应用。2.3网络的数学模型电力网络是保证电力安全生产的重要物质基础，电网结构是系统的重要资源，充分利用现有电网结构可以降低供电成本。电网网架结构对系统的稳定有关键作用，如网架结构的紧密程度、系统送受端容量的大小及比例、送电距离等。系统的暂态稳定水平依赖电网结构，同时改变电网结构可提高系统暂态稳定水平。暂态稳定计算时，网络方程用代数方程描述，通常用直角坐标下的节点电压方程表示，可以表示为功率平衡形式或电流平衡形式。1).功率平衡形式，写成展开形式为: (2-2)2).电流平衡形式。写成展开形式为: (2-3)2.6电流源、负荷与网络方程联立求解在电力系统中，电力网络将系统中看起来相互独立的所有暂态元件联系在一起，在暂态过程中的任一时刻，各暂态元件注入网络的电流不但由其本身的特性决定，而且整个电力网络必须满足基尔霍夫定律，其中前者由各暂态元件自身的代数方程描述，后者反映在电力网络方程中。因此，为了求解网络方便，需要列出各暂态元件自身的代数方程，并对其进行处理，从而可以和网络方程联立求解。对于各种发电机模型，d-q坐标系下得定子电压都可统一为 （2-4）式中: 、。分别表示同步电机d轴和q轴的电势和电抗，随着同步电机所用模型的不同而不同，它们的取值可通过对比式(2-10)与原始定子电压方程而得到。把坐标变化式用于式(2.12)，可得到x一坐标系下的定子电压过程: (2-5)将上式加以整理，即为发电机节点注入电流的表达式 (2-6)式中：， ， ， ， 将式(2.14)得到注入电流表达式代入网络方程式，经整理后可以看出，台发电机接于系统相当于在网络中的相应母线上接入一个电流源。 (2-7)这个电流也称为发电机的虚拟电流，并且相当于在网络导纳矩阵中的相应对角块上应加上矩阵。由此可见，发电机接入系统后，在暂态过程中的任何时刻网络方程仍为线性方程，但其中的发电机虚拟注入电流及相应的导纳矩阵是发电机本身的状态量 、 、的函数，因此这个线性方程是时变的。在同步电机采用简化模型时，相应的网络方程能够得到简化，方程可保持为n阶的复型方程。除非网络方程发生故障或操作，否则网络方程的系数矩阵为定常矩阵。这样在整个暂态过程中，为了求解网络方程，只需在网络发生故障或操作时对方程的系数矩阵重新进行求解。在暂态稳定计算中，需要特别指出以下几点:1)当发生网络故障或操作时，需要修改导纳矩阵，以反映故障或操作。2)由于在调节系统中存在各种限制环节，在计算过程中当有关变量超出下界或上界时，它们将被限制在其下界或上界处，直至变量重新回到其上、下界范围内为止。3)由于忽略网络中的电磁暂态过程，各节点的电压、电流以及发电机和负荷的功率，在网络故障或操作瞬间将发生突变，但微分方程中的状态变量是连续变化的。为此，在发生故障或操作后，需要根据故障或操作瞬间状态变量的取值重新求解网络方程。首先将其等值导纳并入电力网络，从而得到考虑了恒定阻抗负荷后的电力网络方程，它在整个暂态过程中是保持不变的。而负荷的等值导纳根据潮流计算的结果可以很容易的得到： 其中，是负荷所吸收的功率，是各负荷节点电压第3章 故障类型3.1 电力系统故障类型电力系统可能发生的故障类型比较多，其中短路故障是电力系统中最常见危害最严重的故障。所谓短路，是指一切不正常的相与相之间或相与地(对于中性点接地的系统)发生通路的情况。产生短路的原因很多，主要有如下几个方面:（1）元件损坏，例如绝缘材料的自然老化，设计、安装及维护不良所带来的设备缺陷发展成短路等;（2）气象条件恶化，例如雷击造成的闪络放电或避雷器动作，架空线路由于大风或导线覆冰引起电杆倒塌等;(3)人为事故，例如运行人员带负荷拉刀闸，线路或设备检修后未拆除接地线就加上电压等;(4)其它，例如挖沟损伤电缆，鸟兽跨接在裸露的载流部分等。在三相系统中，可能发生的短路有:三相短路、两相短路、两相短路接地和单相接地短路。三相短路也称为对称短路，系统各相与正常运行时一样仍处于对称状态。其他类型的短路都是不对称短路。电力系统的运行经验表明，在各种类型的短路中，单相短路占大多数。两相短路较少，三相短路的机会最少。三相短路虽然很少发生，但情况较严重，应给以足够的重视。况且，从短路计算方法来看，一切不对称短路的计算，在采用对称分量法后，都归结为对称短路的计算。因此，对三相短路的研究是有其重要意义的。3.2 三相对称短路计算假定系统中的节点f经过渡阻抗Zf发生短路。这个过渡阻抗Zf均不参与形成网络的节点导纳(或阻抗)短阵。图3-1中方框内的有源网络代表系统正常状态的等值网络。现在我们保持故障处的边界条件不变，把网络的原有部分同故障支路分开(见图3-1)。容易看出，对于正常状态的网络而言，发生短路相当于在故障节点f增加了一个注入电流-If(短路电流以流出故障点为正，节点电流则以注入为正)。因此，网络中任一节点i的电压可表示为： （3-1）式中，G为网络内有源节点的集合。 图3-1 对称短路分析由式(3-1)可见，任一节点i的电压都由两项叠加而成。第一项是符号下的总和，它表示当If=0时由网络内所有电源在节点i产生的电压，也就是短路前瞬间正常运行状态下的节点电压，这是节点电压的正常分量，记为 Vi(0)。第二项是当网络中所有电流源都断开，电势源都短接时，仅仅由短路电流If在节点i产生的电压，这就是节点电压的故障分量。上述两个分量的登加，就等于发生短路后节点i的实际电压，即 (3-2)公式(3-2)也适用于故障节点f，于是有, （3-3）方程式(3-3)也可以根据戴维南定理直接写出，含有两个未知量Vf和If，需要根据故障点的边界条件再写出一个方程才能求解。这个条件是 （3-4）由方程式(3-3)和（3-4）可解出: (3-5)如果网络在正常状态下的节点电压为已知，为了进行短路计算，只须利用节点阻抗矩阵中与故障点f对应的一列元素。因此，尽管是采用了阻抗型的节点方程，但是并不需要作出全部阻抗矩阵。在短路的实际计算中，一般只需形成网络的节点导纳矩阵，并根据具体要求，用第四章所讲的方法求出阻抗矩阵的某一列或某几列元素即可。在应用节点阻抗矩阵进行短路计算时，我们都将采用这种算法。在不要求精确计算的场合，可以不计负荷电流的影响。在形成节点导纳矩阵时，所有节点的负荷都略去不计，短路前网络处于空载状态，各节点电压的正常分量的标么值都取作等于l。图3-2给出了对称短路简化计算的原理框图。图3-2 对称短路计算原理框图3.3 导纳矩阵的修改在现代电力系统分析中，往往需要研究不同接线方式情况下的运行状态，例如某台变压器或某条输电线路的投入或切除，对某些元件的参数进行修改等。由于改变一条支路的开合状态只影响该支路两端节点的自导纳及其互导纳，因此在这种情况下不必重新形成导纳矩阵，仅仅需要在原有导纳矩阵的基础上进行必要的修改就可以得到所要求的导纳矩阵。下面分几种情况进行介绍。(1)从原有网络引出一条新的支路，同时增加一个新的节点，见图3-3(a) (a) (b) (c) (d) 图3-3 电力网络接线变更示意图设i为原有网络N中任意一节点，j为新增加节点，Zij为新增加的支路阻抗。由于增加了一个新的节点，因而导纳矩阵相应增加一阶。因为j点只有一条支路，所以并增加非对角元素在这种情况下，原有网络i节点的自导纳应有如下的增量: (2)在原有节点i和j间增加一条支路，见图3-3(b)。在这种情况下，虽增加了支路，但并不增加节点数，导纳矩阵的阶数不变。但是，与节点i,j有关的元素应作以下修正: （3-6）(3)在原有网络节点i和j间切除一条阻抗Zij为的支路。 在这种情况下，相当于在节点i和j间追加一条阻抗为一Zij的支路，见图3-3(c)。因此导纳矩阵有关元素应作以下修改: （3-7）(4)原有网络节点i和j之间支路阻抗由Zij改变为Zij。 在这种情况下，可以看作首先在节点i和j间切除阻抗为Zij的支路，然后再在节点i和j间追加阻抗为Zij的支路，如图3-3(d)所示。根据式(3-6)及式(3-7),求出在此情况下导纳矩阵有关元素的修正量。3.4 快速切除故障对暂态稳定性的影响快速切除故障对于提高系统的暂态稳定性有着决定性的作用，因为快速切除故障减小了加速面积，增加了减速面积，提高了发电机之间并列运行的稳定性。另一方面，快速切除故障也可使负荷中的电动机端电压迅速回升，减小了电动机失速的危险。第四章 简单模型下的暂态稳定分析在电力系统暂态稳定分析的过程中，我们首先要通过潮流计算来确定故障或操作发生前一时刻系统的状态，即各个变量的初始值，然后再根据这一系列的初始值，运用直接法求解代数方程来确定故障或操作发生以及以后的各个时刻各个发电机的电磁功率，然后运用基于显式积分的龙格-库塔法来求解各发电机的转子角度，描绘其摇摆曲线，从而判断系统的运行状态是否稳定。本章的主要内容就是具体的对简单模型下的暂态稳定分析过程进行说明，主要包括初值计算、微分方程的求解、转子摇摆曲线的绘制以及临界切除时间的确定。在本次毕业设计中采用的模型和计算方法如下：发电机：用发电机暂态电势保持恒定来模拟。负荷：负荷采用恒定阻抗来模拟。电力网络：用导纳矩阵描述。微分方程：用四阶龙格-库塔法进行求解。网络方程：用直接法求解。4.1 初值计算在进行暂态计算前，首先应对故障前正常运行的电力系统进行潮流计算确定微分方程求解所需的初值。本文采用了牛顿-拉夫逊法潮流计算。牛顿-拉夫逊法（以下简称牛顿法）作为解非线性方程式非常有效的方法，被广泛用于电力系统潮流计算中。其核心部分就是将非线性方程的求解过程变成反复对相应的线性方程的求解过程，通常称为逐次线性化过程。 图 4-1 牛顿法的几何解释下面以单个变量为例具体说明一下牛顿法的计算过程。试求非线性方程式： （4-1）的解。设为该方程式的初值，真正解在它附近： （4-2）式中为初值的修正量。如果求得，则由式（4-2）就可得到真正解，将式 按泰劳级数展开： （4-3）很小时，式（4-3）可以简化为： （4-4）这是对于变量的线性方程式，即修正方程式，用它可以求出修正量，用对修正以后得到的： （4-5）再以作为初值得 （4-6）解式（4-6）得更趋近于真正解的： （4-7）如此反复的计算，当第t次迭代时的修正方程式为： （4-8）当时，此式满足原式（4-1），因而成为该方程的解。牛顿法不仅用于求解单变量方程，它也是求解多变量非线性代数方程的有效方法。利用牛顿法进行潮流计算的一般过程是这样的：(1)给定各节点电压初值，；(2)将电压初值，代入功率方程，求修正方程式的常数项，；(3)利用电压初值，求雅可比矩阵中各元素；(4)解修正方程式，求修正量，；(5)修正各节点电压向量： （4-9）(6)以，代入式功率方程中求，；(7)校验是否收敛；如收敛则输出结果，否则再以，为初值转回第三步进行下一次迭代。如前所述，牛顿法求解潮流问题的过程，实际上是不断形成并求解修正方程的过程。 首先，我们对发电机的初值进行求解。由潮流计算可得到干扰前各发电机的端电压和各发电机注入网络的功率，然后就可以计算出发电机注入网络的电流： （4-10）这样，就可以求出虚构电势，即 （4-11）依次就可以确定发电机转子角度的初值 （4-12）在系统稳态运行时，发电机转子以同步转速旋转，于是有 （4-13）利用坐标变换可以求出发电机定子电压和电流的d、q分量 (4-14）然后再算出暂态电势的值 （4-15）另外，稳态运行时发电机的电磁功率等于原动机的机械功率，即有 （4-16） 负荷初值的计算比较简单。由潮流计算结果可知干扰前各负荷节点电压和负荷所吸收的功率，这样，我们可以得到负荷的等值导纳 （4-17）由于它在整个暂态过程中保持恒定，所以，可以将它包括在网络导纳矩阵里。4.2 用直接法求解网络方程 在这种方法中，网络方程采用实数形式描述，如式（4-18）。其导纳矩阵包括了恒定阻抗模拟的负荷，它在整个暂态过程中是保持恒定的。 （4-18） 设发电机接在网络中的节点，再简单模型中，得到发电机节点注入电流的表达式 （4-19）其中的元素： （4-20）这样，我们得到的新的网络方程只是对原网络方程的简单修改：导纳矩阵的相应对角块发生变化，电流向量除了发电机节点有虚拟注入电流，其余节点的注入电流为零，即： 导纳矩阵的第个对角块变为 （4-21）发电机节点的虚拟注入电流为： （4-22）这样，在每个积分步得到的线性方程足可以用Matlab直接进行求解，从而求出此时刻网络各节点电压的实部和虚部、。在得到发电机节点的电压后，即可按式（4-19）算出发电机节点的注入电流、。4.3 用龙格-库塔法求解微分方程4.3.1龙格-库塔法的基本原理龙格-库塔法近似于泰勒级数求解法，然而不同于正规的泰勒级数求解，龙格-库塔法并不需要显式计算高于一阶的微分。高阶微分的作用包括在n次一阶微分的求值中。取决于泰勒级数保留的有效项的数目，因而有不同阶数的龙格-库塔法。本次毕业设计采用的是四阶龙格-库塔法。这种方法用，区间四个点的导数去推算，从而拟合了泰劳级数的前五项： （4-23）因此，它的局部截断误差是阶的，其全局截断误差是阶的。而改进欧拉法用，区间两点的导数（或斜率）推算，拟合了积分函数泰勒级数的前三项，局部阶段误差达到了阶。对于一阶微分方程式： （4-24）给定初值为、。当利用四阶龙格库塔法求解时，可按下式进行计算： （4-25）其中： （4-26）用同样的方法可以求出， 。一般四阶龙格库塔法的递推公式为： （4-27）其中： （4-28）应用龙格库塔法也可以求解一阶微分方程组。例如： （4-29）递推公式可进行如下计算： （4-30）其中： （4-31） （4-32）多个方程多个变量的一阶微分方程组，可以仿照上面的方法得到递推的公式。4.3.2利用四阶龙格-库塔法求解转子微分方程 在本次毕业设计中，我们利用了四阶龙格-库塔法求解了转子的运动微分方程。在简单模型下的暂态稳定分析中，转子运动方程如下： （4-33）假设电力系统的机电暂态过程已经计算到时刻，现在计算时刻系统运行状态的过程。首先，我们要对此时刻系统是否发生操作或故障进行判断，若无故障或操作发生，则直接以时刻系统的状态作为初值求解微分方程；否则需要计算出故障后电力系统的运行参数，用它和时刻的状态变量作为初值求解微分方程。龙格-库塔法求解微分方程的步骤如下： （1）由时刻各发电机的计算出系统所有节点的电压、和发电机节点的注入电流、。 （2）根据龙格-库塔法的递推公式求出、。然后修正、，即算出，由于节点电压和发电机节点的注入电流都为的函数，所以要修正网络方程的导纳矩阵，重新求解网络方程，得出节点的电压、和发电机节点的注入电流、，计算出电磁功率，求出。其中，电磁功率的按下式计算： （4-34） （3）重复第二步，求出、。最后，根据式（4-30）计算出时刻的发电机转子角度以及转速。通过以上的三个步骤，我们就可以求出各时刻发电机的转子角度，然后就可以绘制出发电机转子的相对摇摆曲线，确定故障临界的切除时间。 四阶龙格-库塔法具有计算精度高，数值稳定好，模型能力强的优点，但是它的计算量比较大。第五章 程序设计与仿真测试电力系统暂态稳定计算也是电力系统不正常运行方式的一种计算。它的任务是已知电力系统某一正常运行状态和受到某种扰动，计算电力系统中所有发电机能否保持同步运行。电力系统受到扰动后是一个复杂的机电暂态过程。但是不管怎样复杂，电力系统各元件的暂态过程都可以用一组微分方程和代数方程描述。电力系统暂态稳定计算就是在一定的初值条件（正常运行状态和扰动方式）下，联合求解这些微分方程和代数方程，从而得到发电机转子位置角的摇摆曲线，来判别系统在这种运行状态和扰动方式下能否保持稳定。5.1 程序设计首先，采用牛顿-拉夫逊法进行潮流计算，得出电力系统在正常运行方式下的电压、原动机功率以及发电机的转子相位角，即暂态稳定计算的初始值。然后，根据所选故障类型，修改导纳矩阵，并将导纳矩阵写成实数形式，用龙格-库塔法对微分方程进行求解，求解过程中，每一个积分步都要进行代数方程和微分方程的交替求解，得出发电机的转子角。最后，根据发电机的相对转子角度画出摇摆曲线，对摇摆曲线进行分析，判别系统的暂态稳定性，找出系统故障的临界切除时间。程序的具体实现过程如图5-1所示。图中，框至为准备工作，其中进行牛顿拉夫逊法的潮流计算。在框中，用潮流计算出的结果计算发电机的内电势的幅值和相角；同时把负荷换算成等值固定阻抗，并把这个阻抗归并到潮流的导纳矩阵中去，一起形成新的暂态计算用的导纳矩阵。框中，把初始参数值传入相应的变量，把计算时间初始化为0时刻。中判断是否有操作或者线路故障，如果有，就用框的程序修改导纳矩阵。、框为龙格库塔的核心部分，这部分里将对给定的时刻的值计算的值。框中判定给定的仿真时间是否已到，一般给定5秒。如果仿真时间到了，就描绘出发电机相对角度的摇摆曲线。从曲线的摇摆规律我们就可以判断系统的稳定性了。初始化输入信息及原始数据用潮流计算扰动前的正常运行状态计算系统初值，形成导纳矩阵 有操作或故障修改导纳矩阵 用龙格库塔法求解转子运动方程，得时刻的和仿真总时间到？画摇摆曲线停机有无是否图51 基于龙格库塔法的暂态稳定计算程序流程图5.2 仿真测试5.2.1基于9节点标准系统的仿真测试（一）假设干扰在线路5-7靠近母线7处发生三相接地短路，故障由断开线路5-7而被消除。（1）步长为0.001s，总时段为2000，第70时段切除故障时发电机转子相对角度摇摆曲线如下图所示：图5-2 在第70时段切除故障时发电机转子摇摆曲线由图（5-2）可知，当系统在70*0.001=0.07s时刻切除故障时，系统是暂态稳定的。（2）此时，可选较大的时段作为切除时段来绘制摇摆曲线，在第90时段切除故障时发电机转子相对角度摇摆曲线如下图所示： 图5-3 在第90时段切除故障时发电机转子摇摆曲线由图（5-3）可知，当系统在90*0.001=0.09s时刻切断故障时，系统是暂态不稳定的。（3）在第70时段和第90时段的发电机转子相对摇摆曲线可知，系统暂态稳定临界时间在0.07s到0.09s之间，在第80时段切除故障时，发电机转子相对角度摇摆曲线如下图：图5-4 在第80时段切除故障时发电机转子摇摆曲线由图（5-4）可知，系统在80*0.001=0.08s时还是暂态稳定的，所以系统的临界时间在0.08s和0.09s之间，取0.085s的发电机转子相对角度摇摆曲线如下图所示：图5-5 在第85时段切除故障时发电机转子摇摆曲线与以上分析同理，此时系统是暂态不稳定的，取0.084s做出发电机转子相对角度摇摆曲线如下图所示： 图5-6 在第84时段切除故障时发电机转子摇摆曲线由图可知，系统仍然是暂态不稳定的，所以继续减小切除时刻，取0.083s,作出发电机相对角度摇摆曲线如下图所示： 图5-7 在第83时段切除故障时发电机转子摇摆曲线由图可知，系统在此时是暂态稳定的，通过上述的分析也可以得到，0.083s是系统的暂态稳定临界时间。即在0.083s之前切除故障，系统都是暂态稳定的。能够在暂态