一种基于最大指数绝对值目标函数的抗差状态估计方法

摘要

本发明提出一种基于最大指数绝对值目标函数的抗差状态估计方法，包括步骤：提供基于最大指数绝对值目标函数的抗差状态估计基本模型；对所述基于最大指数绝对值目标函数的抗差状态估计基本模型引进辅助变量，得到了基于最大指数绝对值目标函数的抗差状态估计等价模型；以及利用原-对偶内点算法，对所述基于最大指数绝对值目标函数的抗差状态估计等价模型求解。算例分析表明，本发明具有很强的抗差性和很高的计算效率，具有良好的工程应用前景。

说明书

技术领域

本发明属于电力系统调度自动化领域，具体涉及一种基于最大指数绝对值目标函数的抗 差状态估计方法。

背景技术

电力系统状态估计是能量管理系统的基础和核心。现在几乎每一个大型调度中心都安装 了状态估计器，状态估计已成为电网安全运行的基石。自1970国外学者首次提出状态估计 以来，人们对状态估计的研究和应用已经有40多年的历史了，这期间涌现出了各种各样的 状态估计方法。

目前，在国内外应用最为广泛的状态估计是加权最小二乘法(Weighted least squares, WLS)。WLS模型简洁，求解容易，但是其抗差性很差。为了增强抗差性，一般有两种方法。 第一种是在WLS估计之后加入不良数据辨识环节，例如最大正则化残差检验法(LNR)或估 计辨识方法等；另一种是采用抗差状态估计方法(Robust state estimation)。目前，国内外学者 已经提出的抗差状态估计方法包括加权最小绝对值估计(Weighted least absolute value, WLAV)、非二次准则法(QL、QC等)、以合格率最大为目标的状态估计(Maximum normal  measurement rate,MNMR)以及指数型目标函数状态估计(Maximum exponential square,MES) 等。但是这些抗差状态估计方法普遍存在计算效率不够高的特点，从而在一定程度上影响了 它们在实际系统中的应用。

发明内容

本发明旨在至少在一定程度上解决上述技术问题之一或至少提供一种有用的商业选择。 为此，本发明的一个目的在于提出一种抗差性好、计算效率高的基于最大指数绝对值目标函 数的抗差状态估计方法(Maximum exponential absolute value state estimation,MEAV)。

根据本发明实施例的基于最大指数绝对值目标函数的抗差状态估计方法，包括步骤：步 骤A.提供基于最大指数绝对值目标函数的抗差状态估计基本模型；步骤B.对所述基于最大 指数绝对值目标函数的抗差状态估计基本模型引进辅助变量，变换得到基于最大指数绝对值 目标函数的抗差状态估计等价模型；以及步骤C.利用原-对偶内点算法，对所述基于最大指 数绝对值目标函数的抗差状态估计等价模型求解。

在本发明的一个实施例中，所述基于最大指数绝对值目标函数的抗差状态估计基本模型 为：s.t.g(x)=0，r=z-h(x)，其中：z∈R^m为量测矢量，包括节点 注入有功和无功、支路有功和无功以及节点电压幅值量测；x∈Rⁿ为状态矢量，包括节点电 压幅值和平衡节点除外的其他各个节点相角；h:Rⁿ→R^m为由状态矢量到量测矢量的非线性映 射；r_i为残差矢量r的第i个元素；g(x):Rⁿ→R^c为零注入功率等式约束；w_i为第i个量测量的 权重，同是为窗宽参数。

在本发明的一个实施例中，所述步骤B包括：引进非负松弛变量u,v∈R^m，变换得到的 所述基于最大指数绝对值目标函数的抗差状态估计等价模型为：s.t.g(x)=0，z-h(x)-u+v＝0，u,v≥0。

