CN102801162B - 一种两阶段线性加权最小二乘电力***状态估计方法 - Google Patents

Abstract

本发明提出一种两阶段线性加权最小二乘电力***状态估计方法，其特征在于，包括以下步骤：形成网络模型，计算节点导纳矩阵及支路-节点关联矩阵；对量测矢量和状态矢量进行变换；形成精确线性化的量测方程；进行第一阶段的线性加权最小二乘估计，得到变换后的状态矢量的估计值；进行逆变换，并进行第二阶段的线性加权最小二乘估计，得到所有节点的电压幅值和相角的估计值；以及进行不良数据辨识。根据本发明提出的两阶段线性加权最小二乘状态估计方法可以得到更为科学的状态估计结果，而且计算效率更高，具有很好的工程应用前景。

Description

一种两阶段线性加权最小二乘电力***状态估计方法

技术领域

本发明涉及电力***调度自动化领域，特别涉及一种两阶段线性加权最小二乘电力***状态估计方法。

背景技术

电力***状态估计是能量管理***的基础和核心，其作用是对数据采集与监控***(SCADA)提供的实时信息进行滤波，从而得到全网状态变量(电压幅值和相角)的估计值，进而可以得到支路功率、节点注入功率等的估计值。

自从国外学者在1970年提出状态估计的第一个模型以来，国内外学者和工程人员对状态估计进行了大量、深入的研究和实践，现在全世界每一个调度中心几乎都部署了状态估计器，状态估计在电网安全运行中的基础性地位得到了广泛的认可。

到目前为止，应用最为广泛的状态估计方法是非线性加权最小二乘法(WLS)，该法模型简单、计算方便，在没有不良数据时可以得到状态变量的无偏估计值。为了提高非线性WLS的计算效率，人们又提出了快速解耦法状态估计方法。因为非线性WLS不具有鲁棒性，人们又提出了基于非线性WLS的最大标准化残差(LNR)检验方法以及不良数据的估计辨识方法(EI)等来辨识不良数据。

与此同时，各种各样的鲁棒性状态估计方法也不断被提出。鲁棒性状态估计能够在估计的过程中自动抑制不良数据的影响，从而得到正确的状态变量估计值。在各种鲁棒性状态估计方法中，M-估计的应用最为普遍。M-估计一般包括加权最小绝对值估计(WLAV)、非二次准则状态估计(如QL估计和QC估计等)等，最近几年刚被提出的以合格率最大为目标的状态估计以及含指数型目标函数的状态估计也属于M-估计的范畴。鲁棒性状态估计面临的问题：1）全局寻优困难；2）计算效率低于非线性WLS。以上原因影响了鲁棒性状态估计在实际中的应用。

综合起来看，已有的各种状态估计方法毫无例外都是一个非线性最优化问题，一般利用牛顿法或者内点法进行求解。从数学角度来看，这种非线性优化模型及其对应的求解方法具有几个不可避免的缺点：1）容易陷入局部最优点，即很难获得全局最优解，特别是一些鲁棒性状态估计；2）为了获得一个可以接受的解须要进行多次迭代，这种非线性最优化问题有时收敛困难，甚至不收敛，并且求解比较费时；3）雅可比矩阵以及残差灵敏度矩阵不是常数矩阵，因此在利用LNR进行不良数据辨识时需要进行多次删去不良数据和重新运行状态估计的循环操作，而如果利用EI方法辨识不良数据，估计出多个量测误差时，同样需要重新运行非线性状态估计。学术界对以上三个问题进行了一些研究，例如国外学者曾提出利用信赖域法来增强状态估计的收敛性，但是到目前为止，以上三个问题并没有得到很好的解决，其根本原因在于传统状态估计模型的非线性特性所造成的。如果能够建立精确线性化的状态估计模型，则可以彻底解决以上三个问题，从而在理论上保证得到更为科学的状态估计结果，促进状态估计的进一步实用化。

发明内容

本发明旨在彻底解决上述技术问题或至少提供一种有用的商业选择。为此，本发明的一个目的在于提出一种结果更为科学且计算效率更高的两阶段线性加权最小二乘电力***状态估计方法（two-stage linear WLS,TLWLS）。

根据本发明实施例的两阶段线性加权最小二乘电力***状态估计方法，包括以下步骤：A.形成网络模型，计算节点导纳矩阵及支路-节点关联矩阵；B.对量测矢量和状态矢量进行变换；C.形成精确线性化的量测方程；D.进行第一阶段的线性加权最小二乘估计，得到变换后的状态矢量的估计值；E进行逆变换，并进行第二阶段的线性加权最小二乘估计，得到所有节点的电压幅值和相角的估计值；以及F.进行不良数据辨识。

在本发明的一个实施例中，所述步骤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分别为节点导纳矩阵的实部和虚部；以及计算支路-节点关联矩阵为A={a_ij}(1≤i≤b,1≤j≤N-1)，其中各个元素定义为：

