一种计及网络攻击和参数不确定性的动态状态估计方法与流程

本发明属于电力系统分析和监测技术领域，尤其涉及一种计及网络攻击和参数不确定性的动态状态估计方法。

背景技术：

近年来，随着能源资源大范围优化配置格局的初步形成、电力市场化改革的稳步推进、新能源开发步伐的加快、“建设坚强智能电网”举措的提出，中国电网结构日益庞大，运行方式日趋复杂，保障电网的安全经济运行意义重大，任务艰巨。电力系统调度中心依靠静态状态估计可以掌握电力系统实时运行状态，而分析和预测系统的运行趋势，对运行中发生的各种问题提出对策，则需要依靠兼备预测功能的动态状态估计。

目前，电力系统动态状态估计主要以ekf及其改进方法为主，如计入非线性卡尔曼滤波，自适应预报动态状态估计，光滑增平面模糊控制动态状态估计等。上述这些方法在一定程度上改善了状态估计的结果。但是，需要指出的是传统基于ekf框架的动态状态估计方法无法计及网络攻击造成量测数据缺失的影响，且对模型的精度要求较高。但在实际电力系统应用中，不仅pmu量测单元易受到网络攻击，而且系统模型精确的参数和系统噪声统计特性往往比较难获取，其无疑会严重影响动态状态估计的结果，大大降低状态估计精度。

技术实现要素：

发明目的：本发明目的是针对现有技术存在的不足，设计一种计及网络攻击和模型参数不确定性的电力系统动态状态估计方法。

技术方案：本发明包括以下步骤：

(1)建立电力系统动态状态估计模型；

(2)初始化hekf-gm状态估计方法参数值；

(3)基于hekf预测步，计算t时刻的状态预测值与预测误差协方差pt|t-1；

(4)建立线性批处理回归模型，增加状态估计的量测冗余量；

(5)采用鲁棒投影统计方法，求取数据点在所有可能向量的投影值，检测步骤(4)线性回归模型中的网络攻击值；

(6)对步骤(4)中的线性批处理回归模型进行白噪化处理；

(7)计算迭代加权最小二乘法的初始权矩阵；

(8)利用迭代加权最小二乘法求解步骤(6)，获取状态估计值；

(9)计算t时刻的估计误差协方差；

(10)按照步骤(3)-(9)依据时间序列对系统状态进行动态估计，直至t+1＞n时迭代停止，输出状态估计结果。

所述步骤(1)中的状态估计模型包含系统方程和量测方程，分别表示为：

式中xt表示t时刻的状态变量，xt＝[δt，ωt]^t由发电机运行功角和电角速度组成，f(·)为发电机系统函数，yt∈r^m为t时刻的量测变量，h为量测输出矩阵，wt-1∈rⁿ，et∈r^m分别是系统噪声和量测噪声值，二者为高斯白噪声序列。

所述步骤(2)中的状态估计方法参数值包含估计初始值估计误差协方差p0|0，系统和量测噪声协方差矩阵分别为w0和r0及最大估计时刻n。

所述步骤(3)中t时刻的状态预测值与预测误差协方差pt|t-1的计算方法如下：

式中表示t-1时刻状态估计值，ft-1表示函数f(·)在处的雅可比矩阵，(·)^t表示矩阵转置运算，wt为t时刻系统噪声协方差矩阵。

所述步骤(4)中的线性批处理回归模型的具体形式如下：

式中i为单位矩阵，xt表示t时刻状态真实值，δt|t-1为状态真值与预测值xt差值，该表达式可进一步表示为

式中则所满足的协方差矩阵为

其中，lt可以通过柯列斯基分解获取，rt表示t时刻et满足的协方差矩阵。

所述步骤(5)中的鲁棒投影统计方法的原理如下

式中psi表示第i行对应的投影值，上标(·)^t表示矩阵转置，medt(·)为求取中值的运算。

所述步骤(6)中的白噪化处理方法为：两端同时乘以即

进一步整理，表示为

yt＝atxt+ηt

式中

