一种基于残差归一化的权函数最小二乘状态估计方法与流程

本发明属于电力系统调度自动化领域，具体涉及一种基于残差归一化的权函数最小二乘状态估计方法。

背景技术：

电力系统状态估计是能量管理系统的重要组成部分，是运行电力系统其它应用软件的基础，其结果直接影响电网调度的智能化分析与决策。自1970年Schweppe首次将状态估计引入电力系统以来，国内外学者和工程人员对状态估计进行了大量、深入的研究和实践，这期间出现了各种各样的状态估计方法。

目前，加权最小二乘估计法(Weighted least squares,WLS)作为在状态估计中应用最为广泛的方法之一，其模型简单、易于求解，但不具有抗干扰性，因此出现了抗差状态估计。现有的抗差状态估计方法包括指数型目标函数抗差状态估计法(Maximum exponential square,MES)、基于标准化残差倒数的权函数最小二乘状态估计、以及指数加权最小二乘抗差状态估计法(Exponential function weighted least squares,EFWLS)等。但是这些方法在计算中由于标准化残差的计算会使每次迭代消耗大量时间，从而影响了它们在实际系统中的应用。

基于此，曾有提出的方法可在一定程度上提高计算效率，但求解过程复杂，对一些大型的电力系统计算量太大，因此不易得到广泛运用。也有人提出了指数权函数的状态估计算法(Exponential least absolute value,E-LAV)，提高了计算效率，但抗不良杠杆量测的能力较低。因此，需要为电力系统状态估计提供一种既能满足计算效率要求，又有良好抗差性能的估计方法，以满足实际工程运用的需要。

技术实现要素：

本发明的目的是针对现有状态估计研究方法的不足，提出一种计算效率高、收敛性好、且能有效抑制不良杠杠量测的基于残差归一化的权函数最小二乘状态估计方法(Residuals normalized weighted least squares,RNWLS)。

实现本发明的技术方案是：一种基于残差归一化的权函数最小二乘状态估计方法，利用计算机，通过程序，首先输入数据终端采集到的任一时间断面的SCADA数据、网络结构及参数信息并初始化，然后计算网络的节点导纳矩阵，形成零注入等式约束方程，接着综合考虑电压幅值量测方程、注入功率量测方程、支路功率量测方程，以节点电压的幅值与相角为状态变量，计算相应残差、雅克比矩阵以及权函数，最后更新状态变量，进行收敛性判断，来实现电网的状态估计。所述方法的具体步骤如下：

(1)输入基础数据及初始化

1)输入基础数据

首先输入任一时间断面下电网的数据采集与监视控制系统即SCADA数据、网络结构及参数信息；所述电网的SCADA数据包括节点电压幅值、节点注入功率以及支路功率；所述网络结构及参数信息包括网络元件的电阻、电抗、电纳、额定电压以及功率基准；所述网络元件包括线路以及变压器。

2)参数初始化

设置m阶单位矩阵R^-1，所述单位矩阵R^-1对角线元素全为1，非对角线元素全为0，m为状态估计中实际的量测变量个数；设置状态估计时节点电压幅值的初值均为标幺值1，相角均为0；所述网络元件的阻抗和电纳均归算为标幺值；节点注入功率以及支路功率也归算至标幺值；初始化最大迭代次数Tmax为40～60；收敛精度ε为10^-3～10^-5，以及权重计算时小残差的检测门槛值r_min＝0.002，并设置迭代次数time＝1。

(2)计算节点导纳矩阵

第(1)步完成后，利用公式(1)计算节点导纳矩阵Y：

式中，节点导纳矩阵的对角线元素表示为Y_ii，节点导纳矩阵的非对角线元素表示为Y_ij，均由上式计算得到；y_ij为节点i和节点j之间的支路阻抗z_ij的倒数；符号j∈i表示节点j和节点i直接相连，且当节点i有接地支路时还应包括j＝0时的情况。

(3)形成零注入等式约束方程

第(2)步完成后，选取量测量中既不是发电机节点也不是负荷节点的节点，将所述节点的注入功率作为零注入等式约束方程，表示为：

c(x)＝0 (2)

式中，c(x)是零注入节点的量测方程，其表示的节点注入有功功率和无功功率值可由公式(3)计算得到；x是n维的状态量，n是状态估计中实际的状态变量个数。c(x)的计算公式为：

