CN104252571B - 基于多预测‑校正内点法的wlav抗差状态估计方法 - Google Patents

Abstract

基于多预测‑校正内点法的WLAV抗差状态估计方法，涉及一种电力***运行和控制方法。状态估计采用加权最小二乘状态估计，估计精度难于进一步提高。本发明包括以下步骤：1）获取电力***参数；2）获取检测数据；3）初始化；4）申请各量测量海森矩阵内存空间并求解；5）计算对偶间隙C_Gap=α^Tl+β^Tu，判断是否满足C_Gap＜ε或K＜K_max；6）预测步：设置扰动因子μ＝0，根据公式进行预测，求出仿射方向；7）校正步；8）判断校正次数计数器t＜4;9）对方程求解，得到Δλ_co，并校正；10）计算新的迭代步长若大于原步长，按公式Δλ_new＝Δλ_af+ωΔλ_co进行更新，t＝t+1，并执行步骤8)；否则，迭代计算器K＝K+1，并执行步骤5)。本技术方案减少迭代次数，提高处理速度；进一步提高算法的收敛特性。

Description

基于多预测-校正内点法的WLAV抗差状态估计方法

技术领域

本发明涉及一种电力***运行和控制方法。

背景技术

随着西电东送的实施、电力市场的迅速推行，特高压、远距离、交直流混合输电技术在我国电网中的发展迅速，电力***调度中心自动化水平也需要逐步提高。状态估计作为现代能量管理***(EMS)的核心与基石，通过对数据采集监控(SCADA)***传送的量测数据进行实时处理，可大幅度提高电力***数据精度。传统的基于最小二乘法及由最小二乘法派生出的快速解耦法状态估计程序在现场已有多年实际运行经验。然而由于最小二乘法是对服从高斯分布样本的最优估计，而实际SCADA中可能存在误差较大的不良数据，这些不良数据不服从高斯分布。因此，传统状态采用加权最小二乘状态估计精度难于进一步提高。

针对状态估计中如何抑制不良数据影响的问题，国内外学者提出了大量解决方法，主要包括以下2个方面：一是寻求新的不良数据辨识方法，将其从有效量测***中剔除；二是提出新的抗差估计器。当SCADA中存在不良杠杆量测或多个强相关的不良数据时，很难完全辨识出来，而抗差估计无需不良数据检测，能够有效抑制不良数据的影响。其中，WLAV状态估计由于能够保证多个量测残差为0，从而有效利用正确的量测值摒弃不良数据，近年来受到多者学者的关注。Irving及Owen最早将WLAV引入电力***状态估计，他们将WLAV转化为线性规划进行求解。Kotiuga与Vidyasagar认为，WLAV估计的本质是量测集的插值，故具有一定的排除不良数据能力。韦化提出基于原-对偶内点法(primal-dual interiorpoint method,PDIPM)的WLAV状态估计，该方法具有较好的数值稳定性，但存在迭代次数偏多的缺点。

发明内容

本发明要解决的技术问题和提出的技术任务是对现有技术方案进行完善与改进，提供基于多预测-校正内点法的WLAV抗差状态估计方法，达到抑制状态估计中不良数据影响、并提高处理速度的目的。为此，本发明采取以下技术方案。

基于多预测-校正内点法的WLAV抗差状态估计方法，其特征在于包括以下步骤：

1）获取电力***参数，包括：母线编号、基于多预测-校正内点法的WLAV抗差状态估计方法、补偿电容、输电线路的支路号、首端节点和末端节点编号、串联电阻、串联电抗、并联电导、并联电纳、变压器变比和阻抗；

2）获取检测数据，包括电压幅值、发电机有功功率、发电机无功功率、负荷有功功率、负荷无功功率、线路首端有功功率、线路首端无功功率、线路末端有功功率以及线路末端无功功率；

3）初始化，包括：对状态量设置初值(直角坐标系)、对拉格朗日乘子和罚因子设置初值、节点次序优化、形成节点导纳矩阵、恢复迭代计算器K＝1，设置最大迭代次数K_max，设置收敛精度要求ε；