所述步骤(7)中的初始权矩阵为q¹＝diag{q(rsi)}，其中q(rsi)＝ψ(rsi)/rsi，ψ(·)表示的函数为

式中c＝1.5为阈值，参数rsi的计算方法为

式中yt(i)表示t时刻量测值的第i行，ai为输出矩阵at第i行。

所述步骤(8)中的求解方法为

式中为t时刻第v次寻优迭代结果，q^(v)为第v次迭代的权矩阵。

所述步骤(9)中的估计误差协方差pt|t的计算公式如下

式中re，t的计算公式为

式中i为对应维度的单位矩阵，γ为参数不确定性约束上界。

有益效果：本发明提出的方法可有效抑制网络攻击和模型参数不确定性所引起的估计偏差，甚至发散等问题；有效提升了状态估计精度，具有较强的鲁棒性；为电力系统动态监测与分析提供了坚实的数据信息。

附图说明

图1为本发明的流程图；

图2为本发明采用不同方法对发电机功角的动态估计结果对比；

图3为本发明采用不同方法对发电机电角速度的动态估计结果对比。

具体实施方式

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

如图1所示，本发明包括以下步骤：

(1)建立电力系统动态状态估计模型

状态估计模型包含系统方程和量测方程，其形式可表示为：

式中xt表示t时刻的状态变量，xt＝[δt，ωt]^t由发电机运行功角和电角速度组成，f(·)为发电机系统函数，yt∈r^m为t时刻的量测变量，h为量测输出矩阵，wt-1∈rⁿ，et∈r^m分别是系统噪声和量测噪声值，二者为高斯白噪声序列。

(2)初始化hekf-gm(hinfinityextendedkalmanfilter，hekf-generalizedmaximumlikelihood，gm)状态估计方法参数值，即包含估计初始值估计误差协方差p0|0，系统和量测噪声协方差矩阵分别为w0和r0，以及最大估计时刻n。

(3)基于hekf预测步，计算t时刻的状态预测值与预测误差协方差pt|t-1，方法如下

式中表示t-1时刻状态估计值，ft-1表示函数f(·)在处的雅可比矩阵，(·)^t表示矩阵转置运算，wt为t时刻系统噪声协方差矩阵。

(4)结合状态预测值和量测值zt，建立线性批处理回归模型，增加状态估计的量测冗余量，具体形式如下：

式中i为单位矩阵，xt表示t时刻状态真实值，δt|t-1为状态真值与预测值xt差值，该表达式可进一步表示为

式中则所满足的协方差矩阵为

其中，lt可以通过柯列斯基分解获取，rt表示t时刻et满足的协方差矩阵。

(5)采用鲁棒投影统计方法，求取数据点h在所有可能向量u的投影值，用以检测步骤(4)线性回归模型中的网络攻击值，其原理如下

式中psi表示第i行对应的投影值，上标(·)^t表示矩阵转置，medt(·)为求取中值的运算。设定判定阈值d＝1.5，若则判定该行为奇异值，为了克服其对状态估计的影响，须降低其对应权重即

(6)对步骤(4)中线性批处理回归模型进行白噪化处理，两端同时乘以即

进一步整理，表示为

yt＝atxt+ηt

式中

(7)计算迭代加权最小二乘法的初始权矩阵q¹＝diag{q(rsi)}，其中q(rsi)＝ψ(rsi)/rsi，ψ(·)表示的函数为

式中c＝1.5为阈值(一般取1.5)，参数rsi的计算方法为

s＝1.4826·medi|rt(i)|，

式中yt(i)表示t时刻量测值的第i行，ai为输出矩阵at第i行。

(8)利用迭代加权最小二乘法求解步骤(6)中的方程，获取状态估计值，其计算方法为

式中为t时刻第v次寻优迭代结果，q^(v)为第v次迭代的权矩阵。

(9)计算t时刻的估计误差协方差pt|t，计算公式如下

式中re，t的计算公式为

式中i为对应维度的单位矩阵，γ为参数不确定性约束上界。

(10)按照步骤(3)-(9)依据时间序列对系统状态进行动态估计，直至t+1＞n时迭代停止，输出状态估计结果。