在本发明的一个实施例中，所述步骤B包括：将状态矢量变换为 $X={[v_1^2,v_2^2,...,v_N^2,{v_l}_i{v_l}_j\cos {\mathrm{\θ}}_{l_i{l}_j}(1\≤l\≤b),{v_l}_i{v}_{l_j}\sin {\mathrm{\θ}}_{l_i{l}_j}(1\≤l\≤b)]}^T,$ 其中，N为网络中所有节点的总数目，b为网络中所有支路的数目，l为支路编号，l_i和l_j为支路l的两端节点号，和分别是节点l_i和l_j的电压幅值，和分别是节点l_i和l_j的相角，为相角差，代表所有的b条支路对状态矢量X的贡献，也代表所有的b条支路对状态矢量X的贡献，X∈R^N+2b为状态矢量；以及将量测矢量变换为y∈R^m，包括节点电压幅值的平方、支路有功、支路无功、注入有功、注入无功，支路电流幅值的平方，其中，m为量测量的总个数，当用变换后的状态矢量X表示时，节点电压幅值的平方为v_i为节点i的电压，从节点i到节点j的支路有功为从节点i到节点j的支路无功为 $Q_{\mathrm{ij}}=-v_i^2(b_{\mathrm{si}}+b_{\mathrm{ij}})+v_i{v}_j{b}_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}}-v_i{v}_j{g}_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}},$ 节点i的注入有功为 $P_i=v_i\underset{j\∈N_i}{\mathrm{\Σ}}{v}_j(G_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}}),$ 节点i的注入无功为 $Q_i=v_i\underset{j\∈N_i}{\mathrm{\Σ}}{v}_j(G_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}}-B_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}}),$ G_ij+jB_ij为节点导纳矩阵中的对应元素，支路电流幅值的平方为其中，I_ij为π型支路ij的电流幅值，A＝(g_si+g_ij)²+(b_si+b_ij)²， D=-g_sib_ij+b_sig_ij。

在本发明的一个实施例中，所述步骤C包括：设J∈R^m×(N+2b)为雅可比矩阵，其中，节点电压幅值量测的平方对应的雅可比矩阵元素为 $\frac{{\∂v}_i^2}{{\∂v}_i^2}=1,$ $\frac{{\∂v}_i^2}{{\∂v}_j^2}=0,$ $\frac{{\∂v}_i^2}{{\∂v}_i{v}_j\cos {\mathrm{\θ}}_{\mathrm{ij}}}=0,$ $\frac{{\∂v}_i^2}{{\∂v}_i{v}_j\sin {\mathrm{\θ}}_{\mathrm{ij}}}=0,$ 支路功率量测对应的雅可比矩阵元素为 $\frac{{\∂P}_{\mathrm{ij}}}{{\∂v}_i^2}=g_{\mathrm{si}}+g_{\mathrm{ij}},$ $\frac{{\∂P}_{\mathrm{ij}}}{{\∂v}_j^2}=0,$ $\frac{{\∂P}_{\mathrm{ij}}}{{\∂v}_i{v}_j\cos {\mathrm{\θ}}_{\mathrm{ij}}}=-g_{\mathrm{ij}},$ $\frac{{\∂P}_{\mathrm{ij}}}{{\∂v}_i{v}_j\sin {\mathrm{\θ}}_{\mathrm{ij}}}=-b_{\mathrm{ij}},$ $\frac{{\∂Q}_{\mathrm{ij}}}{{\∂v}_i^2}=-(b_{\mathrm{si}}+b_{\mathrm{ij}}),$ $\frac{{\∂Q}_{\mathrm{ij}}}{{\∂v}_j^2}=0,$ $\frac{{\∂Q}_{\mathrm{ij}}}{{\∂v}_i{v}_j\cos {\mathrm{\θ}}_{\mathrm{ij}}}=b_{\mathrm{ij}},$ $\frac{{\∂Q}_{\mathrm{ij}}}{{\∂v}_i{v}_j\sin {\mathrm{\θ}}_{\mathrm{ij}}}=-g_{\mathrm{ij}},$ 注入功率量测对应的雅可比矩阵元素为 $\frac{{\∂P}_i}{{\∂v}_i^2}=G_{\mathrm{ii}},$ $\frac{{\∂P}_i}{{\∂v}_j^2}=0,$ $\frac{{\∂P}_i}{{\∂v}_i{v}_j\cos {\mathrm{\θ}}_{\mathrm{ij}}}=G_{\mathrm{ij}},$ $\frac{{\∂P}_i}{{\∂v}_i{v}_j\sin {\mathrm{\θ}}_{\mathrm{ij}}}=B_{\mathrm{ij}},$ $\frac{{\∂Q}_i}{{\∂v}_i^2}=-B_{\mathrm{ii}},$ $\frac{{\∂Q}_i}{{\∂v}_j^2}=0,$ $\frac{{\∂Q}_i}{{\∂v}_i{v}_j\cos {\mathrm{\θ}}_{\mathrm{ij}}}=-B_{\mathrm{ij}},$ 支路电流幅值量测的平方对应的雅可比矩阵元素为 以及根据步骤B得到的变换后的量测矢量和状态矢量，得到精确线性化的量测方程为：y=JX+τ，其中，τ∈R^m为量测误差矢量，J∈R^m×(N+2b)为常数雅可比矩阵。

在本发明的一个实施例中，所述步骤D包括：构造 $\mathrm{MinJ}\left(X\right)={\mathrm{\Σ}}_{i=1}^m{(y_i-J_{i}X)}^2/R_{\mathrm{ii}}={(y-\mathrm{JX})}^{T}W(y-\mathrm{JX})$ 的线性加权最小二乘问题，其中，W=R^-1为权重矩阵，其中最优解应满足条件记G＝J^TWJ为信息矩阵，求解得到变换的状态矢量X的估计值X＝G^-1J^TWy。