4）申请各量测量海森矩阵内存空间并求解；

5）计算对偶间隙C_Gap=α^Tl+β^Tu，判断是否满足C_Gap＜ε或K＜K_max，若是，则输出计算结果，退出循环；若否，则执行步骤6)；

6）预测步：设置扰动因子μ＝0，根据以下公式进行预测，求出仿射方向：

式中：x为状态量，包含节点电压幅值和相角；l、u为松弛变量，即原变量；y、α、β为拉格朗日乘子，即对偶变量；ΔxΔx、Δy、Δw、Δl、Δu、Δα、Δβ分别为x、y、w、l、u、α、β的修正量;▽_xh(x)为h(x)的雅可比矩阵，为h(x)的海森矩阵,μ为扰动因子，L＝diag(l₁,…,l_m)、U＝diag(u₁,…,u_m)，A＝diag(a₁,…,a_m)、B＝diag(β₁,…,β_m)，e＝[1,…,1]^T，利用阻尼牛顿法对求解，由式解出Δλ，并对原、对偶变量进行修正：λ^(k+1)＝λ^(k)+αΔλ，α为迭代步长，其大小确定如下：

计算仿射方向的互补间隙：

C_gap ^af＝(α+α^*Δα_af)^T(l+α^*Δl_af)+(β+α^*Δβ_af)^T(u+α^*Δu_af)

动态估计中心参数：

μ＝min{(C_gap ^af/C_gap)³,0.1}C_gap/2m

7）校正步，设置校正次数计数器t＝1；

8）判断校正次数计数器t＜4，若是，则执行步骤9);若否，则执行步骤5);

9）对方程进行求解，得到Δλ_co，并采用动态选择校正方向在总的牛顿方向中所占的权重，即Δλ_new＝Δλ_af+ωΔλ_co；采用2阶段法对权重进行搜索，即：第1阶段，在[α_pα_d,1]中进行线性搜索，并允许在原、对偶空间寻找不同的最优权重；在确定高效的最优权值搜索子区间后，第2阶段在该子区间里找到最优的权值，与阶段1类似，在子区间进行线性搜索，最终找到最优的和

10）计算新的迭代步长若大于原步长，按公式Δλ_new＝Δλ_af+ωΔλ_co进行更新，校正次数计数器t＝t+1，并执行步骤8)；否则，迭代计算器K＝K+1，并执行步骤5)。

本技术方案保留非线性项的高阶信息，并利用多次校正的方式增大迭代步长，减少迭代次数，提高处理速度。同时，对校正方向进行加权处理，寻找出在总的牛顿方向中的最优比重，从而保证迭代点向中心轨迹靠拢，进一步提高算法的收敛特性。

作为对上述技术方案的进一步完善和补充，本发明还包括以下附加技术特征。

t_max=4。

对偶间隙C_Gap计算公式为α^Tl+β^Tu。

在第6步预测步时，首先对L(x,y,l,u,α,β)＝c^T(l+u)-y^T[z-h(x)+l-u]-α^Tl-β^Tu

按照KKT条件进行求导，得：

有益效果：本技术方案能够有效抑制不良数据的影响，有效利用正确的量测值摒弃不良数据，具有较好的数值稳定性。且保留非线性项的高阶信息，并利用多次校正的方式增大迭代步长，进一步减少迭代次数，提高处理速度。同时，对校正方向进行加权处理，寻找出在总的牛顿方向中的最优比重，从而保证迭代点向中心轨迹靠拢，进一步提高收敛特性。

附图说明

图1是本发明方法流程图。

图2（a）是本发明采用的线路Π形等值电路图。

图2（b）是本发明采用的变压器Π形等值电路图。

图3是本发明所应用的IEEE-14标准算例***。

图4是算例二中IEEE-57节点***不同内点法的迭代步长比较。

图5是算例二中IEEE-57节点***不同内点法的对偶间隙收敛过程比较。

图6是算例二中IEEE-57节点***某次迭代阶段1搜索中ω与α_p,α_d的关系。

具体实施方式

以下结合说明书附图对本发明的技术方案做进一步的详细说明。

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

