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

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

Info

Publication number
CN102664397A
CN102664397A CN2012100794539A CN201210079453A CN102664397A CN 102664397 A CN102664397 A CN 102664397A CN 2012100794539 A CN2012100794539 A CN 2012100794539A CN 201210079453 A CN201210079453 A CN 201210079453A CN 102664397 A CN102664397 A CN 102664397A
Authority
CN
China
Prior art keywords
state
generator
transient
value
constantly
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
CN2012100794539A
Other languages
English (en)
Other versions
CN102664397B (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 CN201210079453.9A priority Critical patent/CN102664397B/zh
Publication of CN102664397A publication Critical patent/CN102664397A/zh
Application granted granted Critical
Publication of CN102664397B publication Critical patent/CN102664397B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Supply And Distribution Of Alternating Current (AREA)

Abstract

本发明公开了一种基于精细积分算法的电力系统暂态稳定仿真方法。与已有的电力系统暂态稳定数值积分方法相比，本发明将描述电力系统暂态过程的非线性微分方程组，表示为线性部分和非线性部分。通过合理地选取线性部分的系统矩阵和用一个可积函数逼近非线性项函数，推导出了隐式精细单步积分公式。该方法通过精细化地求解线性部分系统矩阵所对应的状态转移矩阵，使得积分公式有较高的精度，从而减少了暂态稳定仿真的计算量。

Description

一种基于隐式精细数值积分的电力系统暂态稳定仿真方法
技术领域
本发明属于电力系统自动化技术领域，特别涉及了一种基于隐式精细数值积分的电力系统暂态稳定仿真方法。
背景技术
电力系统暂态稳定分析是电力系统分析计算中最基础、最核心的内容之一。随着电力系统运行控制技术的发展，在线动态安全分析、安全稳定紧急控制、预防控制、智能调度等先进技术已逐步在电力系统中推广应用。而实现这些先进技术的前提条件则是能够对大规模电力系统暂态稳定计算做到快速、准确、可靠。
数值积分法和直接法，以及将数值积分和直接法相结合的混合分析方法是电力系统暂态稳定分析当中常用的主要方法。在这些方法当中，数值积分是电力系统暂态稳定计算最准确、最可靠的方法。数值积分法的最大缺点是计算量大。尽管计算机运算速度较之过去已经有了提高，但对于大规模电力系统，计算速度要想满足在线动态安全分析、预防控制、紧急控制等计算的要求，仍然是远远不够的。
电力系统的各种元件都可以由数学模型表示出来，其暂态过程可用如下形式的微分-代数方程组描述
Figure 2012100794539100002DEST_PATH_IMAGE002
                                      (1)
式中，
Figure 2012100794539100002DEST_PATH_IMAGE004
表示微分方程组中描述系统动态特性的状态变量，通常向量
Figure 401760DEST_PATH_IMAGE004
包含发电机功角和转速等动态元件的状态量；
Figure 2012100794539100002DEST_PATH_IMAGE006
表示代数方程组中系统的运行变量, 通常包含与网络相关的运行变量，如节点电压的幅值和相位等。
用数值积分法求解电力系统暂态过程的一般流程如图1所示。其核心步骤是框⑧所示的利用
Figure 2012100794539100002DEST_PATH_IMAGE008
、
Figure 2012100794539100002DEST_PATH_IMAGE010
求解（1）式所表示的微分-代数方程组得到下一积分步的解
Figure 2012100794539100002DEST_PATH_IMAGE012
和
Figure 2012100794539100002DEST_PATH_IMAGE014
。目前，在电力系统数值仿真领域求解（1）式中微分方程组的常用方法有隐式梯形积分法、改进欧拉法、龙格-库塔法等。隐式梯形积分法数值稳定性好，但需要多次迭代求解，计算量大，目前采用的这种积分方法的电力系统商业计算程序有PSD-BPA、PSASP。龙格-库塔法，高阶Taylor展开法为典型的显式积分方法，虽然无需迭代，但为了达到一定的计算精度，还需在每一积分步求解多次微分-代数方程组来保证其精度，而且其数值稳定性较差，容易导致计算失败。因此，一般情况下显式积分算法要根据算法的截断误差，通过选择合理的积分步长，来保证算法的收敛性。
为了保证算法的稳定性和仿真精度，所取积分步长要与算法的截断误差成反比，即数值积分算法的截断误差越小，在相同精度要求下，积分步长
Figure 2012100794539100002DEST_PATH_IMAGE016
可取得大一些，反之积分步长
Figure 194267DEST_PATH_IMAGE016
要取得小一些。通常每一积分步的截断误差越小，计算量也越大。如欧拉法的局部截断误差为，每一积分步只需计算一次微分代数方程；改进欧拉法的局部截断误差为
Figure 2012100794539100002DEST_PATH_IMAGE020
，每一积分步需计算两次微分代数方程；四阶显式龙格－库塔法的局部截断误差为
Figure 2012100794539100002DEST_PATH_IMAGE022
，每一积分步需计算四次微分代数方程。而隐式梯形积分法的局部截断误差为
Figure 659490DEST_PATH_IMAGE020
，则需经过多次迭代求解微分-代数方程，才能得到满足精度要求的解。若能在提高算法截断误差的同时，不增加算法的计算量，则能减少整个暂态仿真的计算量，加快计算速度。
目前，在电力系统暂态稳定数值积分仿真方法中所采用的数值积分方法，均直接采用计算方法理论中的通用算法，如隐式梯形积分法、改进欧拉法、龙格-库塔法等，并没有根据描述电力系统暂态过程的微分方程的特点对算法进行改进。在交替求解微分代数方程时，通常是求出全部状态变量、再求取运行变量。本发明中的数值积分法基于发电机模型的特点，将描述电力系统暂态过程的非线性微分方程组，表示为线性部分和非线性部分。通过合理地选取线性部分的系统矩阵和用一个可积函数逼近非线性项函数，推导出了隐式精细单步积分公式。该算法通过精细化地求解线性部分系统矩阵所对应的状态转移矩阵，使得积分公式有较高的精度，从而减少了暂态稳定仿真的计算量。
发明内容
本发明目的是为了解决电力系统暂态稳定仿真计算中，现有的数值积分方法计算量大，计算速度不能满足电力系统在线计算要求的缺点，基于非线性方程的多哈米尔积分和精细积分算法，提出了一种新的暂态稳定数值积分仿真方法。该方法充分利用了发电机模型及其他设备可以表示为线性和非线性部分并能写成传输函数的特点，推导出了基于隐式精细积分法的单步积分公式。该积分方法计算量小于局部截断误差为的隐式梯形法；该方法还充分利用了发电机转子角方程可以同其他系统元件解耦的特点，提高了算法的通用性。
本发明目的是通过以下技术方案实现的：一种基于隐式精细积分的电力系统暂态稳定仿真方法，包括以下步骤：
步骤1：输入系统的原始参数和信息，进行潮流计算得到稳态工况下的运行变量值
Figure 2012100794539100002DEST_PATH_IMAGE024
，包括各台发电机节点的电压，注入网络的电流
Figure 2012100794539100002DEST_PATH_IMAGE028
，发电机电磁功率
Figure 2012100794539100002DEST_PATH_IMAGE030
，其中
Figure 2012100794539100002DEST_PATH_IMAGE032
（为发电机台数）。
步骤2：计算状态变量初值包括发电机功角
Figure 2012100794539100002DEST_PATH_IMAGE036
、角频率
Figure 2012100794539100002DEST_PATH_IMAGE038
和发电机的暂态和次暂态电势、励磁及调速系统各动态环节状态变量初值。
步骤3：形成描述系统暂态过程的微分方程和网络代数方程，并且进行网络代数方程因子表分解。
步骤4：置暂态稳定计算初值时刻。
步骤5：判断是否有故障或操作发生。若无，则转向步骤8；若有，则执行步骤6。
步骤6：依据故障或操作情况，修改微分方程和网络代数方程及其因子表。
步骤7：求解网络代数方程，得到
Figure 2012100794539100002DEST_PATH_IMAGE042
时刻的运行变量。
步骤8：计算
Figure 2012100794539100002DEST_PATH_IMAGE044
时刻的系统的状态变量值包括各台发电机功角
Figure 2012100794539100002DEST_PATH_IMAGE046
、角频率
Figure 2012100794539100002DEST_PATH_IMAGE048
和发电机的暂态和次暂态电势、励磁及调速系统各动态环节状态变量值，及运行变量值包括发电机节点的电压
Figure 2012100794539100002DEST_PATH_IMAGE050
，注入网络的电流
Figure 2012100794539100002DEST_PATH_IMAGE052
和电磁功率，本步骤具体过程如下：
步骤8.1：判断本步的步长
Figure 706182DEST_PATH_IMAGE016
是否与前一时刻步长相同，若相同则跳转到步骤8.3，若不同，则进行步骤8.2的操作。
步骤8.2：根据微分方程得到状态矩阵
Figure 2012100794539100002DEST_PATH_IMAGE056
和非线性项
Figure 2012100794539100002DEST_PATH_IMAGE058
的线性化表达式，利用步长
Figure 661631DEST_PATH_IMAGE016
，计算出状态转移矩阵
Figure 2012100794539100002DEST_PATH_IMAGE060
、非线性项状态线性化后的常数项状态转移矩阵的
Figure 2012100794539100002DEST_PATH_IMAGE062
和非线性项线性化后的一次项状态转移矩阵的。
步骤8.3：给定
Figure 2012100794539100002DEST_PATH_IMAGE066
时刻的状态变量值及运行变量值，置迭代次数
Figure 2012100794539100002DEST_PATH_IMAGE068
。
步骤8.4：预估
Figure 2012100794539100002DEST_PATH_IMAGE070
时刻的状态变量值。
步骤8.5：根据网络代数方程以及
Figure 676510DEST_PATH_IMAGE070
状态变量值计算运行变量值。
步骤8.6：根据计算求得的
Figure 407706DEST_PATH_IMAGE070
的运行状态变量值，按下列公式求解时刻的状态变量值
Figure 2012100794539100002DEST_PATH_IMAGE072
，
其中，式中的
Figure 2012100794539100002DEST_PATH_IMAGE074
表示
Figure 843815DEST_PATH_IMAGE070
时刻的各台发电机的发电机功角
Figure 643143DEST_PATH_IMAGE046
、角频率
Figure 727381DEST_PATH_IMAGE048
、暂态和次暂态电势、励磁及调速系统各动态环节状态变量组成的状态向量子向量，
Figure 2012100794539100002DEST_PATH_IMAGE076
为
Figure 219542DEST_PATH_IMAGE066
时刻各台发电机的发电机功角
Figure 2012100794539100002DEST_PATH_IMAGE078
、角频率、暂态和次暂态电势、励磁及调速系统各动态环节非线性项线性化后的常数项组成的向量，
Figure 2012100794539100002DEST_PATH_IMAGE082
为
Figure 211900DEST_PATH_IMAGE066
和
Figure DEST_PATH_IMAGE084
时刻各台发电机的发电机功角
Figure DEST_PATH_IMAGE086
、角频率
Figure DEST_PATH_IMAGE088
、暂态和次暂态电势、励磁及调速系统各动态环节非线性项线性化后的线性部分组成的向量。
步骤8.7：检查两次迭代各发电机最大电磁功率偏差值，若偏差大于给定精度，令
Figure DEST_PATH_IMAGE092
，返回步骤8.5继续迭代；否则，执行步骤9。
步骤9：判断系统是否稳定，即任意两台发电机的最大相对摇摆功角是否大于某一给定值，若是执行步骤12；否则，执行步骤10。
步骤10：将仿真时间推进一个步长，令，并以此时的
Figure DEST_PATH_IMAGE096
值作为下一时刻的计算起点值，即
Figure 353601DEST_PATH_IMAGE066
值。
步骤11：判断是否到达事先给定的仿真时间
Figure DEST_PATH_IMAGE098
。若
Figure DEST_PATH_IMAGE100
则执行步骤12，否则返回步骤5。
步骤12：输出计算结果并结束计算。
在暂态稳定仿真计算步骤8.2 中计算状态转移矩阵和非线性项线性化后的常数项的状态转移矩阵
Figure 20653DEST_PATH_IMAGE062
和非线性项线性化后的一次项的状态转移矩阵使用的方法是精细积分法。
本发明的有益效果：该方法充分利用了精细积分算法的特点，推导出了精细积分法单步积分公式，该积分公式的局部截断误差低，而计算量则小于局部截断误差为
Figure 372744DEST_PATH_IMAGE020
的隐式梯形积分法。
 
附图说明
图1为暂态稳定数值解法的一般流程图；
图2为暂态稳定计算每一积分步的计算流程；
图3为IEEE39节点系统最大功角摇摆曲线；
图4为隐式梯形积分法最大功角计算误差；
图5为隐式精细积分法最大功角计算误差。
具体实施方式
以下结合附图对本发明作进一步说明。   
本发明基于非线性方程的多哈米尔积分，提出了一种新的高精度的数值仿真方法。本发明方法与传统的数值积分仿真方法的区别在于，将描述电力系统暂态过程的非线性微分方程组，表示为线性和非线性部分。通过精细化的求解线性部分的状态转移矩阵和合理的选取逼近非线性部分的可积函数，得到较高精度的数值积分公式。本方法对于模型的适应性较好，可以方面的添加新建模型，而不需要繁琐的步骤。
以下结合附图1对本发明作进一步说明。
本发明提出的一种基于隐式精细数值积分的电力系统暂态稳定仿真方法，包括以下步骤：
步骤1：输入系统的原始参数和信息，进行潮流计算得到稳态工况下的运行变量值，包括各台发电机节点的电压
Figure 74170DEST_PATH_IMAGE026
，注入网络的电流
Figure 845816DEST_PATH_IMAGE028
，电磁功率
Figure 858772DEST_PATH_IMAGE030
，其中
Figure 305059DEST_PATH_IMAGE032
Figure 719860DEST_PATH_IMAGE034
（
Figure 724725DEST_PATH_IMAGE034
为发电机台数）。
步骤2：计算状态变量初值包括各台发电机功角、角频率
Figure 645593DEST_PATH_IMAGE038
和发电机的暂态和次暂态电势、励磁及调速系统各动态环节状态变量初值。
步骤3：形成描述系统暂态过程的微分方程和网络代数方程，并形成因子表。
步骤4：置暂态稳定计算初值时刻
Figure 684874DEST_PATH_IMAGE040
。
步骤5：判断是否有故障或操作发生。若无，则转向步骤8；若有，则执行步骤6。
步骤6：依据故障或操作情况，修改微分方程和网络代数方程及其因子表。
步骤7：求解网络代数方程，得到
Figure 860641DEST_PATH_IMAGE042
时刻的运行变量。
步骤8：计算
Figure 848188DEST_PATH_IMAGE044
时刻的系统的状态变量值包括各台发电机功角
Figure 134813DEST_PATH_IMAGE046
、角频率
Figure 196310DEST_PATH_IMAGE048
和发电机的暂态和次暂态电势、励磁及调速系统各动态环节状态变量值，及运行变量值包括发电机节点的电压
Figure 44443DEST_PATH_IMAGE050
，注入网络的电流
Figure 519286DEST_PATH_IMAGE052
和电磁功率
Figure 609602DEST_PATH_IMAGE054
，本步骤具体过程如下：
步骤8.1：判断本步的步长
Figure 853502DEST_PATH_IMAGE016
是否与前一时刻步长相同，若相同则跳转到步骤8.3，若不同，则进行步骤8.2的操作。
步骤8.2：根据微分方程得到状态矩阵
Figure 43175DEST_PATH_IMAGE056
和非线性项
Figure 503850DEST_PATH_IMAGE058
的线性化表达式，利用步长，计算出状态转移矩阵
Figure 496262DEST_PATH_IMAGE060
 、非线性项状态线性化后的常数项状态转移矩阵的
Figure 919153DEST_PATH_IMAGE062
和非线性项线性化后的一次项状态转移矩阵的
Figure 870054DEST_PATH_IMAGE064
。
计算状态转移矩阵、非线性项状态线性化后的常数项状态转移矩阵的
Figure 192768DEST_PATH_IMAGE062
和非线性项线性化后的一次项状态转移矩阵的
Figure 52140DEST_PATH_IMAGE064
的方法使用的是精细积分算法。
步骤8.3：给定时刻的状态变量值及运行变量值，置迭代次数
Figure 475215DEST_PATH_IMAGE068
。
步骤8.4：预估
Figure 485896DEST_PATH_IMAGE044
时刻的状态变量值。
步骤8.5：根据网络代数方程以及状态变量值计算运行变量值。
按照求解，其中
Figure DEST_PATH_IMAGE106
为电力网络导纳矩阵，首先求解虚拟注入电流
Figure DEST_PATH_IMAGE108
，得到
Figure 300717DEST_PATH_IMAGE070
时刻系统的运行变量，进而得到各发电机的电磁功率
Figure DEST_PATH_IMAGE112
。
步骤8.6根据计算求得的
Figure 199271DEST_PATH_IMAGE044
的运行状态变量值，按下列公式求解
Figure 625311DEST_PATH_IMAGE044
时刻的各台发电机状态变量值
Figure 764169DEST_PATH_IMAGE072
，
其中，式中的
Figure 409914DEST_PATH_IMAGE074
表示
Figure 721946DEST_PATH_IMAGE070
时刻的各台发电机的发电机功角
Figure 769537DEST_PATH_IMAGE046
、角频率
Figure 643077DEST_PATH_IMAGE048
、暂态和次暂态电势、励磁及调速系统各动态环节状态变量组成的状态向量子向量，
Figure 776118DEST_PATH_IMAGE076
为时刻各台发电机的发电机功角
Figure 731622DEST_PATH_IMAGE078
、角频率
Figure 274598DEST_PATH_IMAGE080
、暂态和次暂态电势、励磁及调速系统各动态环节非线性项线性化后的常数项组成的向量，
Figure 133751DEST_PATH_IMAGE082
为
Figure 53165DEST_PATH_IMAGE066
和
Figure 747452DEST_PATH_IMAGE084
时刻各台发电机的发电机功角
Figure 461330DEST_PATH_IMAGE086
、角频率
Figure 568963DEST_PATH_IMAGE088
、暂态和次暂态电势、励磁及调速系统各动态环节非线性项线性化后的一次项部分组成的向量。
步骤8.7：检查两次迭代各发电机最大电磁功率偏差值，若偏差大于给定精度
Figure 527954DEST_PATH_IMAGE090
，令
Figure 139064DEST_PATH_IMAGE092
，返回步骤8.5继续迭代；否则，执行步骤9。
步骤9：判断系统是否稳定，即任意两台发电机的最大相对摇摆功角是否大于某一给定值，若是执行步骤12；否则，执行步骤10。
步骤10：将仿真时间推进一个步长，令
Figure 227106DEST_PATH_IMAGE094
，并以此时的
Figure 556456DEST_PATH_IMAGE096
值作为下一时刻的计算起点值即
Figure 817673DEST_PATH_IMAGE066
。
步骤11：判断是否到达事先给定的仿真时间
Figure 781825DEST_PATH_IMAGE098
。若
Figure 103085DEST_PATH_IMAGE100
则执行步骤12，否则返回步骤5。
步骤12：输出计算结果并结束计算。
以下详细介绍本发明方法的具体过程。
微分方程组（1）主要包括描述发电机组和其它动态装置动态特性的微分方程，其中各台发电机组的微分方程可表示为：
   
Figure DEST_PATH_IMAGE114
                 （2）
式中，分别表示各台发电机的发电机功角、角频率、原动机机械功率、电磁功率及惯性时间常数，
Figure DEST_PATH_IMAGE118
为系统同步电角速度，表示发电机阻尼系数。
Figure DEST_PATH_IMAGE122
为各台发电机的暂态和次暂态电势、励磁及调速系统各动态环节状态变量组成的状态向量子向量，
Figure DEST_PATH_IMAGE124
为与状态向量子向量对应的微分方程右侧的函数向量。这样，每一台发电机组的状态向量
Figure DEST_PATH_IMAGE126
可表示为：。
各台发电机微分方程组(2)就可进一步表示为：
Figure DEST_PATH_IMAGE130
                （3）
式中，
Figure DEST_PATH_IMAGE132
为非线性微分方程组线性部分的系统矩阵，
Figure DEST_PATH_IMAGE134
为发电机的暂态和次暂态电势、励磁及调速系统各动态环节微分方程组线性部分的系统矩阵；
为非线性方程组，非线性部分的函数向量。
非线性微分方程组线性部分的系统矩阵当中的子矩阵的选取规则是：根据每一台发电机组的状态向量和微分方程组，提取出方程组当中发电机组状态向量
Figure 891744DEST_PATH_IMAGE126
对应的各元素的线性项的系数，并列写成矩阵形式。
对非线性函数向量
Figure DEST_PATH_IMAGE138
，用一个近似的线性函数表示：
Figure DEST_PATH_IMAGE140
（4）  
其中
Figure DEST_PATH_IMAGE142
为线性化后的常数项，为线性化后的一次项，具体表示为：
Figure DEST_PATH_IMAGE146
，      
Figure DEST_PATH_IMAGE148
，
其中，
Figure 744031DEST_PATH_IMAGE076
为时刻各台发电机的发电机功角、角频率
Figure 638541DEST_PATH_IMAGE080
、暂态和次暂态电势、励磁及调速系统各动态环节非线性项线性化后的常数项组成的向量，
Figure 284066DEST_PATH_IMAGE082
为
Figure 782044DEST_PATH_IMAGE066
和
Figure 350428DEST_PATH_IMAGE084
时刻各台发电机的发电机功角
Figure 628963DEST_PATH_IMAGE086
、角频率
Figure 573785DEST_PATH_IMAGE088
、暂态和次暂态电势、励磁及调速系统各动态环节非线性项线性化后的一次项部分组成的向量。
Figure DEST_PATH_IMAGE150
为
Figure 552368DEST_PATH_IMAGE066
时刻发电机的暂态和次暂态电势、励磁及调速系统各动态环节非线性项线性化后的常数项组成的向量，
Figure DEST_PATH_IMAGE152
为
Figure 557233DEST_PATH_IMAGE070
时刻发电机的暂态和次暂态电势、励磁及调速系统各动态环节非线性项线性化后的常数项组成的向量, 
Figure DEST_PATH_IMAGE154
为发电机惯性时间常数，
Figure 618336DEST_PATH_IMAGE120
表示各台发电机阻尼系数，为发电机原动机机械功率, 
Figure 101270DEST_PATH_IMAGE118
为系统同步电角速度。
对于以上建立的微分方程如何求解，是本发明的核心所在。对于微分方程
Figure DEST_PATH_IMAGE158
，本发明的精细积分算法主要通过下式：
Figure DEST_PATH_IMAGE160
求得。代入线性化后的
Figure DEST_PATH_IMAGE162
进一步整理可得
Figure DEST_PATH_IMAGE164
，式中
Figure DEST_PATH_IMAGE166
表示线性化后的常数项，
Figure DEST_PATH_IMAGE168
表示线性化后的一次项。
将发电机微分方程代入精细积分计算后可得下列公式：
，
其中，式中的
Figure 875507DEST_PATH_IMAGE074
表示时刻的各台发电机的发电机功角
Figure 87362DEST_PATH_IMAGE046
、角频率
Figure 476755DEST_PATH_IMAGE048
、暂态和次暂态电势、励磁及调速系统各动态环节状态变量组成的状态向量子向量，
Figure 823423DEST_PATH_IMAGE076
为
Figure 799732DEST_PATH_IMAGE066
时刻各台发电机的发电机功角
Figure 624468DEST_PATH_IMAGE078
、角频率
Figure 868368DEST_PATH_IMAGE080
、暂态和次暂态电势、励磁及调速系统各动态环节非线性项线性化后的常数项组成的向量，
Figure 323620DEST_PATH_IMAGE082
为
Figure 285760DEST_PATH_IMAGE066
和
Figure 412722DEST_PATH_IMAGE084
时刻各台发电机的发电机功角
Figure 511128DEST_PATH_IMAGE086
、角频率
Figure 199598DEST_PATH_IMAGE088
、暂态和次暂态电势、励磁及调速系统各动态环节非线性项线性化后的一次项部分组成的向量。这里
Figure DEST_PATH_IMAGE170
对应微分方程中的，
Figure DEST_PATH_IMAGE174
对应
Figure DEST_PATH_IMAGE176
，
Figure DEST_PATH_IMAGE178
对应
Figure DEST_PATH_IMAGE180
为此，根据本发明，一种基于隐式精细积分的电力系统暂态稳定仿真方法每一积分步的计算步骤如下：
1：判断本步的步长
Figure 966478DEST_PATH_IMAGE016
是否与前一时刻步长相同，若相同则跳转到3，若不同，则进行步骤2的操作。
2：根据微分方程得到状态矩阵，利用步长变化后步长
Figure 133017DEST_PATH_IMAGE016
，计算出状态转移矩阵
Figure 351509DEST_PATH_IMAGE060
 、非线性项线性化后的常数项状态转移矩阵的
Figure 210881DEST_PATH_IMAGE062
和非线性项线性化后的线性项状态转移矩阵的
Figure 383498DEST_PATH_IMAGE064
。
3：给定
Figure 619308DEST_PATH_IMAGE066
时刻的状态变量值及运行变量值，置迭代次数
Figure 692306DEST_PATH_IMAGE068
。
4：预估
Figure 660262DEST_PATH_IMAGE044
时刻的状态变量值。
5：根据网络代数方程以及
Figure 818711DEST_PATH_IMAGE044
状态变量值计算运行变量值。
6：根据计算求得的
Figure 825588DEST_PATH_IMAGE044
的运行状态变量值，按公式（5）、（6）和（7）求解
Figure 18671DEST_PATH_IMAGE044
时刻的状态变量值。
7：检查两次迭代各发电机最大电磁功率偏差值，若偏差大于给定精度，令
Figure 803274DEST_PATH_IMAGE092
，返回步骤5继续迭代；否则，本积分步迭代过程结束;
其计算流程如附图2所示。
以下是本发明方法的一个实施例，以IEEE39节点系统进行仿真实验作实施例，进一步说明如下：
仿真算例所采用的发电机模型为发电机三阶模型即
Figure DEST_PATH_IMAGE182
变化模型，发电机附属设备为励磁调节器，负荷采用恒阻抗模型。故障为在34号线路及节点28与节点29之间的线路的始端上在0时刻发生三相短路故障，故障在0.1s时清除。整个故障持续时间为0.1s。仿真整个时间长度为3.0s，系统仿真结果为发生失稳现象。
图3为IEEE39节点系统最大相对功角摇摆曲线，步长为0.01s，采用隐式梯形积分算法，其结果作为标准数据值。图4为隐式梯形积分法最大功角计算误差，步长分别为0.02s和0.05s，图5为隐式精细积分法最大功角计算误差，步长分别为0.02s和0.05s。
由图4和图5可见当积分步长取0.02-0.05s时本发明的最大功角误差都不超过10度。而且在大步长下本发明的误差精度要远小于隐式梯形积分算法。而在计算量上，本发明在步长取0.05s时求解网络代数方程的次数为539次，而隐式梯形积分求解网络代数方程的次数为577次。本发明方法节约了大约6-7%的计算量。

Claims (2)

1. 一种基于隐式精细数值积分的电力系统暂态稳定仿真方法，其特征在于该方法包括如下步骤：
步骤1：输入系统的原始参数和信息，进行潮流计算，得到稳态工况下的运行变量值                                               
Figure 2012100794539100001DEST_PATH_IMAGE002
，包括各台发电机节点的电压，注入网络的电流
Figure 2012100794539100001DEST_PATH_IMAGE006
，电磁功率
Figure 2012100794539100001DEST_PATH_IMAGE008
，其中
Figure DEST_PATH_IMAGE010
，
Figure 7274DEST_PATH_IMAGE012
为发电机台数；
步骤2：计算状态变量初值，包括发电机功角
Figure DEST_PATH_IMAGE014
、角频率
Figure DEST_PATH_IMAGE016
和发电机的暂态和次暂态电势、励磁及调速系统各动态环节状态变量初值；
步骤3：形成描述系统暂态过程的微分方程和网络代数方程，并且进行网络代数方程因子表分解；
步骤4：置暂态稳定计算初值时刻
Figure DEST_PATH_IMAGE018
；
步骤5：判断是否有故障或操作发生；若无，则转向步骤8；若有，则执行步骤6；
步骤6：依据故障或操作情况，修改微分方程和网络代数方程及其因子表；
步骤7：求解网络代数方程，得到
Figure DEST_PATH_IMAGE020
时刻的运行变量；
步骤8：计算
Figure DEST_PATH_IMAGE022
时刻的系统的状态变量值包括发电机功角
Figure DEST_PATH_IMAGE024
、角频率
Figure DEST_PATH_IMAGE026
和发电机的暂态和次暂态电势、励磁及调速系统各动态环节状态变量值，及运行变量值包括发电机节点的电压
Figure DEST_PATH_IMAGE028
，注入网络的电流
Figure DEST_PATH_IMAGE030
和电磁功率
Figure DEST_PATH_IMAGE032
，具体过程如下：
步骤8.1：判断本步的步长
Figure DEST_PATH_IMAGE034
是否与前一时刻步长相同，若相同，则跳转到步骤8.3，若不同，则进行步骤8.2的操作；
步骤8.2：根据微分方程得到状态矩阵
Figure DEST_PATH_IMAGE036
和非线性项
Figure DEST_PATH_IMAGE038
的线性化表达式，利用步长
Figure 372308DEST_PATH_IMAGE034
，计算出状态转移矩阵
Figure DEST_PATH_IMAGE040
、非线性项状态线性化后的常数项状态转移矩阵的
Figure DEST_PATH_IMAGE042
和非线性项线性化后的线性项状态转移矩阵的；
步骤8.3：给定
Figure DEST_PATH_IMAGE046
时刻的状态变量值及运行变量值，置迭代次数
Figure DEST_PATH_IMAGE048
； 
步骤8.4：预估
Figure 55968DEST_PATH_IMAGE022
时刻的状态变量值；
步骤8.5：根据网络代数方程以及
Figure DEST_PATH_IMAGE050
状态变量值计算运行变量值；
步骤8.6：根据计算求得的
Figure 251326DEST_PATH_IMAGE022
的运行状态变量值，按下式求解
Figure 304733DEST_PATH_IMAGE050
时刻的状态变量值
Figure DEST_PATH_IMAGE052
，
其中，式中的
Figure DEST_PATH_IMAGE054
表示
Figure 824140DEST_PATH_IMAGE050
时刻的各台发电机的发电机功角、角频率
Figure 223339DEST_PATH_IMAGE026
、暂态和次暂态电势、励磁及调速系统各动态环节状态变量组成的状态向量子向量，
Figure DEST_PATH_IMAGE056
为
Figure 142754DEST_PATH_IMAGE046
时刻各台发电机的发电机功角
Figure DEST_PATH_IMAGE058
、角频率
Figure DEST_PATH_IMAGE060
、暂态和次暂态电势、励磁及调速系统各动态环节非线性项线性化后的常数项组成的向量，
Figure DEST_PATH_IMAGE062
为
Figure 584843DEST_PATH_IMAGE046
和时刻各台发电机的发电机功角
Figure DEST_PATH_IMAGE066
、角频率
Figure DEST_PATH_IMAGE068
、暂态和次暂态电势、励磁及调速系统各动态环节非线性项线性化后的一次项部分组成的向量；
步骤8.7：检查两次迭代各发电机最大电磁功率偏差值，若偏差大于给定精度
Figure DEST_PATH_IMAGE070
，令
Figure DEST_PATH_IMAGE072
，返回步骤8.5继续迭代；否则，执行步骤9；
步骤9：判断系统是否稳定，即任意两台发电机的最大相对摇摆功角是否大于某一给定值，若是执行步骤12；否则，执行步骤10；
步骤10：将仿真时间推进一个步长，令
Figure DEST_PATH_IMAGE074
，并以此时的
Figure DEST_PATH_IMAGE076
值作为下一时刻的计算起点值，即值；
步骤11：判断是否到达事先给定的仿真时间
Figure DEST_PATH_IMAGE078
；若
Figure DEST_PATH_IMAGE080
则执行步骤12，否则返回步骤5；
步骤12：输出计算结果并结束计算。
2.根据权利要求1所述的一种基于隐式精细数值积分的电力系统暂态稳定仿真方法，其特征在于：暂态稳定计算步骤8.2 中计算状态转移矩阵
Figure 848432DEST_PATH_IMAGE040
 和非线性项线性化后的常数项的状态转移矩阵
Figure 305958DEST_PATH_IMAGE042
和非线性项线性化后的线性项的状态转移矩阵
Figure DEST_PATH_IMAGE082
使用的方法是精细积分法。
CN201210079453.9A 2012-03-23 2012-03-23 一种基于隐式精细数值积分的电力系统暂态稳定仿真方法 Expired - Fee Related CN102664397B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201210079453.9A CN102664397B (zh) 2012-03-23 2012-03-23 一种基于隐式精细数值积分的电力系统暂态稳定仿真方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201210079453.9A CN102664397B (zh) 2012-03-23 2012-03-23 一种基于隐式精细数值积分的电力系统暂态稳定仿真方法

Publications (2)

Publication Number Publication Date
CN102664397A true CN102664397A (zh) 2012-09-12
CN102664397B CN102664397B (zh) 2014-01-29

Family

ID=46773839

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201210079453.9A Expired - Fee Related CN102664397B (zh) 2012-03-23 2012-03-23 一种基于隐式精细数值积分的电力系统暂态稳定仿真方法

Country Status (1)

Country Link
CN (1) CN102664397B (zh)

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105867360A (zh) * 2016-06-14 2016-08-17 江南大学 一种机电控制系统的初值预估迭代学习故障诊断算法
CN107706948A (zh) * 2017-11-03 2018-02-16 华北电力大学 多维阶数控制的多步Taylor级数暂态稳定分析方法
CN105468864B (zh) * 2015-12-14 2018-05-15 三峡大学 基于増维精细积分的高压输电线路电磁暂态数值计算方法
CN110119523A (zh) * 2019-01-29 2019-08-13 中国电力科学研究院有限公司 一种基于功角等效法的网络代数方程求解预处理方法及系统
CN110135031A (zh) * 2019-04-30 2019-08-16 东南大学 基于半隐式龙格库塔法的电力系统暂态稳定计算方法
CN111404825A (zh) * 2020-03-13 2020-07-10 重庆第二师范学院 数据传输方法、系统、计算机设备及存储介质
CN113221298A (zh) * 2021-04-21 2021-08-06 南方电网科学研究院有限责任公司 一种机电暂态过程的仿真方法及系统

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2007288878A (ja) * 2006-04-14 2007-11-01 Hitachi Ltd 電力系統安定度判定方法及び装置
CN101699448A (zh) * 2009-10-26 2010-04-28 清华大学 一种电力系统暂态稳定的分布式仿真方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2007288878A (ja) * 2006-04-14 2007-11-01 Hitachi Ltd 電力系統安定度判定方法及び装置
CN101699448A (zh) * 2009-10-26 2010-04-28 清华大学 一种电力系统暂态稳定的分布式仿真方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
王宇宾等: "基于隐式Taylor级数法的电力系统暂态稳定计算", 《华北电力大学学报》 *

Cited By (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105468864B (zh) * 2015-12-14 2018-05-15 三峡大学 基于増维精细积分的高压输电线路电磁暂态数值计算方法
CN105867360A (zh) * 2016-06-14 2016-08-17 江南大学 一种机电控制系统的初值预估迭代学习故障诊断算法
CN105867360B (zh) * 2016-06-14 2018-05-08 江南大学 一种机电控制系统的初值预估迭代学习故障诊断算法
CN107706948A (zh) * 2017-11-03 2018-02-16 华北电力大学 多维阶数控制的多步Taylor级数暂态稳定分析方法
CN107706948B (zh) * 2017-11-03 2020-12-01 华北电力大学 多维阶数控制的多步Taylor级数暂态稳定分析方法
CN110119523A (zh) * 2019-01-29 2019-08-13 中国电力科学研究院有限公司 一种基于功角等效法的网络代数方程求解预处理方法及系统
CN110135031A (zh) * 2019-04-30 2019-08-16 东南大学 基于半隐式龙格库塔法的电力系统暂态稳定计算方法
CN111404825A (zh) * 2020-03-13 2020-07-10 重庆第二师范学院 数据传输方法、系统、计算机设备及存储介质
CN113221298A (zh) * 2021-04-21 2021-08-06 南方电网科学研究院有限责任公司 一种机电暂态过程的仿真方法及系统
CN113221298B (zh) * 2021-04-21 2023-02-24 南方电网科学研究院有限责任公司 一种机电暂态过程的仿真方法及系统

Also Published As

Publication number Publication date
CN102664397B (zh) 2014-01-29

Similar Documents

Publication Publication Date Title
CN102664397B (zh) 一种基于隐式精细数值积分的电力系统暂态稳定仿真方法
CN102609575B (zh) 一种基于隐式数值积分的电力系统暂态稳定仿真方法
Ghandhari et al. Control Lyapunov functions for controllable series devices
Singh et al. IEEE PES task force on benchmark systems for stability controls report on the 68-bus 16-machine 5-area system
Yao et al. Efficient and robust dynamic simulation of power systems with holomorphic embedding
CN104375876B (zh) 一种输入量突变情况下的0+误差免疫电磁暂态仿真方法
CN106099922B (zh) 基于启发式能量函数法的电力系统暂态电压稳定性判断方法
CN102545263B (zh) 一种基于显式数值积分的电力系统暂态稳定仿真方法
CN105512367A (zh) 并网火电机组一次调频转速不等率的临界稳定值确定方法
WO2018102720A1 (en) System and method for a fast power network simulator
CN112380670A (zh) 基于虚拟动子的分段供电直线感应电机建模方法、系统
Ajala et al. A library of second-order models for synchronous machines
Pandey et al. Unified power system analyses and models using equivalent circuit formulation
CN106208879A (zh) 一种基于自抗扰控制器控制异步电动机变频调速系统的方法
Song et al. Identification of PMSM based on EKF and elman neural network
CN103117692B (zh) 多种外部干扰下的带有机械弹性储能的永磁电动机组控制方法
Al-Hinai Dynamic stability enhancement using genetic algorithm power system stabilizer
Basic et al. Euler-Lagrange models with complex currents of three-phase electrical machines and observability issues
CN102609576B (zh) 预估-校正数值积分的电力系统暂态稳定仿真方法
CN102354332B (zh) 一种用于简化柔性交直流输电系统中rga计算的方法
Ramanujam Power system dynamics: analysis and simulation
Karaagac et al. An efficient voltage-behind-reactance formulation-based synchronous machine model for electromagnetic transients
CN105140957B (zh) 基于风电场和光伏电站聚合模型的机电振荡模式估算方法
Yan et al. Disturbance Observer‐Based Backstepping Control of PMSM for the Mine Traction Electric Locomotive
Hossain et al. Enhancement of transient stability limit and voltage regulation with dynamic loads using robust excitation control

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: 20140129

Termination date: 20210323

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