具体实施方法为：

(a)模型建立

同步发电机的经典二阶模型具体形式如下：

式中δ为发电机转子功角，rad；ω，ω0分别为发电机转子电角速度和同步转速，pu；pm和pe分别为发电机的机械功率和电磁功率，pu；tj和d分别为发电机参数中的惯性时间常数和阻尼系数。

对电力系统动态变量进行动态估计时，选取发电机的状态变量为x＝(δ，ω)^t，把发电机的机械功率和电磁功率作为已知的输入量，记为u＝(pm，pe)^t，此时发电机转子的运动方程将和外部网络解耦。则二阶模型对应的状态方程形式如下

式中δ为单位为度。

另一方面，随着同步相量量测单元(pmu)的快速推广和应用，使得发电机功角和电角速度的直接量测成为可能，所以此处的量测方程设置为

式中y为量测变量。

(b)实施例分析

为验证本发明所提出的hekf-gm状态估计方法的有效性和实用性，本发明选取某大区域电网中一台实际参数机组的扰动过程进行了仿真验证，发电机惯性时间常数tj取值为527.861729，阻尼因子d为2，故障设置在第40周波时，该发电机一出线回路发生三相短路故障，第43周波时短路故障消失。运用bpa软件模拟pmu设备进行量测数据采集，获取发电机运行真实值，量测数据值由真实值叠加随机噪声形成。本发明在进行仿真实验时取前300周波(1周波为0.02s)量测值进行算法验证，即n为300。

在算法进行验证时，取发电机功角和电角速度作为状态估计变量，并将调速器的作用考虑在内，发电机采用经典二阶模型。估计时状态变量初值取上一时刻的静态值，初始协方差矩阵p0|0取对应维度的单位矩阵。过程噪声和量测噪声所满足的协方差矩阵q，r真值如下

q＝diag(10^-6，10^-6)，r＝diag(10^-6，10^-6)

而在进行状态估计时，假设二者均存在不确定性，其取值均取为10^-4。

另外，设定在60-63时刻，系统受到网络攻击，造成量测序列y2(t)，t＝60，…63量测数据丢失，其值为0。

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

不同方法对发电机功角与电角速度的动态估计结果对比如图2和图3所示，从仿真结果可以看出，在0-40周波，系统处于稳定运行状态，此时三种方法均能高效追踪发电机运行动态变量，而本发明所提方法精度更高，因为其计及了模型参数不确定性所引起的状态估计偏差；但在40周波三相短路故障发生后，ekf和hekf方法仅能跟踪电角速度状态变量的大致趋势，虽然hekf方法较ekf方法有所改善，但误差依然较大，且迟滞现象严重，而本发明所提的hekf-gm方法依然能够准确跟踪状态变化。

另外，当系统受到网络攻击，造成量测序列y2(t)，t＝60，…63量测数据丢失时，ekf和hekf已经无法追踪发电机电角速度的变化趋势，而本发明所提出的hekf-gm方法能够较好的抑制网络攻击引起的状态精度下降，实现状态变量的准确估计，显示出hekf-gm方法具有更强的鲁棒性。

为了对不同算法状态估计结果做进一步对比分析，本发明采用平均相对估计误差和最大绝对误差xm作为指标进行算法间性能对比。

式中为k时刻第i个状态量的滤波值(i＝1，2)，为k时刻第i个状态量的真实值(bpa数据)，为平均相对估计误差，xm为最大绝对估计误差，n为总的采样周期数。

表1给出了不同算法对实施例中的系统动态估计结果的性能指标数据。从表中性能数据可以看出，在计及网络攻击和模型参数不确定情形下，本发明提出的hekf-gm方法，其各项性能指标均优于ekf和hekf方法，凸显了所提方法的优越性和实用性。

本发明提出的hekf-gm电力系统动态状态估计方法具有更好的鲁棒性，可有效抑制网络攻击和模型参数不确定性所引起的估计偏差，能够为电力系统动态监测与分析提供了坚实的数据信息，保障电力系统安全稳定运行。

表1不同算法动态估计结果指标