CN105403834B - 一种发电机动态状态估计方法 - Google Patents

Abstract

本发明公开了一种发电机动态状态估计方法，该方法将时变多维观测噪声尺度因子引入到CKF动态状态估计中，根据量测新息对量测误差进行在线调整，使其更加逼近真实噪声。再利用调整后的误差计算滤波增益，使其能够在量测量含有坏数据的情况下对状态量预报值进行准确修正，得到精确的状态量估计值；同时，还给出了时变多维观测噪声尺度因子的具体构造方法；并针对滤波增益求逆发生奇异的问题，提出解决方案；通过采用本发明公开的方法，在保证实时性需求的前提下，能够有效抑制量测坏数据对发电机动态状态估计的影响。

Description

一种发电机动态状态估计方法

技术领域

本发明涉及电力***动态状态估计技术领域，尤其涉及一种发电机动态状态估计方法。

背景技术

同步相量测量单元(phasor measurement unit，PMU)的出现为电力***暂态稳定分析与控制提供了新的技术手段。然而，在干扰、测量设备故障、同步信号丢失等情况下，往往导致坏数据的出现。坏数据可能使得PMU在应用过程中导致错误的分析结果和控制策略。状态估计能够剔除量测量中存在的坏数据，因此，研究基于PMU的电力***机电暂态过程中动态状态估计至关重要。

近年来，基于PMU的机电暂态过程发电机动态状态估计问题是电力***动态状态估计领域的热点。针对具有非线性的发电机动态方程，其动态状态估计问题是一个典型的非线性问题，采用以卡尔曼滤波算法为基础方案的是一种普遍的解决思路，如基于扩展卡尔曼滤波(EKF)、无迹卡尔曼滤波(UKF)的动态状态估计。EKF的线性化过程导致截断误差过大，UKF需要确定参数值，灵活性差，应用不便。针对该问题，加拿大学者Simon Haykin于2009年提出容积卡尔曼滤波(CKF)算法。然而，无论UKF还是CKF，当量测量存在坏数据时，估计精度都会在一定程度上受到影响，甚至导致估计结果严重偏离实际值，动态状态估计失败。

发明内容

本发明的目的是提供一种发电机动态状态估计方法，该方法使其能够在PMU量测量含有坏数据的情况下对发电机状态量预报值进行准确修正，得到准确的状态量估计值。

本发明的目的是通过以下技术方案实现的：

一种发电机动态状态估计方法，包括：

根据***方程对k-1时刻的发电机动态状态估计结果进行计算，获得第k时刻发电机动态状态估计结果；

利用CKF滤波算法对第k时刻发电机动态状态估计结果进行预报与滤波处理；其中，通过预报处理获得k+1时刻的发电机动态状态预报结果；在滤波处理时，引入时变多维观测噪声尺度因子并结合PMU的实时量测值对k+1时刻的发电机动态状态预报结果进行滤波，从而实现发电机动态状态结果的准确估计。

进一步的，所述***方程为：

式中，下标k+1与k均相应的表示k+1时刻与k时刻；F和H分别为状态方程函数和量测方程函数，x、u和z分别为状态量、控制量和量测量；v和w分别为***噪声和量测噪声，误差方差阵分别为Q和R的正态分布；

其中，状态量x＝[δωE′_qE′_d]^T，控制量状态方程的具体形式为：

式中，δ为发电机转子绝对功角标幺值，ω为电角速度标幺值；和是对时间微分的简化写法，字母上方的·表示d/dt微分算子；T_J为发电机惯性时间常数；T_m为机械转矩；U_t和分别为发电机出口电压相量的幅值和相角；X′_q和X′_d分别为q轴和d轴瞬变电抗；E′_q和E′_d分别为q轴和d轴瞬态电动势；D为阻尼系数；T′_q0和T′_d0分别为q轴和d轴开路瞬变时间常数；E_f为定子励磁电动势；X_q和X_d分别为q轴和d轴同步电抗；

量测量z＝[δ_z ω_z P_e]^T，量测方程的具体形式为：

