CN102609576B - 预估-校正数值积分的电力系统暂态稳定仿真方法 - Google Patents

Abstract

本发明公开了一种适用发电机采用二阶模型，基于预估-校正数值积分的电力系统暂态稳定仿真方法。与已有的电力系统暂态稳定数值积分方法相比，该方法将一显式二阶积分用于状态变量的预估计算，将一四阶精度的功角积分用作校正计算。该方法为单步算法，且考虑了角频率变化对转矩的影响。组成的预估-校正算法在较大积分步长(0.1～0.12秒)下分析结果仍有足够精度，从而显著地加快了暂态稳定的仿真速度。

Description

预估-校正数值积分的电力系统暂态稳定仿真方法

技术领域

本发明属于电力系统自动化，特别涉及了一种电力系统暂态稳定仿真方法。

背景技术

电力系统暂态稳定分析是电力系统分析计算中最基础、最核心的内容之一。在线动态安全分析，安全稳定紧急控制、预防控制，智能调度等先进技术已逐步在电力系统中推广应用。实现这些先进技术的前提条件是能够快速、准确、可靠地对大规模电力系统进行的暂态稳定仿真计算。

电力系统暂态稳定分析的常用方法主要有数值积分法，直接法，以及将数值积分和直接法相结合的混合分析方法。数值积分是电力系统暂态稳定计算最准确、最可靠的方法，其它方法是否准确均以数值积分法为标准结果来评判。数值积分法的最大缺点是计算量大，尽管计算机速度已经有了飞速提高，但对于大规模电力系统，计算时间仍难以满足在线动态安全分析、预防控制、紧急控制等计算的要求。

电力系统的暂态过程可用如下形式的微分-代数方程组描述

$\frac{\mathrm{dx}}{\mathrm{dt}}=f(x,y)---\left(1\right)$

0＝g(x，y)                                (2)

式中，x表示微分方程组中描述系统动态特性的状态变量；y表示代数方程组中系统的运行变量。通常向量x包含发电机功角和转速等描述系统中各动态环节的状态变量，而向量y通常包含与网络相关的运行变量，如节点电压的幅值和相位等。

用数值积分法求解电力系统暂态过程的一般流程如图1所示(夏道止电力系统分析下册水利电力出版社1995年11月)。其核心步骤是框⑧所示的在每一积分步根据t时刻的状态变量和运行变量x(t)、y(t)求解(1)、(2)式所表示的微分-代数方程组，得到t+h时刻的状态变量和运行变量x(t+h)、y(t+h)。根据求解微分方程所采用的数值积分方法的不同，以及对微分-代数方程组是联立求解，还是交替求解，构成了不同的暂态稳定数值积分算法。目前，在电力系统数值仿真领域求解(1)式中微分方程组的常用方法有隐式梯形积分法、改进欧拉法、龙格-库塔等。隐式梯形积分数值稳定性好，但需要多次迭代求解，计算量大，目前电力系统商业计算程序BPA、PSASP采用的就是这种积分方法。改进欧拉法和龙格-库塔法为显式积分方法，无需迭代，计算量小，但数值稳定性较差。显式积分算法要根据算法的截断误差，通过选择合理的积分步长，来保证算法的数值稳定性，如在电力系统中广泛应用的PSS/E程序采用的就是改进欧拉法。

为了保证算法的稳定性和仿真精度，所取积分步长要与算法的截断误差成反比，即数值积分算法的截断误差越小，在相同精度要求下，积分步长h可取的大一些，反之积分步长h要取的小一些。通常每一积分步的截断误差越小，计算量也越大。如欧拉法的局部截断误差为O(h²)，每一积分步只需计算一次微分-代数方程；改进欧拉法的局部截断误差为O(h³)，每一积分步需计算两次微分-代数方程；四阶显式龙格-库塔法的局部截断误差为O(h⁵)，每一积分步需计算四次微分-代数方程。而隐式梯形积分法的局部截断误差为O(h³)，则需经过多次迭代求解微分-代数方程，才能得到满足精度要求的解。若能在减少算法截断误差的同时，不增加算法的计算量，则能减少整个暂态仿真的计算量，加快计算速度。

目前，在电力系统暂态稳定数值积分方法中所采用的积分方法，均直接采用计算方法理论中的通用算法，如隐式梯形积分法、改进欧拉法、龙格-库塔以及其他方法，并没有根据描述电力系统暂态过程的微分方程的特点对算法进行改进。

发明内容

本发明目的是为了解决电力系统暂态稳定分析计算中，现有的数值积分方法计算量大，计算速度不能满足电力系统在线计算要求的这个问题，提出了一种基于预估-校正数值积分的暂态稳定仿真方法。

本发明目的是通过以下技术方案实现的：预估-校正数值积分的电力系统暂态稳定仿真方法，包括以下步骤：

