本发明属于电力系统调度自动化
技术领域：
：，具体涉及一种反比型最大指数绝对值抗差状态估计方法(inverseproportionalmaximumexponentialabsolutevaluestateestimation,inv-meav)。
背景技术：
：：电力系统状态估计是能量管理系统的基础和核心。现在几乎每一个大型调度中心都安装了状态估计器，状态估计已成为电网安全运行的基石。自1970年国外学者首次提出状态估计以来，人们对状态估计的研究和应用已经有40多年的历史了，这期间涌现出了各种各样的状态估计方法。目前，在国内外应用最为广泛的状态估计是加权最小二乘法(weightedleastsquares,wls)。wls模型简洁，求解容易，但是其抗差性很差。为了增强抗差性，一般有两种方法。第一种是在wls估计之后加入不良数据辨识环节，例如最大正则化残差检验法(lnr)或估计辨识方法等；另一种是采用抗差状态估计方法。目前，国内外学者已经提出的抗差状态估计方法(robuststateestimation)包括加权最小绝对值估计(weightedleastabsolutevalue,wlav)、非二次准则法(ql、qc等)、以合格率最大为目标的状态估计(maximumnormalmeasurementrate,mnmr)以及指数型目标函数状态估计(maximumexponentialsquare,mes)等。但是这些抗差状态估计方法的估计性能仍有待提高。技术实现要素：本发明旨在至少在一定程度上解决上述技术问题之一或至少提供一种有用的商业选择。为此，本发明的一个目的在于提出一种抗差性好、计算效率高的反比型最大指数绝对值抗差状态估计(inverseproportionalmaximumexponentialabsolutevaluestateestimation,inv-meav)。本发明的技术方案是一种反比型最大指数绝对值抗差状态估计方法，包括以下步骤：步骤a.提供反比型最大指数绝对值抗差状态估计基本模型；步骤b.对所述反比型最大指数绝对值抗差状态估计基本模型引进辅助变量，变换得到反比型最大指数绝对值抗差状态估计等价模型；以及步骤c.利用原-对偶内点算法，对所述反比型最大指数绝对值抗差状态估计等价模型求解。在本发明的一个实施例中，所述反比型最大指数绝对值抗差状态估计基本模型为：其中：z∈rm为量测矢量，包括节点注入有功和无功、支路有功和无功以及节点电压幅值量测；x∈rn为状态矢量，包括节点电压幅值和平衡节点除外的其他各个节点相角；h:rn→rm为由状态矢量到量测矢量的非线性映射；ri为残差矢量r的第i个元素；g(x):rn→rc为零注入功率等式约束；wi为第i个量测量的权重，σ0和σ1为窗宽参数。在本发明的一个实施例中，所述步骤b包括：引进非负松弛变量u,v∈rm，变换得到的所述反比型最大指数绝对值抗差状态估计等价模型为：在本发明的一个实施例中，所述步骤c包括：步骤c1：令x为平启动状态变量；选择λ(0)＝π(0)＝0及u(0),v(0),α(0),β(0)>0；令中心参数ρ∈(0,1)及收敛判据ε＝10-3，置迭代计数器k＝0；步骤c2：计算对偶间隙gap＝αtv+βtu，判断是否收敛，若gap<ε,则转步骤c7，否则进入步骤c3；步骤c3：求解修正方程，以完成对原变量和对偶变量的修正,得到[dxtdλtdπt]t，dv,du，dα和dβ；步骤c4：计算原问题和对偶问题的修正步长θp和θd，其中：步骤c5：分别修正原问题和对偶问题的变量为：步骤c6：令迭代计数器k＝k+1，进入步骤c2；以及步骤c7：输出最优解，结束。在本发明的一个实施例中，所述步骤c3包括：步骤c31：计算扰动参数μ＝ρ·gap/2m；步骤c32：形成量测方程以及零注入功率约束对应的雅克比矩阵及形成量测方程以及零注入功率约束对应的海森矩阵及其中h(x)为状态矢量到量测矢量的映射，即为量测估计值，g(x)＝0为零注入功率约束；步骤c33：计算lx＝gtλ-htπ，lλ＝g(x)，lπ＝z-h(x)-u+v，步骤c34:计算γ＝z-h(x)-u+v+aa-bb，其中aa,bb∈rm，步骤c35：求解方程得到[dxtdλtdπt]t；步骤c36：求解dvi＝k1idπi+aai及dui＝k2idπi+bbi，其中k1i＝aivi-biui，k2i＝civi-diui；以及步骤c37：求解求解dαi＝pidui+pidvi-dπi+lvi,及dβi＝pidui+pidvi+dπi+lui。本发明反比型最大指数绝对值抗差状态估计方法(inverseproportionalmaximumexponentialabsolutevaluestateestimation,inv-meav)在估计过程中可有效抑制包括一致性不良数据在内的多个不良数据，显示了良好的抗差性，并具有很高的计算效率，非常适宜于实际工程应用。本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显。附图说明图1为本发明反比型最大指数绝对值抗差状态估计方法流程图。具体实施方式下面详细描述本发明的实施例，所述实施例旨在用于解释本发明，而不能理解为对本发明的限制。如图1所示，本发明反比型最大指数绝对值抗差状态估计方法(inverseproportionalmaximumexponentialabsolutevaluestateestimation,inv-meav)包括下列步骤：步骤a：提供反比型最大指数绝对值抗差状态估计(inverseproportionalmaximumexponentialabsolutevaluestateestimation,inv-meav)基本模型。具体地，本发明提出的inv-meav的基本模型如下所示s.t.g(x)＝0(2)r＝z-h(x)(3)式中：z∈rm为量测矢量，常包括节点注入有功和无功、支路有功和无功以及节点电压幅值量测等；x∈rn为包括节点电压幅值和相角的状态矢量(平衡节点相角除外)；h:rn→rm为由状态矢量到量测矢量的非线性映射；ri是残差矢量r的第i个元素；g(x):rn→rc为零注入功率等式约束；wi为第i个量测量的权重，σ0和σ1为窗宽参数。步骤b：对反比型最大指数绝对值抗差状态估计基本模型引进辅助变量，变换得到反比型最大指数绝对值抗差状态估计等价模型。具体地，inv-meav基本模型的目标函数虽然处处连续，但是在0处不可导，因而直接求解比较困难。可将模型(1)～(3)转化为一个处处可导的等价模型。引进变量ξ∈rm，使其满足|ri|≤ξii＝1,…,m(4)由式(1)和(4)可得，最大等价为最大。引进非负松弛变量l,k∈rm，将不等式(4)转化为两个等式约束为ri-li＝-ξii＝1,…,m(5)ri+ki＝ξii＝1,…,m(6)引进非负松弛变量u,v∈rm，使其满足ui＝li/2i＝1,…,m(7)vi＝ki/2i＝1,…,m(8)由式(5)～(8)，可得ri＝ui-vii＝1,…,m(9)ξi＝ui+vii＝1,…,m(10)将式(9)带入式(3)，可得等价量测约束条件为z-h(x)-u+v＝0(11)则式(1)～(3)给出的inv-meav基本模型的等价模型为s.t.g(x)＝0(13)z-h(x)-u+v＝0(14)u,v≥0(15)模型(12)～(15)即为inv-meav等价模型，该模型处处连续可导，可用基于梯度的方法来求解。步骤c：利用原-对偶内点算法，对所述反比型最大指数绝对值抗差状态估计等价模型求解。(1)inv-meav等价模型的求解方法注意到inv-meav的等价模型(12)～(15)是一个含有等式约束和不等式约束的最优化问题，适宜用原-对偶内点算法进行求解。为使本领域技术人员更好地理解本发明，首先给出详细的推导过程如下：引入拉格朗日函数式中：λ∈rc及π,α,β∈rm为拉格朗日乘子矢量。为取得最优值，根据kkt条件，可得式中：为有效解决以上问题，现代内点法引入扰动参数μ>0对式(22)、(23)进行松弛，从而得以上方程由牛顿法求解可得gdx＝-lλ(27)-hdx-du+dv＝-lπ(28)其中，由式(29)和(30)，可得将式(33)、(34)带入(31)、(32)可得令由式(35)和(36)可得dvi＝k1idπi+aai(37)dui＝k2idπi+bbi(38)式中：k1i＝aivi-biui,k2i＝civi-diui，将式(37)、(38)带入(28)可得hdx+qdπ＝γ(39)式中：q为rm×m的对角阵，其对角元素为qii＝-k1i+k2i；γ＝z-h(x)-u+v+aa-bb，aa,bb∈rm，aai,bbi与式(37)、(38)中同。根据式(39)、(26)及(27)，可得修正方程为求解式(40)可得dx,dλ和dπ；由式(37)、(38)可得dv和du；将所得结果带入式(33)、(34)可得dα和dβ，则迭代即可持续进行。(2)inv-meav等价模型的求解步骤在介绍inv-meav等价模型的求解推导过程之后，将求解步骤归纳如下：步骤c1:进行初始化，令x为平启动状态变量；选择λ(0)＝π(0)＝0及u(0),v(0),α(0),β(0)>0；令中心参数ρ∈(0,1)及确定收敛判据值，以及置迭代计数器为零。具体地，令x(0)∈rn代表由所有节点电压幅值和相角组成的的平启动状态变量(参考节点相角除外)；选择λ(0)＝π(0)＝0及u(0),v(0),α(0),β(0)>0，其中λ∈rc及π,α,β∈rm为拉格朗日乘子矢量，m为量测量的个数，而c为零注入功率约束的个数；令中心参数ρ∈(0,1)及收敛判据ε＝10-3，置迭代计数器k＝0。步骤c2：计算对偶间隙gap＝αtv+βtu，判断是否收敛。具体地，若gap<ε，则认为收敛，可直接进入步骤c7；否则为不收敛，进入步骤c3。步骤c3：求解修正方程，以完成对原变量和对偶变量的修正，得到[dxtdλtdπt]t，dv,du，dα和dβ。具体地，首先计算扰动参数μ＝ρ·gap/2m，然后求解式(40)得[dxtdλtdπt]t；求解式(37)、(38)得dv,du；求解(33)、(34)得dα,dβ。步骤c4：计算原问题和对偶问题的修正步长θp和θd，其中：步骤c5：分别修正原问题和对偶问题的变量为：步骤c6：令迭代计数器k＝k+1，进入步骤c2；以及步骤c7：输出最优解，结束。为使本领域技术人员更好地理解本发明以及了解本发明相对现有技术的优点，申请人结合具体实施例进行进一步的阐释。设定利用ieee标准系统检验基于原-对偶内点算法的inv-meav的性能。试验采用全量测，量测值通过在潮流计算的结果上叠加白噪声(均值为0，标准差为τ)来获得。对于电压量测，取τv＝0.005p.u.；对于功率量测，取τpq＝1mw/mvar。cpu为intel(r)core(tm)i54210m、主频为2.60ghz、内存4.00gb。1.抗差性能的比较将本发明的inv-meav与其他状态估计器进行比较，测试inv-meav的抗差性。在ieee-14系统上设置4个一致性不良数据(p1-2、q1-2、p1、q1)。所设置的不良量测值以及量测量的正确值如表1所示。表1cosh-meav对ieee14系统一致性不良数据的辨识table1estimationofconformingbaddatafortheieee14-bussystembycosh-meav作为对比，首先用广为应用的wls进行估计，并用lnr进行不良数据的辨识(简记为wls+lnr)。首次辨识的结果为：10个量测量的标准化残差大于门槛值(3.0)，这10个量测量被认为是可疑数据；其中标准化残差最大的量测量为p2-1，删去该量测后重新运行wls；此时发现p2的标准化残差最大。以上过程循环4次，4个良好的量测量被lnr误认为是可疑数据而被删去，但真正的不良数据仍然存在。可见，wls+lnr不能辨识一致性不良数据。应用inv-meav方法的估计结果如表1所示。可以发现，即使量测量中存在一致性不良数据，inv-meav的估计值与真值也可很好地吻合。在ieee其他系统的多次试验也表明inv-meav在估计的过程中可以自动抑制不良数据，具有良好的抗差性。尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在不脱离本发明的原理和宗旨的情况下在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。当前第1页12当前第1页12