在本发明的一个实施例中，所述步骤E包括：利用所述步骤D得到的状态矢量X，根据公式 $\left\{\begin{array}{c}{v}_i=\sqrt{v_i^2(1\≤i\≤N)}\\ {\mathrm{\θ}}_{l_i{l}_j}=\arcsin (v_{l_i}{v}_{l_j}\sin {\mathrm{\θ}}_{l_i{l}_j}/\left(v_{l_i}{v}_{l_j}\right))(1\≤l\≤b)\\ {\mathrm{\θ}}_{l_i{l}_j}=\arccos (v_{l_i}{v}_{l_j}\cos {\mathrm{\θ}}_{l_i{l}_j}/\left(v_{l_i}{v}_{l_j}\right))(1\≤l\≤b)\end{array}\right.$ 进行逆变换，得到所有节点电压的幅值以及所有支路两端相角差的估计值θ_2b∈R^2b，即 ${\mathrm{\θ}}_{2b}={[\arcsin (v_{l_i}{v}_{l_j}\sin {\mathrm{\θ}}_{l_i{l}_j}/\left(v_{l_i}{v}_{l_j}\right)),\arccos (v_{l_i}{v}_{l_j}\cos {\mathrm{\θ}}_{l_i{l}_j}/\left(v_{l_i}{v}_{l_j}\right))(1\≤l\≤b)]}^T;$ 以及利用所述所有支路两端相角差的估计值θ_2b构造Min J(θ)＝(θ_2b-A₂θ)^TW_θ(θ_2b-A₂θ)的线性加权最小二乘问题，其中W_θ为权重矩阵，取值为单位矩阵，A₂＝[A^T A^T]^T，A为步骤A得到的支路-节点关联矩阵，θ∈R^N-1为除参考节点外的所有节点的相角，其中最优解应满足条件可得将A代入A₂，并求解可得θ＝(A^T A)^-1A^Tθ_b，其中，

${\mathrm{\θ}}_b={\left[\right[\arcsin (v_{l_i}{v}_{l_j}\sin {\mathrm{\θ}}_{l_i{l}_j}/\left(v_{l_i}{v}_{l_j}\right))+\arccos (v_{l_i}{v}_{l_j}\cos {\mathrm{\θ}}_{l_i{l}_j}/\left(v_{l_i}{v}_{l_j}\right))]/2(1\≤l\≤b)]}^T.$

在本发明的一个实施例中，所述步骤F包括：计算第i个量测的正则化残差值其中，r＝[I-J[J^TWJ]^-1J^TW]τ＝Sτ，S＝I-J(J^TWJ)^-1J^TW为残差灵敏度矩阵，是常数矩阵，r~N(0，Ω)，其中Ω＝SR；以及若某个量测的正则化残差值大于预定值，则去掉该量测，重新运行所述D步骤，直至所有量测的正则化残差都小于所述预定值。

本发明的两阶段线性加权最小二乘电力***状态估计方法（TLWLS）的优点至少包括：TLWLS模型本质是一个二次规划问题，从数学上保证可求得唯一的全局最优解；TLWLS仅须求解线性方程组，无需迭代，不存在收敛性问题，估计效率远高于常规的非线性加权最小二乘状态估计方法（WLS）；TLWLS的雅可比矩阵和残差灵敏度矩阵是常数矩阵，这给不良数据的辨识提供了便利，可提高不良数据的辨识效率。综上，本发明提出的两阶段线性加权最小二乘状态估计方法可以得到更为科学的状态估计结果，而且计算效率更高，具有很好的工程应用前景。

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

附图说明

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

图1为本发明的两阶段线性加权最小二乘电力***状态估计方法的流程图；

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

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

图4为某个三节点***单线图及量测配置图。

具体实施方式

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

如图1所示，本发明的电力***状态估计量测方程的精确线性化方法包括如下步骤：

步骤S101，形成网络模型，计算节点导纳矩阵及支路-节点关联矩阵。

具体地，将网络中的三绕组变压器等效为三个两绕组变压器，则网络中所有的线路和变压器可以用统一的π型支路表示，如图2所示。图2中，y_s＝1/(r_ij+jx_ij)＝g_s+jb_s为支路ij的串联电纳；r_ij+jx_ij为串联阻抗值；b_c为支路的接地电纳。其中，若π型支路ij为变压器支路，则b_c=0且k为理想变压器的变比；若π型支路ij为普通线路，则k＝1，并联的多条支路等效为一条支路。需要说明的是，并联的多条支路等效为一条支路。

图2的π型支路的等值电路如图3所示。图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分别为节点导纳矩阵的实部和虚部，以及形成支路-节点关联矩阵为A={a_ij}(1≤i≤b,1≤j≤N-1)，其各元素定义为：

步骤S102，对量测矢量和状态矢量进行变换。

首先，将状态矢量变换为 $X={[v_1^2,v_2^2,...,v_N^2,{v_l}_i{v_l}_j\cos {\mathrm{\θ}}_{l_i{l}_j}(1\≤l\≤b),{v_l}_i{v}_{l_j}\sin {\mathrm{\θ}}_{l_i{l}_j}(1\≤l\≤b)]}^T.$

其中，N为网络中所有节点的总数目；b为网络中所有支路的数目(并联的多条支路等效为一条支路)；l为支路编号，l_i和l_j为支路l的两端节点号，和分别是节点l_i和l_j的电压幅值，和分别是节点l_i和l_j的相角，为相角差；代表所有的b条支路对状态矢量X的贡献，也代表所有的b条支路对状态矢量X的贡献；X∈R^N+2b为状态矢量。

其次，将量测矢量变换为y∈R^m，包括的量测类型有：节点电压幅值的平方、支路有功、支路无功、注入有功、注入无功，支路电流幅值的平方；m为量测量的总个数。则量测量可以用状态矢量X进行表达。

1）电压幅值量测的平方

$v_i^2=v_i^2---\left(1\right)$

其中，v_i为节点i的电压。

2）从节点i到节点j的支路有功和无功量测

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

$Q_{\mathrm{ij}}=-v_i^2(b_{\mathrm{si}}+b_{\mathrm{ij}})+v_i{v}_j{b}_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}}-v_i{v}_j{g}_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}}---\left(3\right)$

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

