CN102609575A - 一种基于隐式数值积分的电力系统暂态稳定仿真方法 - Google Patents

一种基于隐式数值积分的电力系统暂态稳定仿真方法 Download PDF

Info

Publication number
CN102609575A
CN102609575A CN2012100175933A CN201210017593A CN102609575A CN 102609575 A CN102609575 A CN 102609575A CN 2012100175933 A CN2012100175933 A CN 2012100175933A CN 201210017593 A CN201210017593 A CN 201210017593A CN 102609575 A CN102609575 A CN 102609575A
Authority
CN
China
Prior art keywords
integration
generator
implicit
transient stability
transient
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN2012100175933A
Other languages
English (en)
Other versions
CN102609575B (zh
Inventor
王建全
李梦骄
赵志奇
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Zhejiang University ZJU
Original Assignee
Zhejiang University ZJU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Zhejiang University ZJU filed Critical Zhejiang University ZJU
Priority to CN2012100175933A priority Critical patent/CN102609575B/zh
Publication of CN102609575A publication Critical patent/CN102609575A/zh
Application granted granted Critical
Publication of CN102609575B publication Critical patent/CN102609575B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
    • Y02E60/00—Enabling technologies; Technologies with a potential or indirect contribution to GHG emissions mitigation

Abstract

本发明公开了一种基于隐式数值积分的电力系统暂态稳定仿真方法。与已有的电力系统暂态稳定数值仿真隐式梯形积分方法相比，采用了局部截断误差更小的功角积分公式，即将描述电力系统暂态过程的非线性微分方程组，表示为线性部分和非线性部分。通过合理将线性部分的系统矩阵选为奇异阵，以得到状态转移矩阵的准确解析表达式，并用线性可积函数去近似微分方程组的非线性部分，以得到一组隐式积分公式。其中，发电机功角隐式积分公式的局部截断误差为，高于局部截断误差为
Figure 31330DEST_PATH_IMAGE004
的隐式梯形积分，而每一次积分的计算量则与隐式梯形积分相当。由于采用了高精度的隐式积分公式，使得在相同迭代精度下每一积分步迭代次数减少，从而显著地减少了仿真的计算量。

Description

一种基于隐式数值积分的电力系统暂态稳定仿真方法
技术领域
本发明属于电力系统自动化专业领域，涉及到了一种用于计算电力系统暂态稳定的数值积分方法。
背景技术
电力系统暂态稳定分析是电力系统分析计算中最核心、最基础的内容之一。由于现代电力系统的规模不断扩大，在线动态安全分析，安全稳定紧急控制、预防控制，智能调度等控制技术已逐步在电力系统中推广使用。实现这些先进技术的前提条件是能够对大规模电力系统进行快速准确可靠的暂态稳定仿真计算。
用于电力系统暂态稳定计算的分析方法主要有数值积分法，直接求解法，以及将数值积分和直接求解法相结合的混合分析方法。其中，数值积分法是电力系统暂态稳定计算方法中最准确、最可靠的方法。数值积分法最大的缺点是计算量大，尽管计算机计算速度已经有了飞速提高，但对于大规模电力系统，计算的速度难以满足在线动态安全分析、预防控制、紧急控制的要求。
电力系统的暂态过程可用如下形式的微分-代数方程组描述
Figure 619345DEST_PATH_IMAGE001
                         （1）               
Figure 946421DEST_PATH_IMAGE002
                          （2） 
式中，
Figure 40891DEST_PATH_IMAGE003
表示微分方程组中描述系统动态特性的状态变量；表示代数方程组中系统的运行变量。通常向量
Figure 278155DEST_PATH_IMAGE003
包含发电机功角和转速等描述系统中各动态环节的状态变量，而向量通常包含与网络相关的运行变量，如节点电压的幅值和相位等。
用数值积分法求解电力系统暂态过程的一般流程如图1所示。其核心步骤是框⑧所示的在每一积分步根据
Figure 728039DEST_PATH_IMAGE005
、
Figure 690178DEST_PATH_IMAGE006
求解（1），（2）式所表示的微分-代数方程组，得到
Figure 256289DEST_PATH_IMAGE007
和。目前，在电力系统数值仿真领域求解（1）式微分方程组的常用方法有隐式梯形积分法、改进欧拉法、龙格-库塔法等。隐式梯形积分数值稳定性好，但需要进行多次迭代求解，计算量大，目前电力系统商业计算程序BPA、PSASP采用的就是这种积分方法。改进欧拉法和龙格-库塔法为显式积分方法，无需迭代，计算量小，但数值稳定性较差。另外，显式积分算法要根据算法的截断误差，通过选择合理的积分步长，来保证算法的数值稳定性，如在电力系统中广泛应用的PSS/E程序采用的就是改进欧拉法。
为了同时保证算法的仿真精度和稳定性，计算时所取积分步长要与算法的截断误差成反比，即要使数值积分算法的截断误差越小，在同样精度要求下，积分步长可以取大一些，反之积分步长要取小一些。通常每一个积分步的截断误差越小，计算量也越大。如欧拉法的局部截断误差为
Figure 285239DEST_PATH_IMAGE010
，每一个积分步只需计算一次微分代数方程组；改进欧拉法的局部截断误差为
Figure 254463DEST_PATH_IMAGE011
，每一个积分步需计算两次微分代数方程组；四阶显式龙格－库塔法的局部截断误差为
Figure 51518DEST_PATH_IMAGE012
，每一个积分步需计算四次微分代数方程组。而隐式梯形积分法的局部截断误差为
Figure 722671DEST_PATH_IMAGE011
，则需经过多次迭代求解微分-代数方程，才能得到满足精度要求的解。以此看来，若能在提高算法截断误差的同时，不增加算法的计算量，则能减少整个暂态仿真的计算量，加快计算速度。 
目前，电力系统暂态稳定数值积分方法中均直接采用计算方法理论中的通用算法，如隐式梯形积分法、改进欧拉法、龙格-库塔法以及其他方法，并没有根据描述电力系统暂态过程的微分方程的特点对算法进行改进。
发明内容
本发明的目的是为了解决电力系统暂态稳定仿真方法中，现有的数值积分方法计算量大，计算速度不能满足电力系统在线计算要求的这个问题，提出了一种基于隐式数值积分的暂态稳定数值仿真方法。
本发明目的是通过以下技术方案实现的：一种基于隐式数值积分的电力系统暂态稳定仿真方法，包括以下步骤：
步骤1：输入系统的原始参数和信息，进行潮流计算得到稳态工况下的运行变量值
Figure 712142DEST_PATH_IMAGE013
，包括发电机节点的电压
Figure 722823DEST_PATH_IMAGE014
，注入网络的电流
Figure 753096DEST_PATH_IMAGE015
及发电机电磁功率；
步骤2：计算状态变量功角的初值
Figure 108302DEST_PATH_IMAGE017
、角频率的初值
Figure 301386DEST_PATH_IMAGE018
、发电机的暂态和次暂态电势、励磁及调速系统各动态环节状态变量组成的状态向量子向量的初值；
步骤3：形成描述系统暂态过程的微分方程和网络代数方程,并且进行网络代数方程因子表分解；
步骤4：置暂态稳定计算初值时刻
Figure 836721DEST_PATH_IMAGE020
，确定暂态稳定计算采用的积分步长
Figure 86437DEST_PATH_IMAGE009
，进行暂态稳定仿真计算；
步骤5：判断是否有故障或操作发生。若无，则转向步骤8；若有则执行步骤6；
步骤6：依据故障或操作情况，修改网络代数方程的因子表；
步骤7：求解网络代数方程，得到
Figure 134027DEST_PATH_IMAGE022
时刻的运行变量；
步骤8：计算
Figure 443786DEST_PATH_IMAGE023
时刻的系统的状态变量，运行变量，本步骤具体过程如下： 
步骤8.1：根据
Figure 327559DEST_PATH_IMAGE022
时刻系统的状态变量包括各台发电机功角、角频率
Figure 220746DEST_PATH_IMAGE025
和各发电机的暂态和次暂态电势、励磁及调速系统各动态环节状态变量组成的状态向量子向量
Figure 763723DEST_PATH_IMAGE026
，其中下标
Figure 56164DEST_PATH_IMAGE028
表示第
Figure 723381DEST_PATH_IMAGE028
台发电机。运行变量
Figure 479984DEST_PATH_IMAGE006
包括各节点电压
Figure 131546DEST_PATH_IMAGE029
和各节点的注入电流
Figure 989911DEST_PATH_IMAGE030
。 
用如下显式积分公式预测各发电机功角在时刻的初值：
                  
Figure 996231DEST_PATH_IMAGE031
其中，
Figure 84272DEST_PATH_IMAGE032
为系统同步角频率， 
Figure 164355DEST_PATH_IMAGE033
为第
Figure 97676DEST_PATH_IMAGE028
台发电机惯性时间常数，
Figure 828871DEST_PATH_IMAGE034
为第
Figure 87814DEST_PATH_IMAGE028
台原动机机械功率；
用显式欧拉公式预测各发电机角频率
Figure 389614DEST_PATH_IMAGE035
的初值
Figure 392205DEST_PATH_IMAGE036
，和其它各状态变量
Figure 977907DEST_PATH_IMAGE037
的初值
Figure 142172DEST_PATH_IMAGE038
。
 取运行变量初值
Figure 188058DEST_PATH_IMAGE039
。
步骤8.2：置迭代次数
Figure 728761DEST_PATH_IMAGE040
。
步骤8.3：按如下积分公式
Figure 168969DEST_PATH_IMAGE041
Figure 769715DEST_PATH_IMAGE042
    
求出发电机功角
Figure 46107DEST_PATH_IMAGE043
和角频率；
按一般隐式梯形积分公式：
Figure 950795DEST_PATH_IMAGE045
求出其它状态变量
Figure 456862DEST_PATH_IMAGE046
。
步骤8.4：求解网络代数方程
                     
       首先求解等式右边的虚拟电流
Figure 368635DEST_PATH_IMAGE048
，从而得到
Figure 783436DEST_PATH_IMAGE049
时刻系统的运行变量
Figure 460405DEST_PATH_IMAGE050
，进而得到各发电机的电磁功率
Figure 711388DEST_PATH_IMAGE051
。
步骤8.5：检查两次迭代各发电机最大电磁功率偏差值，若偏差大于给定精度
Figure 397585DEST_PATH_IMAGE052
，令，返回步骤8.3继续迭代；否则，执行步骤9；
步骤9：判断系统是否稳定，即任意两台发电机的最大相对摇摆功角是否大于某一给定值，若是执行步骤12；否则，执行步骤10；
步骤10：将仿真时间推进一个步长，令
Figure 780342DEST_PATH_IMAGE054
；
步骤11：判断是否到达事先给定的仿真时间
Figure 515692DEST_PATH_IMAGE055
。若
Figure 740000DEST_PATH_IMAGE056
则执行步骤12，否则返回步骤5；
步骤12：输出计算结果并结束计算。
在暂态稳定仿真过程中，每一积分步功角
Figure 67076DEST_PATH_IMAGE057
是采用如下所示局部截断误差为的隐式单步积分公式计算得到：
Figure 560691DEST_PATH_IMAGE059
功角预测是采用如下所示的局部截断误差为
Figure 401739DEST_PATH_IMAGE060
的单步显式数值积分： 
Figure 583322DEST_PATH_IMAGE061
本发明的有益效果：本发明将描述电力系统暂态过程的非线性微分方程组，表示为线性部分和非线性部分。通过合理将线性部分的系统矩阵选为奇异阵，就以得到状态转移矩阵的准确解析表达式；通过用一个线性可积函数去近似微分方程组的非线性部分，就可以求出一组隐式积分公式。其中，发电机功角隐式积分公式的局部截断误差为，高于局部截断误差为
Figure 714DEST_PATH_IMAGE060
的隐式梯形积分，而每一次积分的计算量则与隐式梯形积分相当。由于采用了高精度的隐式积分公式，使得在相同迭代精度下每一积分步迭代次数减少，从而显著地减少了仿真的计算量。
附图说明
图1为暂态稳定数值解法的一般流程图；
图2为隐式积分暂态稳定计算每一积分步的计算流程图；
图3为步长h=0.01时相对摇摆角最大的两台发电机间的转角差；
图4为本发明方法在不同步长下最大相对摇摆角偏差曲线；
图5为隐式梯形法在不同步长下最大相对摇摆角偏差曲线。
 
具体实施方式
以下结合附图对本发明作进一步说明。
本发明提出的一种基于隐式数值积分的电力系统暂态稳定仿真方法，包括以下步骤：
步骤1：输入系统的原始参数和信息，进行潮流计算得到稳态工况下的运行变量值
Figure 379874DEST_PATH_IMAGE013
，包括发电机节点的电压
Figure 415963DEST_PATH_IMAGE014
，注入网络的电流及发电机电磁功率
Figure 225973DEST_PATH_IMAGE016
；
步骤2：计算状态变量功角的初值
Figure 408824DEST_PATH_IMAGE017
、角频率的初值
Figure 564999DEST_PATH_IMAGE018
、发电机的暂态和次暂态电势、励磁及调速系统各动态环节状态变量组成的状态向量子向量的初值
Figure 424370DEST_PATH_IMAGE019
；
步骤3：形成描述系统暂态过程的微分方程和网络代数方程,并且进行网络代数方程因子表分解；
步骤4：置暂态稳定计算初值时刻
Figure 33206DEST_PATH_IMAGE020
，确定暂态稳定计算采用的积分步长
Figure 757098DEST_PATH_IMAGE009
，进行暂态稳定仿真计算；
步骤5：判断是否有故障或操作发生。若无，则转向步骤8；若有则执行步骤6；
步骤6：依据故障或操作情况，修改网络代数方程的因子表；
步骤7：求解网络代数方程，得到时刻的运行变量；
步骤8：计算
Figure 798052DEST_PATH_IMAGE023
时刻的系统的状态变量，运行变量，本步骤具体过程如下： 
步骤8.1：根据
Figure 894184DEST_PATH_IMAGE022
时刻系统的状态变量包括各台发电机功角
Figure 418838DEST_PATH_IMAGE024
、角频率和各发电机的暂态和次暂态电势、励磁及调速系统各动态环节状态变量组成的状态向量子向量
Figure 750779DEST_PATH_IMAGE026
，其中下标
Figure 68628DEST_PATH_IMAGE028
表示第
Figure 583923DEST_PATH_IMAGE028
台发电机。运行变量
Figure 444563DEST_PATH_IMAGE006
包括各节点电压
Figure 567371DEST_PATH_IMAGE029
和各节点的注入电流
Figure 372515DEST_PATH_IMAGE030
。 
用如下显式积分公式预测各发电机功角在
Figure 488239DEST_PATH_IMAGE023
时刻的初值：
                  
其中，
Figure 618799DEST_PATH_IMAGE032
为系统同步角频率， 
Figure 176819DEST_PATH_IMAGE062
为第
Figure 830654DEST_PATH_IMAGE028
台发电机惯性时间常数，为第
Figure 255130DEST_PATH_IMAGE028
台原动机机械功率；
用显式欧拉公式预测各发电机角频率
Figure 34868DEST_PATH_IMAGE035
的初值
Figure 492394DEST_PATH_IMAGE036
，和其它各状态变量
Figure 306766DEST_PATH_IMAGE037
的初值
Figure 942278DEST_PATH_IMAGE038
。
 取运行变量初值
Figure 209311DEST_PATH_IMAGE039
。
步骤8.2：置迭代次数
Figure 470528DEST_PATH_IMAGE040
。
步骤8.3：按如下积分公式
Figure 952456DEST_PATH_IMAGE041
Figure 945820DEST_PATH_IMAGE042
    
求出发电机功角
Figure 590164DEST_PATH_IMAGE043
和角频率
Figure 592755DEST_PATH_IMAGE044
；
按一般隐式梯形积分公式：
Figure 178457DEST_PATH_IMAGE045
求出其它状态变量
Figure 155772DEST_PATH_IMAGE046
。
步骤8.4：求解网络代数方程
                     
Figure 397397DEST_PATH_IMAGE047
   首先求解等式右边的虚拟电流
Figure 417DEST_PATH_IMAGE048
，从而得到
Figure 378308DEST_PATH_IMAGE049
时刻系统的运行变量
Figure 792103DEST_PATH_IMAGE050
，进而得到各发电机的电磁功率
Figure 255446DEST_PATH_IMAGE051
。
步骤8.5：检查两次迭代各发电机最大电磁功率偏差值，若偏差大于给定精度
Figure 662156DEST_PATH_IMAGE052
，令
Figure 894555DEST_PATH_IMAGE053
，返回步骤8.3继续迭代；否则，执行步骤9；
步骤9：判断系统是否稳定，即任意两台发电机的最大相对摇摆功角是否大于某一给定值，若是执行步骤12；否则，执行步骤10；
步骤10：将仿真时间推进一个步长，令
Figure 479251DEST_PATH_IMAGE054
；
步骤11：判断是否到达事先给定的仿真时间
Figure 429889DEST_PATH_IMAGE055
。若则执行步骤12，否则返回步骤5；
步骤12：输出计算结果并结束计算。
在暂态稳定仿真过程中，每一积分步功角
Figure 727195DEST_PATH_IMAGE057
是采用如下所示局部截断误差为
Figure 479863DEST_PATH_IMAGE058
的隐式单步积分公式计算得到：
功角预测是采用如下所示的局部截断误差为
Figure 666311DEST_PATH_IMAGE060
的单步显式数值积分： 
Figure 873301DEST_PATH_IMAGE061
以下详细介绍本发明方法的具体过程。
微分方程组（1）主要包括描述发电机组、感应电动机和其它动态装置动态特性的微分方程，其中每一台发电机组的微分方程可表示为：
                          　   （3）
式中，
Figure 725031DEST_PATH_IMAGE064
分别表示第
Figure 11656DEST_PATH_IMAGE028
台发电机的功角、角频率、机械功率、电磁功率及惯性时间常数，
Figure 338732DEST_PATH_IMAGE032
为系统同步角频率。
Figure 357504DEST_PATH_IMAGE065
为
Figure 583080DEST_PATH_IMAGE028
台第发电机的暂态和次暂态电势、励磁及调速系统各动态环节状态变量组成的状态向量子向量，
Figure 611079DEST_PATH_IMAGE066
为与状态向量子向量
Figure 854978DEST_PATH_IMAGE065
对应的微分方程右侧的函数向量。这样，每一台发电机组的状态向量可表示为：
Figure 757523DEST_PATH_IMAGE068
。
微分方程组(3)就可进一步表示为：
Figure 323634DEST_PATH_IMAGE069
      　　（4）
式中，
Figure 687619DEST_PATH_IMAGE070
为非线性微分方程组线性部分的系统矩阵，
Figure 313772DEST_PATH_IMAGE071
为常数；
Figure 251291DEST_PATH_IMAGE072
为非线性方程组，非线性部分的函数向量。
线性系统矩阵
Figure 621093DEST_PATH_IMAGE073
是一奇异矩阵。显然有：
Figure 839584DEST_PATH_IMAGE074
这样，可得到其所对应的状态转移矩阵准确的解析表达式为：
Figure 636639DEST_PATH_IMAGE075
                   　        　 （5）
对非线性函数向量
Figure 58524DEST_PATH_IMAGE076
，用一个近似的线性函数表示：
Figure 966437DEST_PATH_IMAGE077
    （6）
将（5）、（6）式代入Duhamel积分方程 ： 
Figure 39436DEST_PATH_IMAGE078
当取时，可得：
Figure 882758DEST_PATH_IMAGE080
（7）                                                                           
式（7）中除了功角
Figure DEST_PATH_IMAGE081
的积分公式外，均为局部截断误差为
Figure 791939DEST_PATH_IMAGE011
的隐式梯形积分，而功角积分公式的局部截断误差为
Figure 503543DEST_PATH_IMAGE012
。
当取时，可得：
Figure DEST_PATH_IMAGE083
                    （8）
为此，根据本发明，基于隐式数值积分的电力系统暂态稳定数值仿真方法每一积分步的计算步骤如下：
1. 根据
Figure 380025DEST_PATH_IMAGE022
时刻系统的状态变量包括各台发电机功角
Figure 963453DEST_PATH_IMAGE024
、角频率
Figure 541065DEST_PATH_IMAGE025
和各发电机的暂态和次暂态电势、励磁及调速系统各动态环节状态变量组成的状态向量子向量
Figure 260759DEST_PATH_IMAGE026
，其中下标
Figure 383567DEST_PATH_IMAGE028
表示第
Figure 454291DEST_PATH_IMAGE028
台发电机。运行变量包括各节点电压
Figure 222844DEST_PATH_IMAGE029
和各节点的注入电流
Figure 765821DEST_PATH_IMAGE030
。 用(8)式预测各发电机功角在时刻的初值
Figure 728409DEST_PATH_IMAGE084
,用显式欧拉公式预测各发电机角频率
Figure 422695DEST_PATH_IMAGE035
的初值
Figure 402153DEST_PATH_IMAGE036
，和其它各状态变量
Figure 181890DEST_PATH_IMAGE037
的初值
Figure 381359DEST_PATH_IMAGE038
。取运行变量初值。 
2.置迭代次数
Figure DEST_PATH_IMAGE085
。
3.按积分公式（7）求出发电机功角
Figure 80511DEST_PATH_IMAGE043
、角频率
Figure 160594DEST_PATH_IMAGE086
和其它状态变量
Figure 359494DEST_PATH_IMAGE046
；
4.求解网络代数方程
 
Figure DEST_PATH_IMAGE087
首先求解等式右边的虚拟电流
Figure 90689DEST_PATH_IMAGE048
，从而得到
Figure 897102DEST_PATH_IMAGE049
时刻系统的运行变量
Figure 651432DEST_PATH_IMAGE050
，进而得到各发电机的电磁功率
Figure 450761DEST_PATH_IMAGE051
。
5. 检查两次迭代各发电机最大电磁功率偏差值，若偏差大于给定精度，令
Figure 217039DEST_PATH_IMAGE053
，返回3继续迭代；否则，本积分步迭代过程结束。
其计算流程如图2所示。
将本发明所提出的电力系统暂态稳定计算每一积分步的计算流程（图2）嵌入暂态稳定计算一般流程（图1）的框⑧，就可实现本发明所提出的隐式积分暂态稳定仿真方法。
以下是本发明方法的一个实施例，以IEEE145节点系统进行仿真实验作实施例，进一步说明如下：
IEEE145系统中所有发电机均采用恒定模型，负荷采用恒定阻抗模型。在节点131与节点130线路的始端
Figure DEST_PATH_IMAGE089
发生三相短路故障，
Figure 606225DEST_PATH_IMAGE090
切除故障线路。图3为积分步长取，用隐式梯形积分法得出的发电机间最大相对摇摆角曲线，并以此作为准确结果。
图4、图5分别给出了本发明所述方法和隐式梯形数值积分法在不同积分步长下计算结果与图3所示标准结果的偏差。由图4和图5可见，积分步长取0.02～0.06秒时本发明所述方法最大相对摇摆角的误差均不超过1.5度，均小于隐式梯形积分方法的误差。
表1对本发明方法(用DI表示)与隐式梯形积分方法(用TI表示)在不同积分步长下整个方程过程求解微分-代数方程组的次数以及计算CPU时间进行了比较：
表1. 本发明方法与隐式梯形积分法计算量比较
Figure DEST_PATH_IMAGE093
从表一可以看出本发明方法在整个暂态稳定仿真过程中求解网络方程的次数和CPU计算时间约比隐式梯形方法少10%左右。当积分步长增加到0.06秒时，隐式梯形积分方法求解微分代数方程组的次数反而大于积分步长为0.05时的情况，说明此时隐式梯形法趋于发散，而本发明方法则计算量继续减少。
从实施的仿真算例可知，无论是计算量还是计算精度均优于传统的隐式梯形积分仿真方法。

Claims (3)

1.一种基于隐式数值积分的电力系统暂态稳定仿真方法，其特征在于该方法包括以下步骤：
步骤1：输入系统的原始参数和信息，进行潮流计算得到稳态工况下的运行变量值                                               
Figure DEST_PATH_IMAGE002
，包括发电机节点的电压
Figure DEST_PATH_IMAGE004
，注入网络的电流
Figure DEST_PATH_IMAGE006
及发电机电磁功率
Figure DEST_PATH_IMAGE008
；
步骤2：计算状态变量功角的初值
Figure DEST_PATH_IMAGE010
、角频率的初值
Figure DEST_PATH_IMAGE012
、发电机的暂态和次暂态电势、励磁及调速系统各动态环节状态变量组成的状态向量子向量的初值
Figure 157899DEST_PATH_IMAGE014
；
步骤3：形成描述系统暂态过程的微分方程和网络代数方程,并且进行网络代数方程因子表分解；
步骤4：置暂态稳定计算初值时刻
Figure 555382DEST_PATH_IMAGE016
，确定暂态稳定计算采用的积分步长
Figure 284304DEST_PATH_IMAGE018
，进行暂态稳定仿真计算；
步骤5：判断是否有故障或操作发生；若无，则转向步骤8；若有则执行步骤6；
步骤6：依据故障或操作情况，修改网络代数方程的因子表；
步骤7：求解网络代数方程，得到时刻的运行变量；
步骤8：计算
Figure 674145DEST_PATH_IMAGE022
时刻的系统的状态变量，运行变量，本步骤具体过程如下： 
步骤8.1：根据
Figure 445792DEST_PATH_IMAGE020
时刻系统的状态变量包括各台发电机功角、角频率
Figure 341253DEST_PATH_IMAGE026
和各发电机的暂态和次暂态电势、励磁及调速系统各动态环节状态变量组成的状态向量子向量
Figure 506786DEST_PATH_IMAGE028
，其中下标表示第
Figure 949586DEST_PATH_IMAGE030
台发电机；运行变量包括各节点电压
Figure 655822DEST_PATH_IMAGE034
和各节点的注入电流
Figure 831588DEST_PATH_IMAGE036
； 
用如下显式积分公式预测各发电机功角在时刻的初值：
                   
Figure 797106DEST_PATH_IMAGE038
其中，为系统同步角频率， 
Figure 18320DEST_PATH_IMAGE042
为第
Figure 430847DEST_PATH_IMAGE030
台发电机惯性时间常数，
Figure 521162DEST_PATH_IMAGE044
为第
Figure 702745DEST_PATH_IMAGE030
台原动机机械功率；
用显式欧拉公式预测各发电机角频率
Figure 705467DEST_PATH_IMAGE046
的初值
Figure 605290DEST_PATH_IMAGE048
，和其它各状态变量的初值
Figure 535386DEST_PATH_IMAGE052
；
 取运行变量初值
Figure 709009DEST_PATH_IMAGE054
；
步骤8.2：置迭代次数；
步骤8.3：按如下积分公式
Figure DEST_PATH_IMAGE058
Figure DEST_PATH_IMAGE060
    
求出发电机功角
Figure DEST_PATH_IMAGE062
和角频率
Figure DEST_PATH_IMAGE064
；
按一般隐式梯形积分公式：
求出其它状态变量
Figure DEST_PATH_IMAGE068
；
步骤8.4：求解网络代数方程
          
Figure DEST_PATH_IMAGE070
      式中，
Figure DEST_PATH_IMAGE072
为节点导纳矩阵；
   首先求解等式右边的虚拟电流
Figure DEST_PATH_IMAGE074
，从而得到
Figure DEST_PATH_IMAGE076
时刻系统的运行变量，进而得到各发电机的电磁功率
Figure DEST_PATH_IMAGE080
；
步骤8.5：检查两次迭代各发电机最大电磁功率偏差值，若偏差大于给定精度
Figure DEST_PATH_IMAGE082
，令
Figure DEST_PATH_IMAGE084
，返回步骤8.3继续迭代；否则，执行步骤9；
步骤9：判断系统是否稳定，即任意两台发电机的最大相对摇摆功角是否大于某一给定值，若是执行步骤12；否则，执行步骤10；
步骤10：将仿真时间推进一个步长，令
Figure DEST_PATH_IMAGE086
；
步骤11：判断是否到达事先给定的仿真时间
Figure DEST_PATH_IMAGE088
；若
Figure DEST_PATH_IMAGE090
则执行步骤12，否则返回步骤5；
步骤12：输出计算结果并结束计算。
2.根据权利1要求所述的基于隐式数值积分的电力系统暂态稳定数值仿真方法，其特征在于：在步骤8.3中，即每一积分步状态变量功角
Figure DEST_PATH_IMAGE092
是采用如下所示局部截断误差为的隐式单步积分公式计算得到：
Figure DEST_PATH_IMAGE096
。
3.根据权利1要求所述的基于隐式数值积分的电力系统暂态稳定数值仿真方法，其特征在于：在每一积分步功角预测是采用如下所示的局部截断误差为
Figure DEST_PATH_IMAGE098
的单步显式数值积分： 
Figure DEST_PATH_IMAGE100
。
CN2012100175933A 2012-01-19 2012-01-19 一种基于隐式数值积分的电力系统暂态稳定仿真方法 Expired - Fee Related CN102609575B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN2012100175933A CN102609575B (zh) 2012-01-19 2012-01-19 一种基于隐式数值积分的电力系统暂态稳定仿真方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN2012100175933A CN102609575B (zh) 2012-01-19 2012-01-19 一种基于隐式数值积分的电力系统暂态稳定仿真方法

Publications (2)

Publication Number Publication Date
CN102609575A true CN102609575A (zh) 2012-07-25
CN102609575B CN102609575B (zh) 2013-12-04

Family

ID=46526943

Family Applications (1)

Application Number Title Priority Date Filing Date
CN2012100175933A Expired - Fee Related CN102609575B (zh) 2012-01-19 2012-01-19 一种基于隐式数值积分的电力系统暂态稳定仿真方法

Country Status (1)

Country Link
CN (1) CN102609575B (zh)

Cited By (18)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102928014A (zh) * 2012-10-23 2013-02-13 保定市三川电气有限责任公司 电力系统数字测量或遥测处理的方法及装置
CN104079228A (zh) * 2014-07-09 2014-10-01 国家电网公司 双馈异步电机的带阻尼的隐式梯形积分电磁暂态建模方法
CN104156542A (zh) * 2014-08-26 2014-11-19 天津大学 一种基于隐式投影的有源配电系统稳定性仿真方法
CN104701857A (zh) * 2015-03-30 2015-06-10 河南行知专利服务有限公司 一种电力系统稳态运行控制方法
CN106295001A (zh) * 2016-08-10 2017-01-04 华北电力大学 适用于电力系统中长时间尺度的准稳态变步长仿真方法
CN106855909A (zh) * 2017-01-13 2017-06-16 国网北京经济技术研究院 一种适于有源配电网随机动态仿真的显隐混合积分方法
CN107247686A (zh) * 2017-05-22 2017-10-13 电子科技大学 一种基于并行算法的fetd仿真模拟方法
CN107526856A (zh) * 2016-06-22 2017-12-29 南京理工大学 可并行的显隐式混合不连续伽辽金时域有限元法
CN107706948A (zh) * 2017-11-03 2018-02-16 华北电力大学 多维阶数控制的多步Taylor级数暂态稳定分析方法
CN108599198A (zh) * 2018-05-15 2018-09-28 中国电力科学研究院有限公司 一种基于模型预测控制的电力系统频率和电压协调校正方法及系统
CN108763790A (zh) * 2018-06-01 2018-11-06 三峡大学 一种基于扩展临界阻尼调整法的电力系统电磁暂态仿真方法
CN110119523A (zh) * 2019-01-29 2019-08-13 中国电力科学研究院有限公司 一种基于功角等效法的网络代数方程求解预处理方法及系统
CN110135031A (zh) * 2019-04-30 2019-08-16 东南大学 基于半隐式龙格库塔法的电力系统暂态稳定计算方法
CN110781443A (zh) * 2019-10-14 2020-02-11 南京理工大学 一种多尺度量子电磁耦合的含时计算方法
CN110991034A (zh) * 2019-11-29 2020-04-10 浙江大学 一种基于全并行嵌套对角加边形式的电力系统暂态稳定仿真并行计算方法
CN111900718A (zh) * 2020-06-17 2020-11-06 东南大学 基于多级优化追赶变分迭代法的有源配电网动态仿真方法
CN113221298A (zh) * 2021-04-21 2021-08-06 南方电网科学研究院有限责任公司 一种机电暂态过程的仿真方法及系统
CN116244894A (zh) * 2022-12-09 2023-06-09 山东大学 一种基于大步长的电力系统暂态仿真方法及系统

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101673948A (zh) * 2009-09-28 2010-03-17 国网电力科学研究院 电力系统暂态稳定分析的n-2组合故障快速排序方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101673948A (zh) * 2009-09-28 2010-03-17 国网电力科学研究院 电力系统暂态稳定分析的n-2组合故障快速排序方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
汤涌: "电力系统稳定计算隐式积分交替求解", 《电网技术》 *
王宇宾 等: "基于隐式Taylor级数法的电力系统暂态稳定计算", 《华北电力大学学报》 *

Cited By (30)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102928014B (zh) * 2012-10-23 2015-05-13 保定市三川电气有限责任公司 电力系统数字测量或遥测处理的方法及装置
CN102928014A (zh) * 2012-10-23 2013-02-13 保定市三川电气有限责任公司 电力系统数字测量或遥测处理的方法及装置
CN104079228A (zh) * 2014-07-09 2014-10-01 国家电网公司 双馈异步电机的带阻尼的隐式梯形积分电磁暂态建模方法
CN104079228B (zh) * 2014-07-09 2016-06-29 国家电网公司 双馈异步电机的带阻尼的隐式梯形积分电磁暂态建模方法
CN104156542B (zh) * 2014-08-26 2018-01-02 天津大学 一种基于隐式投影的有源配电系统稳定性仿真方法
CN104156542A (zh) * 2014-08-26 2014-11-19 天津大学 一种基于隐式投影的有源配电系统稳定性仿真方法
CN104701857A (zh) * 2015-03-30 2015-06-10 河南行知专利服务有限公司 一种电力系统稳态运行控制方法
CN107526856A (zh) * 2016-06-22 2017-12-29 南京理工大学 可并行的显隐式混合不连续伽辽金时域有限元法
CN107526856B (zh) * 2016-06-22 2022-01-14 南京理工大学 可并行的显隐式混合不连续伽辽金时域有限元法
CN106295001B (zh) * 2016-08-10 2020-02-18 华北电力大学 适用于电力系统中长时间尺度的准稳态变步长仿真方法
CN106295001A (zh) * 2016-08-10 2017-01-04 华北电力大学 适用于电力系统中长时间尺度的准稳态变步长仿真方法
CN106855909A (zh) * 2017-01-13 2017-06-16 国网北京经济技术研究院 一种适于有源配电网随机动态仿真的显隐混合积分方法
CN106855909B (zh) * 2017-01-13 2020-03-20 国网北京经济技术研究院 一种适于有源配电网随机动态仿真的显隐混合积分方法
CN107247686B (zh) * 2017-05-22 2020-11-17 电子科技大学 一种基于并行算法的fetd仿真模拟方法
CN107247686A (zh) * 2017-05-22 2017-10-13 电子科技大学 一种基于并行算法的fetd仿真模拟方法
CN107706948A (zh) * 2017-11-03 2018-02-16 华北电力大学 多维阶数控制的多步Taylor级数暂态稳定分析方法
CN107706948B (zh) * 2017-11-03 2020-12-01 华北电力大学 多维阶数控制的多步Taylor级数暂态稳定分析方法
CN108599198A (zh) * 2018-05-15 2018-09-28 中国电力科学研究院有限公司 一种基于模型预测控制的电力系统频率和电压协调校正方法及系统
CN108599198B (zh) * 2018-05-15 2021-06-04 中国电力科学研究院有限公司 一种基于模型预测控制的电力系统频率和电压协调校正方法及系统
CN108763790A (zh) * 2018-06-01 2018-11-06 三峡大学 一种基于扩展临界阻尼调整法的电力系统电磁暂态仿真方法
CN110119523A (zh) * 2019-01-29 2019-08-13 中国电力科学研究院有限公司 一种基于功角等效法的网络代数方程求解预处理方法及系统
CN110135031A (zh) * 2019-04-30 2019-08-16 东南大学 基于半隐式龙格库塔法的电力系统暂态稳定计算方法
CN110781443A (zh) * 2019-10-14 2020-02-11 南京理工大学 一种多尺度量子电磁耦合的含时计算方法
CN110781443B (zh) * 2019-10-14 2023-09-26 南京理工大学 一种多尺度量子电磁耦合的含时计算方法
CN110991034A (zh) * 2019-11-29 2020-04-10 浙江大学 一种基于全并行嵌套对角加边形式的电力系统暂态稳定仿真并行计算方法
CN111900718A (zh) * 2020-06-17 2020-11-06 东南大学 基于多级优化追赶变分迭代法的有源配电网动态仿真方法
CN111900718B (zh) * 2020-06-17 2022-08-09 东南大学 基于多级优化追赶变分迭代法的有源配电网动态仿真方法
CN113221298A (zh) * 2021-04-21 2021-08-06 南方电网科学研究院有限责任公司 一种机电暂态过程的仿真方法及系统
CN116244894A (zh) * 2022-12-09 2023-06-09 山东大学 一种基于大步长的电力系统暂态仿真方法及系统
CN116244894B (zh) * 2022-12-09 2023-09-15 山东大学 一种基于大步长的电力系统暂态仿真方法及系统

Also Published As

Publication number Publication date
CN102609575B (zh) 2013-12-04

Similar Documents

Publication Publication Date Title
CN102609575B (zh) 一种基于隐式数值积分的电力系统暂态稳定仿真方法
CN103995734A (zh) 基于rtds的电力系统混合实时仿真系统及仿真方法
CN103700036B (zh) 一种适于电力系统多时间尺度的暂态稳定性投影积分方法
CN102664397B (zh) 一种基于隐式精细数值积分的电力系统暂态稳定仿真方法
CN103500245B (zh) 一种基于多回路法的场路瞬态-瞬态耦合仿真方法
CN102545263B (zh) 一种基于显式数值积分的电力系统暂态稳定仿真方法
CN104732033A (zh) 一种电力系统机电暂态过程的仿真分析方法
Di Leonardo et al. Transient analysis of PM synchronous motor drives by finite element model co-simulation
US11569682B2 (en) System and method for a fast power network simulator
Therrien et al. Constant-parameter interfacing of induction machine models considering main flux saturation in EMTP-type programs
CN115021638A (zh) 一种高效同步机的电磁暂态建模方法、系统及设备
CN108959671B (zh) 半桥和全桥型模块化多电平换流器的实时仿真建模方法
Pandey et al. Unified power system analyses and models using equivalent circuit formulation
WO2006054784A1 (ja) モータモデルの演算方法、モータシミュレーション方法、モータシミュレーション装置、モータモデル演算プログラム、シミュレーション方法及びシミュレーションプログラム
Xia et al. An efficient phase domain synchronous machine model with constant equivalent admittance matrix
Gao et al. Shifted frequency‐based electromagnetic transient simulation for AC power systems in symmetrical component domain
CN109379003A (zh) 高速永磁同步电机控制的方法、装置及计算机存储介质
CN109150049A (zh) 一种静止坐标系电机分布式参数模型建立方法
CN104778317A (zh) 一种基于套接字技术的电力系统机电-电磁混合仿真方法
CN106992732A (zh) 一种电机磁共能模型建立系统
CN102609576B (zh) 预估-校正数值积分的电力系统暂态稳定仿真方法
CN105720874B (zh) 基于分布参数的电机气隙磁场建模方法及其应用
CN102545216A (zh) 用于电力系统暂态稳定仿真的发电机节点电压预测方法
CN114977165A (zh) 基于饱和效应的同步机电磁暂态建模方法、系统及设备
CN113326673B (zh) 一种同步电机的vbr模型电磁暂态仿真方法及系统

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20131204

Termination date: 20210119

CF01 Termination of patent right due to non-payment of annual fee