式中，P_i′为节点i注入的有功功率；Q_i′为节点i注入的无功功率；u_i为节点i的电压幅值，u_j为节点j的电压幅值；θ_i为节点i的电压相角，θ_j为节点j的电压相角，且θ_ij＝θ_i-θ_j，表示节点i和节点j的电压相角差；G_ij、B_ij分别为节点导纳矩阵中对应节点i和j之间元素的实部和虚部；N为电网的节点个数。

(4)计算残差、雅克比矩阵及权函数

第(3)步完成后，以电网中节点电压幅值、除去零注入节点后的节点注入功率以及支路功率为量测量，计算电网中量测量的残差、雅克比矩阵及权函数。具体步骤内容如下：

I)计算残差

基于公式(4)计算状态估计中各量测量的参数：

r＝z-h(x) (4)

式中，z是m维的量测量，m是状态估计中实际的量测变量个数；x是n维的状态量，n是状态估计中实际的状态变量个数。h(x)是量测方程，包括注入功率对应的量测方程、线路支路功率对应的量测方程、变压器支路功率对应的量测方程和节点电压对应的量测方程；h(x)可由公式(5)-公式(8)计算得到；r是量测残差。

基于公式(5)得到除去零注入节点后的节点注入功率对应的量测方程：

式中，P_i为节点i注入的有功功率；Q_i为节点i注入的无功功率；u_i为节点i的电压幅值，u_j为节点j的电压幅值；θ_i为节点i的电压相角，θ_j为节点j的电压相角，且θ_ij＝θ_i-θ_j，表示节点i和节点j的电压相角差；G_ij、B_ij分别为节点导纳矩阵中对应节点i和j之间元素的实部和虚部。

基于公式(6)得到线路支路功率对应的量测方程：

式中，P_ij、Q_ij分别为线路支路节点i侧的有功功率、无功功率，P_ji、Q_ji为支路节点j侧的有功功率、无功功率。u_i为节点i的电压幅值，u_j为节点j的电压幅值；θ_i为节点i的电压相角，θ_j为节点j的电压相角，且θ_ij＝θ_i-θ_j，表示节点i和节点j的电压相角差。g为线路电导，b为线路电纳，y_c为线路对地电纳。

基于公式(7)得到变压器支路功率对应的量测方程：

式中，P_ijk、Q_ijk分别为变压器支路节点i侧的有功功率、无功功率，P_jik、Q_jik分别为变压器支路节点j侧的有功功率、无功功率。u_i为节点i的电压幅值，u_j为节点j的电压幅值。θ_i为节点i的电压相角，θ_j为节点j的电压相角，且θ_ij＝θ_i-θ_j，表示节点i和节点j的电压相角差。K为变压器非标准变比：j为标准侧，变比为1，i为非标准侧，变比为K；b_T为变压器标准侧的电纳。

基于公式(8)得到节点电压对应的量测方程：

U_i＝u_i (8)

式中，U_i、u_i均表示为节点i的电压幅值。

II)计算雅克比矩阵

基于公式(9)-公式(19)形成雅克比矩阵H和C。其中H为量测量的雅克比矩阵，包括注入功率、线路支路功率、变压器支路功率以及节点电压幅值所形成的雅克比矩阵元素。C为零注入等式约束的雅克比矩阵，由零注入节点的注入功率所形成的雅克比矩阵元素组成。

对于注入功率，基于公式(9)和(10)形成其雅克比矩阵元素：

式中，P_i为节点i注入的有功功率；Q_i为节点i注入的无功功率。u_i为节点i的电压幅值，u_j为节点j的电压幅值。θ_i为节点i的电压相角，θ_j为节点j的电压相角，且θ_ij＝θ_i-θ_j，表示节点i和节点j的电压相角差。G_ij、B_ij分别为节点导纳矩阵中对应节点i和j之间元素的实部和虚部。G_ii、B_ii分别为节点导纳矩阵中对应节点i处主对角线上的元素的实部和虚部。

对于线路i侧支路功率，基于公式(11)和(12)形成其雅克比矩阵元素：