在本发明的一个实施例中，所述步骤C包括：步骤C1：令x为平启动状态变量；选 择λ⁽⁰⁾=π⁽⁰⁾=0及u⁽⁰⁾,v⁽⁰⁾,α⁽⁰⁾,β⁽⁰⁾>0；令中心参数ρ∈(0,1)及收敛判据ε=10^-3，置迭代计数器 k=0；步骤C2：计算对偶间隙Gap=α^Tv+β^Tu，判断是否收敛，若Gap＜ε,则转步骤C7，否则 进入步骤C3；步骤C3：求解修正方程，以完成对原变量和对偶变量的修正,得到 [dx^T dλ^T dπ^T]^T，dv,du，dα和dβ；步骤C4：计算原问题和对偶问题的修正步长θ_P和θ_D， 其中：${\theta }_P=0.9995\min \{\underset{i}{\min }(-\frac{v_i}{d{v}_i}:d{v}_i<0;-\frac{u_i}{d{u}_i}:d{u}_i<0),1\},$${\theta }_D=0.9995\min \{\underset{i}{\min }(-\frac{{\alpha }_i}{d{\alpha }_i}:d{\alpha }_i<0;-\frac{{\beta }_i}{d{\beta }_i}:d{\beta }_i<0),1\};$步骤C5：分别修正原问题和对偶问题的变量为：$\left(\begin{array}{c}{x}^{(k+1)}\\ v^{(k+1)}\\ u^{(k+1)}\end{array}\right)=\left(\begin{array}{c}{x}^{\left(k\right)}\\ v^{\left(k\right)}\\ u^{\left(k\right)}\end{array}\right)+{\theta }_P\left(\begin{array}{c}\mathrm{dx}\\ \mathrm{dv}\\ \mathrm{du}\end{array}\right),$$\left(\begin{array}{c}{\lambda }^{(k+1)}\\ {\pi }^{(k+1)}\\ {\alpha }^{(k+1)}\\ {\beta }^{(k+1)}\end{array}\right)=\left(\begin{array}{c}{\lambda }^{\left(k\right)}\\ {\pi }^{\left(k\right)}\\ {\alpha }^{\left(k\right)}\\ {\beta }^{\left(k\right)}\end{array}\right)+{\theta }_D\left(\begin{array}{c}\mathrm{d\lambda }\\ \mathrm{d\pi }\\ \mathrm{d\alpha }\\ \mathrm{d\beta }\end{array}\right);$步骤C6：令迭代计数器k=k+1，进入步骤C2；以及步骤C7：输出最优解，结束。

在本发明的一个实施例中，所述步骤C3包括：步骤C31：计算扰动参数μ=ρ·Gap/2m； 步骤C32：形成量测方程以及零注入功率约束对应的雅克比矩阵及形成量测方程以及零注入功率约束对应的海森矩阵及其中h(x)为状态矢量到 量测矢量的映射，即为量测估计值，g(x)=0为零注入功率约束；步骤C33：计算L_x=G^Tλ-H^Tπ， L_λ=g(x)，L_π=z-h(x)-u+v，$L_{v_i}=-{\omega }_i-{\pi }_i-{\alpha }_i,$$L_{u_i}=-{\omega }_i-{\pi }_i-{\beta }_i,$$L_{{\alpha }_i}^{\mu }={\alpha }_i{v}_i-\mu$及$L_{{\beta }_i}^{\mu }={\beta }_i{u}_i-\mu ,$其 中ω_i=exp(-(u_i+v_i)/w_i)，w_i为第i个量测量对应的权重；步骤C34:计算 γ=z-h(x)-u+v+AA-BB  ，其中AA,BB∈R^m，${\mathrm{AA}}_i=-a_i(v_i{L}_{v_i}+L_{{\alpha }_i}^{\mu })-b_i(u_i{L}_{u_i}+L_{{\beta }_i}^{\mu }),$${\mathrm{BB}}_i=-c_i(v_i{L}_{v_i}+L_{{\alpha }_i}^{\mu })-d_i(u_i{L}_{u_i}+L_{{\beta }_i}^{\mu }),$$\left(\begin{array}{cc}{a}_i& b_i\\ c_i& d_i\end{array}\right)={\left(\begin{array}{cc}{v}_i{\omega }_i{w}_i^{-1}+{\alpha }_i& v_i{\omega }_i{w}_i^{-1}\\ u_i{\omega }_i{w}_i^{-1}& u_i{\omega }_i{w}_i^{-1}+{\beta }_i\end{array}\right)}^{-1},$z∈R^m为量测矢量；步骤 C35：求解方程$\left(\begin{array}{ccc}{▿}^{2}g\left(x\right)\lambda -{▿}^{2}h\left(x\right)\pi & G^T& -H^T\\ H& 0& Q\\ G& 0& 0\end{array}\right)\left(\begin{array}{c}\mathrm{dx}\\ \mathrm{d\lambda }\\ \mathrm{d\pi }\end{array}\right)=\left(\begin{array}{c}-L_x\\ \gamma \\ -L_{\lambda }\end{array}\right)$得到[dx^T dλ^T dπ^T]^T；步骤C36：求 解dv_i=k_1idπ_i+AA_i及du_i=k_2idπ_i+BB_i，其中k_1i=a_iv_i-b_iu_i，k_2i=c_iv_i-d_iu_i；以及步骤C37：求解 $d{\alpha }_i={\omega }_i{w}_i^{-1}d{u}_i+{\omega }_i{w}_i^{-1}d{v}_i-d{\pi }_i+L_{v_i}$及$d{\beta }_i={\omega }_i{w}_i^{-1}d{u}_i+{\omega }_i{w}_i^{-1}d{v}_i+d{\pi }_i+L_{u_i}.$