式中，δ_z、ω_z和P_e分别为转子绝对功角、电角速度和输出电磁功率的PMU量测值；

根据上述状态方程与量测方程的具体形式，以及k时刻***噪声方差矩阵Q_k和k+1时刻量测噪声方差阵R_k+1，则能够实现发电机动态状态估计；

其中的k+1时刻量测噪声方差阵R_k+1根据PMU实际量测误差进行取值，k时刻***噪声方差矩阵Q_k表示为：

式中，和分别为k时刻δ、ω、E′_q和E′_d的***噪声方差；

通过误差传递公式进行计算：

式中：σ为噪声方差；为发电机出口电压相量的幅值PMU量测误差标准差，为发电机出口电压相量的相角PMU量测误差标准差；

则有：

式中，Δt为步长。

进一步的，所述通过预报处理获得k+1时刻的发电机动态状态预报结果的步骤包括：

所述的第k时刻发电机动态状态估计结果包括：第k时刻发电机动态状态量的估计值与估计误差方差阵P_k|k；

对P_k|k进行Cholesky分解，得到k时刻估计误差方差阵的平方根矩阵S_k|k：

利用下式对第k时刻发电机动态状态量的估计值生成一组等权值的容积点X_i,k|k：

式中，参数n为状态量维数；

利用下式对每一个容积点进行变换，得到所有容积点的预报值

对所有发电机状态量容积点的预报值进行加权求和，得到状态量预报值

并通过下式获得预报误差协方差阵P_k+1|k：

进一步的，所述在滤波处理时，引入时变多维观测噪声尺度因子并结合PMU的实时量测值对k+1时刻的发电机动态状态预报结果进行滤波，从而实现发电机动态状态结果的准确估计包括：

所述的k+1时刻的发电机动态状态预报结果包括：k+1时刻的发电机动态状态量预报值与发电机状态量预报误差协方差阵P_k+1|k；

对预报误差协方差阵P_k+1|k进行Cholesky分解，得到k+1时刻发电机状态量预报误差协方差阵的平方根矩阵S_k+1|k：

利用下式对k+1时刻的发电机动态状态量预报值生成一组等权值的状态量预报值容积点X_i,k+1|k：

利用下式对每一个状态量预报值容积点进行变换，得到PMU量测量预报值的容积点Z_i,k+1|k：

Z_i,k+1|k＝H(X_i,k+1|k,u_k)；

对所有PMU量测量预报值的容积点进行加权求和，进而得到PMU量测量预报值

计算PMU量测量预报误差方差阵P_vv,k+1：

利用得到PMU的实时量测值z_k+1＝[δ_zk+1 ω_zk+1 P_ek+1]^T计算新息向量e_k+1，进而得到时变多维观测噪声尺度因子γ_k+1：

式中，M为开窗估计法的窗长；

再利用下式计算对角阵γ′_k+1：

γ′_k+1＝diag(γ′₁,γ′₂,…,γ′_m)；

式中，对角元素γ′₁的取值为：γ′_i＝max{1,γ_k+1,ii}，i＝1,2,…,n；γ_k+1,ii为γ_k+1的第i个对角元素；

根据下式计算发电机状态量预报值和PMU量测量预报值之间的互协方差矩阵P_xz,k+1|k：

再计算卡尔曼滤波增益W_k+1：

W_k+1＝P_xz,k+1|k(P_vv,k+1+γ′_k+1R_k+1)^-1；

利用k+1时刻新息向量e_k+1，并通过卡尔曼滤波增益W_k+1对k+1时刻的发电机动态状态量预报值进行滤波，获得k+1时刻发电机动态状态量的估计值

并通过下式计算获得发电机动态状态量估计误差方差阵P_k+1|k+1：

由上述本发明提供的技术方案可以看出，该估计方法针对PMU量测量中坏数据导致发电机动态状态估计结果偏离真实值的问题，将时变多维观测噪声尺度因子引入到CKF滤波算法中，不仅提高了该方法的鲁棒性，还可得到准确的状态量估计值。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域的普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图。