式中，P_ij、Q_ij为线路支路节点i侧的有功功率、无功功率。u_i为节点i的电压幅值，u_j为节点j的电压幅值。θ_i为节点i的电压相角，θ_j为节点j的电压相角，且θ_ij＝θ_i-θ_j，表示节点i和节点j的电压相角差。g为线路电导，b为线路电纳，y_c为线路对地电纳。

对于线路j侧支路功率，基于公式(13)和(14)形成其雅克比矩阵元素：

式中，P_ji、Q_ji为线路支路节点j侧的有功功率、无功功率。u_i为节点i的电压幅值，u_j为节点j的电压幅值。θ_i为节点i的电压相角，θ_j为节点j的电压相角，且θ_ij＝θ_i-θ_j，表示节点i和节点j的电压相角差。g为线路电导，b为线路电纳，y_c为线路对地电纳。

对于变压器i侧支路功率，基于公式(15)和(16)形成其雅克比矩阵元素：

式中，P_ijk、Q_ijk为变压器支路节点i侧的有功功率、无功功率。u_i为节点i的电压幅值，u_j为节点j的电压幅值。θ_i为节点i的电压相角，θ_j为节点j的电压相角，θ_ij＝θ_i-θ_j，表示节点i和节点j的电压相角差。K为变压器非标准变比：j为标准侧，变比为1，i为非标准侧，变比为K；b_T为变压器标准侧的电纳。

对于变压器j侧支路功率，基于公式(17)和(18)形成其雅克比矩阵元素：

式中，P_jik、Q_jik为变压器支路节点j侧的有功功率、无功功率。u_i为节点i的电压幅值，u_j为节点j的电压幅值。θ_i为节点i的电压相角，θ_j为节点j的电压相角，θ_ij＝θ_i-θ_j，表示节点i和节点j的电压相角差。K为变压器非标准变比：j为标准侧，变比为1，i为非标准侧，变比为K；b_T为变压器标准侧的电纳。

对于节点电压，基于公式(19)形成其雅克比矩阵元素：

式中，U_i、u_i均为节点i的电压幅值，u_j为节点j的电压幅值，θ_i为节点i的电压相角，θ_j为节点j的电压相角。

零注入等式约束的雅可比矩阵为其中c(x)是零注入节点的量测方程，其表示节点的注入有功功率和无功功率值；x是n维的状态量，n是状态估计中实际的状态变量个数。

III)计算权函数

基于公式(20)计算权函数W的对角阵元素，计算公式为：

式中，r_min为小残差的检测门槛值，取0.002，w_i^*为量测i的权函数，为量测i的固定权重；r_i为量测i的残差。

(5)状态变量更新和收敛性判断

i)状态变量更新

第(4)步完成后，根据公式(21)计算状态变量的修正量Δx^(time)，然后更新状态变量，得到状态变量新值x^(time+1)＝x^(time)+Δx^(time)，time＝time+1。

式中，time为计算迭代次数；x^(time)为第time次迭代时的状态量；W为权函数对角阵，其对角元素等于权函数，即W_ii＝w_i^*。为量测量的雅可比矩阵，H^T为其转置。c(x^(time))为迭代值是x^(time)时的零注入等式约束，为零注入等式约束的雅可比矩阵，C^T为其转置。z-h(x^(time))表示迭代值为x^(time)时的残差；λ^(time)为第time次迭代时的拉格朗日乘子向量。

ii)收敛性判断

当状态变量的修正量Δx^(time)满足max(|Δx^(time)|)＜ε，则结束迭代计算，输出结果；当max(|Δx^(time)|)≥ε且迭代次数time≥Tmax，则停止迭代，输出“不收敛！”。

当max(|Δx^(time)|)≥ε且迭代次数time＜Tmax，使迭代次数time增加1，返回第(3)步，进行重新迭代计算。

本发明采用上述技术方案后，主要有以下效果：

与现有技术的权系数赋值方式比较，本发明采用基于残差归一化的权函数，在迭代过程中只需对权重进行改变，便可以自动消除坏数据的影响。抗差性能良好，且能有效抑制不良杠杆量测对估计结果的不良影响。

本发明在权函数计算时引入门槛值判据，即对于残差小于门槛值的量测其权重均取为1，而残差大于门槛值的量测其权重均小于1。从而充分利用了量测量中的有效信息，避免了因量测量权重的波动太大而导致状态估计不收敛，因此极大改善了坏数据与量测误差共存时的收敛性和数值稳定性。

