一种电力系统动态状态估计多时段滚动优化方法与流程

本发明实施例属于电网优化领域，具体而言，涉及一种电力系统动态状态估计多时段滚动优化方法。

背景技术：

电力系统静态状态估计是ems的核心功能之一，一直受学术界和电力工业界的普遍重视，然而，随着pmu的广泛应用，电力系统动态状态估计逐渐得到关注。目前，电力系统动态状态估计主要是基于扩展卡尔曼滤波原理的各类递推估计算法，其特点是：须假设误差满足高斯分布。但在电力生产实际中，量测误差分布往往不确定且不服从这一分布，因而导致基于扩展卡尔曼滤波原理的估计方法结果不准确，甚至发散。因此，摒弃量测误差分布假设条件，研究对量测误差具有强鲁棒性的电力系统动态状态估计方法是十分必要的。

技术实现要素：

本发明针对现有技术中的不足，从多时段优化的角度，提出一种电力系统动态状态估计的滚动优化方法。

为实现上述发明目的，本发明采用了以下技术方案：

步骤1，初始化系统变量x，时刻权重wt，以及时段长度k，向前读取k个时刻量测数据，基于该k个时段量测数据，建立动态状态估计滚动优化模型。本文中所述的量测数据均为pmu量测数据。

步骤2，分析当前网络拓扑结构，读取系统参数，确定各量测量对应的量测方程，删除最早时刻量测数据与变量，用当前时刻量测数据和变量代替，使步骤1建立的模型随系统动态变化而实时更新；

模型随系统动态变化而实时更新的更新特征如下：

对任意时刻t原模型为：

上式中，下标k代表当前时刻数据，也代表了估计时段长度，vz和vc为残差向量，z为量测向量，h(x)为与z对应的量测函数，g为差分方程，为i时刻量测误差权重；w^ci为i时刻差分误差权重，做如下更新：

系统变量：

量测约束：

z^k＝[z¹z²…z^k]→z^k+1＝[z²z³…z^k+1]

更新内容为，xk和z^k中第一列为删除量，xk+1和z^k+1中最后一列为添加量，n表示发电机台数，xk+1和z^k+1分别表示t+1时刻的系统变量和量测量。

步骤3，将量测权重与时刻权重相乘得到计算权重，方法为：

量测矩阵：z^t＝[z¹z²…z^k]；

量测权重：wz＝[w¹w²…w^k]；

时刻权重：wt＝[0.1*10.1*2…0.1*k]；

计算权重：wx＝wz*diag(wt)；

上式中，z^k为当前时刻量测向量，wz为z^t对应的量测权重矩阵，wt为按时刻排列的权重向量，称为时刻权重，0.1为权重系数参考值可按需调整，wx为计算权重矩阵，定义为时刻权重和量测权重乘积。

步骤4，执行动态估计计算，求出该时段系统状态最优解。

步骤5，判断是否有不良数据，若有，则删除不良数据，返回步骤4，否则，输出最终结果，并更新数据库。

步骤6，新的采样数据上传以后进入下一时刻，回到步骤2。

相比于现有技术，本发明的该方法依据k时段量测量建立优化模型求解；k控制模型规模和精度，k值越大模型规模越大且结果精度越高，但计算速度降低；反之k越小模型规模越小，结果精度越低，但计算速度提高；时段内各时刻设定时刻权重，距当前时刻越近的时刻权重值越大；新时刻数据自动替代最早时刻数据，使模型完成系统状态跟踪的同时k值不发生变化。

综上所述，本发明的方法具有以下几大优势：

1.从优化的角度建立电力系统动态状态估计滚动优化模型，估计结果为系统k时段的最优解。

2.用k控制模型中变量和约束规模，可根据实际工况对计算精度和速度的需要，动态调整k值。

3.本方法实施无需假设量测误差满足某种分布，更符合电力生产实际。

附图说明

图1为电力系统动态状态估计多时段滚动优化方法算法框架。

图2为k值对状态估计模、计算速度和精度的影响示意图。

具体实施方式

为了能够更清楚地理解本发明的上述目的、特征和优点，下面结合附图和具体实施方式对本发明技术方案作进一步非限制性的详细描述。

见图1，一种电力系统动态状态估计多时段滚动优化方法，其步骤如下所示：