3）节点i的注入有功和注入无功量测

$P_i=v_i\underset{j\∈N_i}{\mathrm{\Σ}}{v}_j(G_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}})---\left(4\right)$

$Q_i=v_i\underset{j\∈N_i}{\mathrm{\Σ}}{v}_j(G_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}}-B_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}})---\left(5\right)$

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

4）支路电流幅值量测的平方

$I_{\mathrm{ij}}^2={\mathrm{Av}}_i^2+{\mathrm{Bv}}_j^2-{2v}_i{v}_j(C\cos {\mathrm{\θ}}_{\mathrm{ij}}-D\sin {\mathrm{\θ}}_{\mathrm{ij}})---\left(6\right)$

其中，I_ij为支路ij的电流幅值；A＝(g_si+g_ij)²+(b_si+b_ij)²； D=-g_sib_ij+b_sig_ij。

步骤S103，形成精确线性化的量测方程。

根据以上的变换后的量测矢量和新的状态矢量，即可得到新的精确线性化的量测方程为：

y=JX+τ                    （7）

其中，y∈R^m为变换后的量测矢量，包括节点电压幅值的平方、支路有功、支路无功、注入有功、注入无功，支路电流幅值的平方；X∈R^N+2b为变换后的状态矢量；τ~N(0,R)为量测误差矢量，其中其中为τ_i的方差；J∈R^m×(N+2b)为常数雅可比矩阵，其各个部分的元素的表达式如下所示。

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

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

$\frac{{\∂v}_i^2}{{\∂v}_i^2}=1,$ $\frac{{\∂v}_i^2}{{\∂v}_j^2}=0,$ $\frac{{\∂v}_i^2}{{\∂v}_i{v}_j\cos {\mathrm{\θ}}_{\mathrm{ij}}}=0,$ $\frac{{\∂v}_i^2}{{\∂v}_i{v}_j\sin {\mathrm{\θ}}_{\mathrm{ij}}}=0.$

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

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

$\frac{{\∂P}_{\mathrm{ij}}}{{\∂v}_i^2}=g_{\mathrm{si}}+g_{\mathrm{ij}},$ $\frac{{\∂P}_{\mathrm{ij}}}{{\∂v}_j^2}=0,$ $\frac{{\∂P}_{\mathrm{ij}}}{{\∂v}_i{v}_j\cos {\mathrm{\θ}}_{\mathrm{ij}}}=-g_{\mathrm{ij}},$ $\frac{{\∂P}_{\mathrm{ij}}}{{\∂v}_i{v}_j\sin {\mathrm{\θ}}_{\mathrm{ij}}}=-b_{\mathrm{ij}},$

$\frac{{\∂Q}_{\mathrm{ij}}}{{\∂v}_i^2}=-(b_{\mathrm{si}}+b_{\mathrm{ij}}),$ $\frac{{\∂Q}_{\mathrm{ij}}}{{\∂v}_j^2}=0,$ $\frac{{\∂Q}_{\mathrm{ij}}}{{\∂v}_i{v}_j\cos {\mathrm{\θ}}_{\mathrm{ij}}}=b_{\mathrm{ij}},$ $\frac{{\∂Q}_{\mathrm{ij}}}{{\∂v}_i{v}_j\sin {\mathrm{\θ}}_{\mathrm{ij}}}=-g_{\mathrm{ij}}.$

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

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

$\frac{{\∂P}_i}{{\∂v}_i^2}=G_{\mathrm{ii}},$ $\frac{{\∂P}_i}{{\∂v}_j^2}=0,$ $\frac{{\∂P}_i}{{\∂v}_i{v}_j\cos {\mathrm{\θ}}_{\mathrm{ij}}}=G_{\mathrm{ij}},$ $\frac{{\∂P}_i}{{\∂v}_i{v}_j\sin {\mathrm{\θ}}_{\mathrm{ij}}}=B_{\mathrm{ij}},$

$\frac{{\∂Q}_i}{{\∂v}_i^2}=-B_{\mathrm{ii}},$ $\frac{{\∂Q}_i}{{\∂v}_j^2}=0,$ $\frac{{\∂Q}_i}{{\∂v}_i{v}_j\cos {\mathrm{\θ}}_{\mathrm{ij}}}=-B_{\mathrm{ij}},$ $\frac{{\∂Q}_i}{{\∂v}_i{v}_j\sin {\mathrm{\θ}}_{\mathrm{ij}}}=G_{\mathrm{ij}}.$

4）支路电流幅值量测的平方对应的雅可比矩阵元素

对于支路电流幅值量测的平方，其对应的雅可比矩阵元素为

$\frac{{\∂I}_{\mathrm{ij}}^2}{{\∂v}_i^2}=A,$ $\frac{{\∂I}_{\mathrm{ij}}^2}{{\∂v}_j^2}=B,$ $\frac{{\∂I}_{\mathrm{ij}}^2}{{\∂v}_i{v}_j\cos {\mathrm{\θ}}_{\mathrm{ij}}}=-2C,$ $\frac{{\∂I}_{\mathrm{ij}}^2}{{\∂v}_i{v}_j\sin {\mathrm{\θ}}_{\mathrm{ij}}}=-2D.$