本发明实施例的基于最大指数绝对值目标函数的抗差状态估计方法(Maximum  exponential absolute value state estimation，MEAV)在估计过程中可有效抑制包括一致性不良 数据在内的多个不良数据，显示了良好的抗差性，并具有很高的计算效率，非常适宜于实际 工程应用。

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

附图说明

本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和 容易理解，其中：

图1是本发明的基于最大指数绝对值目标函数的抗差状态估计方法的流程示意图；

图2是IEEE 300系统节点电压幅值最大估计误差曲线；

图3是IEEE 300系统节点电压相角最大估计误差曲线：和

图4是IEEE 300系统四种状态估计器的合格率比较。

具体实施方式

下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或 类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的 实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。

如图1所示，本发明实施例的基于最大指数绝对值目标函数的抗差状态估计方法包括下 列步骤：

步骤A：提供基于最大指数绝对值目标函数的抗差状态估计(Maximum exponential  absolute value state estimation，MEAV)基本模型。

具体地，本发明提出的MEAV的基本模型如下所示

$\underset{x}{\max }J\left(x\right)=\underset{i=1}{\overset{m}{\Sigma }}{w}_i\exp (-\frac{\left|r_i\right|}{w_i})---\left(1\right)$

s.t.g(x)=0                                       (2)

r=z-h(x)                               (3)

式中：z∈R^m为量测矢量，常包括节点注入有功和无功、支路有功和无功以及节点电压 幅值量测等；x∈Rⁿ为包括节点电压幅值和相角的状态矢量(平衡节点相角除外)；h:Rⁿ→R^m为 由状态矢量到量测矢量的非线性映射；r_i是残差矢量r的第i个元素；g(x):Rⁿ→R^c为零注入 功率等式约束；w_i为第i个量测量的权重，同是为窗宽参数。

步骤B：对基于最大指数绝对值目标函数的抗差状态估计基本模型引进辅助变量，变换 得到基于最大指数绝对值目标函数的抗差状态估计等价模型。

具体地，MEAV基本模型的目标函数虽然处处连续，但是在0处不可导，因而直接求解 比较困难。可将模型(1)~(3)转化为一个处处可导的等价模型。

引进变量ξ∈R^m，使其满足

|r_i|≤ξ_i  i＝1,…,m                   (4)

由式(1)和(4)可得，最大等价为最大。

引进非负松弛变量l,k∈R^m，将不等式(4)转化为两个等式约束为

r_i-l_i=-ξ_i i=1,…,m                    (5)