步骤1：输入系统的原始参数和信息，进行潮流计算得到稳态工况下的运行变量值y(0)，包括发电机节点电压V(0)，注入网络的电流I(0)及各发电机电磁功率P_ei(0)，其中i＝1，2，…N_G(N_G为发电机台数)；

步骤2：计算状态变量功角的初值δ_i(0)、角频率的初值ω_i(0)；

步骤3：形成描述系统暂态过程的微分方程和网络代数方程，并对网络代数方程进行因子表分解；

步骤4：置暂态稳定计算初值时刻t＝0，确定暂态稳定计算采用的积分步长h，进行暂态稳定仿真计算；

步骤5：判断是否有故障或操作发生。若无，则转向步骤8；若有则执行步骤6；

步骤6：依据故障或操作情况，修改网络代数方程的因子表；

步骤7：求解网络代数方程，得到t时刻的运行变量；

步骤8：计算t+h时刻的系统的状态变量值，运行变量值，本步骤具体过程如下：

步骤8.1：根据t时刻各发电机的电磁功率P_ei(t)和状态变量功角δ_i(t)、角频率ω_i(t)，按下式得到

t+h时刻的功角的预估值

${\mathrm{\δ}}_i^{\left(0\right)}(t+h)={\mathrm{\δ}}_i\left(t\right)+h({\mathrm{\ω}}_i\left(t\right)-1){\mathrm{\ω}}_s+\frac{{\mathrm{\ω}}_s{h}^2}{2T_{\mathrm{Jj}}}(\frac{P_{\mathrm{mi}}}{{\mathrm{\ω}}_i\left(t\right)}-\frac{P_{\mathrm{ei}}\left(t\right)}{{\mathrm{\ω}}_i\left(t\right)})$

其中，ω_s为系统同步角速度，D_i、T_Ji、P_mi分别为各台发电机阻尼系数、惯性时间常数和各台发电机机械功率；

步骤8.2：求解网络代数方程YV⁽⁰⁾(t+h)＝I′(δ⁽⁰⁾(t+h)，V⁽⁰⁾(t+h))，并计算各发电机电磁功率在t+h时刻的预估值

其中I′(δ⁽⁰⁾(t+h)，V⁽⁰⁾(t+h))为t+h时刻节点虚拟注入电流预估值，Y为节点导纳矩阵，V⁽⁰⁾(t+h)为t+h时刻节点电压预估值。

步骤8.3：按如下式计算得到t+h时刻各发电机角频率的预估值

${\mathrm{\ω}}_i^{\left(0\right)}(t+h)=\frac{1}{2}({\mathrm{\ω}}_i\left(t\right)+\frac{h(P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right))}{2T_{\mathrm{Jj}}{\mathrm{\ω}}_i\left(t\right)}+\sqrt{{({\mathrm{\ω}}_i\left(t\right)+\frac{h(P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right))}{2T_{\mathrm{Jj}}{\mathrm{\ω}}_i\left(t\right)})}^2+\frac{2h}{T_{\mathrm{Jj}}}(P_{\mathrm{mi}}-P_{\mathrm{ei}}^{\left(0\right)}(t+h))})$

此时得到了t+h时刻各发电机运行变量和状态变量的预估值

和 ${\mathrm{\ω}}_i^{\left(0\right)}(t+h);$

步骤8.4：按如下积分公式

${\mathrm{\δ}}_i(t+h)={\mathrm{\δ}}_i\left(t\right)+\frac{{\mathrm{\ω}}_s}{2}h({\mathrm{\ω}}_i\left(t\right)+{\mathrm{\ω}}_i^{\left(0\right)}(t+h)-2)+\frac{1}{12T_{\mathrm{Jj}}}{h}^2{\mathrm{\ω}}_s(\frac{P_{\mathrm{mi}}-P_e\left(t\right)}{{\mathrm{\ω}}_i\left(t\right)}-\frac{P_{\mathrm{mi}}-P_{\mathrm{ei}}^{\left(0\right)}(t+h)}{{\mathrm{\ω}}_i^{\left(0\right)}(t+h)})$

求出发电机功角δ_i(t+h)；

步骤8.5：求解网络方程

YV(t+h)＝I′(δ(t+h)，V(t+h))

并计算各发电机在t+h时刻的电磁功率P_ei(t+h)；其中I′(δ(t+h)，V(t+h))为t+h时刻节点虚拟注入电流，Y为节点导纳矩阵，V(t+h)为t+h时刻节点电压：

步骤8.6：按如下式计算得到t+h时刻各发电机的角频率ω_i(t+h)

${\mathrm{\ω}}_i(t+h)=\frac{1}{2}({\mathrm{\ω}}_i\left(t\right)+\frac{h(P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right))}{2T_{\mathrm{Jj}}{\mathrm{\ω}}_i\left(t\right)}+\sqrt{{({\mathrm{\ω}}_i\left(t\right)+\frac{h(P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right))}{2T_{\mathrm{Jj}}{\mathrm{\ω}}_i\left(t\right)})}^2+\frac{2h}{T_{\mathrm{Jj}}}(P_{\mathrm{mi}}-P_{\mathrm{ei}}^{\left(0\right)}(t+h))});$

