结构风险最小化的加权最小二乘电力系统状态估计方法

摘要

本发明公开了一种基于结构风险最小化的加权最小二乘电力系统状态估计方法，该方法针对电力系统状态估计的量测数目有限的特点，从统计学习理论出发，提出基于结构风险最小化的加权最小二乘估计模型，可在最小化残差的范数的同时最小化状态变量的置信区间；并给出了该方法的详细求解过程。该方法符合统计学习理论中的结构风险最小化思想，在有限量测条件下可得到更接近于状态变量真值的估计结果，具有良好的工程应用前景。

说明书

技术领域

本发明属于电力系统调度自动化领域，具体涉及一种基于结构风险最小化的加权最 小二乘(Structural Risk Minimization Based Weighted Least Squares,SWLS)电力系统状态估计 方法。

背景技术

电力系统状态估计是能量管理系统(Energy Management System，EMS)的基础和核心 组成部分。高性能的状态估计(State Estimator，SE)可为EMS提供准确可靠的实时运行数 据，是EMS各项高级应用正常运行的保证。现在几乎每一个大型的电网调度中心都装设 了状态估计软件，SE已成为保证电网安全、可靠运行必不可少的环节。

目前，在国内外应用最为广泛的状态估计方法是加权最小二乘法(Weighted Least  Squares,WLS)。WLS模型简洁，求解容易，收敛性能好。WLS的理论基础是传统统计学 的经验风险最小化思想(Empirical Risk Minimization,ERM)，即当样本(量测量)数目趋近于 无穷大时，状态变量估计值以概率逼近于状态变量真值。但是在实际工程中，量测量的数 目毕竟是有限的，此时基于ERM思想的WLS方法从理论上并不能保证估计性能。

发明内容

本发明旨在至少解决现有技术中存在的技术问题。统计学习理论是有限样本时的最佳 学习理论，统计学习理论认为有限样本时更为合理的建模方法是结构风险最小化思想 (Structural Risk Minimization，SRM)。将SRM思想运用于解决状态估计问题在理论更具合 理性。基于SRM思想进行状态估计建模就是要在最小化残差表示的某种范数的同时，最小 化状态变量的置信区间。

为此，本发明的目的在于提出一种有限量测数据条件下更为合理的电力系统状态估计 方法，即基于结构风险最小的电力系统状态估计方法。

为了实现上述目的，根据本发明的实施例的结构风险最小化的加权最小二乘电力系统 状态估计方法，可以包括以下步骤：A.在电力系统中形成网络模型，并提出基于结构风险 最小化的加权最小二乘估计模型为 Min S(x)=k_p[z-h(x)]^TR^-1[z-h(x)]+tr([H^T(x)R^-1H(x)]^-1)，其中：S(x)代表状态估计学 习模型的结构风险，x代表状态变量，k_p>0为常数参数，z∈R^m为量测矢量，包括节点 电压幅值量测、支路有功和无功量测、节点注入有功和无功量测，h为量测表达式，tr(·) 代表矩阵的迹，为雅可比矩阵；B.计算所述雅可比矩阵H(x^(k))及增益矩 阵G(x^(k))=H^T(x^(k))R^-1H(x^(k))，其中k为迭代计数；C.对所述增益矩阵G(x^(k))进行因子分 解，并计算矩阵A=[A₁,A₂,…,A_n]^T∈Rⁿ，A_i=tr(G^-1(x)B_iG^-1(x))∈R(i=1,2,…,n)， 所述矩阵A为求解所述基于结构风险最小化的加权最小二乘估计模 型过程中的中间运算物理量；D.根据所述矩阵A，计算H^T(x^(k))R^-1[z-h(x^(k))]+A/k_p；E. 根据步骤C得到的所述增益矩阵G(x^(k))的因子式，计算Δx^(k)；F.判断Δx^(k)是否收敛，如果 收敛则输出状态变量估计值结果，如果不收敛则更新x^(k+1)=x^(k)+Δx^(k)，k=k+1然后转至 步骤B。

