CN110112770A - 一种基于自适应h∞容积卡尔曼滤波的发电机动态状态估计方法 - Google Patents

Abstract

本发明公开了一种基于自适应H_∞容积卡尔曼滤波的动态状态估计方法，用于增强模型参数不确定情形下动态状态估计器的鲁棒性，实现电力系统发电机状态的准确估计。该方法在H_∞容积卡尔曼滤波(H_∞cubature Kalman filter,HCKF)基础上，通过引入自适应方法，建立了自适应H_∞容积卡尔曼滤波动态估计方法，该方法不仅可以界定系统模型参数不确定性对状态估计结果的影响，而且避免了HCKF方法不确定上界参数难选取的问题。本发明方法能够有效提升状态估计精度，增强状态估计器的鲁棒性，且实现流程清晰，具有较高的工程应用价值。

Description

一种基于自适应H∞容积卡尔曼滤波的发电机动态状态估计 方法

技术领域

本发明涉及一种基于自适应H_∞容积卡尔曼滤波的发电机动态状态估计方法，属于电力监测与控制领域。

背景技术

近年来，广域量测系统(WAMS)的同步相量量测单元(PMU)应用日趋广泛，PMU可提供带有时标的高频率系统信息采样值，能够应用于电力系统机电暂态的过程分析与控制。然而，WAMS作为一个量测系统，在量测过程中会不可避免受到随机干扰等因素的影响，造成量测数据的污染。因此，由PMU获取的量测信息无法直接用于电力系统机电暂态分析。动态状态估计不仅可以有效的滤除量测数据中的误差和噪声值，而且，借助其预测功能可以为系统未来可能的变化制定相应的控制策略。所以，提高电力系统中发电机动态状态估计跟踪精度对于电网安全稳定运行有重要意义。

目前，针对电力系统发电机进行动态状态估计，基于确定的状态估计模型，研究人员开发卡尔曼滤波，扩展卡尔曼滤波(EKF)、粒子滤波、无迹变换卡尔曼滤波等状态估计方法。但是，在实际的电力系统分析时，系统的模型往往存在不确定性，诸如系统噪声和量测噪声所满足的统计规律无法准确获知，而且模型的一些参数辨识结果也存在一定的程度的不确定性，这些不确定性因素将会严重影响上述状态估计器的性能，甚至导致状态估计结果无法收敛。

发明内容

本发明针对上述电力系统动态状态估计中存在的困难和亟需解决的问题，设计了一种基于自适应H_∞容积卡尔曼滤波的发电机动态状态估计器。

本发明公开了一种基于自适应H∞容积卡尔曼滤波的发电机动态状态估计方法，包括以下步骤：

步骤一：获取自适应H_∞容积卡尔曼滤波的滤波参数初始值；

步骤二：建立发电机动态状态估计模型；

步骤三：基于自适应H_∞容积卡尔曼滤波的预测步，计算k时刻的状态预测值

步骤四：计算k时刻预测误差协方差矩阵

步骤五：生成围绕状态预测值的等权值Cubature点集X_i,k，通过状态估计模型对状态预测值Cubature点进行变换，获取对应量测预测值的Cubature点Z_i,k；

步骤六：基于量测预测值的Cubature点Z_i,k，求解k时刻量测预测值

步骤七：计算k时刻量测预测误差协方差矩阵P_zz,k和状态预测值与量测预测值之间的交互协方差矩阵P_xz,k；

步骤八：计算k时刻的滤波增益K_k，并采用量测值z_k对状态预测值进行修正更新；

步骤九：根据外界情况变化，自适应计算并更新k时刻估计误差协方差矩阵

步骤十：按照步骤三至步骤九依据量测信息时间序列进行发电机动态状态估计，直至k+1＞N时状态估计程序结束，输出发电机动态估计结果。

可选的，所述状态估计模型，表示为：

式中，f(·)表示发电机状态方程，h(·)表示量测方程，x，u，z分别对应表示状态变量，控制变量及量测向量；下标k和k+1表示时刻，w表示系统噪声，v为量测噪声，二者均为高斯白噪声，所满足的均值均为0，对应的协方差矩阵分别为Q与R，w与v相互独立且和状态变量无关。

可选的，所述步骤三中，k时刻的状态预测值求解如下：

式中，X_i,k-1和分别表示k-1时刻发电机状态量的第i个Cubature点及其k时刻预测值；为k-1时刻估计误差协方差，表示k-1时刻的发电机状态估计值；n表示发电机状态估计变量个数，表示对矩阵进行Cholesky分解运算，u_k-1表示k-1时刻控制向量；ξ_i为容积点集合{ξ_i},i＝1,2,…2n的第i个容积点，容积点集合形式为：