步骤9：判断系统是否稳定，即任意两台发电机的最大相对摇摆功角是否大于某一给定值，若是，执行步骤12；否则，执行步骤10；

步骤10：将仿真时间推进一个步长，令t＝t+h；

步骤11：判断是否到达事先给定的仿真时间T。若t≥T则执行步骤12，否则返回步骤5；

步骤12：输出计算结果并结束计算。

在暂态稳定仿真过程中，在预估步用下式计算状态变量预估值：

${\mathrm{\δ}}_i^{\left(0\right)}(t+h)={\mathrm{\δ}}_i\left(t\right)+h({\mathrm{\ω}}_i\left(t\right)-1){\mathrm{\ω}}_s+\frac{{\mathrm{\ω}}_s{h}^2}{2T_{\mathrm{Jj}}}(\frac{P_{\mathrm{mi}}}{{\mathrm{\ω}}_i\left(t\right)}-\frac{P_{\mathrm{ei}}\left(t\right)}{{\mathrm{\ω}}_i\left(t\right)})$

${\mathrm{\ω}}_i^{\left(0\right)}(t+h)=\frac{1}{2}({\mathrm{\ω}}_i\left(t\right)+\frac{h(P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right))}{2T_{\mathrm{Jj}}{\mathrm{\ω}}_i\left(t\right)}+\sqrt{{({\mathrm{\ω}}_i\left(t\right)+\frac{h(P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right))}{2T_{\mathrm{Jj}}{\mathrm{\ω}}_i\left(t\right)})}^2+\frac{2h}{T_{\mathrm{Jj}}}(P_{\mathrm{mi}}-P_{\mathrm{ei}}^{\left(0\right)}(t+h))})$

在校正步用下式计算状态变量：

${\mathrm{\δ}}_i(t+h)={\mathrm{\δ}}_i\left(t\right)+\frac{{\mathrm{\ω}}_s}{2}h({\mathrm{\ω}}_i\left(t\right)+{\mathrm{\ω}}_i^{\left(0\right)}(t+h)-2)+\frac{1}{12T_{\mathrm{Jj}}}{h}^2{\mathrm{\ω}}_s(\frac{P_{\mathrm{mi}}-P_e\left(t\right)}{{\mathrm{\ω}}_i\left(t\right)}-\frac{P_{\mathrm{mi}}-P_{\mathrm{ei}}^{\left(0\right)}(t+h)}{{\mathrm{\ω}}_i^{\left(0\right)}(t+h)})$

${\mathrm{\ω}}_i(t+h)=\frac{1}{2}({\mathrm{\ω}}_i\left(t\right)+\frac{h(P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right))}{2T_{\mathrm{Jj}}{\mathrm{\ω}}_i\left(t\right)}+\sqrt{{({\mathrm{\ω}}_i\left(t\right)+\frac{h(P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right))}{2T_{\mathrm{Jj}}{\mathrm{\ω}}_i\left(t\right)})}^2+\frac{2h}{T_{\mathrm{Jj}}}(P_{\mathrm{mi}}-P_{\mathrm{ei}}(t+h))})$

本发明的有益效果：本发明方法充分利用了发电机二阶模型的特点，用高精度的积分公式组成预估-校正积分，可在较大积分步长(0.1～0.12秒)下分析结果仍有足够精度，从而显著地加快了暂态稳定的仿真速度。

附图说明

图1为暂态稳定数值仿真的一般流程图；

图2为预估-校正数值积分暂态稳定计算每一积分步的计算流程；

图3为步长h＝0.01秒时相对摇摆角最大的两台发电机间的转角差；

图4为本发明方法在不同步长下最大相对摇摆角误差曲线；

图5为改进欧拉法在不同步长下最大相对摇摆角误差曲线。

具体实施方式

以下结合附图对本发明作进一步说明。

本发明提出的一种基于预估-校正数值积分的电力系统暂态稳定仿真方法，包括以下步骤：

步骤1：输入系统的原始参数和信息，进行潮流计算得到稳态工况下的运行变量值y(0)，包括发电机节点电压V(0)，注入网络的电流I(0)及各发电机电磁功率P_ei(0)，其中i＝1，2，…N_G(N_G为发电机台数)；

步骤2：计算状态变量功角的初值δ_i(0)、角频率的初值ω_i(0)；

步骤3：形成描述系统暂态过程的微分方程和网络代数方程，并对网络代数方程进行因子表分解；

步骤4：置暂态稳定计算初值时刻t＝0，确定暂态稳定计算采用的积分步长h，进行暂态稳定仿真计算；

步骤5：判断是否有故障或操作发生。若无，则转向步骤8；若有则执行步骤6；