根据本发明实施例的结构风险最小化的加权最小二乘电力系统状态估计方法，符合统 计学习理论中的结构风险最小化思想，在有限量测条件下可得到更接近于状态变量真值的 估计结果，具有良好的工程应用前景。

另外，根据本发明实施例的结构风险最小化的加权最小二乘电力系统状态估计方法还 可以具有如下附加技术特征：

在本发明的一个实施例中，所述步骤A具体包括：将网络中所有的线路和变压器等效 为π型支路ij，记y_s=1/(r_ij+jx_ij)=g_s+jb_s为π型支路ij的串联电纳，r_ij+jx_ij为π型支路ij的串联 阻抗值，b_c为π型支路ij的接地电纳，其中，若π型支路ij为变压器支路，则b_c=0且k为理 想变压器的变比，若π型支路ij为普通线路，则k=1，并联的多条支路等效为一条支路；在 等效后的电路中，记g_ij=g_s/k，b_ij=b_s/k，g_si=(1-k)g_s/k²，b_si=(1-k)b_s/k²+b_c/2，g_sj=(k-1)g_s/k， b_sj=(k-1)b_s/k+b_c/2；计算节点导纳矩阵Y=G+jB,G和B分别为节点导纳矩阵的实部和虚部。

在本发明的一个实施例中，所述雅可比矩阵元素包括：电压幅值量测对应的雅可比矩 阵元素、支路功率量测对应的雅可比矩阵元素以及注入功率量测对应的雅可比矩阵元素。

在本发明的一个实施例中，所述步骤C具体包括以下步骤：C1.采用数值方法求G(x^(k)) 的逆矩阵G^-1(x^(k))；C2.计算矩阵元素，其中 ；C3.计算A_i=tr(G^-1(x)B_iG^-1(x))(i=1,2,…,n)以得到 A=[A₁,A₂,…,A_n]^T(i=1,2,…,n)。

在本发明的一个实施例中，所述步骤D中的公式H^T(x^(k))R^-1[z-h(x^(k))]+A/k_p中，z 为量测矢量，h(x)的表达式如下：电压幅值量测对应h(x)：第i个电压幅值量测对应的 h(x)元素为v_i=v_i；支路功率量测对应h(x)：支路ij的有功和无功量测对应的h(x)元素 为$P_{\mathrm{ij}}=v_i^2(g_{\mathrm{si}}+g_{\mathrm{ij}})-v_i{v}_j{g}_{\mathrm{ij}}\cos {\theta }_{\mathrm{ij}}-v_i{v}_j{b}_{\mathrm{ij}}\sin {\theta }_{\mathrm{ij}},Q_{\mathrm{ij}}=-v_i^2(b_{\mathrm{si}}+b_{\mathrm{ij}})+v_i{v}_j{b}_{\mathrm{ij}}\cos {\theta }_{\mathrm{ij}}-v_i{v}_j{g}_{\mathrm{ij}}\sin {\theta }_{\mathrm{ij}},$其中，P_ij和 Q_ij分别为节点i流向节点j的支路有功和无功；注入功率量测对应h(x)：节点i的注入有 功和注入无功量测对应的h(x)元素为

其中，P_i和Q_i分别为节点i的注入有功和注入无功，G_ij+jB_ij为 节点导纳矩阵中的对应元素。

在本发明的一个实施例中，根据步骤C得到的所述增益矩阵G(x^(k))的因子式，求解下 式以得到Δx^(k)：$\left(\begin{array}{c}{x}^{(k+1)}=x^{\left(k\right)}-K^{-1}\left(x^{\left(k\right)}\right)g\left(x^{\left(k\right)}\right)\\ =x^{\left(k\right)}+[k_{p}G\left(x^{\left(k\right)}\right)-\partial A/\partial x{|}_{x=x^{\left(k\right)}}{]}^{-1}g\left(x^{\left(k\right)}\right)\\ \approx x^{\left(k\right)}+[k_{p}G\left(x^{\left(k\right)}\right){]}^{-1}g\left(x^{\left(k\right)}\right)\\ =x^{\left(k\right)}+[G\left(x^{\left(k\right)}\right){]}^{-1}\left(H^T\left(x^k\right)R^{-1}\right[z-h\left(x^{\left(k\right)}\right)]+A/k_p)\end{array}\right).$

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