本发明在权函数计算时直接以量测量的残差作为变量，避免了标准化残差的计算，因此大大提高了计算速度，节省了计算时间，非常适宜于工程实际的应用。

本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。

附图说明

图1为本发明的基于残差归一化的权函数最小二乘状态估计方法的流程示意图；

图2为3节点系统的网络参数和量测分布；

图3为IEEE-9节点系统的接线图。

具体实施方式

下面结合附图和实施例对本发明作进一步说明，但不应该理解为本发明上述主题范围仅限于下述实施例。在不脱离本发明上述技术思想的情况下，根据本领域普通技术知识和惯用手段，做出各种替换和变更，均应包括在本发明的保护范围内。

如1所示，一种基于残差归一化的权函数最小二乘状态估计方法的具体步骤如下：

(1)输入基础数据及初始化

1)输入基础数据

首先输入任一时间断面下电网的SCADA数据、网络结构及参数信息。即输入电网在任一时间断面下的节点电压幅值、节点注入功率以及支路功率，网络结构及参数信息为网络元件(包括线路、变压器)的电阻、电抗和电纳参数、以及额定电压、功率基准。

2)参数初始化

根据Anderson P M与Fouad A A所著“Power system control and stability”一文关于IEEE-9节点系统的标准数据，输入网络结构、负荷数据及相关参数信息。以潮流计算所得节点电压幅值、节点注入功率、支路功率为量测量，假设量测中无不良数据，且正常量测均带有一定的量测误差，该误差符合正态分布。设置m阶单位矩阵R^-1(即m阶单位矩阵R^-1的对角线元素全为1，非对角线元素全为0)，m为状态估计中实际的量测变量个数，为57；置状态估计时节点电压幅值的初值均为标幺值1，相角均为0；网络元件的阻抗、电纳均归算为标幺值；节点注入功率以及支路功率也归算至标幺值，基准功率设为100MVA；以及权重计算时小残差的检测门槛值r_min＝0.002，并设置迭代次数time＝1。对于初始化最大迭代次数和收敛精度的选取采用，本实施例优选Tmax＝40、ε＝10^-4。

(2)计算节点导纳矩阵

第(1)步完成后，计算该网络的节点导纳矩阵，计算公式为技术方案中的公式(1)。

根据IEEE-9节点系统的网络结构和参数信息，按照技术方案中的公式(1)，计算得到该网络的节点导纳矩阵Y为：

(3)形成零注入等式约束方程

第(2)步完成后，选取量测量中既不是发电机节点也不是负荷节点的节点，其注入功率作为零注入等式约束方程，按照技术方案中的公式(2)，以第1次迭代计算的结果举例，计算得到零注入功率c(x)为：

(4)计算残差、雅克比矩阵及权函数

第(3)步完成后，以电网中节点电压幅值、除去零注入节点后的节点注入功率、以及支路功率为量测量，计算电网中量测量的残差、雅可比矩阵及权函数，计算公式为技术方案中的公式(4)-公式(20)。

I)计算残差

根据残差及量测方程的定义，按照技术方案中的公式(4)-公式(8)，以第1次迭代计算的结果举例，计算得到量测量的残差r为：

II)计算雅克比矩阵

根据雅克比矩阵元素的定义，按照技术方案中的公式(9)-公式(19)计算雅可比矩阵H和C，以第1次迭代计算的结果举例，计算得到雅可比矩阵H和C为：

III)计算权函数

根据权函数的定义，按照技术方案中的公式(20)，以第1次迭代计算的结果举例，计算得到权函数W为：

(5)状态变量更新和收敛性判断

i)状态变量更新

第(4)步完成后，按照技术方案中的公式(21)计算状态变量的修正量Δx^(time)，然后更新状态变量，得到状态变量新值，即：x^(time+1)＝x^(time)+Δx^(time)，time＝time+1。

以第1次迭代计算的结果举例，按技术方案中的公式(21)，计算得到状态变量的修正量Δx⁽¹⁾为：

ii)收敛性判断