需要说明的是，上述式（7）即为精确线性化的量测方程。基于该量测方程，可以采用已有的任何一种状态估计方法求解得到变换后的状态矢量X，得到X后，即可进一步得到所有的支路功率、节点注入功率以及支路电流幅值等的估计值，这样也就获得了对全网状态的精确感知。

步骤S104，进行第一阶段的线性加权最小二乘估计，得到变换后状态矢量X的估计值。

具体地，求解下面的线性加权最小二乘问题，获得状态矢量X的估计值。

$\mathrm{MinJ}\left(X\right)={\mathrm{\Σ}}_{i=1}^m{(y_i-J_{i}X)}^2/R_{\mathrm{ii}}={(y-\mathrm{JX})}^{T}W(y-\mathrm{JX})---\left(8\right)$

其中W=R^-1为权重矩阵。为了获得最优值，必须满足最优条件：即

GX=J^TWy                    （9）

其中，G=J^TWJ为信息矩阵。对信息矩阵G进行正交分解，即可求解得到状态矢量X的估计值，即X＝G^-1J^TWy。

步骤S105，进行逆变换，并进行第二阶段的线性加权最小二乘估计，得到所有节点的电压幅值和相角的估计值。

根据变换后的状态矢量X的定义，利用反变换可以得到所有节点电压幅值以及所有支路两端相角差的估计值，即