附图说明

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

图1为本发明实施例的基于结构风险最小化的加权最小二乘电力系统状态估计方法的 流程示意图；

图2为π型支路的示意图；

图3为π型支路等值电路的示意图；

图4为典型的两节点电力系统状态估计算例的示意图；

图5为1000次试验由EWLS和SWLS得到的解的散布图；

图6为10000次试验由EWLS和SWLS得到的解的散布图。

具体实施方式

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

本发明涉及一种基于SRM准则的WLS模型(SWLS)，该模型符合统计学习理论中的 SRM思想，在有限量测条件下可得到更接近于状态变量真值的估计结果。

以下首先给出传统的基于ERM准则的WLS模型(Empirical Risk Minimization Based  Weighted Least Squares,EWLS)，即传统的WLS状态估计模型；然后给出SWLS的模型；

最后详细推导了SWLS的求解方法，并给出了求解步骤。

（1）模型

现有技术中基于ERM准则的EWLS模型如下所示

Min J(x)=[z-h(x)]^TR^-1[z-h(x)]          (1) 式中：z∈R^m为量测矢量，一般包括支路功率量测、节点注入功率量测和节点电压幅值量 测等，PMU量测中还包括相角量测，m为量测量的总个数；x∈Rⁿ为包括所有节点电压幅 值和相角（参考节点相角除外）的状态变量，n=2N-1，N为网络中节点的总数目； h:Rⁿ→R^m为状态变量到量测矢量的非线性映射，即量测表达式。

EWLS模型(1)的迭代求解公式如下式所示

G(x^(k))Δx^(k+1)=H^T(x^(k))R^-1[z-h(x^(k))]          (2)

式中：，G(x^(k))=H^T(x^(k))R^-1H(x^(k))为增益矩阵（信息矩阵）； Δx^(k+1)=x^(k+1)-x^(k)。

EWLS模型在收敛时，可得状态估计的误差方差矩阵为

$(x-\hat{x}){(x-\hat{x})}^T{G}^{-1}\left(x\right)---\left(3\right)$

式中：x为状态变量的真值；为EWLS得到的状态变量的估计值，即式（2）收敛时得 到的状态变量；G(x)=H^T(x)R^-1H(x)。

状态估计误差方程阵[H^T(x)R^-1H(x)]^-1的对角元表示状态变量的误差方差。

基于SRM准则的状态估计一般模型如式(4)所示

$\mathrm{MinS}\left(x\right)=\frac{1}{m}\underset{i=1}{\overset{l}{\Sigma }}L\left(\right|z_i-h_i\left(x\right)\left|\right)+\mathrm{k\Phi }(m/h)---\left(4\right)$

式中：S(x)代表状态统计学习模型的结构风险；x代表状态变量；Φ(m/h)代表状态变量 的置信区间；k用来衡量置信区间与经验风险的比例；z_i和h_i分别为z和h的第i个分量。

将模型(5)应用于WLS时，基于此，可得如下模型

Min S(x)=k_p[z-h(x)]^TR^-1[z-h(x)]+tr([H^T(x)R^-1H(x)]^-1)     (5)

式中：tr(·)代表矩阵的迹；k_p>0为某一常数参数。

模型（5）即为基于SRM准则的WLS模型，即SWLS，S(x)中的第一项用残差的加 权平方和代表经验风险，与常规WLS(EWLS)相同；第二项用所有变量的误差方差之和代 表状态变量的置信区间。k_p代表经验风险占结构风险的比例，k_p随着量测量数目的增大而 增大；k_p趋近于无穷大时，模型（5）等价于基于EWLS；k_p趋近于0时，模型（5）的意 义是仅仅最小化状态变量的置信区间。