r_i+k_i＝ξ_i i=1，…,m                   (6)

引进非负松弛变量u,v∈R^m，使其满足

u_i＝l_i/2  i=1,…,m                     (7)

v_i＝k_i/2  i=1,…,m                     (8)

由式(5)~(8)，可得

r_i＝u_i-v_i i=1,…,m                    (9)

ξ_i=u_i+v_i i=1,…,m                    (10)

将式(9)带入式(3)，可得等价量测约束条件为

z-h(x)-u+v=0                          (11)

则式(1)~(3)给出的MEAV基本模型的等价模型为

$\underset{x}{\max }J\left(x\right)=\underset{i=1}{\overset{m}{\Sigma }}{w}_i\exp (-\frac{u_i+v_i}{w_i})---\left(12\right)$

s.t.g(x)=0                            (13)

z-h(x)-u+v=0                          (14)

u,v ≥0                               (15)

模型(12)~(15)即为MEAV等价模型，该模型处处连续可导，可用基于梯度的方法来求解。 步骤C：利用原-对偶内点算法，对所述基于最大指数绝对值目标函数的抗差状态估计等 价模型求解。

（1）MEAV等价模型的求解方法

注意到MEAV的等价模型(12)~(15)是一个含有等式约束和不等式约束的最优化问题，适 宜用原-对偶内点算法进行求解。为使本领域技术人员更好地理解本发明，首先给出详细的 推导过程如下：

引入拉格朗日函数

$L\equiv \underset{i=1}{\overset{m}{\Sigma }}\left[w_i\exp (-\frac{u_i+v_i}{w_i})\right]-{\lambda }^{T}g\left(x\right)-{\pi }^T(z-h\left(x\right)-u+v)-{\alpha }^{T}v-{\beta }^{T}u---\left(16\right)$

式中：λ∈R^c及π,α,β∈R^m为拉格朗日乘子矢量。

为取得最优值，根据KKT条件，可得

$L_x\equiv \partial L/\partial x=G^T\lambda -H^T\pi =0---\left(17\right)$

$L_{\lambda }\equiv \partial L/\partial \lambda =g\left(x\right)=0---\left(18\right)$

$L_{\pi }\equiv \partial L/\partial \pi =z-h\left(x\right)-u+v=0---\left(19\right)$

$L_{v_i}\equiv \partial L/\partial v_i=-{\omega }_i-{\pi }_i-{\alpha }_i=0---\left(20\right)$

$L_{u_i}\equiv \partial L/\partial u_i=-{\omega }_i+{\pi }_i-{\beta }_i=0---\left(21\right)$

$L_{{\alpha }_i}\equiv \partial L/\partial {\alpha }_i={\alpha }_i{v}_i=0---\left(22\right)$

$L_{{\beta }_i}\equiv \partial L/\partial {\beta }_i={\beta }_i{u}_i=0---\left(23\right)$

式中：ω_i=exp(-(u_i+v_i)/w_i)，$H=\partial h\left(x\right)/\partial x,$$G=\partial g\left(x\right)/\partial x.$

为有效解决以上问题，现代内点法引入扰动参数μ>0对式(22)、(23)进行松弛，从而得

$L_{{\alpha }_i}^{\mu }\equiv {\alpha }_i{v}_i-\mu =0---\left(24\right)$

$L_{{\beta }_i}^{\mu }\equiv {\beta }_i{u}_i-\mu =0---\left(25\right)$

以上方程由牛顿法求解可得

$[{▿}^{2}g\left(x\right)\lambda -{▿}^{2}h\left(x\right)\pi ]\mathrm{dx}+G^T\mathrm{d\lambda }-H^T\mathrm{d\pi }=-L_x---\left(26\right)$

Gdx=-L_λ                              (27)

-Hdx-du+dv=-L_π                       (28)

${\omega }_i{w}_i^{-1}d{u}_i+{\omega }_i{w}_i^{-1}d{v}_i-d{\pi }_i-d{\alpha }_i=-L_{v_i}---\left(29\right)$