步骤6：依据故障或操作情况，修改网络代数方程的因子表；

步骤7：求解网络代数方程，得到t时刻的运行变量；

步骤8：计算t+h时刻的系统的状态变量值，运行变量值，本步骤具体过程如下：

步骤8.1：根据t时刻各发电机的电磁功率P_ei(t)和状态变量功角δ_i(t)、角频率ω_i(t)，按下式得到

t+h时刻的功角的预估值

${\mathrm{\δ}}_i^{\left(0\right)}(t+h)={\mathrm{\δ}}_i\left(t\right)+h({\mathrm{\ω}}_i\left(t\right)-1){\mathrm{\ω}}_s+\frac{{\mathrm{\ω}}_s{h}^2}{2T_{\mathrm{Jj}}}(\frac{P_{\mathrm{mi}}}{{\mathrm{\ω}}_i\left(t\right)}-\frac{P_{\mathrm{ei}}\left(t\right)}{{\mathrm{\ω}}_i\left(t\right)})$

其中，ω_s为系统同步角速度，D_i、T_Ji、P_mi分别为各台发电机阻尼系数、惯性时间常数和各台发电机机械功率；

步骤8.2：求解网络代数方程YV⁽⁰⁾(t+h)＝I′(δ⁽⁰⁾(t+h)，V⁽⁰⁾(t+h))，并计算各发电机电磁功率在t+h时刻的预估值

其中I′(δ⁽⁰⁾(t+h)，V⁽⁰⁾(t+h))为t+h时刻节点虚拟注入电流预估值，Y为节点导纳矩阵，V⁽⁰⁾(t+h)为t+h时刻节点电压预估值。

步骤8.3：按如下式计算得到t+h时刻各发电机角频率的预估值

${\mathrm{\ω}}_i^{\left(0\right)}(t+h)=\frac{1}{2}({\mathrm{\ω}}_i\left(t\right)+\frac{h(P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right))}{2T_{\mathrm{Jj}}{\mathrm{\ω}}_i\left(t\right)}+\sqrt{{({\mathrm{\ω}}_i\left(t\right)+\frac{h(P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right))}{2T_{\mathrm{Jj}}{\mathrm{\ω}}_i\left(t\right)})}^2+\frac{2h}{T_{\mathrm{Jj}}}(P_{\mathrm{mi}}-P_{\mathrm{ei}}^{\left(0\right)}(t+h))})$

此时得到了t+h时刻各发电机运行变量和状态变量的预估值

和 ${\mathrm{\ω}}_i^{\left(0\right)}(t+h);$

步骤8.4：按如下积分公式

${\mathrm{\δ}}_i(t+h)={\mathrm{\δ}}_i\left(t\right)+\frac{{\mathrm{\ω}}_s}{2}h({\mathrm{\ω}}_i\left(t\right)+{\mathrm{\ω}}_i^{\left(0\right)}(t+h)-2)+\frac{1}{12T_{\mathrm{Jj}}}{h}^2{\mathrm{\ω}}_s(\frac{P_{\mathrm{mi}}-P_e\left(t\right)}{{\mathrm{\ω}}_i\left(t\right)}-\frac{P_{\mathrm{mi}}-P_{\mathrm{ei}}^{\left(0\right)}(t+h)}{{\mathrm{\ω}}_i^{\left(0\right)}(t+h)})$

求出发电机功角δ_i(t+h)；

步骤8.5：求解网络方程

YV(t+h)＝I′(δ(t+h)，V(t+h))

并计算各发电机在t+h时刻的电磁功率P_ei(t+h)；其中I′(δ(t+h)，V(t+h))为t+h时刻节点虚拟注入电流，Y为节点导纳矩阵，V(t+h)为t+h时刻节点电压：

步骤8.6按如下式计算得到t+h时刻各发电机的角频率ω_i(t+h)

${\mathrm{\ω}}_i(t+h)=\frac{1}{2}({\mathrm{\ω}}_i\left(t\right)+\frac{h(P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right))}{2T_{\mathrm{Jj}}{\mathrm{\ω}}_i\left(t\right)}+\sqrt{{({\mathrm{\ω}}_i\left(t\right)+\frac{h(P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right))}{2T_{\mathrm{Jj}}{\mathrm{\ω}}_i\left(t\right)})}^2+\frac{2h}{T_{\mathrm{Jj}}}(P_{\mathrm{mi}}-P_{\mathrm{ei}}(t+h))})$

步骤9：判断系统是否稳定，即任意两台发电机的最大相对摇摆功角是否大于某一给定值，若是，执行步骤12；否则，执行步骤10；

步骤10：将仿真时间推进一个步长，令t＝t+h；

步骤11：判断是否到达事先给定的仿真时间T。若t≥T则执行步骤12，否则返回步骤5；

步骤12：输出计算结果并结束计算。在暂态稳定仿真过程中，在预估步用下式计算状态变量预估值：