可选的，步骤四中，k时刻预测误差协方差矩阵求解计算如下：

式中，Q_k-1分别表示k-1时刻系统噪声满足的协方差矩阵。

可选的，步骤五中，等权值Cubature点集X_i,k和Cubature点Z_i,k，计算公式如下：

可选的，步骤六中，k时刻量测预测值的计算公式如下：

可选的，步骤七中，k时刻量测预测误差协方差矩阵P_zz,k和状态预测值与量测预测值之间的交互协方差矩阵P_xz,k，计算公式如下：

式中，R_k表示k时刻量测误差协方差矩阵。

可选的，步骤八中，k时刻的滤波增益K_k表达式为：

采用量测值z_k对状态预测值进行修正更新，计算公式为：

式中，上标(·)^-1表示矩阵的求逆运算，为k时刻状态估计值，z_k表示量测值。

可选的，步骤九中，更新k时刻估计误差协方差矩阵表达式如下：

式中，α是一个待设定的正常数，用于调节动态过程中估计误差协方差自适应变换的阈值，其中R_e,k的计算方法如下：

式中，表示k时刻的信息矩阵，为k时刻信息序列，ρ为遗忘因子，γ为界定发电机不确定性引起状态估计误差的正标量参数，L_k矩阵计算公式为：

式中表示对矩阵进行Cholesky分解运算，I为适当维度的单位矩阵，为估计误差协方差最大上界设置参数。

有益效果：本发明的一种基于自适应H_∞容积卡尔曼滤波的发电机动态状态估计方法，有效结合H_∞滤波理论与容积卡尔曼滤波，通过创新性的引入自适应技术动态更新调整估计误差协方差矩阵，不仅可以有效界定模型不确定性引起的估计误差，增强估计器的鲁棒性，提升状态估计精度，且避免了传统H_∞容积卡尔曼滤波估计误差上界难选取的问题。

附图说明

图1：为本发明的方法流程图；

图2：为IEEE 10机39节点系统结构图；

图3：为利用CKF、HCKF和本发明方法的发电机功角和角速度估计结果对比；

图4：为利用CKF、HCKF和本发明方法的发电机暂态电动势估计结果对比；

图5：为利用CKF、HCKF和本发明方法的发电机状态估计结果误差对比。

具体实施方式

下面结合附图和实施例进一步阐述本发明。

如图1所示，一种基于自适应H_∞容积卡尔曼滤波的发电机动态状态估计方法，其包含如下步骤：

(1)状态估计模型建立

发电机动态状态估计的状态方程与量测方程，其一般形式可以表示为：

式中，f(·)表示发电机状态方程，h(·)表示量测方程，x，u，z分别对应表示状态变量，控制变量及量测向量；下标k和k+1表示时刻，w表示系统噪声，v为量测噪声，二者均为高斯白噪声，所满足的均值均为0，对应的协方差矩阵分别为Q与R，w与v相互独立且和状态变量无关。

(2)设定自适应H_∞容积卡尔曼滤波的滤波参数初始值，包含设定初始时刻的状态变量值状态估计误差协方差矩阵控制变量值u₀取为稳态运行值；设定系统噪声和量测噪声所满足的初始协方差矩阵Q₀，R₀，以及最大估计时刻N；

(3)基于自适应H_∞容积卡尔曼滤波的预测步，计算k时刻的状态预测值求解计算公式如下：

式中X_i,k-1和分别表示k-1时刻发电机状态量的第i个Cubature点及其k时刻预测值；为k-1时刻估计误差协方差，表示k-1时刻的发电机状态估计值；n表示发电机状态估计变量个数，表示对矩阵进行Cholesky分解运算，u_k-1表示k-1时刻控制向量；ξ_i为容积点集合{ξ_i},i＝1,2,…2n的第i个容积点，容积点集合形式为

(4)计算k时刻预测误差协方差矩阵求解计算公式如下

式中上标T表示矩阵的转置运算，Q_k-1分别表示k-1时刻系统噪声满足的协方差矩阵。

(5)生成围绕状态预测值的等权值Cubature点集X_i,k；之后，通过状态估计模型的量测方程对状态预测值Cubature点进行变换，获取对应量测预测值的Cubature点Z_i,k，计算公式如下

(6)基于量测预测值的Cubature点Z_i,k，求解k时刻量测预测值计算公式如下

(7)计算k时刻量测预测误差协方差矩阵P_zz,k和状态预测值与量测预测值之间的交互协方差矩阵P_xz,k，计算公式如下

式中R_k表示k时刻量测误差协方差矩阵，上标T表示矩阵的转置运算。