$\left\{\begin{array}{c}{v}_i=\sqrt{v_i^2(1\≤i\≤N)}\\ {\mathrm{\θ}}_{l_i{l}_j}=\arcsin (v_{l_i}{v}_{l_j}\sin {\mathrm{\θ}}_{l_i{l}_j}/\left(v_{l_i}{v}_{l_j}\right))(1\≤l\≤b)\\ {\mathrm{\θ}}_{l_i{l}_j}=\arccos (v_{l_i}{v}_{l_j}\cos {\mathrm{\θ}}_{l_i{l}_j}/\left(v_{l_i}{v}_{l_j}\right))(1\≤l\≤b)\end{array}\right.---\left(10\right)$

在式（10）中，arccos的符号应该与arcsin的符号保持一致。

至此，已经得到所有节点电压幅值的估计值，但是所有节点的相角还是未知的，显然可以利用所有支路两端的相角差进行估计。

令 ${\mathrm{\θ}}_{2b}={\left[\arcsin (v_{l_i}{v}_{l_j}\sin {\mathrm{\θ}}_{l_i{l}_j}/\left(v_{l_i}{v}_{l_j}\right))(1\≤l\≤b),\arccos (v_{l_i}{v}_{l_j}\cos {\mathrm{\θ}}_{l_i{l}_j}/\left(v_{l_i}{v}_{l_j}\right))(1\≤l\≤b)\right]}^T,$ 显然，可通过式（10）得到θ_2b。以下利用θ_2b估计得到除参考节点外的所有节点的相角。

令θ＝[θ₂，…,θ_N]^T代表所有节点的相角(参考节点除外)，则θ_2b与θ的关系是

θ_2b=A₂θ                                 （11）

其中，A₂＝[A^T A^T]^T，A₂为2b×(N-1)的矩阵。

因为θ_2b是由第一阶段的状态估计结果推导得到的，所以可将θ_2b视为含有噪声的量测量。而θ为待估计值，则式（11）可以改写为下面的关于θ的量测方程

θ_2b=A₂θ+τ                              （12）

则可构造如下的线性加权最小二乘问题，对θ进行估计

Min J(θ)=(θ_2b-A₂θ)^TW_θ(θ_2b-A₂θ)      （13）

其中，W_θ为权重矩阵，不失一般性，可取W_θ为单位矩阵。当式（13）获得最优值，须满足

$A_2^T{A}_2\mathrm{\θ}=A_2^T{\mathrm{\θ}}_{2b}---\left(14\right)$

将A代入A₂，可得

A^TAθ＝A^Tθ_b                    （15）

其中， ${\mathrm{\θ}}_b={\left[\right[\arcsin (v_{l_i}{v}_{l_j}\sin {\mathrm{\θ}}_{l_i{l}_j}/\left(v_{l_i}{v}_{l_j}\right))+\arccos (v_{l_i}{v}_{l_j}\cos {\mathrm{\θ}}_{l_i{l}_j}/\left(v_{l_i}{v}_{l_j}\right))]/2(1\≤l\≤b)]}^T$ 为b维矢量，显然可由θ_2b得到θ_b。

由式（15）可得θ的估计值为θ＝(A^T A)^-1A^Tθ_b。

至此，已经得到所有节点的电压幅值和相角的估计值。

步骤S106，进行不良数据辨识。

在本发明中，仅仅需要在两阶段线性加权最小二乘状态估计的第一阶段进行不良数据的辨识即可。

第一阶段线性加权最小二乘状态估计的残差为

r=y-JX                            （16）

其中，X为第一阶段状态估计的结果，由式（9）求解得到。

由式（9）、（16）可以推导得到残差与量测误差的关系为

r＝[I-J[J^TWJ]^-1J^TW]τ＝Sτ        （17）

其中，S＝I-J(J^TWJ)^-1J^TW为残差灵敏度矩阵。由于雅可比矩阵J是常数矩阵，因此残差灵敏度矩阵S亦是常数矩阵。这意味着残差灵敏度矩阵不受不良数据的影响。

容易证明，r~N(0,Ω)，其中Ω＝SR，则可得到正则化残差（normalized residual）r^N为

${r_i}^N=\frac{\left|r_i\right|}{\sqrt{{\mathrm{\Ω}}_{\mathrm{ii}}}}---\left(18\right)$

其中，|r_i|是残差r_i的绝对值。

如果则第k个量测即被认为是不良数据，可将其从量测集中去掉，重新运行第一阶段的线性加权最小二乘状态估计。重复以上步骤，直到所有的正则化残差都小于3，即可完成对所有不良数据的辨识。由于雅可比矩阵是常数矩阵，显然基于线性加权最小二乘状态估计的不良数据辨识效率是非常高的。

为使本领域技术人员更好地理解，申请人以一个三节点***为例，来说明本发明的方法。该三节点***的单线图及量测配置图如图4所示，其网络参数如表1所示。选取节点1为参考节点。

表1三节点***网络参数

Table 1 The network data of the 3-bus system

（1）常规负荷时非线性WLS及本发明提出的TLWLS的估计结果常规负荷时的量测量如表2所示。

表2三节点***量测量

Table 2 Measurements of the 3-bus system

若采用传统的非线性WLS，当收敛精度为10^-6时，经过4次迭代非线性WLS收敛，所得到的所有节点的电压幅值和相角的估计值如表3所示。

表3非线性WLS的估计结果

Table 3 The Estimated State Vector by Nonlinear WLS

以下采用本发明提出的TLWLS进行估计。

将状态矢量变换为 $X={[v_1^2,v_2^2,v_3^2,v_1{v}_2\cos {\mathrm{\θ}}_{12},v_1{v}_3\cos {\mathrm{\θ}}_{13},v_2{v}_3\cos {\mathrm{\θ}}_{23},v_1{v}_2\sin {\mathrm{\θ}}_{12},v_1{v}_3\sin {\mathrm{\θ}}_{13},v_2{v}_3\sin {\mathrm{\θ}}_{23}]}^T,$ 将量测矢量变换为 $y={[v_1^2,v_2^2,Q_{12},Q_{21},Q_{13},Q_{31},Q_2,P_{12},P_{21},P_{13},P_{31},P_2]}^T,$ 则可得到精确线性化的量测方程为y＝JX+τ。其中常数雅可比可由S103给出的公式求得。通过求解式（9）可变换后的状态矢量X的估计值，进一步由式（10）可得所有节点电压幅值以及所有支路两端相角差的估计值，如表4所示。在表4中，和分别是X_i、v_i和θ_2bi的估计值。

表4两阶段线性加权最小二乘状态估计第一阶段的估计结果

Table 4 The estimation value by the first stage linear WLS of the TLWLS

根据该网路的接线图，可得该网络的支路-节点关联矩阵（节点1是参考节点）是

$A=\left[\begin{array}{cc}-1& 0\\ 0& -1\\ 1& -1\end{array}\right]$

通过求解式（15），可得所有节点相角的估计值为θ=[θ₂,θ₃]^T=[-0.0213,-0.0479]^T。从而由本发明提出的两阶段线性加权最小二乘状态估计得到的所有节点的电压幅值和相角的估计值如表5所示。

表5两阶段线性加权最小二乘状态估计的估计结果

Table 5 The Estimated State Vector by TLWLS

对比表3和表5，可以发现，在常规负荷时由本发明提出的两阶段线性加权最小二乘状态估计（TLWLS）得到的状态变量的估计值与由传统非线性WLS得到的状态变量的估计值吻合的很好。但是TLWLS无需迭代，而传统的非线性WLS需要4次迭代才能收敛，在本算例中，TLWLS的计算效率是非线性WLS的3倍多。

（2）非线性WLS陷入局部最优解的算例

以下仍以该3节点***为例，给出非线性WLS陷入局部最优解的例子。当***负荷持续增大时，可以得到另一组量测量，如表6所示。显然，此时发生了电压崩溃。

表6重负荷下三节点***量测量

Table 6 Measurements of the 3-bus system at hevy load

用非线性WLS对该算例进行估计时，经过35次收敛，收敛时非线性WLS的目标函数数值为0.1088。对于此算例，状态变量的真值、非线性WLS的估计值以及两阶段线性加权最小二乘状态估计(TLWLS)的估计值如表7所示。

表7重负荷时状态变量的真值、非线性WLS及TLWLS的估计结果

Table 7 The True value,Estimate Value by WLS and TLWLS at hevy load

由表7可见，在重负荷特别是在电压崩溃时，传统的非线性WLS陷入局部最优解，其估计值出现了负的电压幅值和很大的相角值。由此可见，在重负荷特别是在电压崩溃时，传统非线性WLS的估计结果不可靠，而此时电力***却要求状态估计具有尽可能好的性能。在重负荷特别是电压崩溃时，之所以会产生这种尖锐的矛盾是由传统非线性WLS的固有特性所造成。

由表7同样可见，即使在电压崩溃时，本发明提出的两阶段线性加权最小二乘状态估计方法（TLWLS）的估计结果也是非常精确的，即在任何时候本发明提出的TLWLS均可得到全局最优解，且不需迭代。因而本发明提出的TLWLS可以得到更为科学的状态估计结果，且计算效率更高。

(3)辨识不良数据

在表2所示的常规负荷时的量测量P₁₂上迭代20%的噪声，从而P₁₂变为不良数据。由式（18）可得P₁₂的正则化残差，发现其大于3，将其从量测集中去掉，重新运行第一阶段的线性加权最小二乘状态估计，发现所有的正则化残差都小于3，不良数据辨识结束。可见，正则化残差检验法可以有效进行第一阶段的TLWLS的不良数据辨识，而且发现其辨识效率是传统基于非线性WLS的正则化残差检验法辨识效率的3倍多。

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

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

Claims (1)

1.一种两阶段线性加权最小二乘电力***状态估计方法，其特征在于，包括以下步骤：

A.形成网络模型，计算节点导纳矩阵及支路-节点关联矩阵；

B.对量测矢量和状态矢量进行变换；

C.形成精确线性化的量测方程；

D.进行第一阶段的线性加权最小二乘估计，得到变换后的状态矢量的估计值；

E.进行逆变换，并进行第二阶段的线性加权最小二乘估计，得到所有节点的电压幅值和相角的估计值；以及

F.进行不良数据辨识，

其中，所述步骤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分别为节点导纳矩阵的实部和虚部；以及

计算支路-节点关联矩阵A={ai_j}(1≤i≤b,1≤j≤N-1)，所述N为网络中所有节点的总数目，b为网络中所有支路的总数目，其中各个元素定义为：

其中，所述步骤B包括：

将状态矢量变换为 $X={[v_1^2,v_2^2,...,v_N^2,v_{l_i}{v}_{l_j}\cos {\mathrm{\θ}}_{l_i{l}_j}(1\≤l\≤b),v_{l_i}{v}_{l_j}\sin {\mathrm{\θ}}_{l_i{l}_j}(1\≤l\≤b)]}^T,$ 其中，l为支路编号，l_i和l_j为支路l的两端节点号，和分别是节点l_i和l_j的电压幅值，和分别是节点l_i和l_j的相角，为相角差，(1≤l≤b)代表所有的b条支路对状态矢量X的贡献，(1≤l≤b)也代表所有的b条支路对状态矢量X的贡献，X∈R^N+2b为状态矢量；以及

将量测矢量变换为y∈R^m，包括节点电压幅值的平方、支路有功、支路无功、注入有功、注入无功，支路电流幅值的平方，其中，m为量测量的总个数，当用变换后的状态矢量X表示时，节点电压幅值的平方为v_i为节点i的电压，从节点i到节点j的支路有功为 $P_{\mathrm{ij}}=v_i^2(g_{\mathrm{si}}+g_{\mathrm{ij}})-v_i{v}_j{g}_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}}-v_i{v}_j{b}_{\mathrm{ij}}{\sin \mathrm{\θ}}_{\mathrm{ij}},$ 从节点i到节点j的支路无功为 $Q_{\mathrm{ij}}=-v_i^2(b_{\mathrm{si}}+b_{\mathrm{ij}})+v_i{v}_j{b}_{\mathrm{ij}}{\cos \mathrm{\θ}}_{\mathrm{ij}}-v_i{v}_j{g}_{\mathrm{ij}}{\sin \mathrm{\θ}}_{\mathrm{ij}},$ 节点i的注入有功为 $P_i=v_i\underset{j\∈N_i}{\mathrm{\Σ}}{v}_j(G_{\mathrm{ij}}{\cos \mathrm{\θ}}_{\mathrm{ij}}+B_{\mathrm{ij}}{\sin \mathrm{\θ}}_{\mathrm{ij}}),$ 节点i的注入无功为G_ij+jB_ij为节点导纳矩阵中的对应元素，支路电流幅值的平方为 $I_{\mathrm{ij}}^2=A^{\′}{v}_i^2+B^{\′}{v}_j^2-{2v}_i{v}_j({C\cos \mathrm{\θ}}_{\mathrm{ij}}-{D\sin \mathrm{\θ}}_{\mathrm{ij}}),$ 其中，I_ij为π型支路ij的电流幅值，A'=(g_si+g_ij)²+(b_si+b_ij)2， $B^{\′}=g_{\mathrm{ij}}^2+b_{\mathrm{ij}}^2,$ $C=g_{\mathrm{ij}}^2+b_{\mathrm{ij}}^2+g_{\mathrm{si}}{g}_{\mathrm{ij}}+b_{\mathrm{si}}{b}_{\mathrm{ij}},$ D=-g_sib_ij+b_sig_ij，