（2）SWLS模型的求解方法

为了求得最优值，模型(5)必须满足一阶最优条件，其中涉及对矩阵的逆进行求导的运 算。

令A(x)∈R^n×n为可逆矩阵，x为标量，则有则模型（5）中的第二 项对状态变量x（x∈Rⁿ）的第i(i=1,2,…,n)个分量x_i的导数为

$\frac{d\left(\mathrm{tr}\left(\right[H^T\left(x\right)R^{-1}H\left(x\right){]}^{-1})\right)}{d{x}_i}=-\mathrm{tr}\left(G^{-1}\left(x\right)\right[B_i^T+B_i\left]G^{-1}\left(x\right)\right)---\left(6\right)$

式中：$B_i=H^T\left(x\right)R^{-1}\frac{\mathrm{dH}\left(x\right)}{d{x}_i}.$

根据基本公式则式（6）等号右边的两项满足

$\mathrm{tr}\left(G^{-1}\left(x\right)B_i^T{G}^{-1}\left(x\right)\right)=\mathrm{tr}\left({\left(G^{-1}\left(x\right)\right)}^2{B}_i^T\right)---\left(7\right)$

tr(G^-1(x)B_iG^-1(x))=tr((G^-1(x))²B_i)         (8)

由于G^-1(x)=[G^-1(x)]^T，因此(G^-1(x))²=[(G^-1(x))²]^T，结合式（7）和式（8），可得

$\mathrm{tr}\left(G^{-1}\left(x\right)B_i^T{G}^{-1}\left(x\right)\right)=\mathrm{tr}\left({G^{-1}\left(x\right)B}_i{G}^{-1}\left(x\right)\right)---\left(9\right)$

将式（9）代入式（6），可得

$\frac{d\left(\mathrm{tr}\left(\right[H^T\left(x\right)R^{-1}H\left(x\right){]}^{-1})\right)}{d{x}_i}=2\mathrm{tr}\left(G^{-1}\left(x\right)B_i{G}^{-1}\left(x\right)\right)---\left(10\right)$

令A_i=tr(G^-1(x)B_iG^-1(x))∈R(i=1,2,…,n)，A=[A₁,A₂,…,A_n]^T∈Rⁿ，则有

$\frac{d\left(\mathrm{tr}\left(\right[H^T\left(x\right)R^{-1}H\left(x\right){]}^{-1})\right)}{\mathrm{dx}}=-2A---\left(11\right)$

至此，对模型（5）中等号右边第2项的求导完成。为取得最优值，模型（5）须满足 下面的一阶最优条件

$\frac{\partial S\left(x\right)}{\mathrm{dx}}=-2k_p{H}^T\left(x\right)R^{-1}[z-h\left(x\right)]-2A=0---\left(12\right)$

式（12）可改写为下面的等价形式

k_pH^T(x)R^-1[z-h(x)]+A=0                   (13)

令g(x)=k_pH^T(x)R^-1[z-h(x)]+A=0，将g(x)进行泰勒级数展开，可得

g(x)=g(x^(k))+K(x^(k))(x-x^(k))+…=0              (14)

式中：$K\left(x^{\left(k\right)}\right)=\partial g\left(x^{\left(k\right)}\right)/\partial x=-k_{p}G\left(x^{\left(k\right)}\right)+\partial A/\partial x{|}_{x=x^{\left(k\right)}}.$

忽略式（14）中的高阶项，可得如下迭代求解公式