1）获取电力***参数，包括：母线编号、基于多预测-校正内点法的WLAV抗差状态估计方法、补偿电容、输电线路的支路号、首端节点和末端节点编号、串联电阻、串联电抗、并联电导、并联电纳、变压器变比和阻抗；

2）获取检测数据，包括电压幅值、发电机有功功率、发电机无功功率、负荷有功功率、负荷无功功率、线路首端有功功率、线路首端无功功率、线路末端有功功率以及线路末端无功功率；

3）初始化，包括：对状态量设置初值(直角坐标系)、对拉格朗日乘子和罚因子设置初值、节点次序优化、形成节点导纳矩阵、恢复迭代计算器K＝1，设置最大迭代次数K_max，设置最大校正次数t_max，设置收敛精度要求ε；

4）申请各量测量海森矩阵内存空间并求解；

5）计算对偶间隙C_Gap，判断是否满足C_Gap＜ε，若是，则输出计算结果，退出循环；若否，则执行步骤6)；

6）预测步：设置扰动因子μ＝0，根据以下公式进行预测，求出仿射方向：

式中：x为状态量，包含节点电压幅值和相角；l、u为松弛变量，即原变量；y、α、β为拉格朗日乘子，即对偶变量；ΔxΔx、Δy、Δw、Δl、Δu、Δα、Δβ分别为x、y、w、l、u、α、β的修正量;▽_xh(x)为h(x)的雅可比矩阵，的海森矩阵,μ为扰动因子，L＝diag(l₁,…,l_m)、U＝diag(u₁,…,u_m)，A＝diag(a₁,…,a_m)、B＝diag(β₁,…,β_m)，e＝[1,…,1]^T，利用阻尼牛顿法对求解，由式解出Δλ，并对原、对偶变量进行修正：λ^(k+1)＝λ^(k)+αΔλ，α为迭代步长，其大小确定如下：

计算仿射方向的互补间隙：

C_gap ^af＝(α+α^*Δα_af)^T(l+α^*Δl_af)+(β+α^*Δβ_af)^T(u+α^*Δu_af)

动态估计中心参数：

μ＝min{(C_gap ^af/C_gap)³,0.1}C_gap/2m

7）校正步，设置校正次数计数器t＝1；

8）判断校正次数计数器t<t_tmax，若是，则执行步骤9);若否，则执行步骤5);

9）对方程进行求解，得到Δλ_co，并采用动态选择校正方向在总的牛顿方向中所占的权重，即Δλ_new＝Δλ_af+ωΔλ_co；采用2阶段法对权重进行搜索，即：第1阶段，在[α_pα_d,1]中进行线性搜索，并允许在原、对偶空间寻找不同的最优权重；在确定高效的最优权值搜索子区间后，第2阶段在该子区间里找到最优的权值，与阶段1类似，在子区间进行线性搜索，最终找到最优的和

10）计算新的迭代步长若大于原步长，按公式Δλ_new＝Δλ_af+ωΔλ_co进行更新，校正次数计数器t＝t+1，并执行步骤8)；否则，迭代计算器K＝K+1，并执行步骤5)。

在给定网络接线、支路参数和量测***的条件下，非线性量测方程可表示为：

z＝h(x)+r

式中，z为量测值矢量，绝大多数是通过遥测得到的，也有一小部分是人工设定的；h(x)为由基尔霍夫等基本电路定律所建立的量测函数；x为***状态变量；r为量测随机误差。设***有n个节点，以节点复电压实部与虚部作为状态变量，平衡节点不参与迭代。

如图2（a），图2（b）所示。在电力***状态估计中，量测量配置的类型要比常规潮流多，不仅包括了各节点的注入功率量测P_i、Q_i，还可以包括支路的功率量测P_ij、Q_ij、P_ji、Q_ji以及节点的电压幅值量测V_i，量测方程如下式所示：

节点i电压：

节点i注入功率：

线路i-j上始端功率：

线路i-j上末端功率：

变压器线路i-j上始端功率：

变压器线路i-j上末端功率：