${\omega }_i{w}_i^{-1}d{u}_i+{\omega }_i{w}_i^{-1}d{v}_i+d{\pi }_i-d{\beta }_i=-L_{u_i}---\left(30\right)$

${\alpha }_{i}d{v}_i+v_{i}d{\alpha }_i=-L_{{\alpha }_i}^{\mu }---\left(31\right)$

${\beta }_{i}d{u}_i+u_{i}d{\beta }_i=-L_{{\beta }_i}^{\mu }---\left(32\right)$

由式(29)和(30)，可得

$d{\alpha }_i={\omega }_i{w}_i^{-1}d{u}_i+{\omega }_i{w}_i^{-1}d{v}_i-d{\pi }_i+L_{v_i}---\left(33\right)$

$d{\beta }_i={\omega }_i{w}_i^{-1}d{u}_i+{\omega }_i{w}_i^{-1}d{v}_i+d{\pi }_i+L_{u_i}---\left(34\right)$

将式(33)、(34)带入(31)、(32)可得

$[v_i{\omega }_i{w}_i^{-1}+{\alpha }_i]d{v}_i+v_i{\omega }_i{w}_i^{-1}d{u}_i=v_{i}d{\pi }_i-v_i{L}_{v_i}-L_{{\alpha }_i}^{\mu }---\left(35\right)$

$u_i{\omega }_i{w}_i^{-1}d{v}_i+[u_i{\omega }_i{w}_i^{-1}+{\beta }_i]d{u}_i=-u_{i}d{\pi }_i-u_i{L}_{u_i}-L_{{\beta }_i}^{\mu }---\left(36\right)$

令$\left(\begin{array}{cc}{a}_i& b_i\\ c_i& d_i\end{array}\right)={\left(\begin{array}{cc}{v}_i{\omega }_i{w}_i^{-1}+{\alpha }_i& v_i{\omega }_i{w}_i^{-1}\\ u_i{\omega }_i{w}_i^{-1}& u_i{\omega }_i{w}_i^{-1}+{\beta }_i\end{array}\right)}^{-1},$由式(35)和(36)可得

dv_i＝k_1idπ_i+AA_i                          (37)

du_i=k_2idπ_i+BB_i                           (38)

式中：k_1i=a_iv_i-b_iu_i,k_2i=c_iv_i-d_iu_i， ${\mathrm{AA}}_i=-a_i(v_i{L}_{v_i}+L_{{\alpha }_i}^{\mu })-b_i(u_i{L}_{u_i}+L_{{\beta }_i}^{\mu }),{\mathrm{BB}}_i=-c_i(v_i{L}_{v_i}+L_{{\alpha }_i}^{\mu })-d_i(u_i{L}_{u_i}+L_{{\beta }_i}^{\mu }).$

将式(37)、(38)带入(28)可得

Hdx+Qdπ＝γ                              (39)

式中：Q为R^m×m的对角阵，其对角元素为Q_ii=-k_1i+k_2i；γ=z-h(x)-u+v+AA-BB， AA,BB ∈R^m，AA_i,BB_i与式(37)、(38)中同。

根据式(39)、(26)及(27)，可得修正方程为

$\left(\begin{array}{ccc}{▿}^{2}g\left(x\right)\lambda -{▿}^{2}h\left(x\right)\pi & G^T& -H^T\\ H& 0& Q\\ G& 0& 0\end{array}\right)\left(\begin{array}{c}\mathrm{dx}\\ \mathrm{d\lambda }\\ \mathrm{d\pi }\end{array}\right)=\left(\begin{array}{c}-L_x\\ \gamma \\ -L_{\lambda }\end{array}\right)---\left(40\right)$

求解式(40)可得dx,dλ和dπ；由式(37)、(38)可得dv和du；将所得结果带入式(33)、(34) 可得dα和dβ，则迭代即可持续进行。

（2）MEAV等价模型的求解步骤

在介绍MEAV等价模型的求解推导过程之后，发明人将求解步骤归纳如下：