$\left(\begin{array}{c}{x}^{(k+1)}=x^{\left(k\right)}-K^{-1}\left(x^{\left(k\right)}\right)g\left(x^{\left(k\right)}\right)\\ =x^{\left(k\right)}+[k_{p}G\left(x^{\left(k\right)}\right)-\partial A/\partial x{|}_{x=x^{\left(k\right)}}{]}^{-1}g\left(x^{\left(k\right)}\right)\\ \approx x^{\left(k\right)}+[k_{p}G\left(x^{\left(k\right)}\right){]}^{-1}g\left(x^{\left(k\right)}\right)\\ =x^{\left(k\right)}+[G\left(x^{\left(k\right)}\right){]}^{-1}\left(H^T\left(x^k\right)R^{-1}\right[z-h\left(x^{\left(k\right)}\right)]+A/k_p)\end{array}\right)---\left(15\right)$

在迭代公式（15）中忽略了这一项，其原因在于一方面的求解比较困难； 另一方面是忽略后不影响收敛时的状态变量估计值。

式（15）为SWLS模型（5）的迭代求解公式。注意到当样本数目趋近于无穷大时，k_p等于无穷大，其迭代等同于常规WLS（EWLS）的迭代求解公式。

（3）模型的求解步骤

图1是本发明的方法的流程示意图。如图1所示，可以包括以下步骤：

A.在电力系统中形成网络模型，并提出基于结构风险最小化的加权最小二乘估计模型 为Min S(x)=k_p[z-h(x)]^TR^-1[z-h(x)]+tr([H^T(x)R^-1H(x)]^-1)，其中式中：S(x)代表状态 估计学习模型的结构风险，x代表状态变量，k_p>0为常数参数，z∈R^m为量测矢量，包括 节点电压幅值量测、支路有功和无功量测、节点注入有功和无功量测，h为量测表达式，tr(·) 代表矩阵的迹，为雅可比矩阵。

B.计算雅可比矩阵H(x^(k))及增益矩阵G(x^(k))=H^T(x^(k))R^-1H(x^(k))，其中k为迭代计 数。

C.对增益矩阵G(x^(k))进行因子分解，并计算矩阵A=[A₁,A₂,…,A_n]^T∈Rⁿ， A_i=tr(G^-1(x)B_iG^-1(x))∈R(i=1,2,…,n)，矩阵A为求解基于结构 风险最小化的加权最小二乘估计模型过程中的中间运算物理量。

D.根据矩阵A，计算H^T(x^(k))R^-1[z-h(x^(k))]+A/k_p。

E.根据步骤C得到的增益矩阵G(x^(k))的因子式，计算Δx^(k)。

F.判断Δx^(k)是否收敛，如果收敛则输出状态变量估计值结果，如果不收敛则更新 x^(k+1)=x^(k)+Δx^(k)，k=k+1然后转至步骤B。

根据本发明实施例的基于结构风险最小化的加权最小二乘电力系统状态估计方法，符 合统计学习理论中的结构风险最小化思想，在有限量测条件下可得到更接近于状态变量真 值的估计结果，具有良好的工程应用前景。

为使本领域技术人员更好地理解本发明，对SWLS模型（5）的求解步骤阐述如下：

Step1:电力系统中的状态变量包括除参考节点外的所有节点的电压幅值和相角，记为x^(k)， 即其中N为网络中节点的总数目，节点1为参考节点； 置迭代计数器k=0；令x^(k)为平启动状态变量，即假设所有节点的电压幅值为1，所有节点 的相角为0弧度；

Step2:计算雅可比矩阵H(x^(k))及增益矩阵G(x^(k))；

Step2包括：A形成网络模型；B计算雅可比矩阵及增益矩阵；

其中，A包括：网络中任一支路(线路或变压器支路)的π型支路如图2所示，图3中是 图2的等效电路。在图2和图3中记y_s=1/(r_ij+jx_ij)=g_s+jb_s为π型支路ij的串联电纳，r_ij+jx_ij为π型支路ij的串联阻抗值，b_c为π型支路ij的接地电纳，其中，若π型支路ij为变 压器支路，则b_c=0且k为理想变压器的变比，若π型支路ij为普通线路，则k=1；在图3 中，记g_ij=g_s/k，b_ij=b_s/k，g_si=(1-k)g_s/k²，b_si=(1-k)b_s/k²+b_c/2，g_sj=(k-1)g_s/k， b_sj=(k-1)b_s/k+b_c/2；计算节点导纳矩阵Y=G+jB,G和B分别为节点导纳矩阵的实部和虚 部；