其中，所述步骤C包括：

设J∈R^m×(N+2b)为雅可比矩阵，其中，节点电压幅值量测的平方对应的雅可比矩阵元素为 $\frac{{\∂v}_i^2}{{\∂v}_i^2}=1,\frac{{\∂v}_i^2}{{\∂v}_j^2}=0,\frac{{\∂v}_i^2}{{\∂v}_i{v}_j{\cos \mathrm{\θ}}_{\mathrm{ij}}}=0,\frac{{\∂v}_i^2}{{\∂v}_i{v}_j{\sin \mathrm{\θ}}_{\mathrm{ij}}}=0,$ 支路功率量测对应的雅可比矩阵元素为

$\frac{{\∂P}_{\mathrm{ij}}}{{\∂v}_i^2}=g_{\mathrm{si}}+g_{\mathrm{ij}},\frac{{\∂P}_{\mathrm{ij}}}{{\∂v}_j^2}=0,\frac{{\∂P}_{\mathrm{ij}}}{{\∂v}_i{v}_j{\cos \mathrm{\θ}}_{\mathrm{ij}}}=-g_{\mathrm{ij}},\frac{{\∂P}_{\mathrm{ij}}}{{\∂v}_i{v}_j{\sin \mathrm{\θ}}_{\mathrm{ij}}}=-b_{\mathrm{ij}},\frac{{\∂Q}_{\mathrm{ij}}}{{\∂v}_i^2}=-(b_{\mathrm{si}}+b_{\mathrm{ij}}),\frac{{\∂Q}_{\mathrm{ij}}}{{\∂v}_j^2}=0,\frac{{\∂Q}_{\mathrm{ij}}}{{\∂v}_i{v}_j{\cos \mathrm{\θ}}_{\mathrm{ij}}}=b_{\mathrm{ij}},$ $\frac{{\∂Q}_{\mathrm{ij}}}{{\∂v}_i{v}_{j}s{\mathrm{in\θ}}_{\mathrm{ij}}}=-g_{\mathrm{ij}},$ 注入功率量测对应的雅可比矩阵元素为 $\frac{{\∂P}_i}{{\∂v}_i^2}=G_{\mathrm{ii}},\frac{{\∂P}_i}{{\∂v}_j^2}=0,\frac{{\∂P}_i}{{\∂v}_i{v}_j{\cos \mathrm{\θ}}_{\mathrm{ij}}}=G_{\mathrm{ij}},$ $\frac{{\∂P}_i}{{\∂v}_i{v}_j{\sin \mathrm{\θ}}_{\mathrm{ij}}}=B_{\mathrm{ij}},\frac{{\∂Q}_i}{{\∂v}_i^2}=-B_{\mathrm{ii}},\frac{{\∂Q}_i}{{\∂v}_j^2}=0,\frac{{\∂Q}_i}{{\∂v}_i{v}_j{\cos \mathrm{\θ}}_{\mathrm{ij}}}=-B_{\mathrm{ij}},\frac{{\∂Q}_i}{{\∂v}_i{v}_j{\sin \mathrm{\θ}}_{\mathrm{ij}}}=G_{\mathrm{ij}},$ 支路电流幅值量测的平方对应的雅可比矩阵元素为 $\frac{{\∂I}_{\mathrm{ij}}^2}{{\∂v}_i^2}=A^{\′},\frac{{\∂I}_{\mathrm{ij}}^2}{{\∂v}_j^2}=B^{\′},\frac{{\∂I}_{\mathrm{ij}}^2}{{\∂v}_i{v}_j{\cos \mathrm{\θ}}_{\mathrm{ij}}}=-2C,\frac{{\∂I}_{\mathrm{ij}}^2}{{\∂v}_i{v}_j{\sin \mathrm{\θ}}_{\mathrm{ij}}}=2D;$ 以及