${\mathrm{\δ}}_i^{\left(0\right)}(t+h)={\mathrm{\δ}}_i\left(t\right)+h({\mathrm{\ω}}_i\left(t\right)-1){\mathrm{\ω}}_s+\frac{{\mathrm{\ω}}_s{h}^2}{2T_{\mathrm{Jj}}}(\frac{P_{\mathrm{mi}}}{{\mathrm{\ω}}_i\left(t\right)}-\frac{P_{\mathrm{ei}}\left(t\right)}{{\mathrm{\ω}}_i\left(t\right)})$

${\mathrm{\ω}}_i^{\left(0\right)}(t+h)=\frac{1}{2}({\mathrm{\ω}}_i\left(t\right)+\frac{h(P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right))}{2T_{\mathrm{Jj}}{\mathrm{\ω}}_i\left(t\right)}+\sqrt{{({\mathrm{\ω}}_i\left(t\right)+\frac{h(P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right))}{2T_{\mathrm{Jj}}{\mathrm{\ω}}_i\left(t\right)})}^2+\frac{2h}{T_{\mathrm{Jj}}}(P_{\mathrm{mi}}-P_{\mathrm{ei}}^{\left(0\right)}(t+h))})$

在校正步用下式计算状态变量：

${\mathrm{\δ}}_i(t+h)={\mathrm{\δ}}_i\left(t\right)+\frac{{\mathrm{\ω}}_s}{2}h({\mathrm{\ω}}_i\left(t\right)+{\mathrm{\ω}}_i^{\left(0\right)}(t+h)-2)+\frac{1}{12T_{\mathrm{Jj}}}{h}^2{\mathrm{\ω}}_s(\frac{P_{\mathrm{mi}}-P_e\left(t\right)}{{\mathrm{\ω}}_i\left(t\right)}-\frac{P_{\mathrm{mi}}-P_{\mathrm{ei}}^{\left(0\right)}(t+h)}{{\mathrm{\ω}}_i^{\left(0\right)}(t+h)})$

${\mathrm{\ω}}_i(t+h)=\frac{1}{2}({\mathrm{\ω}}_i\left(t\right)+\frac{h(P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right))}{2T_{\mathrm{Jj}}{\mathrm{\ω}}_i\left(t\right)}+\sqrt{{({\mathrm{\ω}}_i\left(t\right)+\frac{h(P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right))}{2T_{\mathrm{Jj}}{\mathrm{\ω}}_i\left(t\right)})}^2+\frac{2h}{T_{\mathrm{Jj}}}(P_{\mathrm{mi}}-P_{\mathrm{ei}}(t+h))})$

以下详细介绍本发明方法的具体过程。

当发电机采用二阶模型时，方程(1)中每一台发电机组的微分方程可如下表示：

$\left[\begin{array}{c}\frac{d{\mathrm{\δ}}_i}{\mathrm{dt}}\\ \frac{d{\mathrm{\ω}}_i}{\mathrm{dt}}\end{array}\right]=\left[\begin{array}{c}{f}_{i1}\left(t\right)\\ f_{i2}\left(t\right)\end{array}\right]=\left[\begin{array}{c}{\mathrm{\ω}}_s({\mathrm{\ω}}_i\left(t\right)-1)\\ \frac{1}{T_{J_i}}(\frac{P_{\mathrm{mi}}\left(t\right)}{{\mathrm{\ω}}_i\left(t\right)}-\frac{P_{\mathrm{ei}}\left(t\right)}{{\mathrm{\ω}}_i\left(t\right)})\end{array}\right]---\left(3\right)$

式中，δ_i、ω_i、P_mi(t)、P_ei(t)、T_Ji分别表示第i台发电机的功角、角频率、机械功率、电磁功率及惯性时间常数，ω_s为系统同步角频率。

微分方程组(3)可进一步表示为：

$\left[\begin{array}{c}\frac{d{\mathrm{\δ}}_i}{\mathrm{dt}}\\ \frac{d{\mathrm{\ω}}_i}{\mathrm{dt}}\end{array}\right]=\left[\begin{array}{cc}0& c{\mathrm{\ω}}_s\\ 0& 0\end{array}\right]\left[\begin{array}{c}{\mathrm{\δ}}_i\\ {\mathrm{\ω}}_i\end{array}\right]+\left[\begin{array}{c}-{\mathrm{\ω}}_s+{\mathrm{\ω}}_s(1-c){\mathrm{\ω}}_i\left(t\right)\\ \frac{1}{T_{J_i}}(\frac{P_{\mathrm{mi}}\left(t\right)}{{\mathrm{\ω}}_i\left(t\right)}-\frac{P_{\mathrm{ei}}\left(t\right)}{{\mathrm{\ω}}_i\left(t\right)})\end{array}\right]=H{x}_i+F_i(x_i\left(t\right),y\left(t\right))---\left(4\right)$