B电力系统SE中的量测量主要包括节点电压幅值量测v_i，支路功率量测P_ij，Q_ij，节 点注入功率量测P_i和Q_i，所有的量测量记为z；量测量的真值记为h(x)，雅可比矩阵定义 为，雅可比矩阵的各元素如下(在第k次迭代中，仅选用x^(k)代替x即可)：

1）电压幅值量测对应的雅可比矩阵元素

对于电压幅值量测，其对应的雅可比矩阵元素为

$\frac{\partial v_i}{\partial v_i}=1,\frac{\partial v_i}{\partial v_j}=0,\frac{\partial v_i}{\partial {\theta }_i}=0,\frac{\partial v_i}{\partial {\theta }_j}=0.$

2）支路功率量测对应的雅可比矩阵元素

对于支路功率量测，其对应的雅可比矩阵元素为

$\left(\begin{array}{c}\frac{\partial P_{\mathrm{ij}}}{\partial v_i}=2v_i(g_{\mathrm{si}}+g_{\mathrm{ij}})-v_j{g}_{\mathrm{ij}}\cos {\theta }_{\mathrm{ij}}-v_j{b}_{\mathrm{ij}}\sin {\theta }_{\mathrm{ij}},\frac{\partial P_{\mathrm{ij}}}{\partial v_j}=-v_i{g}_{\mathrm{ij}}\cos {\theta }_{\mathrm{ij}}-v_i{b}_{\mathrm{ij}}\sin {\theta }_{\mathrm{ij}},\\ \frac{\partial P_{\mathrm{ij}}}{\partial {\theta }_i}{=v}_i{v}_j{g}_{\mathrm{ij}}\sin {\theta }_{\mathrm{ij}}-v_i{v}_j{b}_{\mathrm{ij}}\cos {\theta }_{\mathrm{ij}},\frac{\partial P_{\mathrm{ij}}}{\partial {\theta }_j}=-v_i{v}_j{g}_{\mathrm{ij}}\sin {\theta }_{\mathrm{ij}}+v_i{v}_j{b}_{\mathrm{ij}}\cos {\theta }_{\mathrm{ij}}.\\ \frac{\partial Q_{\mathrm{ij}}}{\partial v_i}=-2v_i(b_{\mathrm{si}}+b_{\mathrm{ij}})-v_j{b}_{\mathrm{ij}}\cos {\theta }_{\mathrm{ij}}-v_j{g}_{\mathrm{ij}}\sin {\theta }_{\mathrm{ij}},\frac{\partial Q_{\mathrm{ij}}}{\partial v_j}=v_i{b}_{\mathrm{ij}}\cos {\theta }_{\mathrm{ij}}-v_i{g}_{\mathrm{ij}}\sin {\theta }_{\mathrm{ij}},\\ \frac{\partial Q_{\mathrm{ij}}}{\partial {\theta }_i}{=-v}_i{v}_j{b}_{\mathrm{ij}}\sin {\theta }_{\mathrm{ij}}-v_i{v}_j{g}_{\mathrm{ij}}\cos {\theta }_{\mathrm{ij}},\frac{\partial Q_{\mathrm{ij}}}{\partial {\theta }_j}=v_i{v}_j{b}_{\mathrm{ij}}\sin {\theta }_{\mathrm{ij}}+v_i{v}_j{g}_{\mathrm{ij}}\cos {\theta }_{\mathrm{ij}}.\end{array}\right)$

3）注入功率量测对应的雅可比矩阵元素

对于注入功率量测，其对应的雅可比矩阵元素为