图1为本发明实施例提供的一种发电机动态状态估计方法的简要流程图；

图2为本发明实施例提供的IEEE9节点测试***示意图；

图3为本发明实施例提供的IEEE9节点仿真测试中发电机G1动态状态估计结果图；

图4为本发明实施例提供的IEEE9节点仿真测试中发电机G2动态状态估计结果图；

图5为本发明实施例提供的IEEE9节点仿真测试中发电机G3动态状态估计结果图；

图6为本发明实施例提供的IEEE9节点仿真测试中多维时变观测噪声尺度因子变化情况图。

具体实施方式

下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明的保护范围。

本发明实施例提供一种发电机动态状态估计方法，其主要包括如下步骤：

根据***方程对k-1时刻的发电机动态状态估计结果进行计算，获得第k时刻发电机动态状态估计结果；

利用CKF滤波算法对第k时刻发电机动态状态估计结果进行预报与滤波处理；其中，通过预报处理获得k+1时刻的发电机动态状态预报结果；在滤波处理时，引入时变多维观测噪声尺度因子并结合PMU的实时量测值对k+1时刻的发电机动态状态预报结果进行滤波，从而实现发电机动态状态结果的准确估计。

发电机动态状态估计的简要流程图如图1所示。发电机端口电压可直接由PMU直接量测，与外部网络解耦，不必考虑网络拓扑结构，从而进行分布式发电机动态状态估计。发电机转子的内在巨大惯性使得发电机转子功角和电角速度在机电暂态过程中不会发生突变，并且满足发电机运动方程的约束条件。PMU能够对发电机功角和角速度进行直接测量，利用GPS可以得到电磁功率的同步量测值；前一时刻状态量结合发电机端口电压相量，通过***方程获得下一时刻发电机状态量的估计结果；然后，通过预报获得状态量的预报值，再结合PMU对状态量的量测值，通过改进的具有鲁棒性CKF算法，最终得到下一时刻发电机的动态状态估计结果。

为了便于理解，下面对本发明做详细的说明。

发电机动态状态估计问题是非线性滤波问题，非线性动态***的***方程可以表示为：

式中，下标k+1与k均相应的表示k+1时刻与k时刻；F和H分别为状态方程函数和量测方程函数，x、u和z分别为状态量、控制量和量测量；v和w分别为***噪声和量测噪声，误差方差阵分别为Q和R的正态分布；

由于现有控制备难以对次暂态过程做出及时反应，因此，忽略次暂态，动态方程由六阶降为四阶；

状态量x＝[δωE′_qE′_d]^T，控制量状态方程的具体形式为：

式中，δ为发电机转子绝对功角标幺值，ω为电角速度标幺值；和是对时间微分的简化写法，字母上方的·表示d/dt微分算子；T_J为发电机惯性时间常数；T_m为机械转矩；U_t和分别为发电机出口电压相量的幅值和相角；X′_q和X′_d分别为q轴和d轴瞬变电抗；E′_q和E′_d分别为q轴和d轴瞬态电动势；D为阻尼系数；T′_q0和T′_d0分别为q轴和d轴开路瞬变时间常数；E_f为定子励磁电动势；X_q和X_d分别为q轴和d轴同步电抗；

量测量z＝[δ_z ω_z P_e]^T，量测方程的具体形式为：

式中，δ_z、ω_z和P_e分别为转子绝对功角、电角速度和输出电磁功率的PMU量测值；

根据上述状态方程与量测方程的具体形式，以及k时刻***噪声方差矩阵Q_k和k+1时刻量测噪声方差阵R_k+1，则能够实现发电机动态状态估计。

其中的k+1时刻量测噪声方差阵R_k+1根据PMU实际量测误差进行取值。

针对***噪声，在发电机动态状态估计的***方程中，假设模型参数准确；状态量均采用估计值，精度较高；控制量u中T_m、E_f可以假设为恒定值；U_t和采用PMU量测值，存在量测误差。因此，***噪声主要源于U_t和的量测误差。该量测误差通过***方程进行传递，最终导致***噪声的产生。