式中， $H=\left[\begin{array}{cc}0& c{\mathrm{\ω}}_s\\ 0& 0\end{array}\right]$ 为非线性微分方程组线性部分的系统矩阵；

$F_i(x_i\left(t\right),y\left(t\right))=\left[\begin{array}{c}-{\mathrm{\ω}}_s+{\mathrm{\ω}}_s(1-c){\mathrm{\ω}}_i\left(t\right)\\ \frac{1}{T_{J_i}}(\frac{P_{\mathrm{mi}}\left(t\right)}{{\mathrm{\ω}}_i\left(t\right)}-\frac{P_{\mathrm{ei}}\left(t\right)}{{\mathrm{\ω}}_i\left(t\right)})\end{array}\right]$ 为非线性方程组，非线性部分的函数向量。

线性系统矩阵H是一奇异矩阵。显然有：

H^k＝0，k≥2

这样，可得到其所对应的状态转移矩阵准确的解析表达式为：

$e^{\mathrm{Ht}}=\underset{k=0}{\overset{\∞}{\mathrm{\Σ}}}{\left(\mathrm{Ht}\right)}^k=I+\mathrm{Ht}=\left[\begin{array}{cc}1& c{\mathrm{\ω}}_{s}t\\ 0& 1\end{array}\right]---\left(5\right)$

对非线性函数向量F_i(x_i(t)，y(t))，用一个近似的线性函数表示：

$\begin{array}{cc}{F}_i(x_i\left(t\right),y\left(t\right))\≈F_i(x_i\left(t_k\right),y\left(t_k\right))+\frac{1}{h}[F_i(x_i(t_k+h),y(t_k+h))-F_i(x_i\left(t_k\right),y\left(t_k\right))]t,& t\∈[t_k,t_k+h]\end{array}---\left(6\right)$

将(5)、(6)式代入Duhamel积分方程：

当取c＝1时，可得：

${\mathrm{\δ}}_i(t_k+h)={\mathrm{\δ}}_i\left(t_k\right)+\frac{{\mathrm{\ω}}_s}{2}h({\mathrm{\ω}}_i\left(t_k\right)-2+{\mathrm{\ω}}_i(t_k+h))+\frac{1}{12T_J}{h}^2{\mathrm{\ω}}_s(\frac{P_{\mathrm{mi}}\left(t_k\right)-P_e\left(t_k\right)}{{\mathrm{\ω}}_i\left(t_k\right)}-\frac{P_{\mathrm{mi}}(t_k+h)-P_e(t_k+h)}{{\mathrm{\ω}}_i(t_k+h)})---\left(7\right)$

${\mathrm{\ω}}_i(t+h)=\frac{1}{2}({\mathrm{\ω}}_i\left(t\right)+\frac{h(P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right))}{2T_{\mathrm{Jj}}{\mathrm{\ω}}_i\left(t\right)}+\sqrt{{({\mathrm{\ω}}_i\left(t\right)+\frac{h(P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right))}{2T_{\mathrm{Jj}}{\mathrm{\ω}}_i\left(t\right)})}^2+\frac{2h}{T_{\mathrm{Jj}}}(P_{\mathrm{mi}}-P_{\mathrm{ei}}(t+h))})---\left(8\right)$

式(7)中功角积分公式的局部截断误差为O(h⁵)。

当取c＝3时，可得：

${\mathrm{\δ}}_i(t_k+h)={\mathrm{\δ}}_i\left(t_k\right)+h({\mathrm{\ω}}_i\left(t_k\right)-1){\mathrm{\ω}}_s+\frac{{\mathrm{\ω}}_s{h}^2}{2T_{J_i}}(\frac{P_{\mathrm{mi}}\left(t_k\right)}{{\mathrm{\ω}}_i\left(t_k\right)}-\frac{P_{\mathrm{ei}}\left(t_k\right)}{{\mathrm{\ω}}_i\left(t_k\right)})---\left(9\right)$

为一显式积分公式，其局部截断误差为O(h³)

为此，根据本发明，基于预估-校正数值积分的电力系统暂态稳定数值仿真方法每一积分步的计算步骤如下：

1.根据t时刻各发电机的电磁功率P_ei(t)和状态变量功角δ_i(t)、角频率ω_i(t)，按(9)式得到t+h时刻的功角的预估值δ⁽⁰⁾(t+h)；

2.求解网络代数方程并计算各发电机电磁功率在t+h时刻的预估值

$P_{\mathrm{ei}}^{\left(0\right)}(t+h);$

3.按(8)式计算得到t+h时刻各发电机角频率的预估值

4.按积分公式(7)求出发电机功角δ_i(t+h)；

5.求解网络方程并计算各发电机在t+h时刻的电磁功率P_ei(t+h)；

6.按(8)式计算得到t+h时刻各发电机的角频率ω_i(t+h)；