$\left(\begin{array}{c}\frac{\partial P_i}{\partial v_i}=2v_i{G}_{\mathrm{ii}}+\underset{j=1,j\ne i}{\overset{N}{\Sigma }}{v}_j(G_{\mathrm{ij}}\cos {\theta }_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\theta }_{\mathrm{ij}}),\frac{\partial P_i}{\partial v_j}=v_i(G_{\mathrm{ij}}\cos {\theta }_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\theta }_{\mathrm{ij}}),\\ \frac{\partial P_i}{\partial {\theta }_i}=v_i\underset{j=1,j\ne i}{\overset{N}{\Sigma }}{v}_j(-G_{\mathrm{ij}}\sin {\theta }_{\mathrm{ij}}+B_{\mathrm{ij}}\cos {\theta }_{\mathrm{ij}}),\frac{\partial P_i}{\partial v_j}=v_i{v}_j(G_{\mathrm{ij}}\sin {\theta }_{\mathrm{ij}}-B_{\mathrm{ij}}\cos {\theta }_{\mathrm{ij}}),\\ \frac{\partial Q_i}{\partial v_i}=-2v_i{B}_{\mathrm{ii}}+\underset{j=1,j\ne i}{\overset{N}{\Sigma }}{v}_j(G_{\mathrm{ij}}\sin {\theta }_{\mathrm{ij}}-B_{\mathrm{ij}}\cos {\theta }_{\mathrm{ij}}),\frac{\partial Q_i}{\partial v_j}=v_i(G_{\mathrm{ij}}\sin {\theta }_{\mathrm{ij}}-B_{\mathrm{ij}}\cos {\theta }_{\mathrm{ij}}),\\ \frac{\partial Q_i}{\partial {\theta }_i}=v_i\underset{j=1,j\ne i}{\overset{N}{\Sigma }}{v}_j(G_{\mathrm{ij}}\cos {\theta }_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\theta }_{\mathrm{ij}}),\frac{\partial Q_i}{\partial {\theta }_j}=v_i{v}_j({-G}_{\mathrm{ij}}\cos {\theta }_{\mathrm{ij}}-B_{\mathrm{ij}}\sin {\theta }_{\mathrm{ij}}).\end{array}\right)$

形成H(x^(k))之后，计算G(x^(k))=H^T(x^(k))R^-1H(x^(k))。

Step3:对增益矩阵G(x^(k))进行因子分解，由公式（10）和（11）计算矩阵A；

其中，矩阵A中的各元素的计算步骤如下：

1）采用数值方法求G(x^(k))的逆矩阵G^-1(x^(k))；

2）计算矩阵元素$B_i=H^T\left(x\right)R^{-1}\frac{\mathrm{dH}\left(x\right)}{d{x}_i}(i=\mathrm{1,2},...,n),$其中$\frac{\mathrm{dH}\left(x\right)}{d{x}_i}(i=\mathrm{1,2},...,n)$各元素

的数值如下（相当于海森矩阵元素）：

（1）电压幅值量测对应的雅可比矩阵元素对x的导数

$\frac{{\partial }^2{v}_i}{\partial v_i^2}=0,\frac{{\partial }^2{v}_i}{\partial v_j\partial v_i}=0,\frac{{\partial }^2{v}_i}{\partial {\theta }_i\partial v_i}=0,.\frac{{\partial }^2{v}_i}{\partial {\theta }_j\partial v_i}=0,\frac{{\partial }^2{v}_i}{\partial {\theta }_j\partial v_i}=0,\frac{{\partial }^2{v}_i}{\partial {\theta }_i\partial v_i}=0.$

（2）支路功率量测对应的雅可比矩阵元素对x的导数

3）计算A_i=tr(G^-1(x)B_iG^-1(x))(i=1,2,…,n)，以得到A=[A₁,A₂,…,A_n]^T(i=1,2,…,n)。

以上三步中用到x时需用x(k)代替；

Step4:计算H^T(x^(k))R^-1[z-h(x^(k))]+A/k_p；