以节点注入功率和线路i-j上的始端功率为例，量测矢量的雅克比矩阵H(x)对应的元素为：

节点i注入功率：

(j＝1,2,…,i-1,i+1,…,n)

(j＝1,2,…,i-1,i+1,…,n)

线路i-j上的始端功率

求出各量测方程的海森矩阵，该步只需在迭代前进行计算并存储，迭代过程中无需再进行计算。以线路i-j上的始端功率为例，量测矢量的海森矩阵He(x)对应的元素为：

式中，e_i、e_j分别为节点i、节点j复电压实部；G_ij、B_ij为导纳矩阵的元素；g、b、y_c为线路Π形模型中的参数；K为变压器的非标准变比；b_T为变压器标准侧的电纳。线路和变压器的Π形等值电路如图2所示。

为了验证本发明的有效性，与不同状态估计方法进行比较。其中WLS状态估计称为方法1，WLS+BD称为方法2，预测-校正内点法(predictor-collector PDIPM,PCPDIPM)称为方法3，多预测-校正内点法(multiple PCPDIPM,MPCPDIPM)称为方法4，线性原-对偶内点法称为方法5。为更符合实际物理意义，对各方法状态变量估计结果按极坐标形式给出，变换方法为：

下面介绍本发明的二个实施例：

本发明采用图3所示的IEEE-14节点的标准算例，分3种情况进行计算。1)***含不良数据，且正常量测均取真实值，用于分析不同方法对不良数据的辨识能力；2)***不含不良数据，且所有量测均带有随机误差，用于分析在符合正态分布的正常量测下不同方法的状态估计性能；3)***含不良数据，且量测量也带有随机误差，即情况1)、2)的综合，这也是实际中最常见的一种模式。表1给出了不同估计方法在情况1)下的不良数据辩识结果，表2给出了不同方法状态变量估计结果比较。

表1 不同方法情况1)下状态估计结果

表2 不同方法3）种情况下状态变量误差

由表1、2可知，方法2没有辨别出不良数据5、7号，而方法3和4所有量测量估计误差均小于1%，很好地辨识出了不良数据，满足工程要求。方法5虽然也均能辨识出不良数据，但估计误差明显大于方法3和4。以上试验表明，本发明(方法4)相较线性原-对偶内点法，较大程度地提高了求解精度，具有更强的抗差能力。

算例二：

为了验证本方法的计算效率，对14到3012节点的6个***进行了测试，其中对偶间隙收敛精度均设为10^-6。表3给出了不同***下，利用3种不同内点法的迭代次数，以及与PDIPM相比减少的迭代次数百分比。表4给出了相应方法的计算时间及与PDIPM相比减少的计算时间百分比。

表3不同方法迭代次数比较

注：表中Wp-2383，3012为波兰电网***；黑体部分为表中最好结果，下同

表4不同方法CPU计算时间比较

由表3、4可知，方法3与PDIPM相比，迭代次数至少能够减少13.64%，计算时间可减少13.84%，而方法4在方法3的基础上，通过多次校正，迭代次数可减少20.83%，计算时间可减少16.28%以上，具有更快的收敛速度。

进一步，图4给出了IEEE-57节点***不同内点方法的迭代步长，图5给出了对偶间隙的收敛过程。由图4可知，PDIPM在前3次迭代过程中步长较大，之后逐渐减小，方法3由于增加了校正步，相比PDIPM步长有所提高，但在迭代初期，迭代步长反而变小了，这是由于校正方向并不指向最优方向，而方法4采用加权多预测，并利用2阶段法找出最优校正权值保证校正方向指向最优方向，仅需迭代11次即可收敛，迭代步长始终保持较大值。由图5可知，PDIPM迭代16次后收敛，方法3对偶间隙收敛速度较快，仅需13次即可收敛，相比PDIPM迭代次数减少了18.75%，而采用方法4后，由于在每次迭代过程中，多次进行校正，，从而获得最优迭代步长，对偶间隙速度显著提高，仅迭代11次即可收敛，较PDIPM减少31.25%。

