一种基于自适应容积卡尔曼滤波的动态状态估计方法与流程

本发明属于电力系统分析和监测技术领域，特别涉及一种基于自适应容积卡尔曼滤波的动态状态估计方法。

背景技术：

准确的状态估计对于电力系统的分析与稳定控制具有重要意义。状态估计一般分为两类，一类是静态估计，另一类是动态状态估计。静态状态估计利用某一时刻断面冗余的量测信息实现系统该时刻状态变量估计。虽然静态状态估计精度较高，但是其忽略了电力系统的动态特性。因此，静态状态估计无法应用于电力系统状态的实时在线估计。为了适应电力系统在线监测的需求，具备估计和预测功能的动态状态估计方法近年来得到研究人员的广泛关注。

目前，电力系统动态状态估计方法主要包括卡尔曼滤波、扩展卡尔曼滤波，无迹卡尔曼滤波、粒子滤波及集合卡尔曼滤波等。上述这些方法在一定程度上改善了状态估计的结果。但是，需要指出的是这些方法均假设系统噪声满足的协方差矩阵为常数；而在实际电力系统中，系统噪声的统计特性很难准确获取且是动态变化的，且系统噪声协方差矩阵的设置与状态估计器的性能密切相关。所以，若系统噪声协方差矩阵设置偏离其真实值，会严重降低状态估计精度，甚至导致状态估计器无法收敛于真实值。

技术实现要素：

发明目的：为了克服现有技术中存在的不足，本发明提供一种基于自适应容积卡尔曼滤波(adaptivecubaturekalmanfilter,ackf)的动态状态估计方法，能够提升发电机动态状态估计精度，为电网的安全稳定运行提供坚实的数据信息。

技术方案：为实现上述目的，本发明提供一种基于自适应容积卡尔曼滤波的动态状态估计方法，包括如下步骤：

(1)建立状态估计模型；

(2)初始化自适应容积卡尔曼滤波即ackf滤波初始值；

(3)计算k时刻的状态预测值

(4)计算k时刻的发电机状态预测误差协方差矩阵

(5)计算k时刻的量测预测值

(6)计算k时刻量测预测误差协方差矩阵pzz,k和状态预测值与量测预测值之间的交互协方差矩阵pxz,k；

(7)计算k时刻的滤波增益kk和状态估计值

(8)计算k时刻的状态估计误差协方差矩阵和系统噪声协方差矩阵qk；

(9)按照步骤(3)至(9)依据量测信息时间序列对电力系统发电机状态动态估计，直至k+1＞n时状态估计迭代停止，输出状态估计结果，否则返回步骤(3)继续计算。

进一步的，所述步骤(1)中建立状态估计模型的具体步骤如下：

发电机动态状态估计的状态方程与量测方程，表示为：

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

进一步的，所述步骤(2)中初始化ackf滤波初始值的具体步骤如下：

设定自适应容积卡尔曼滤波的发电机状态估计参数初始值，包含设定初始时刻状态变量值状态估计误差协方差矩阵控制变量值u0取为稳态运行值；设定系统噪声和量测噪声所满足的初始协方差矩阵q0，r0，以及最大估计时刻n。

进一步的，所述步骤(3)中计算k时刻的状态预测值的具体步骤如下：

运用自适应容积卡尔曼滤波的预测步，计算k时刻的状态预测值公式如下

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

进一步的，所述步骤(4)中计算k时刻的发电机状态预测误差协方差矩阵的具体步骤如下：

公式如下

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

进一步的，所述步骤(5)中计算k时刻的量测预测值的具体步骤如下：

基于状态预测值生成等权值cubature点xi,k；并通过状态估计模型的量测函数对状态预测值cubature点进行变换，获取对应点xi,k的量测预测值zi,k，公式如下：

基于量测预测值的cubature点zi,k，求解k时刻量测预测值计算公式如下

进一步的，所述步骤(6)中计算k时刻量测预测误差协方差矩阵pzz,k和状态预测值与量测预测值之间的交互协方差矩阵pxz,k的具体步骤如下：

公式如下

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

进一步的，所述步骤(7)中计算k时刻的滤波增益kk和状态估计值的具体步骤如下：

计算k时刻的滤波增益kk，并基于量测值zk对状态预测值进行修正更新，计算公式分别为

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

进一步的，所述步骤(8)中计算k时刻的状态估计误差协方差矩阵和系统噪声协方差矩阵qk的具体步骤如下：

求解更新k时刻的状态估计误差协方差矩阵，公式如下

式中表示k时刻的状态估计误差协方差矩阵，上标t表示矩阵的转置运算；

利用渐消记忆指数加权sage-husa噪声协方差估计器，动态修正更新k时刻系统噪声协方差矩阵qk，系统噪声协方差估计器的形式如下

dk-1＝(1-b)/(1-b^k)

式中b表示遗忘因子，且0＜b＜1，当系统状态的变化程度越大时b的取值就越大；dk-1为k-1时刻系统噪声协方差矩阵估计器调节参数。

本发明能够解决现有动态状态估计器针对系统噪声设置不当引起的性能下降问题，提升动态状态估计精度。

有益效果：本发明与现有技术相比具有以下优点：