步骤C1:进行初始化，令x为平启动状态变量；选择λ⁽⁰⁾=π⁽⁰⁾=0及u⁽⁰⁾,v⁽⁰⁾,α⁽⁰⁾,β⁽⁰⁾>0； 令中心参数ρ∈(0,1)及确定收敛判据值，以及置迭代计数器为零。

具体地，令x⁽⁰⁾∈Rⁿ代表由所有节点电压幅值和相角组成的的平启动状态变量(参考节点 相角除外)；选择λ⁽⁰⁾=π⁽⁰⁾=0及u⁽⁰⁾,v⁽⁰⁾,α⁽⁰⁾,β⁽⁰⁾>0，其中λ∈R^c及π,α,β∈R^m为拉格朗日乘子矢 量，m为量测量的个数，c为零注入功率约束的个数；令中心参数ρ∈(0,1)及收敛判据 ε=10^-3，置迭代计数器k=0。

步骤C2：计算对偶间隙Gap=α^Tv+β^Tu，判断是否收敛。具体地，若Gap<ε，则认为收 敛，可直接进入步骤C7；否则为不收敛，进入步骤C3。

步骤C3：求解修正方程，以完成对原变量和对偶变量的修正，得到[dx^T dλ^T dπ^T]^T， dv,du，dα和dβ。

具体地，首先计算扰动参数μ=ρ·Gap/2m，然后求解式(40)得[dx^T dλ^T dπ^T]^T;求解式 (37)、(38)得dv,du;求解(33)、(34)得dα,dβ。

步骤C4：计算原问题和对偶问题的修正步长θ_P和θ_D，其中： ${\theta }_P=0.9995\min \{\underset{i}{\min }(-\frac{v_i}{d{v}_i}:d{v}_i<0;-\frac{u_i}{d{u}_i}:d{u}_i<0),1\},$${\theta }_D=0.9995\min \{\underset{i}{\min }(-\frac{{\alpha }_i}{d{\alpha }_i}:d{\alpha }_i<0;-\frac{{\beta }_i}{d{\beta }_i}:d{\beta }_i<0),1\}$

步骤C5：分别修正原问题和对偶问题的变量为：

$\left(\begin{array}{c}{x}^{(k+1)}\\ v^{(k+1)}\\ u^{(k+1)}\end{array}\right)=\left(\begin{array}{c}{x}^{\left(k\right)}\\ v^{\left(k\right)}\\ u^{\left(k\right)}\end{array}\right)+{\theta }_P\left(\begin{array}{c}\mathrm{dx}\\ \mathrm{dv}\\ \mathrm{du}\end{array}\right),\left(\begin{array}{c}{\lambda }^{(k+1)}\\ {\pi }^{(k+1)}\\ {\alpha }^{(k+1)}\\ {\beta }^{(k+1)}\end{array}\right)=\left(\begin{array}{c}{\lambda }^{\left(k\right)}\\ {\pi }^{\left(k\right)}\\ {\alpha }^{\left(k\right)}\\ {\beta }^{\left(k\right)}\end{array}\right)+{\theta }_D\left(\begin{array}{c}\mathrm{d\lambda }\\ \mathrm{d\pi }\\ \mathrm{d\alpha }\\ \mathrm{d\beta }\end{array}\right);$

步骤C6：令迭代计数器k=k+1，进入步骤C2；以及

步骤C7：输出最优解，结束。

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

设定利用IEEE标准系统检验基于原-对偶内点算法的MEAV的性能。试验采用全 量测，量测值通过在潮流计算的结果上叠加白噪声(均值为0，标准差为τ)来获得。对于 电压量测，取τ_V=0.005p.u；对于功率量测，取τ_PQ＝1MW/MVar。测试环境为PC机，CPU 为Intel(R)Core(TM)i3M370、主频为2.40GHz、内存2.00GB。

1．抗差性能的比较

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

（1）IEEE-14系统

在IEEE-14系统上设置4个一致性不良数据(P_1-2、Q_1-2、P₁、Q₁)。所设置的不良量 测值以及量测量的正确值如表1所示。