k时刻***噪声方差矩阵Q_k表示为：

式中，和分别为k时刻δ、ω、E′_q和E′_d的***噪声方差；

通过误差传递公式进行计算：

式中：σ为噪声方差；σ_U为发电机出口电压相量的幅值PMU量测误差标准差，为发电机出口电压相量的相角PMU量测误差标准差；可设置其取值分别为0.2％和0.2°。

按照式(5)，对***方程分别关于发电机出口电压幅值和相角求偏导数，进而推导k时刻状态量ω、E′_q和E′_d过程噪声方差的计算公式。由于发电机功角δ的预报值需要用到电角速度ω的估计值，而ω发估计值误差相对较小，因此，将功角的***噪声误差方差定为一个较小的正数，此处选为10^-6。

则有：

由式(7)-(9)就可以计算Q_k中各元素值。式中，Δt为步长，式(7)-(9)中U_t和取值为k时刻PMU的量测值；E′_q、E′_d和δ的取值为k时刻的估计值。

再利用CKF滤波算法对第k时刻发电机动态状态估计结果进行预报与滤波处理。

1、预报过程。

对第k时刻发电机动态状态估计结果进行预报处理获得k+1时刻的发电机动态状态预报结果；

所述的第k时刻发电机动态状态估计结果包括：第k时刻发电机动态状态量的估计值与估计误差方差阵P_k|k；

对P_k|k进行Cholesky分解，得到k时刻估计误差方差阵的平方根矩阵S_k|k：

利用下式对第k时刻发电机动态状态量的估计值生成一组等权值的容积点X_i,k|k：

式中，参数n为状态量维数；示例性的，若发电机动态状态估计状态量维数n＝4，则容积点个数为8。

利用下式对每一个容积点进行变换，得到所有容积点的预报值

对所有发电机状态量容积点的预报值进行加权求和，得到状态量预报值

示例性的，若发电机动态状态估计状态量维数n＝4，则利用球面-径向规则生成的容积点权重为1/8。

并通过下式获得预报误差协方差阵P_k+1|k：

2、滤波过程。

在滤波处理时，引入时变多维观测噪声尺度因子并结合PMU的实时量测值对k+1时刻的发电机动态状态预报结果进行滤波，从而实现发电机动态状态结果的准确估计。

所述的k+1时刻的发电机动态状态预报结果包括：k+1时刻的发电机动态状态量预报值发电机状态量预报误差协方差阵P_k+1|k；

对预报误差协方差阵P_k+1|k进行Cholesky分解，得到k+1时刻发电机状态量预报误差协方差阵的平方根矩阵S_k+1|k：

利用下式对k+1时刻的发电机动态状态量预报值生成一组等权值的状态量预报值容积点X_i,k+1|k：

利用下式对每一个状态量预报值容积点进行变换，得到PMU量测量预报值的容积点Z_i,k+1|k：

Z_i,k+1|k＝H(X_i,k+1|k,u_k)； (17)

对所有PMU量测量预报值的容积点进行加权求和，进而得到PMU量测量预报值

计算PMU量测量预报误差方差阵P_vv,k+1：

利用得到PMU的实时量测值z_k+1＝[δ_zk+1 ω_zk+1 P_ek+1]^T计算新息向量e_k+1，进而得到时变多维观测噪声尺度因子γ_k+1：

式中，M为开窗估计法的窗长；

再利用下式计算对角阵γ′_k+1：

γ′_k+1＝diag(γ′₁,γ′₂,…,γ′_m)； (22)

式中，对角元素γ′₁的取值为：γ′_i＝max{1,γ_k+1,ii}，i＝1,2,…,n；γ_k+1,ii为γ_k+1的第i个对角元素；

根据下式计算发电机状态量预报值和PMU量测量预报值之间的互协方差矩阵P_xz,k+1|k：

再计算卡尔曼滤波增益W_k+1：