根据步骤B得到的变换后的量测矢量和状态矢量，得到精确线性化的量测方程为：y=JX+τ，其中，τ∈R^m为量测误差矢量，J∈R^m×(N+2b)为常数雅可比矩阵，

其中，所述步骤D包括：

构造 $\mathrm{MinJ}\left(X\right)={\mathrm{\Σ}}_{i=1}^m{(y_i-J_{i}X)}^2/R_{\mathrm{ii}}={(y-\mathrm{JX})}^{T}W(y-\mathrm{JX})$ 的线性加权最小二乘问题，其中，W=R^-1为权重矩阵，其中最优解应满足条件记G=J^TWJ为信息矩阵，求解得到变换的状态矢量X的估计值X=G^-1J^TWy，

其中：

利用所述步骤D得到的状态矢量X，根据公式 $\left\{\begin{array}{c}{v}_i=\sqrt{v_i^2}(1\≤i\≤N)\\ {\mathrm{\θ}}_{l_i{l}_j}=\arcsin (v_{l_i}{v}_{l_j}{\sin \mathrm{\θ}}_{l_i{l}_j}/\left(v_{l_i}{v}_{l_j}\right))(1\≤l\≤b)\\ {\mathrm{\θ}}_{l_i{l}_j}=\arccos (v_{l_i}{v}_{l_j}{\cos \mathrm{\θ}}_{l_i{l}_j}/\left(v_{l_i}{v}_{l_j}\right))(1\≤l\≤b)\end{array}\right.$ 进行逆变换，得到所有节点电压的幅值以及所有支路两端相角差的估计值θ_2b∈R^2b，即

${\mathrm{\θ}}_{2b}={[\arcsin (v_{l_i}{v}_{l_j}{\sin \mathrm{\θ}}_{l_i{l}_j}/\left(v_{l_i}{v}_{l_j}\right)),\arccos (v_{l_i}{v}_{l_j}{\cos \mathrm{\θ}}_{l_i{l}_j}/\left(v_{l_i}{v}_{l_j}\right))(1\≤l\≤b)]}^T;$ 以及

利用所述所有支路两端相角差的估计值θ_2b构造Min J(θ)=(θ_2b-A₂θ)^TW_θ(θ_2b-A₂θ)的线性加权最小二乘问题，其中W_θ为权重矩阵，取值为单位矩阵，A₂=[A^T A^T]^T，A为步骤A得到的支路-节点关联矩阵，θ∈R^N-1为除参考节点外的所有节点的相角，其中最优解应满足条件可得将A代入A₂，并求解可得θ=(A^TA)^-1A^Tθ_b，其中，

${\mathrm{\θ}}_b={\left[\right[\arcsin (v_{l_i}{v}_{l_j}{\sin \mathrm{\θ}}_{l_i{l}_j}/\left(v_{l_i}{v}_{l_j}\right)),+\arccos (v_{l_i}{v}_{l_j}{\cos \mathrm{\θ}}_{l_i{l}_j}/\left(v_{l_i}{v}_{l_j}\right))]/2(1\≤l\≤b)]}^T,$

其中，所述步骤F包括：

计算第i个量测的正则化残差值其中，r=[I-J[J^TWJ]^-1J^TW]τ=Sτ，

S=I-J(J^TWJ)^-1J^TW为残差灵敏度矩阵，是常数矩阵，r～N(0,Ω)，其中Ω=SR，I为m维的单位矩阵；r为m维残差矢量，R为量测误差方差阵，为m维对角阵；以及

若某个量测的正则化残差值大于预定值，则去掉该量测，重新运行所述D步骤，直至所有量测的正则化残差都小于所述预定值。