步骤1，初始化系统变量x，时刻权重wt，以及时段长度k，向前读取k个时刻量测数据，基于该k个时段量测数据，建立动态状态估计滚动优化模型，由于本模型从多时段优化角度建模，无需误差分布信息。

见图2，滚动优化估计时段k可动态调整，其大小用于协调计算速度与估计精度。本模型中，k值大小决定估计时刻长短，即决定模型规模。当k值增加时模型量测约束增加使得量测冗余度增加，因而精度也提高，但会降低计算速度。当k值减小时模型规模减小计算速度提升，但估计精度降低。因此，k的值可按实际需求动态调整。

步骤2，分析当前网络拓扑结构，读取系统参数，确定各量测量对应的量测方程，删除最早时刻量测数据与变量，用当前时刻量测数据和变量代替，使步骤1建立的模型随系统动态变化而实时更新。

模型随系统动态变化而实时更新的更新特征如下：

对任意时刻t原模型为：

上式中，下标k代表当前时刻数据，也代表了估计时段长度；vz和vc为残差向量；z为量测向量，h(x)为与z对应的量测函数；g为差分方程；为i时刻量测误差权重；w^ci为i时刻差分误差权重。做如下更新：

系统变量：

量测约束：

z^k＝[z¹z²…z^k]→z^k+1＝[z²z³…z^k+1]

更新内容为，xk和z^k中第一列为删除量，xk+1和z^k+1中最后一列为添加量，n表示发电机台数，xk+1和z^k+1分别表示t+1时刻的系统变量和量测量。

步骤3，将量测权重与时刻权重相乘得到计算权重，方法为：

量测矩阵：z^t＝[z¹z²…z^k]；

量测权重：wz＝[w¹w²…w^k]；

时刻权重：wt＝[0.1*10.1*2…0.1*k]；

计算权重：wx＝wz*diag(wt)；

上式中，z^k为当前时刻量测向量，wz为z^t对应的量测权重矩阵，wt为按时刻排列的权重向量，称为时刻权重，0.1为权重系数参考值可按需调整，wx为计算权重矩阵，定义为时刻权重和量测权重乘积。

步骤3的权重分配方式合理反映了不同时刻数据对估计结果的影响程度，肯定了当前时刻数据在估计过程中的重要性，也考虑了早期时刻系统状态对当前系统状态的影响。

步骤4，执行动态估计计算，求出该时段系统状态最优解。

估计计算用一范数内点法实现，过程如下：

对任意时刻原模型的等效模型为：

minw^t(l+u)

s.tz-f(x)+l-u＝0

l，u≥0

上式中，f(x)包含了量测方程和差分方程；

构造拉格朗日函数如下：

l＝w^t(l+u)-π^t(z-f(x)+l-u)-α^tl-β^tu

求偏导得到扰动kkt条件如下：

lπ＝z-f(x)+l-u＝0

ll＝w-π-α＝0

lu＝w+π-β＝0

上式中，μ为扰动因子，e为全1列向量；

得到修正方程如下：

dα＝-dπ

dβ＝dπ

上式中

s＝a^-1l+b^-1u

上式中，a，b，l，u分别为向量α，β，l，u的对角阵，和分雅可比和海森矩阵。

以上为修正方程推导过程，执行迭代求解过程如下：

计算互补间隙gap＝α^tl+β^tu。如果gap<ε，则停止计算，通过残差判断是否有不良数据，是则删除该数据并重新启动内点法计算过程，否则输出最终结果；如果gap>ε，则计算扰动因子μ＝σ*gap/2n。求h(x)，g(x)的雅可比矩阵和海森矩阵，得出修正方程，计算修正量[dxdπdαdβdldu]。计算原对偶步长θp和θd：

变量更新如下

然后再次执行此该求解过程。

步骤5，判断是否有不良数据，若有，则删除不良数据，返回步骤4，否则，输出最终结果，并更新数据库。

步骤6，新的采样数据上传以后进入下一时刻，回到步骤2。

需要指出的是，上述求解模型的方法仅为说明本发明的技术构思及特点，其目的在于让熟悉此项技术的人士能够了解本发明的内容并据以实施，并不能以此限制本发明的保护范围。凡根据本发明精神实质所作的等效变化或修饰，都应涵盖在本发明的保护范围之内。