图6给出了某次迭代ω在区间[α_pα_d,1]寻找最优权值的搜索过程，其中α_p=0.5524,α_d=0.7606。为了能够取得原、对偶变量的最优步长，本文允许原、对偶变量取不同的最优权值。图6为第1阶段搜索过程，找出α_p最优子区间[0.47818,0.59414]，α_d的最优子区间为[0.82606,0.94202]。最优步长与权重大多为近似抛物线关系，因此第1阶段搜索时，设置点数不要求太多，只要能够保证得到的子区间为全局最优即可。在第2阶段子区间中设置搜索点数为4，实际上往往误差较大，本实施例将其设置为10，以保证搜索结果更接近真实最优权值。本次最终寻优结果为α_p=0.88857，α_d=0.80510，相比原迭代步长(即ω=1)均得到一定程度的增大。算例二的试验结果表明本方法采用的2阶段线性搜索法寻找最优权重的有效性。

以上图1-6所示的基于多预测-校正内点法的WLAV抗差状态估计方法是本发明的具体实施例，已经体现出本发明实质性特点和进步，可根据实际的使用需要，在本发明的启示下，对其进行形状、结构等方面的等同修改，均在本方案的保护范围之列。

Claims (4)

1.一种基于多预测-校正内点法的WLAV抗差状态估计方法，其特征在于包括以下步骤：

1)获取电力***参数，包括：母线编号、一种基于多预测-校正内点法的WLAV抗差状态估计方法、补偿电容、输电线路的支路号、首端节点和末端节点编号、串联电阻、串联电抗、并联电导、并联电纳、变压器变比和阻抗；

2)获取检测数据，包括电压幅值、发电机有功功率、发电机无功功率、负荷有功功率、负荷无功功率、线路首端有功功率、线路首端无功功率、线路末端有功功率以及线路末端无功功率；

3)初始化，包括：对状态量设置初值(直角坐标系)、对拉格朗日乘子和罚因子设置初值、节点次序优化、形成节点导纳矩阵、恢复迭代计算器K＝1，设置最大迭代次数K_max，设置最大校正次数t_max，设置收敛精度要求ε；

4)申请各量测量海森矩阵内存空间并求解；

5)计算对偶间隙C_Gap，判断是否满足C_Gap＜ε，若是，则输出计算结果，退出循环；若否，则执行步骤6)；

6)预测步：设置扰动因子μ＝0，根据以下公式进行预测，求出仿射方向：

$\left[\begin{array}{cccccc}A& L& 0& 0& 0& 0\\ 0& -I& 0& 0& 0& -I\\ 0& 0& B& U& 0& 0\\ 0& 0& 0& -I& 0& I\\ 0& 0& 0& 0& {\&dtri;}_x^{2}h\left(x\right)y& {\&dtri;}_x^{T}h\left(x\right)\\ 0& 0& 0& 0& {\&dtri;}_{x}h\left(x\right)& H\end{array}\right]\left[\begin{array}{c}\mathrm{\&Delta;}\mathrm{\&alpha;}\\ \mathrm{\&Delta;}l\\ \mathrm{\&Delta;}\mathrm{\&beta;}\\ \mathrm{\&Delta;}u\\ \mathrm{\&Delta;}x\\ \mathrm{\&Delta;}y\end{array}\right]=\left[\begin{array}{c}-ALe\\ -c+y+\mathrm{\&alpha;}\\ -BUe\\ -c-y+\mathrm{\&beta;}\\ -{\&dtri;}_{x}h\left(x\right)y\\ -L_y+M\end{array}\right]+\left[\begin{array}{c}\mathrm{\&mu;}e\\ 0\\ \mathrm{\&mu;}e\\ 0\\ 0\\ 0\end{array}\right]+\left[\begin{array}{c}-\mathrm{\&Delta;}A\mathrm{\&Delta;}Le\\ 0\\ -\mathrm{\&Delta;}B\mathrm{\&Delta;}Ue\\ 0\\ 0\\ 0\end{array}\right]$