W_k+1＝P_xz,k+1|k(P_vv,k+1+γ′_k+1R_k+1)^-1； (24)

利用k+1时刻新息向量e_k+1，并通过卡尔曼滤波增益W_k+1对k+1时刻的发电机动态状态量预报值进行滤波，获得k+1时刻发电机动态状态量的估计值

并通过下式计算获得发电机动态状态量估计误差方差阵P_k+1|k+1：

上述滤波过程中，将时变多维观测噪声尺度因子引入到滤波处理中，根据量测新息对量测误差进行在线调整，使其更加逼近真实噪声；再利用调整后的误差计算滤波增益，使其能够在量测量含有坏数据的情况下对状态量预报值进行准确修正，得到精确的状态量估计值。

本发明实施例对CKF算法进行改进的原理如下：

容积卡尔曼滤波是根据确定的容积点来计算后验概率密度函数，无需计算复杂的雅克比矩阵，比扩展卡尔曼滤波实现更加容易，同时也避免了截断误差。与无迹变换卡尔曼滤波相比，容积卡尔曼滤波的容积点在更低一维的子空间中对称出现，且各容积点权值大小相等，无需像无迹变换卡尔曼滤波那样提前设置参数，在选择方式上更加简单，适应性更强。

作为一种基于卡尔曼滤波构架的非线性状态估计方法，容积卡尔曼滤波需要以准确的***模型和噪声统计特性为基础。然而，实际应用中，一方面获得准确的噪声先验统计特性往往比较困难；另一方面，即使获得了精确的先验噪声统计特性，***在运行过程中受到内部或外部环境不确定因素的影响，可能会引起噪声统计特性的变化。例如，在发电机动态状态估计过程中，PMU量测有可能出现不良数据。一旦量测出现不良数据，就会使得量测误差方差阵R与实际误差不符，通过传统方法来计算会使得量测误差总方差无法正确反映预报值的偏差，传统方法中量测误差总方差P_zz,k+1|k＝P_vv,k+1+R_k+1；从而使得传统方法中卡尔曼滤波增益无法对状态量预报值进行正确修正，使得估计结果不准确，传统方法中卡尔曼滤波增益传统的容积卡尔曼滤波算法对初始误差比较敏感，对于量测噪声统计特性缺乏在线调整的自适应能力。噪声的突然变化、观测值可信度降低等因素都可能影响算法估计精度，甚至导致滤波值发散。

因此，需要通过改变滤波增益以适应噪声统计特性的变化。噪声先验统计特性比较准确时，新息可信程度高，对应的权重比较大，通过一个较大的滤波增益能有效抑制噪声引起的误差；而当PMU出现不良数据时，噪声统计特性偏差较大，新息可信程度降低，需要通过降低滤波增益改变新息权重，对状态更新做出调节，提高滤波精度。因此，可以通过改变滤波增益以适应噪声统计特性的变化。噪声先验统计特性比较准确时，新息可信程度高，对应的权重比较大，通过一个较大的滤波增益能有效抑制噪声引起的误差；而当PMU出现不良数据时，噪声统计特性偏差较大，新息可信程度降低，需要通过降低滤波增益改变新息权重，对状态更新做出调节，提高滤波精度。

本发明实施例中，针对量测出现坏数据导致量测误差方差与实际值不符的问题，将一个时变多维观测噪声尺度因子引入CKF。由式(25)可见，CKF滤波算法需要根据新息向量对状态量的预报值进行修正，进而得到状态量估计值修正程度由新息向量e_k+1和卡尔曼滤波增益W_k+1共同决定。具体来说：

采用传统方法时，当实际量测噪声与给定的量测误差方差阵R相符时，e_k+1和W_k+1能够对预报值进行正确修正，CKF能够得到准确的估计结果。然而，当量测量出现坏数据时，新息向量e_k+1中坏数据对应的元素突然增大，而W_k+1并未随之进行调整，对状态量预报值的修正不准确，导致估计结果准确度下降。