本发明通过在容积卡尔曼滤波中引入渐消记忆指数加权的sage-husa噪声统计估计器，动态估计校正时变系统噪声的均值与方差，可有效避免系统噪声协方差矩阵设置不当对状态估计精度的影响，提升发电机状态变量的估计精度。另外，该方法流程清晰，实现简单方便，能够更好的满足发电机运行状态的监测需求。

附图说明

图1是本发明的方法流程图；

图2是具体实施例中ieee10机39节点系统结构图；

图3是具体实施例中利用传统ckf方法和本发明方法对发电机功角与角速度的动态估计结果对比图；

图4是具体实施例中利用传统ckf方法和本发明方法对发电机暂态电动势的动态估计结果对比图；

图5是具体实施例中利用传统ckf方法和本发明方法对发电机状态估计结果误差对比图。

具体实施方式

下面结合附图对本发明作更进一步的说明。

如图1所示，运用本发明方法对实施例测试系统动态变量进行估计，其包含如下步骤：

(1)状态估计模型建立

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

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

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

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

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

(5)计算k时刻发电机状态预测误差协方差矩阵求解计算公式如下

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

(6)基于状态预测值生成等权值cubature点xi,k；并通过状态估计模型的量测函数对状态预测值cubature点进行变换，获取对应点xi,k的量测预测值zi,k，计算公式如下

zi,k＝h(xi,k,uk),i＝1,…,2n

(7)基于量测预测值的cubature点zi,k，求解k时刻量测预测值计算公式如下

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

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

(9)计算k时刻的滤波增益kk，并基于量测值zk对状态预测值进行修正更新，计算公式分别为

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

(10)求解更新k时刻的状态估计误差协方差矩阵，计算公式如下

式中表示k时刻的状态估计误差协方差矩阵，上标t表示矩阵的转置运算。

(11)利用渐消记忆指数加权sage-husa噪声协方差估计器，动态修正更新k时刻系统噪声协方差矩阵qk，系统噪声协方差估计器的形式如下

dk-1＝(1-b)/(1-b^k)

式中b表示遗忘因子，且0＜b＜1，当系统状态的变化程度越大时b的取值就越大；dk-1为k-1时刻系统噪声协方差矩阵估计器调节参数。

(12)按照(2)-(11)计算步骤依据量测信息时间序列对电力系统发电机状态动态估计，直至k+1＞n时状态估计迭代停止，输出状态估计结果。

实施例：

(a)模型建立

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

式中：δ表示发电机功角，rad；ω和ω0分别为电角速度与同步转速，pu；e'q和e'd分别表示发电机q轴和d轴的暂态电动势；h表示发电机惯性常数，tm和te分别表示发电机机械功率与电磁功率，其中te＝pe/ω；kd表示阻尼因子，efd为定子励磁电压；t′d0和t′q0表示发电机机在d-q坐标系下的开路时间常数；xd和x'd分别表示发电机d轴同步电抗与暂态电抗，xq和x'q分别为发电机q轴同步电抗与暂态电抗；id和iq分别表示发电机d轴和q轴的定子电流。

对电力系统发电机动态变量进行动态估计时，状态估计向量为x＝(δ,ω,e'q,e'd)^t；选取发电机机械功率、定子励磁电压以及定子r轴和i轴的电流ir,ii为控制向量，即u＝(tm,efd,ir,ii)^t；选取发电机绝对功角、发电机角速度及发电机定子r轴和i轴的电压er,ei作为量测值，即量测向量为

z＝(δ,ω,er,ei)^t

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

(b)实施例分析

为了验证本发明所提出基于自适应容积卡尔曼滤波的动态状态估计方法的有效性和实用性，本实施例选取ieee10机39节点系统作为测试系统，系统结构见图2。在算法进行验证时，以系统中发电机g8的状态变量作为估计对象，其中发电机采用四阶模型。发电机g8的惯性时间参数为24.3，阻尼因子为2，d轴和q轴开路瞬变时间常数分别为6.7和0.41。为模拟发电机暂态特征，假定在第20周波时，节点16-节点21支路发生三相金属性短路故障，故障持续6周期(采样周期为0.02s)后消失。

运用bpa软件模拟pmu数据采集，获取发电机运行真实值。量测数据值由真实值叠加随机噪声形成。本发明在进行仿真实验时取前500周波(1周波为0.02s)量测值进行算法验证，即n为500。估计时状态变量的初值选取上一时刻的静态值，系统噪声的初始协防矩阵设置为q0＝10^-5i4×4，其真实值为q0＝10^-6i4×4量测噪声的初始协方差矩阵设置为r0＝10^-6i4×4。

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

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

对上述实施例系统，分别运用传统容积卡尔曼滤波算法(ckf)(其所需的相关参数值和本发明方法的参数初值相同)，以及本发明提出的自适应容积卡尔曼滤波方法(ackf)进行测试。

分别运用ckf方法和本发明方法对发电机g8的状态变量进行估计，其中发电机功角和转速估计结果如图3，图4给出了ckf和本发明方法对发电机g8的暂态电动势估计结果对比；两种方法的状态估计误差对比见图5。这些估计结果证实了本发明所提的自适应容积卡尔曼滤波方法由于能够动态修正调整系统噪声协方差矩阵，可以有效避免因系统噪声协方差矩阵设置不当引起的估计器性能下降，提升发电机动态状态估计精度。因此，本发明方法能够更好的满足电力系统发电机状态在线监测的需求。