式中：x为状态量，包含节点电压幅值和相角；l、u为松弛变量，即原变量；y、α、β为拉格朗日乘子，即对偶变量；Δx、Δy、Δw、Δl、Δu、Δα、Δβ分别为x、y、w、l、u、α、β的修正量；为h(x)的雅可比矩阵，为h(x)的海森矩阵,μ为扰动因子，L＝diag(l₁,…,l_m)、U＝diag(u₁,…,u_m)，A＝diag(α₁,…,α_m)、B＝diag(β₁,…,β_m)，e＝[1,…,1]^T，利用阻尼牛顿法对求解，由式解出Δλ，并对原、对偶变量进行修正：λ^(k+1)＝λ^(k)+αΔλ，α为迭代步长，其大小确定如下：

${\mathrm{\&alpha;}}_p=0.9995min\{\underset{i}{min}(\frac{-l_i}{{\mathrm{\&Delta;l}}_i},{\mathrm{\&Delta;l}}_i<0;\frac{-u_i}{{\mathrm{\&Delta;u}}_i},{\mathrm{\&Delta;u}}_i<0),1\}$

${\mathrm{\&alpha;}}_d=0.9995min\{\underset{i}{min}(\frac{-{\mathrm{\&alpha;}}_i}{{\mathrm{\&Delta;\&alpha;}}_i},{\mathrm{\&Delta;\&alpha;}}_i<0;\frac{-{\mathrm{\&beta;}}_i}{{\mathrm{\&Delta;\&beta;}}_i},{\mathrm{\&Delta;\&beta;}}_i<0),1\}$

计算仿射方向的互补间隙：

C_gap ^af＝(α+α^*Δα_af)^T(l+α^*Δl_af)+(β+α^*Δβ_af)^T(u+α^*Δu_af)

动态估计中心参数：

μ＝min{(C_gap ^af/C_gap)³,0.1}C_gap/2m

7)校正步，设置校正次数计数器t＝1；

8)判断校正次数计数器t<t_tmax，若是，则执行步骤9)；若否，则执行步骤5)；

9)对方程进行求解，得到Δλ_co，并采用动态选择校正方向在总的牛顿方向中所占的权重，即Δλ_new＝Δλ_af+ωΔλ_co；采用2阶段法对权重进行搜索，即：第1阶段，在[α_pα_d,1]中进行线性搜索，并允许在原、对偶空间寻找不同的最优权重；在确定高效的最优权值搜索子区间后，第2阶段在该子区间里找到最优的权值与阶段1类似，在子区间进行线性搜索，最终找到最优的和

10)计算新的迭代步长若大于原步长，按公式Δλ_new＝Δλ_af+ωΔλ_co进行更新，校正次数计数器t＝t+1，并执行步骤8)；否则，迭代计算器K＝K+1，并执行步骤5)。

2.根据权利要求1所述的一种基于多预测-校正内点法的WLAV抗差状态估计方法，其特征在于:t_max＝4。

3.根据权利要求1所述的一种基于多预测-校正内点法的WLAV抗差状态估计方法，其特征在于:对偶间隙C_Gap计算公式为α^Tl+β^Tu。

4.根据权利要求1所述的一种基于多预测-校正内点法的WLAV抗差状态估计方法，其特征在于:在第6步预测步时，首先对L(x,y,l,u,α,β)＝c^T(l+u)-y^T[z-h(x)+l-u]-α^Tl-β^Tu

按照KKT条件进行求导，得：

$\left\{\begin{array}{c}{L}_{\mathrm{\&alpha;}}=-l\&DoubleRightArrow;L_{\mathrm{\&alpha;}}^{\mathrm{\&mu;}}=ALe-\mathrm{\&mu;}e=0\\ L_l=c-y-\mathrm{\&alpha;}=0\\ L_{\mathrm{\&beta;}}=-u\&DoubleRightArrow;L_{\mathrm{\&beta;}}^{\mathrm{\&mu;}}=BUe-\mathrm{\&mu;}e=0\\ L_u=c+y-\mathrm{\&beta;}=0\\ L_x={\&dtri;}_{x}h\left(x\right)y=0\\ L_y=-z+h\left(x\right)-l+u=0\end{array}\right.$

式中，z为量测值矢量。