由前文所述传统方法中的公式：可见，W_k+1是量测噪声方差阵R的函数，为了使W_k+1能够对状态量预报值进行正确修正，本发明实施例中构造时变多维观测噪声尺度因子γ_k+1，对R进行在线调整。当量测量出现坏数据时，γ_k+1的值随之变化，进而调整W_k+1使其能够对状态量预报值进行准确修正。

即，将传统的P_zz,k+1|k＝P_vv,k+1+R_k+1表达式修改为：

P_zz,k+1|k＝P_vv,k+1+γ_k+1R_k+1 (26)

式中，γ_k+1为m维时变多维观测噪声尺度因子，m为量测量维数。由此，传统的表达式变为：

式中，等号右边求逆部分为观测噪声方差，为使W_k+1能够准确修正状态量预报值，需要满足

当量测量没有坏数据时，γ_k+1为m维单位矩阵。当量测量出现坏数据时，新息e_k+1会增大，γ_k+1不再为单位矩阵，具体取值情况需要另行计算。

利用开窗估计法计算新息实时协方差：

式中，P_e,k+1为新息协方差；M为开窗估计法的窗长。当量测出现坏数据时，必然有P_e,k+1＞P_vv,k+1+R_k+1，则需要计算γ_k+1，进而对量测误差方差进行调整。γ_k+1的取值原则是使得新息方差与观测噪声相等，即

由下式(31)求解得到γ_k+1：

使用开窗估计法计算得到的新息协方差可能会导致式(31)计算的时变多维观测噪声尺度因子不是对角阵，从而使得式(27)等号右边矩阵求逆发生奇异。针对该问题，定义对角阵

γ′_k+1＝diag(γ′₁,γ′₂,…,γ′_m) (32)

式(32)中对角元素γ′_i的取值为

γ′_i＝max{1,γ_k+1,ii}，i＝1,2,…,n (33)

式中，γ_k+1,ii为γ_k+1的第i个对角元素。由于γ′_k+1为对角阵且对角元素均大于0，因此，将γ′_k+1作为时变多维观测噪声尺度因子计算式(27)的滤波增益就可以保证求逆不发生奇异。则滤波增益计算公式为

W_k+1＝P_xz,k+1|k(P_vv,k+1+γ′_k+1R_k+1)^-1 (34)

由式(34)可见，滤波增益中引入时变多维观测噪声尺度因子，可以根据量测新息对量测误差方差矩阵进行调整。当量测量出现坏数据后，在观测噪声尺度因子的作用下，滤波增益减小，进而降低坏数据对估计结果的影响。

为了进一步说明本发明，下面再以具体的实例对上述方案进行仿真测试，具体来说：

IEEE9节点仿真测试：

将本发明提供的基于PMU的电力***分布式动态状态估计方法应用于如图2所示的IEEE9节点测试***(***的基准容量为100MVA)，IEEE9节点测试***的相关参数及发电机动态参数和支路参数如表1-表4所示。

表1 IEEE 9节点测试***节点数据表

注：平衡节点为节点1，除平衡节点外，凡未明确给出电压幅值的节点均为PQ节点，明确给出电压幅值的点均为PV节点。

表2 IEEE 9节点测试***线路支路数据表

表3 IEEE 9节点测试***变压器支路数据表

表4 IEEE 9节点测试***发电机动态数据表

本实例中，假设t＝0时刻在节点5和6之间的线路发生三相金属性短路故障，随后断路器动作将故障线路从***中断开。利用BPA对故障过程进行仿真计算，仿真步长为1周波，即0.02s，仿真时间为18s。将仿真结果作为真值，在真值的基础上叠加高斯白噪声作为量测值。发电机转子功角和电角速度的PMU量测误差标准差分别为2°和0.1％，发电机出口电压相量的幅值和相角的PMU量测误差标准差分别为0.1％和0.1°。

图3-图5分别给出了IEEE9节点测试***中发电机G1、G2和G3的动态状态估计结果，其中的“——鲁棒CKF”为本发明所述方案的结果，“……CKF”为传统CKF算法的结果，“—真值”为实际值。