当状态变量的修正量Δx^(time)满足max(|Δx^(time)|)＜ε，则结束迭代计算，输出结果；当max(|Δx^(time)|)≥ε且迭代次数time≥Tmax，则停止迭代，输出“不收敛！”。当max(|Δx^(time)|)≥ε且迭代次数time＜Tmax，使迭代次数time增加1，返回第(3)步，进行重新迭代计算。

以第1次迭代计算的结果举例，此时，ε＝10^-4、Tmax＝40、max(|Δx⁽¹⁾|)＝0.16419＞ε，time＝1＜Tmax。根据收敛性判断，进行如下步骤：time＝time+1＝2，返回第(3)步，重新进行迭代计算。

根据前面的步骤，迭代3次后满足收敛条件，此时max(|Δx⁽³⁾|)＝7.4934E^-05＜ε，状态变量的估计结果如下表所示：

表1 IEEE-9节点系统状态变量估计结果

实验效果比较分析

为使本领域技术人员更好地理解本发明以及了解本发明相对现有技术的优点，申请人结合具体实施例进行进一步的阐释。

(1)抗差性能的比较

发明人将本发明的RNWLS与其他状态估计器进行比较，来测试RNWLS的抗差性。

1)3节点系统

以2所示3节点系统为例，1号量测为杠杆量测，并设置为不良数据，其他量测值取真值(潮流计算值)，其中P、Q分别表示有功和无功功率。用四种状态估计器(RNWLS、EFWLS、E-LAV、WLS)对上述算例进行仿真分析，结果如下所示：

表2 3节点系统中各方法估计结果的比较

其中，粗体字标识的数据为量测坏数据。由上述计算结果可见，在处理不良杠杆量测的能力上，本发明方法RNWLS和EFWLS估计的结果最接近于潮流真值，说明其能准确地估计出系统的真实状态，且1号不良杠杆量测能够被辨识并自动排除；而E-LAV估计的有功功率较接近潮流真值，无功功率则相差较大；WLS估计结果与潮流真值相差最大，表明其抗不良杠杆量测的能力有限，无法估计出系统的真实状态。

2)IEEE-118节点系统

为分析本发明方法RNWLS在规模更大系统上的抗差性能，以IEEE-118节点系统为例，分别对坏数据比例为0％～10％共11种情况进行试验。其中正常量测均带有一定的量测误差，且该误差符合正态分布；坏数据为潮流真值中相应比例电压幅值、有功功率、无功功率值的相反数(其中，有功功率、无功功率均在最大值处设置坏数据)。用四种状态估计器(RNWLS、EFWLS、E-LAV、WLS)对上述算例进行仿真分析，结果如下所示：

表3 IEEE-118节点系统不同估计方法电压幅值误差比较

其中，ΔU_max与ΔU_mean分别表示节点电压幅值的最大估计误差与平均估计误差，考虑到由SCADA所获取的电网潮流信息中节点的电压相角均为0，因此只对节点电压的幅值进行仿真分析。由上述计算结果可见，在不良数据从0％～10％变化的过程中，E-LAV与WLS所得结果误差较大；而RNWLS与EFWLS的估计结果，其误差都接近于0，表明在规模较大系统中，这两种估计方法都能排除不良数据，获得准确的估计结果，抗差性能良好。

(2)计算效率的比较

发明人为了进行效率比较，以IEEE-118节点系统为例，假设量测中无不良数据，且正常量测均带有一定的量测误差，该误差符合正态分布。用四种状态估计器(RNWLS、EFWLS、E-LAV、WLS)对上述算例进行仿真分析，得到平均每次迭代时间和总计算时间如下表所示：

表4 IEEE-118节点系统中各方法的计算效率比较

由上述计算结果可见，本发明方法RNWLS和EFWLS估计收敛性更好，但针对平均每次迭代时间而言，WLS估计最短，但其不具有抗差性；EFWLS估计最长；RNWLS估计与E-LAV估计居中。

综上所述，本发明所提出的RNWLS估计在抗差性能上近似于EFWLS估计，抗差性能良好，且能有效抑制不良杠杆量测对结果的不良影响；而在计算效率上近似于E-LAV估计，计算效率高，收敛性良好。从而有效集中了现有状态估计方法在抗差性能和计算效率上的优点，非常适宜于实际工程的应用。