此时得到了t+h时刻各发电机运行变量和状态变量值P_ei(t+h)、δ_i(t+h)和ω_i(t+h)；将本发明所提出的电力系统暂态稳定计算每一积分步的计算流程(图2)嵌入暂态稳定计算一般流程(图1)的框8，就可实现基于二阶发电机模型的局部截断误差为O(h³)的预估-校正数数值积分的电力系统暂态稳定仿真方法。

以下是本发明方法的一个实施例，以IEEE145节点系统进行仿真实验作实施例，进一步说明如下：

IEEE145系统中所有发电机均采用E′恒定模型，负荷采用恒定阻抗模型。在节点131与节点130线路的始端t＝0s发生三相短路故障，t＝0.1s切除故障线路。图3为积分步长取h＝0.01s，用隐式梯形积分法得出的发电机间最大相对摇摆角曲线，并以此作为准确结果。

本发明提出的预估-校正算法与改进欧拉法在一个积分步的计算量相同，均为求解二次微分-代数方程组。图4、图5分别给出了本发明所述方法和改进欧拉法在不同积分步长下计算结果与图3所示标准结果的偏差。由图4和图5可见，本发明方法在积分步长取0.12秒时仍然保持满意的精度，而改进欧拉法当步长大于0.08秒时即趋于发散，故采用本发明方法进行暂态稳定仿真时可采用较大步长，从而加快了仿真速度。

Claims (3)

1.预估-校正数值积分的电力系统暂态稳定仿真方法，该方法适用于发电机采用二阶模型时的情形，其特征在于该方法包括以下步骤：

步骤1：输入系统的原始参数和信息，进行潮流计算得到稳态工况下的运行变量值y(0)，包括发电机节点电压V(0)，注入网络的电流I(0)及各发电机电磁功率P_ei(0)，其中i=1，2，…N_G，其中N_G为发电机台数；

步骤2：计算状态变量功角的初值δ_i(0)、角频率的初值ω_i(0)；

步骤3：形成描述系统暂态过程的微分方程和网络代数方程，并对网络代数方程进行因子表分解；

步骤4：置暂态稳定计算初值时刻t=0，确定暂态稳定计算采用的积分步长h，进行暂态稳定仿真计算；

步骤5：判断是否有故障或操作发生；若无，则转向步骤8；若有则执行步骤6；

步骤6：依据故障或操作情况，修改网络代数方程的因子表；

步骤7：求解网络代数方程，得到t时刻的运行变量；

步骤8：计算t+h时刻的系统的状态变量值，运行变量值，本步骤具体过程如下：

步骤8.1：根据t时刻各发电机的电磁功率P_ei(t)和状态变量功角δ_i(t)、角频率ω_i(t)，按下式得到

t+h时刻的功角的预估值

${\mathrm{\δ}}_i^{\left(0\right)}(t+h)={\mathrm{\δ}}_i\left(t\right)+h({\mathrm{\ω}}_i\left(t\right)-1){\mathrm{\ω}}_s+\frac{{\mathrm{\ω}}_s{h}^2}{2T_{\mathrm{Ji}}}(\frac{P_{\mathrm{mi}}}{{\mathrm{\ω}}_i\left(t\right)}-\frac{P_{\mathrm{ei}}\left(t\right)}{{\mathrm{\ω}}_i\left(t\right)})$

其中，ω_s为系统同步角速度，T_Ji、P_mi分别为各台发电机惯性时间常数和各台发电机机械功率；

步骤8.2：求解网络代数方程YV⁽⁰⁾(t+h)=I'(δ⁽⁰⁾(t+h),V⁽⁰⁾(t+h))，并计算各发电机电磁功率在t+h时刻的预估值

其中I'(δ⁽⁰⁾(t+h),V⁽⁰⁾(t+h))为t+h时刻节点虚拟注入电流预估值，Y为节点导纳矩阵，V⁽⁰⁾(t+h)为t+h时刻节点电压预估值；

步骤8.3：按如下式计算得到t+h时刻各发电机角频率的预估值

${\mathrm{\ω}}_i^{\left(0\right)}(t+h)=\frac{1}{2}({\mathrm{\ω}}_i\left(t\right)+\frac{h(P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right))}{2T_{\mathrm{Ji}}{\mathrm{\ω}}_i\left(t\right)}+\sqrt{{({\mathrm{\ω}}_i\left(t\right)+\frac{h(P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right))}{2T_{\mathrm{Ji}}{\mathrm{\ω}}_i\left(t\right)})}^2+\frac{2h}{T_{\mathrm{Ji}}}(P_{\mathrm{mi}}-P_{\mathrm{ei}}^{\left(0\right)}(t+h))});$

此时得到了t+h时刻各发电机运行变量

和状态变量的预估值

和

步骤8.4：按如下积分公式