对于发电机G1，分别在t＝4s和t＝8s时发电机绝对功角的PMU量测值中人为叠加10个连续不良数据点。从图3可以看出，由于不良数据的存在，基于容积卡尔曼滤波的发电机动态状态估计出现了较大的波动。这是因为量测误差方差与绝对功角的PMU量测中不良数据的实际误差相差较大，使得容积卡尔曼滤波的滤波增益无法对状态量预报值进行正确修正，最终导致状态量的估计值不准确。而本发明的方案采用时变多维观测噪声尺度因子对量测误差方差进行实时调整，使其能够随量测噪声而改变。由此计算得到的滤波增益能够正确修正发电机状态量的预报值。因此，当发电机功角的PMU量测值出现不良数据时，依然能够得到准确的状态量估计值。

对于发电机G2，分别在t＝6s和t＝12s时向发电机转子电角速度的PMU量测值中人为叠加10个连续坏数据，在t＝14s时叠加单点坏数据。从图4可以看出，由于转子电角速度的量测中存在不良数据，基于容积卡尔曼滤波算法的发电机动态状态估计发生了较大的波动，而本发明的方案依然能够得到较为准确的动态状态估计值。

对于发电机G3，在t＝12s时同时在发电机转子电角速度和转子绝对功角的PMU量测值中人为叠加连续不良数据。从图5可见，即使绝对功角和电角速度的PMU量测同时存在连续多点不良数据，本发明的方案依然能够得到准确的状态估计值。

图6给出了PMU量测量存在不良数据时，量测量的时变多维观测噪声尺度因子的变化情况。可以看出，当PMU量测量出现不良数据时，量测量的时变多维观测噪声尺度因子会突然增大。这就能够有效调节量测误差方差，由此计算出的卡尔曼滤波增益就能够准确修正状态量的预报值，从而提高状态估计精度。

通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到上述实施例可以通过软件实现，也可以借助软件加必要的通用硬件平台的方式来实现。基于这样的理解，上述实施例的技术方案可以以软件产品的形式体现出来，该软件产品可以存储在一个非易失性存储介质(可以是CD-ROM，U盘，移动硬盘等)中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述的方法。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (3)

1.一种发电机动态状态估计方法，其特征在于，包括：

根据***方程对k-1时刻的发电机动态状态估计结果进行计算，获得第k时刻发电机动态状态估计结果；

利用CKF滤波算法对第k时刻发电机动态状态估计结果进行预报与滤波处理；其中，通过预报处理获得k+1时刻的发电机动态状态预报结果；在滤波处理时，利用得到PMU的实时量测值z_k+1＝[δ_zk+1 ω_zk+1 P_ek+1]^T计算新息向量e_k+1，进而得到时变多维观测噪声尺度因子；下标k+1表示k+1时刻，δ_z、ω_z和P_e分别为转子绝对功角、电角速度和输出电磁功率的PMU量测值；引入时变多维观测噪声尺度因子并结合PMU的实时量测值对k+1时刻的发电机动态状态预报结果进行滤波，从而实现发电机动态状态结果的准确估计；

其中，所述***方程为：

式中，下标k+1与k均相应的表示k+1时刻与k时刻；F和H分别为状态方程函数和量测方程函数，x、u和z分别为状态量、控制量和量测量；v和w分别为***噪声和量测噪声，误差方差阵分别为Q和R的正态分布；

其中，状态量x＝[δωE′_qE′_d]^T，控制量状态方程的具体形式为：

式中，δ为发电机转子绝对功角标幺值，ω为电角速度标幺值；和是对时间微分的简化写法，字母上方的·表示d/dt微分算子；T_J为发电机惯性时间常数；T_m为机械转矩；U_t和分别为发电机出口电压相量的幅值和相角；X′_q和X′_d分别为q轴和d轴瞬变电抗；E′_q和E′_d分别为q轴和d轴瞬态电动势；D为阻尼系数；T′_q0和T′_d0分别为q轴和d轴开路瞬变时间常数；E_f为定子励磁电动势；X_q和X_d分别为q轴和d轴同步电抗；

