电力系统机电暂态仿真程序

发电机(经典模型)、线路、负荷;
2)故障只处理对称故障;
3)采用改进欧拉法;

第 1 章 绪论

1.1背景介绍

暂态稳定性分析计算是电力系统三大经典计算之一，它是电力系统规划、设计、调度运行与控制中必不可少的一项计算任务。

1.1.1 电力系统暂态分析的基本概念

在正常的稳态运行修况下，电力系统中的各发电机组输出的电磁转矩和原动机输入的机械转矩平衡，因而所有发电机转子速度保持恒定。但是电力系统不可避免地会受到一些大干扰的冲击，例如，各类短路故障、大容量机组或大负荷的投入和切除等。在遭受大的扰动后电力系统除了经历电磁暂态过程以外，还经历机电暂态过程。事实上，由于系统的结构发生了较大的变化，系统的潮流及各发电机的输出电磁功率也随之发生变化，以而破坏了原动机和发电机之同的功率平衡，在发电机转轴上产生不平衡转矩，从而导致发电机转子加速或减速。通常情况下，大扰动后各发电机组的功率不平衡状况基本上各不相同，加之各发电机转子的转动惯量也有所不同，使得各机组转速变化的情况各不相同。这样，发电机转子之间将产生相对运动，使得转子之间的相对功角发生变化，而转子之间相对功角的变化又反过来影响各发电机的输出功率，从而使各个发电机的功率、转速和转子之间的相对功角继续发生变化。
与此同时，由于发电机端电压和定子电流的变化，引起励磁调节系统的调节过程；由于机组转速的变化，引起调速系统的调节过程；由于电力网络中母线电压的变化，引起负荷功率的交化；网络潮流的变化引起一些其他控制装置(如SVC、TCSC、直流系统中的换流器等）的调节过程等。所有这些变化都直接或间接地影响发电机转轴上的功率平衡状况。
以上各种变化过程互相影响，形成了一个以各发电机转子机械运动和电磁功率变化为主体的机电暂态过程。
电力系统遭受大干扰后所发生的机电暂态过程可能有两个不同的结局。一种结局是各发电机转子之间的相对功角随时间的变化是摇摆(或振荡)状态，且振荡幅值逐渐衰减，各发电机之间的相对运动逐渐平息，从而系统过渡到一个新的稳态运行情况，各发电机仍然保持同步运行，这种情况，称电力系统是暂态稳定的。另一种结局是在暂态过程中某些发电机转子之间始终存在着相对运动，使得转子间的相对功角随时间不断增大，最终导致这些发电机失去同步，这时，称电力系统是暂态不稳定的。
电力系统正常运行的必要条件是所有发电机保持同步，因此，电力系统在大扰动下的稳定性分析，就是分析遭受大干扰后系统中各发电机维持同步运行的能力，常称为电力系统的暂态稳定性分析。需要说明的是电力系统暂态稳定性分析一般只涉及系统在受扰后短时间内(最多10s以内)的暂态行为。有关力系统受扰后的中长期动态行为的分析，是电力系统中、长期动态稳定性分析的内容。从概念上讲，暂态稳定性分析一般不考虑电力系统中的慢过程，例如，锅炉的动态过程、水轮机进水闸门及水管的动态过程、自动发电控制(AGC)等。此外，暂态稳定性分析也不考虑系统受扰后变化极为快速的电磁暂态过程。在暂态稳定性分析中，通常认为电磁暂态过程已经结束，电磁暂态变量已衰减完毕，而中、长期动态过程还没有开始变化。
从理论上讲，电力系统的动态过程是一个连续的过程，上述所谓的电磁暂态过程、机电暂态过程以及中、长期动态过程并不是截然分开的。但是，将电力系统受扰后的动态过程依据其响应速度或其变化的快慢程度划分为3个不同的阶段，这从概念上来讲是清晰、合理的，更有利于对复杂问题的深入理解。
电力系统遭受大干扰是人们所不希望的，但事实上却是无法避免的。系统在遭受大干扰后失去稳定的后果往往非常严重，有时候甚至是灾难性的。事实上电力系统遭受到的各种大干扰，是以一定的概率随机地发生的。因此，系统的设计运行方式的制定总是需要保证系统在合理选择的预想事故下能够保持稳定，但不能要求电力系统能经受所有干扰的冲击。由于各国对系统稳定性的要求不同，所以对预想事故的选择也就有不同的标准。通常情况下，人们希望电力系统能够经受住任一发电机机端发生短时间(故障快速切除时间)短路故障如三相短路故障的冲击。
判断电力系统在预想事故下能否稳定、同步运行，需要对电力系统进行暂态稳定性分析。当系统不稳定时，还需要研究提高电力系统暂态稳定性的措施。迄今为止，有关电力系统暂态稳定性的分析，主要基于两类方法：一类是直接法(direct methods)或称暂态能量函数(transient energy function，TEF)法；另一类是基于状态空间的时域数值计算方法。前者主要是基于动力系统的Lyapunov稳定性定理，其核心思想是构造一个类似于“能量”的暂态能量函数，即Lyapunov函数，可以在不求解微分方程组的情况下，通过检查该函数沿系统积分轨迹随时间变化的性质来判断系统的暂态稳定性。然而，对实际电力系统而言，其暂态能量函数的构造至今尚缺乏系统性的通用方法。因此，直接法的应用目前仅限于简单电力系统。后者又称数值仿真方法，其基本思想是利用解初值问题的数值方法来求解描述电力系统动态特性的微分方程组，并通过发电机组之间的转子功角差来判断系统是否稳定。这种方法直观、精确，易于处理各种复杂模型和复杂故障的影响。

1.1.2 微分代数方程

工程领域中绝大多数随时回演变的现象或过程，主要是采用各种复杂的常微分方程描述。一些典型的常微分方程（线性方程或某些特殊的非线性方程）等可以运用基本方法求出其解析解。但是，大多数常微分非线性方程很难或不可能获得其封闭形式的解析解。从实际应用的角度看，—般只要求得到解在若干个时间点上的近似值或者解的易于计算的近似表达式，这就是所谓的常微分方程的数值解。
高阶常微分方程或方程组通常可化为一阶方程组来研究，因此，主要介绍一阶常微分非线性方程的初值问题。
有关求解常微分方程数值解的方法即常微分方程数值方法的研究，是数值分析计算领域中较古老的一个分支。数值方法的构造首先依赖于方程本身的性质，因为对一个没有解的方程设计其数值方法是没有意义的。
常微分方程组的求解方法主要有线性多步法以及RK方法两大类，但也可以分为显式数值积分方法和隐式数值积分方法两大类。据此分类方法，电力系统暂态稳定性的数值计算，大致可以分为两类：一类是所谓的显式分离求解(partitioned explicit method，PE)；另一类是隐式联立求解法(simultaneous implicit method，SI)。前一类方法的主要特征：采用显式数值积分方法求解微分方程组；微分方程组与代数方程组两者的求解彼此是独立但需要交替进行的。后一类方法的主要特征：利用隐式数值积分方法将微分方程组离散化为非线，然后利用Newton法对差分代数方程性代数方程组(以下简称差分代数方程组)组和网络代数方程组进行联立求解，即微分方程组与代数方程组两者是联立求解的。
迄今为止，在电力系统暂态稳定性计算中，显式RK方法和隐式梯形方法是应用最为普遍的数值积分方法。
显式RK方法，如4级4阶显式RK方该，具有计算过程简单、计算精度较高的优点，电力系统暂态稳定性计算中PE类方法绝大多数采用显式RK方法。
隐式梯形积分方法是2阶、A-稳定的数值方法，它具有很好的数值稳定性，特别适合于刚性问题，因而在电力系统暂态稳定性计算中得到了最为广泛的应用，SI类方法主要采用这一方法。自Dommel和Sato于1972年首次将隐式形积分规则应用于电力系统暂态过程的数值仿真计算后，隐式梯形积分规则在电力系统暂态稳定性分析计算中一直占据主导地位。研究人员普遍认为：隐式梯形积分规则是A-稳定的，具有很好的数值稳定性；它是2阶的单步方法，其计算过程比较简单，而且对电力系统全动态过程——从电磁暂态到长期动态稳定性的分析计算均比较适用。从1972年起，无论在学术界还是在电力系统工程实际应用领域，隐式梯形积分规则仍然是应用最为普遍的数值积分方法。因此，隐式梯形积分规则以其2阶、A-稳定性和对刚性问题的适用性，在电力系统暂态稳定性计算中持续展现出其重要性和广泛应用。自其首次被引入至今，这一方法一直是学术界和电力系统工程领域中最为常用的数值积分方法，为电力系统暂态过程的精确分析和计算提供了可靠的基础和支持。
微分-代数方程在电力系统暂态稳定性分析中扮演着重要的角色。其作为描述系统动态行为的数学模型，为我们提供了理解和解决电力系统暂态过程中复杂动态行为的关键工具。通过对这些方程的求解和分析，我们能够深入了解系统的响应特性、稳定性以及在各种异常情况下的行为。这些研究不仅推动了电力系统稳定性分析的发展，也为提高电力系统的稳定性和可靠性提供了重要的理论基础和方法支持。因此，微分-代数方程的研究不断推动着电力系统暂态稳定性分析领域的进步与发展。

1.2 假设条件

暂态稳定性分析计算的主要目的在于确定系统在给定的大扰动下发电机能否继续保持同步运行。因此，只需研究表征发电机是否同步的转子运动特性，即功角随时间变化特性即可。据此，找出暂态过程中对转子机械运动起主要影响的因素，在分析计算中加以考虑，而对于影响不大的因素，则可以忽略或作近似考虑。具体的近似可以概述如下：
1.只考虑平衡的三相系统和平衡的扰动。因此，只用到正序网络。
2.不考虑频率变化对系统参数的影响。在较短时间内的机电暂态过程中，发电机频率偏离同步频率50Hz很小，忽略直流偏移量和谐波。因此输电线路、变压器和阻抗负荷构成的网络基本上是处于稳态的；电压、电流和功率可以通过代数潮流方程计算得出。
3.忽略发电机定子电流的非周期分量及相对应的转子电流的周期分量。在大扰动特别是发生短路故障时，定子非周期分量电流将在定子回路电阻中产生有功损耗，增加发电机转轴上的电磁功率，从而使发电机产生减速运动。然而，一方面由于定子非周期分量电流衰减时间常数很小，另一方面，定子非周期分量电流产生的磁场在空间上是静止不动的，它与转子绕组直流包括自由电流所产生的转矩以同步频率作周期变化，其平均值很小，由于转子机械惯性较大，所以对转子整体相对运动影响很小。采用这个近似之后，发电机定、转子绕组的电流、电压以及发电机的电磁功率等，在大扰动的瞬间均可以发生突变。同时，这一近似也意味着同时忽略发电机定子以及电力网络中各元件的电磁暂态过程，从而使得发电机定子电压方程和网络方程均为代数方程，且仅包含基频电气分量，因而可以用稳态关系式来描述。

1.3 本文主要工作

本文围绕电力系统数字仿真中的用户自定义建模技术和发电机建模中转速的简化处理对暂态稳定计算的影响进行了研究，主要工作如下：
1.强调了电力系统暂态潮流计算的重要意义；讨论了电力系统暂态分析的基本概念、微分代数方程的发展和求解；介绍了时域仿真法、直接法和机器学习法三种稳定性分析的方法。
2.对后续建模和求解过程中的必要环节和设备进行假设，以使系统处理和操作更加完善，考虑更加全面。
3.介绍电力系统设备的数学模型。对发电机转子运动方程和电压电流方程进行阐释；对考虑不同因素的负荷模型进行模型建立和适用条件的分析；对线路网络模型进行简单介绍；将全电力系统设备进行数学建模，分别介绍考虑保留网络结构的经典模型和基于网络节点收缩的经典模型。
4.基于改进欧拉法对微分-代数方程进行求解，从而实现电力系统进行暂态仿真计算。对建立的数学模型进行整合分析，并介绍数值解法的一般过程；针对数值计算的初值计算、故障/操作处理和基于改进欧拉法的交替迭代计算三部分，进行原理说明、代码编写和过程讲解；最后对主循环和结果输出进行代码阐释，并绘制流程图进行过程说明。
5.应用IEEE14节点系统进行算例仿真与分析。对IEEE14节点系统进行简要介绍；对故障前的稳态进行简要分析，观察各发电机转子角度和角速度的增量；发生三相短路故障后，对各节点故障时最大的功角差和角速度进行统计，对各发电机功角差和各节点电压进行时域仿真分析；对故障切除时间进行简单调研，设故障切除时间分别为0.1s、0.12s和0.2s，对故障切除后的各发电机功角差和各节点电压进行时域仿真分析；应用等面积法则，并确定系统的极限切除角；对暂态稳定性进行分析，并提出提高暂态稳定性的措施，对提高电力系统稳定性具有重要意义。

第 2 章 电力系统设备数学模型

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

2.4 全系统数学模型

图 2 2 暂态分析计算中全系统数学模型框架图
由上图可以看出，电力系统暂态稳定性计算主要包括发电机、与发电机相关的励磁系统、原动机及调速器、负荷、其他调节装置等动态环节的数学模型以及电力网络的数学模型。显然，电力系统不同发电机的动态模型是互相独立的，通过电力网络联系到一起。在同一台发电机中，不同动态环节之间只通过特定的变量发生联系。每个动态环节都有其独特的数学描述方式，通常以微分方程组的形式呈现。这些方程组通过相应的变量联系在一起，形成了一个复杂的动态系统模型。
整个系统的数学模型可用以下形式的微分-代数方程组来描述。
微分方程组主要包括：(1) 描述各同步发电机暂态和次暂态电势变化规律的微分方程。(2) 描述各同步发电机转子运动的摇摆方程。(3) 描述同步发电机组中励磁调节系统动态特性的微分方程。(4) 描述同步发电机组中原动机及其调速系统动态特性的微分方程。(5) 描述各感应电动机和同步电动机负荷动态特性的微分方程。(6) 描述直流系统整流器和逆变器控制行为的微分方程。(7) 描述其他动态装置。
代数方程组主要包括：(1) 电力网络方程，即描述在公共参考坐标系x-y下节点电压与节点注入电流之间关系。(2) 各同步发电机定子电压方程及d-q坐标系与x-y坐标系间联系的坐标变换方程。(3) 各直流线路的电压方程。(4) 负荷的电压静态特性方程等。在具体的暂态稳定性分析计算中，可根据研究问题的性质、对计算结果精度要求不同，依据抓住重点、忽略次要因素的原则采用相应复杂程度的元件数学模型。
因此，本文在暂态计算的通用模型基础上，进一步采用以下简化计算条件：
(1) 对所有的发电机均采用经典模型，即将同步发电机看作带内阻抗( )、恒电动势( )的简单电源。
(2) 不考虑各类调节器及控制装置的作用。
(3) 所有的负荷均为恒阻抗负荷。

第 3 章 电力系统机电暂态计算

电力系统暂态稳定性数值计算，即在时域内求解微分-代数方程组。其中，微分方程组的求解是暂态稳定性计算的核心。
3.1 暂态稳定数值解法
3.1.1 系统模型
电力系统暂态稳定分析的数值解法是用数值求解方法分析求解电力系统全系统的数学模型(方程)。全系统数学模型(方程)是由电力系统各元件模型根据元件间拓扑关系形成的一组微分方程和代数方程。
式中，x表示电力系统有关元件内部的状态参数，y表示电力网络的运行参数。
电力系统基本上是由发电机、励磁系统、原动机及调速器以及网络和负荷组成的。其中发电机分为两部分，即转子运动方程部分和电磁回路方程部分。转子运动方程反映了当发电机输入机械功率 和输出电功率 不平衡时引起发电机转速 和转子角 的变化。发电机转速信号送人调速系统和参考速度比较，其偏差作为调速器的控制输人量，以控制原动机的输出机械功率 。而发电机的输出电磁功率 将影响转子运动的功率平衡及转子速度和角度的变化。网络一般表示为节点导纳矩阵形式，网络除和发电机相连外，还和负荷相连。实际的电网有许多发电机和负荷，通过网络互相联系和互相影响，造成了电力系统暂态稳定分析的复杂性。
暂态稳定分析由于主要研究发电机转子摇摆特性，主要和网络中的工频分量有关，故发电机可忽略定子暂态而采用实用模型，而网络采用准稳态模型，负荷则采用静态模型。为了突出电力系统暂态稳定分析的基本原理和步骤，对发电机采用经典二阶模型，忽略凸极效应，并设暂态电动势 幅值恒定，从而忽略励磁系统的动态，以简化分析。应当指出， 恒定已计及了励磁系统的一定作用，即认为励磁系统足够强，从而能保证 后的暂态电动势 恒定。另外，忽略调速器和原动机动态作用，即认为机械功率 为定常值。

3.1.2 数值解法一般过程

电力系统暂态稳定计算的基本流程如下图所示。为与后文章节相对应，在下面的流程图中，按照循环前的准备、循环的第一个分支、数值解法的核心和其他部分分为四节内容，分别对应于后续的3.2~3.5节。

（1）潮流。在计算电力系统暂态稳定之前，首先应利用潮流计算程序算出扰动前系统的运行状态，求出电力网络的运行参数 。
（2）初值。然后根据潮流计算结果，计算状态变量的初值 ，及 和 的稳态值，采用简化模型时 和 在暂态过程中保持不变。
（3）导纳矩阵。根据网络元件参数和网络拓扑关系形成网络稳态工况下的节点导纳矩阵，也可从潮流输出中读入。将负荷等值导纳 及发电机内部暂态导纳并入导纳矩阵。
（4）开始。将时钟指针 置零。
（5）扰动/操作判断。根据扰动过程参数，判断当前有无扰动发生。
（6）修改故障方程。若有扰动则需要根据扰动参数修改导纳矩阵及微分方程。
（7）计算扰动后参数。设发生故障当前时刻的状态参数不发生突变，据扰动后系统代数方程计算下一时刻的代数量作为下一步的初值。
（8）求解。微分方程和代数方程的求解是暂态稳定数值解法的核心，采用改进欧拉法和交替求解法计算 时刻的 和 。
（9）判稳。判断系统是否稳定。
（10）时间累加。若系统保持稳定，则令 。
（11）循环。若未达到仿真总时间，就进行下一时刻的计算。
（12）结束输出。若系统失稳或达到仿真总时长，则计算结束，输出计算结果。

3.2 初值计算及方程形成

3.2.1 潮流计算
采用“pypower”进行稳态潮流计算。pypower是一个基于Python的开源工具包，用于电力系统的仿真、分析和优化。它建立在MATPOWER的基础上，提供了一系列用于电力系统分析的函数和工具。该工具包允许用户进行电力系统的各种计算和分析，包括潮流计算、稳定性分析、优化等。且支持直接处理潮流计算、最小负荷分布、最优潮流等问题，也能够进行暂态稳定性分析、电压稳定性评估等方面的计算。PyPower提供了一系列用于处理电力系统数据、建立模型、执行计算的函数和工具，使用户能够轻松地进行电力系统的建模和仿真。同时，作为一个开源工具包，PyPower也允许用户根据自己的需求进行定制和扩展。具体代码如下：

from pypower.api import case14， ppoption， runopf， makeYbus 
import numpy as np
def powerflowstable():
    data = case14()
    ppopt = ppoption()
    print(ppopt)
    result = runopf(data， ppopt)
    bus0 = result['bus']
    data['branch'][:，0]=data['branch'][:，0]-1
    data['branch'][:，1] = data['branch'][:，1] - 1
    [Ybus， Yf， Yt] = makeYbus(data['baseMVA']， bus=data['bus']， branch=data['branch'])
    V = np.array(bus0[:，7]  * np.exp(1j * np.radians(bus0[:，8])))
    S = np.array((bus0[:，2] + 1j * bus0[:，3])/data['baseMVA'])
Y = Ybus.todense()  
    Y = np.array(Y)
    return V，S，Y

通过“powerflowstable”函数，返回稳态电压、复功率和节点导纳矩阵。
其部分运算结果如下图

图 3 3 runopf运行结果图