(8)计算k时刻的滤波增益K_k，并采用量测值z_k对状态预测值进行修正更新，计算公式分别为

式中上标(·)^-1表示矩阵的求逆运算，为k时刻状态估计值，z_k表示量测值。

(9)根据外界情况变化，自适应计算并更新k时刻估计误差协方差矩阵计算公式如下：

式中上标(·)^-1表示对矩阵求逆，α是一个待设定的正常数，用于调节动态过程中估计误差协方差自适应变换的阈值，其中R_e,k的计算方法如下

式中表示k时刻的信息矩阵，为k时刻信息序列，ρ＝0.98为遗忘因子，γ为界定发电机不确定性引起状态估计误差的正标量参数，L_k矩阵计算公式为：

式中表示对矩阵进行Cholesky分解运算，I为适当维度的单位矩阵，为估计误差协方差最大上界设置参数。

(10)按照(3)-(9)步骤依据量测信息时间序列进行发电机动态状态估计，直至k+1＞N时状态估计程序结束，输出发电机动态估计结果。

实施例

(a)模型建立

依据发电机四阶动态方程，构建的发电机状态估计方程如下：

式中：δ表示发电机功角，rad；ω和ω₀分别为电角速度与同步转速，pu；e'_q和e'_d分别表示发电机q轴和d轴的暂态电动势；H表示发电机惯性常数，T_m和T_e分别表示发电机机械功率与电磁功率，其中T_e＝P_e/ω；K_D表示阻尼因子，E_fd为定子励磁电压；T’_d0和T’_q0表示发电机机在d-q坐标系下的开路时间常数；x_d和x'_d分别表示发电机d轴同步电抗与暂态电抗，x_q和x'_q分别为发电机q轴同步电抗与暂态电抗；i_d和i_q分别表示发电机d轴和q轴的定子电流。

对电力系统发电机动态变量进行动态估计时，状态估计向量为x＝(δ,ω,e'_q,e'_d)^T；选取发电机机械功率、定子励磁电压以及定子R轴和I轴的电流i_R,i_I为控制向量，即u＝(T_m,E_fd,i_R,i_I)^T；选取发电机绝对功角、发电机角速度及发电机定子R轴和I轴的电压e_R,e_I作为量测值，即量测向量为

z＝(δ,ω,e_R,e_I)^T

其中发电机的绝对功角与角速度均可由PMU量测设备直接量测获取，此情形下系统满足能观性。

(b)实施例分析

为了验证本发明所提出基于鲁棒容积卡尔曼滤波的发电机动态状态估计方法的有效性和实用性，本实施例选取IEEE 10机39节点系统作为测试系统，系统结构见图2。在算法进行验证时，以系统中发电机G2的状态变量作为估计对象，并将调速器的作用考虑在内，其中发电机采用四阶模型。发电机惯性时间参数为30.3，阻尼因子为2，并假定发电机在节点16-节点21支路发生三相金属性短路故障，故障持续6周期(采样周期为0.02s)后消失，仿真时间为10s。运用BPA软件模拟PMU数据采集，获取发电机运行真实值。量测数据值由真实值叠加随机噪声形成。

在发电机动态运行过程中，由于运行工况的变化，系统噪声与量测噪声所满足的统计特性会发生变化，即噪声的统计特性存在不确定性，为验证本发明方法针对此种情形的有效性，设定系统和量测噪声的协方差矩阵存在不确定性，偏离真实值，分别取为Q＝10^{- 5}I_4×4，R＝10^-5I_4×4(真实值分别为Q＝10^-6I_4×4，R＝10^-6I_4×4)，状态变量初值取为上一时刻的静态状态估计值，初始估计误差协方差矩阵模型不确定性界定参数λ＝1.8，状态估计最大时刻N＝500。

为了对不同算法之间的估计结果进行对比分析，本发明采用平均绝对估计误差MAE作为指标进行算法间性能对比。

式中MAE(k)为k时刻发电机各状态变量估计结果的绝对误差之和的平均值，x_i,k为k时刻第i个状态量的真实值(BPA数据)，为其对应估计值，N_s为总的状态变量数目。

对上述实施例系统，分别运用CKF、HCKF算法(其所需的相关参数值和本发明方法的参数初值相同)，以及本发明提出的AHCKF方法对发电机状态进行估计测试。

不同方法对发电机G2功角和角速度估计结果对比如图3所示，图4给出了发电机暂态电动势状态变量的估计结果对比，图5进一步展示了不同方法的状态变量估计平均绝对误差。从状态估计结果的对比图可以看出，本发明的状态估计精度显著高于CKF和HCKF算法，表明方法能够更好的抑制噪声矩阵不确定性对状态估计结果的影响，对模型不确定性具有更强的鲁棒性。