量测量z＝[δ_z ω_z P_e]^T，量测方程的具体形式为：

式中，δ_z、ω_z和P_e分别为转子绝对功角、电角速度和输出电磁功率的PMU量测值；

根据上述状态方程与量测方程的具体形式，以及k时刻***噪声方差矩阵Q_k和k+1时刻量测噪声方差阵R_k+1，则能够实现发电机动态状态估计；

其中的k+1时刻量测噪声方差阵R_k+1根据PMU实际量测误差进行取值，k时刻***噪声方差矩阵Q_k表示为：

式中，和分别为k时刻δ、ω、E′_q和E′_d的***噪声方差；

通过误差传递公式进行计算：

式中：σ为噪声方差；σ_U为发电机出口电压相量的幅值PMU量测误差标准差，为发电机出口电压相量的相角PMU量测误差标准差；

则有：

式中，Δt为步长。

2.根据权利要求1所述的方法，其特征在于，所述通过预报处理获得k+1时刻的发电机动态状态预报结果的步骤包括：

所述的第k时刻发电机动态状态估计结果包括：第k时刻发电机动态状态量的估计值与估计误差方差阵P_k|k；

对P_k|k进行Cholesky分解，得到k时刻估计误差方差阵的平方根矩阵S_k|k：

利用下式对第k时刻发电机动态状态量的估计值生成一组等权值的容积点X_i,k|k：

式中，参数n为状态量维数；

利用下式对每一个容积点进行变换，得到所有容积点的预报值

对所有发电机状态量容积点的预报值进行加权求和，得到状态量预报值

并通过下式获得预报误差协方差阵P_k+1|k：

3.根据权利要求1所述的方法，其特征在于，所述在滤波处理时，引入时变多维观测噪声尺度因子并结合PMU的实时量测值对k+1时刻的发电机动态状态预报结果进行滤波，从而实现发电机动态状态结果的准确估计包括：

所述的k+1时刻的发电机动态状态预报结果包括：k+1时刻的发电机动态状态量预报值与发电机状态量预报误差协方差阵P_k+1|k；

对预报误差协方差阵P_k+1|k进行Cholesky分解，得到k+1时刻发电机状态量预报误差协方差阵的平方根矩阵S_k+1|k：

利用下式对k+1时刻的发电机动态状态量预报值生成一组等权值的状态量预报值容积点X_i,k+1|k：

式中，参数n为状态量维数；

利用下式对每一个状态量预报值容积点进行变换，得到PMU量测量预报值的容积点Z_i,k+1|k：

Z_i,k+1|k＝H(X_i,k+1|k,u_k)；

对所有PMU量测量预报值的容积点进行加权求和，进而得到PMU量测量预报值

计算PMU量测量预报误差方差阵P_vv,k+1：

利用得到PMU的实时量测值z_k+1＝[δ_zk+1 ω_zk+1 P_ek+1]^T计算新息向量e_k+1，进而得到时变多维观测噪声尺度因子γ_k+1：

式中，M为开窗估计法的窗长；

再利用下式计算对角阵γ′_k+1：

γ′_k+1＝diag(γ′₁,γ′₂,…,γ′_m)；

式中，对角元素γ′₁的取值为：γ′_i＝max{1,γ_k+1,ii}，i＝1,2,…,n；γ_k+1,ii为γ_k+1的第i个对角元素；

根据下式计算发电机状态量预报值和PMU量测量预报值之间的互协方差矩阵P_xz,k+1|k：

再计算卡尔曼滤波增益W_k+1：

W_k+1＝P_xz,k+1|k(P_vv,k+1+γ′_k+1R_k+1)^-1；

利用k+1时刻新息向量e_k+1，并通过卡尔曼滤波增益W_k+1对k+1时刻的发电机动态状态量预报值进行滤波，获得k+1时刻发电机动态状态量的估计值

并通过下式计算获得发电机动态状态量估计误差方差阵P_k+1|k+1：