表1MEAV对IEEE14系统一致性不良数据的辨识

Table1 Estimation of conforming bad data for the IEEE 14-bus system by MEAV

作为对比，首先用广为应用的WLS进行估计，并用LNR进行不良数据的辨识(简 记为WLS+LNR)。首次辨识的结果为：10个量测量的标准化残差大于门槛值(3.0)，这 10个量测量被认为是可疑数据；其中标准化残差最大的量测量为P_2-1，删去该量测后重 新运行WLS；此时发现P₂的标准化残差最大。以上过程循环4次，4个良好的量测量 被LNR误认为是可疑数据而被删去，但真正的不良数据仍然存在。可见，WLS+LNR 不能辨识一致性不良数据。

应用MEAV方法的估计结果如表1所示。可以发现，即使量测量中存在一致性不 良数据，MEAV的估计值与真值也可很好地吻合。在IEEE其他系统的多次试验也表明 MEAV在估计的过程中可以自动抑制不良数据，具有良好的抗差性。

（2）IEEE-300系统

为了测试MEAV在规模更大系统上的估计性能，本节基于IEEE 300节点系统，分 别对不良数据比例为0%~10%共11种情况进行试验，每一种比例下均包括50次仿真试 验。不良数据是在正确量测值上迭代50%的误差而产生。图2、3给出了在不同比例的 不良数据下，WLS+LNR及MEAV得到的节点电压幅值和相角估计误差绝对值的最大值 的平均值(50次试验)，其中，|ΔU_max|和|Δθ_max|分别表示节点电压幅值和相角估计误差绝对 值的最大值，η表示不良数据的比例。

由图2、3可见，当没有不良数据时(η=0)，MEAV与WLS+LNR的估计误差均较 小。然而，随着不良数据比例的增大，WLS+LNR的估计误差迅速增大，而MEAV估计 的幅值误差则基本不变，而相角误差只略有增加(即使不良数据比例增大为10%时亦如 此)，体现了良好的抗差性能。

进一步，我们对四种状态估计器（WLS+LNR、WLAV、MNMR及MEAV）在不同 不良数据比例时的合格率进行了测试。50次试验的平均值如图4所示，其中η和Ω分别 表示不良数据的比例及状态估计的合格率。由图4可见，随着不良数据比例的上升， WLS+LNR的合格率迅速下降；而其余三种抗差估计器的合格率下降相对较少。在同一 不良数据比例时，WLS+LNR的合格率最低，而MEAV的合格率最高，MNMR的合格 率次之。由此可见，在以上四种状态估计器中，MEAV的抗差性最好。

2．计算效率的比较

发明人为了进行效率比较，在正常量测条件下分别对四种状态估计器WLS、WLAV、 MNMR以及MEAV进行了测试，其中后三种属于抗差状态估计器。在试验中，WLS采 用牛顿法求解，其他三种状态估计采用内点法求解；且MNMR采用两阶段法，即第一 阶段进行WLS估计，第二阶段将WLS的估计值作为MNMR估计的初值进行计算。

共进行50次仿真试验，状态估计收敛时的迭代次数以及平均计算耗时如表2所示。 由表2可见，在这四种状态估计器中，WLS的计算效率最高；而在后三种抗差状态估 计器中，MEAV的计算效率最高；而且随着系统规模的增大，MEAV的迭代次数以及 计算耗时增长的很缓慢，因而MEAV适用于实际的大规模系统的估计。

表2四种状态估计器的迭代次数以及计算耗时

Table2 Iterations and CPU time of the four estimators

综上所述，本发明提出的MEAV在估计过程中可有效抑制包括一致性不良数据在内的 多个不良数据，显示了良好的抗差性，并具有很高的计算效率，非常适宜于实际工程应用。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、 或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含 于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的 是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或 多个实施例或示例中以合适的方式结合。

尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的， 不能理解为对本发明的限制，本领域的普通技术人员在不脱离本发明的原理和宗旨的情况下 在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。