${\mathrm{\δ}}_i(t+h)={\mathrm{\δ}}_i\left(t\right)+\frac{{\mathrm{\ω}}_s}{2}h({\mathrm{\ω}}_i\left(t\right)+{\mathrm{\ω}}_i^{\left(0\right)}(t+h)-2)+\frac{1}{12T_{\mathrm{Ji}}}{h}^2{\mathrm{\ω}}_s(\frac{P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right)}{{\mathrm{\ω}}_i\left(t\right)}-\frac{P_{\mathrm{mi}}-P_{\mathrm{ei}}^{\left(0\right)}(t+h)}{{\mathrm{\ω}}_i^{\left(0\right)}(t+h)})$

求出发电机功角δ_i(t+h)，其中P_ei(t)表示第i台发电机电磁功率；

步骤8.5：求解网络方程

YV(t+h)=I'(δ(t+h),V(t+h))

并计算各发电机在t+h时刻的电磁功率P_ei(t+h)；其中I'(δ(t+h),V(t+h))为t+h时刻节点虚拟注入电流，Y为节点导纳矩阵，V(t+h)为t+h时刻节点电压；

步骤8.6按如下式计算得到t+h时刻各发电机的角频率ω_i(t+h)

${\mathrm{\ω}}_i(t+h)=\frac{1}{2}({\mathrm{\ω}}_i\left(t\right)+\frac{h(P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right))}{2T_{\mathrm{Ji}}{\mathrm{\ω}}_i\left(t\right)}+\sqrt{{({\mathrm{\ω}}_i\left(t\right)+\frac{h(P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right))}{2T_{\mathrm{Ji}}{\mathrm{\ω}}_i\left(t\right)})}^2+\frac{2h}{T_{\mathrm{Ji}}}(P_{\mathrm{mi}}-P_{\mathrm{ei}}(t+h))});$

步骤9：判断系统是否稳定，即任意两台发电机的最大相对摇摆功角是否大于某一给定值，若是，执行步骤12；否则，执行步骤10；

步骤10：将仿真时间推进一个步长，令t=t+h；

步骤11：判断是否到达事先给定的仿真时间T；若t≥T则执行步骤12，否则返回步骤5；

步骤12：输出计算结果并结束计算。

2.根据权利要求1所述的预估-校正数值积分的电力系统暂态稳定仿真方法，其特征在于：在预估步用下式计算状态变量预估值：

${\mathrm{\δ}}_i^{\left(0\right)}(t+h)={\mathrm{\δ}}_i\left(t\right)+h({\mathrm{\ω}}_i\left(t\right)-1){\mathrm{\ω}}_s+\frac{{\mathrm{\ω}}_s{h}^2}{{2T}_{\mathrm{Ji}}}(\frac{P_{\mathrm{mi}}}{{\mathrm{\ω}}_i\left(t\right)}-\frac{P_{\mathrm{ei}}\left(t\right)}{{\mathrm{\ω}}_i\left(t\right)})$

${\mathrm{\ω}}_i^{\left(0\right)}(t+h)=\frac{1}{2}({\mathrm{\ω}}_i\left(t\right)+\frac{h(P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right))}{2T_{\mathrm{Ji}}{\mathrm{\ω}}_i\left(t\right)}+\sqrt{{({\mathrm{\ω}}_i\left(t\right)+\frac{h(P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right))}{2T_{\mathrm{Ji}}{\mathrm{\ω}}_i\left(t\right)})}^2+\frac{2h}{T_{\mathrm{Ji}}}(P_{\mathrm{mi}}-P_{\mathrm{ei}}^{\left(0\right)}(t+h))}).$

3.根据权利要求1所述的预估-校正数值积分的电力系统暂态稳定仿真方法，其特征在于：在校正步用下式计算状态变量：

${\mathrm{\δ}}_i(t+h)={\mathrm{\δ}}_i\left(t\right)+\frac{{\mathrm{\ω}}_s}{2}h({\mathrm{\ω}}_i\left(t\right)+{\mathrm{\ω}}_i^{\left(0\right)}(t+h)-2)+\frac{1}{12T_{\mathrm{Ji}}}{h}^2{\mathrm{\ω}}_s(\frac{P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right)}{{\mathrm{\ω}}_i\left(t\right)}-\frac{P_{\mathrm{mi}}-P_{\mathrm{ei}}^{\left(0\right)}(t+h)}{{\mathrm{\ω}}_i^{\left(0\right)}(t+h)})$

${\mathrm{\ω}}_i(t+h)=\frac{1}{2}({\mathrm{\ω}}_i\left(t\right)+\frac{h(P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right))}{2T_{\mathrm{Ji}}{\mathrm{\ω}}_i\left(t\right)}+\sqrt{{({\mathrm{\ω}}_i\left(t\right)+\frac{h(P_{\mathrm{mi}}-P_{\mathrm{ei}}\left(t\right))}{2T_{\mathrm{Ji}}{\mathrm{\ω}}_i\left(t\right)})}^2+\frac{2h}{T_{\mathrm{Ji}}}(P_{\mathrm{mi}}-P_{\mathrm{ei}}(t+h))}).$

Images