Claims (9)

1.一种基于自适应H∞容积卡尔曼滤波的发电机动态状态估计方法，其特征在于：包括以下步骤：

步骤一：获取自适应H_∞容积卡尔曼滤波的滤波参数初始值；

步骤二：建立发电机动态状态估计模型；

步骤三：基于自适应H_∞容积卡尔曼滤波的预测步，计算k时刻的状态预测值

步骤四：计算k时刻预测误差协方差矩阵

步骤五：生成围绕状态预测值的等权值Cubature点集X_i,k，通过状态估计模型对状态预测值Cubature点进行变换，获取对应量测预测值的Cubature点Z_i,k；

步骤六：基于量测预测值的Cubature点Z_i,k，求解k时刻量测预测值

步骤七：计算k时刻量测预测误差协方差矩阵P_zz,k和状态预测值与量测预测值之间的交互协方差矩阵P_xz,k；

步骤八：计算k时刻的滤波增益K_k，并采用量测值z_k对状态预测值进行修正更新；

步骤九：根据外界情况变化，自适应计算并更新k时刻估计误差协方差矩阵

步骤十：按照步骤三至步骤九依据量测信息时间序列进行发电机动态状态估计，直至k+1＞N时状态估计程序结束，输出发电机动态估计结果。

2.根据权利要求1所述的一种基于自适应H∞容积卡尔曼滤波的发电机动态状态估计方法，其特征在于：所述状态估计模型，表示为：

式中，f(·)表示发电机状态方程，h(·)表示量测方程，x，u，z分别对应表示状态变量，控制变量及量测向量；下标k和k+1表示时刻，w表示系统噪声，v为量测噪声，二者均为高斯白噪声，所满足的均值均为0，对应的协方差矩阵分别为Q与R，w与v相互独立且和状态变量无关。

3.根据权利要求2所述的一种基于自适应H∞容积卡尔曼滤波的发电机动态状态估计方法，其特征在于：所述步骤三中，k时刻的状态预测值求解如下：

式中，X_i,k-1和分别表示k-1时刻发电机状态量的第i个Cubature点及其k时刻预测值；为k-1时刻估计误差协方差，表示k-1时刻的发电机状态估计值；n表示发电机状态估计变量个数，表示对矩阵进行Cholesky分解运算，u_k-1表示k-1时刻控制向量；ξ_i为容积点集合{ξ_i},i＝1,2,…2n的第i个容积点，容积点集合形式为：

4.根据权利要求3所述的一种基于自适应H∞容积卡尔曼滤波的发电机动态状态估计方法，其特征在于：步骤四中，k时刻预测误差协方差矩阵求解计算如下：

式中，Q_k-1分别表示k-1时刻系统噪声满足的协方差矩阵。

5.根据权利要求4所述的一种基于自适应H∞容积卡尔曼滤波的发电机动态状态估计方法，其特征在于：步骤五中，等权值Cubature点集X_i,k和Cubature点Z_i,k，计算公式如下：

Z_i,k＝h(X_i,k,u_k),i＝1,…,2n (7)。

6.根据权利要求5所述的一种基于自适应H∞容积卡尔曼滤波的发电机动态状态估计方法，其特征在于：步骤六中，k时刻量测预测值的计算公式如下：

7.根据权利要求6所述的一种基于自适应H∞容积卡尔曼滤波的发电机动态状态估计方法，其特征在于：步骤七中，k时刻量测预测误差协方差矩阵P_zz,k和状态预测值与量测预测值之间的交互协方差矩阵P_xz,k，计算公式如下：

式中，R_k表示k时刻量测误差协方差矩阵。

8.根据权利要求7所述的一种基于自适应H∞容积卡尔曼滤波的发电机动态状态估计方法，其特征在于：步骤八中，k时刻的滤波增益K_k表达式为：

采用量测值z_k对状态预测值进行修正更新，计算公式为：

式中，上标(·)^-1表示矩阵的求逆运算，为k时刻状态估计值，z_k表示量测值。

9.根据权利要求8所述的一种基于自适应H∞容积卡尔曼滤波的发电机动态状态估计方法，其特征在于：步骤九中，更新k时刻估计误差协方差矩阵表达式如下：

式中，α是一个待设定的正常数，用于调节动态过程中估计误差协方差自适应变换的阈值，其中R_e,k的计算方法如下：

式中，表示k时刻的信息矩阵，为k时刻信息序列，ρ为遗忘因子，γ为界定发电机不确定性引起状态估计误差的正标量参数，L_k矩阵计算公式为：

式中表示对矩阵进行Cholesky分解运算，I为适当维度的单位矩阵，为估计误差协方差最大上界设置参数。