其中，z为量测矢量，h(x)的表达式如下(在第k次迭代中，仅选用x^(k)代替x即可)：

1）电压幅值量测对应h(x)

第i个电压幅值量测对应的h(x)元素为

v_i=v_i；

2）支路功率量测对应h(x)

支路ij（i、j分别为支路两端节点）的有功和无功量测对应的h(x)元素为

$P_{\mathrm{ij}}=v_i^2(g_{\mathrm{si}}+g_{\mathrm{ij}})-v_i{v}_j{g}_{\mathrm{ij}}\cos {\theta }_{\mathrm{ij}}-v_i{v}_j{b}_{\mathrm{ij}}\sin {\theta }_{\mathrm{ij}};$

$Q_{\mathrm{ij}}={-v}_i^2(b_{\mathrm{si}}+b_{\mathrm{ij}})-v_i{v}_j{b}_{\mathrm{ij}}\cos {\theta }_{\mathrm{ij}}-v_i{v}_j{g}_{\mathrm{ij}}\sin {\theta }_{\mathrm{ij}};$

其中，P_ij和Q_ij分别为节点i流向节点j的支路有功和无功。

3）注入功率量测对应h(x)

节点i的注入有功和注入无功量测对应的h(x)元素为

$P_i=v_i\underset{j\in N_i}{\Sigma }{v}_j(G_{\mathrm{ij}}\cos {\theta }_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\theta }_{\mathrm{ij}});$

$Q_i=v_i\underset{j\in N_i}{\Sigma }{v}_j(G_{\mathrm{ij}}\sin {\theta }_{\mathrm{ij}}-B_{\mathrm{ij}}\cos {\theta }_{\mathrm{ij}});$

其中，P_i和Q_i分别为节点i的注入有功和注入无功，G_ij+jB_ij为节点导纳矩阵中的对应 元素。

Step5:利用Step3得到的G(x^(k))的因子式求解式（15），得到Δx^(k)；

Step6:判断|Δx^(k)|≤ε是否成立，若成立，则转Step8；否则，转Step7；

Step7:更新状态变量为x^(k+1)=x^(k)+Δx^(k)，k=k+1，转Step2；

Step8:输出状态变量估计值，END。

为使本技术领域人员更好地理解本发明，下面结合图4-6列举实施例做进一步说明。

图4为典型的两节点电力系统状态估计算例的示意图。其中节点1(Bus1)是参考节点， 节点2(Bus2)的电压幅值和相角待估计（其状态变量的真值为v₂=1,θ₂=-π/6≈-0.5236）。 测试中在潮流的基础上叠加高斯分布的随机数来模拟量测量。图4中同时给出了支路参数 (r和x分别表示线路的电阻和电抗)、潮流真实值(P_ij表示从节点i流向节点j潮流的真值) 以及其中一组量测值（括号内的数值）。

在测试时，电压幅值量测、支路有功和注入有功量测以及支路无功和注入无功量测中 相对噪声的标准差分别是0.01p.u.、0.02p.u.和0.03p.u.（相对噪声的意思是噪声幅值与对应 量测量真值的比值）。分别独立进行1000次和10000次试验（测试中k_p=10⁴），由EWLS 和SWLS得到的解的散布图分别如图5和图6所示。

在图5和图6中，v₂表示节点2节点幅值的估计值，θ₂表示节点2电压相角的估计值。 由图3和图4可见，当不含不良数据时，SWLS得到的状态变量的散布区间小于EWLS得 到的散布区间，显示了SWLS的优越性。

根据以上测试结果，可对广为应用的EWLS进行改造以改进其估计性能。即首先利用 EWLS进行估计，并利用估计辨识法进行不良数据的辨识；待不良数据被全部辨识出来时， 可用SWLS再次进行估计，即可得到更为可信的估计结果。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、 或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包 含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须 针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一 个或多个实施例或示例中以合适的方式结合。此外，本领域的技术人员可以将本说明书中 描述的不同实施例或示例进行结合和组合。

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