CN105977963A - 一种基于两阶段法的主配网一体化分布式状态估计方法 - Google Patents

Abstract

本发明涉及一种基于两阶段法的主配网一体化分布式状态估计方法，方法包括：建立主配网一体化状态估计模型；将主配网一体化状态估计模型的配网模型解耦为若干子系统；建立解耦前主配网的状态量与解耦后子系统中状态量的线性关系；对解耦后子系统以分布式的方式独立进行状态估计；以解耦前后主配网状态量的线性关系为量测方程，建立状态估计模型进行状态估计，获取主配网的一体化状态。本发明提出的方法降低了计算量，有效且准确地解决了由于配网网络规模庞大造成的主配网一体化状态估计维数灾的问题，而且能够得到主配网一致的、精确的实时状态，为主配网一体化的分析与计算提供基础数据，进而提高了主配网一体化的运行可靠性。

Description

一种基于两阶段法的主配网一体化分布式状态估计方法

技术领域

本发明涉及配电自动化领域，具体涉及一种基于两阶段法的主配网一体化分布式状态估计方法。

背景技术

分布式电源与常规电源相结合的方式将是电力系统发展的趋势，大量分布式电源的接入，改变了配电网的特性，采用基于边界等值实现主网和配网相互独立运行分析控制的方法已不适用，需要从主配网一体化全局考虑，对主配网进行一体化分析。

分布式状态估计利用状态估计算法本身可分区域计算的特性，在各分区进行本区域的状态估计，再通过与其他区域的信息交换，从而使所有分区的状态估计结果达到或接近整体计算的效果，相对于集中式状态估计，不存在因增益矩阵的条件数过大、数值稳定性差等问题导致的整体状态估计失败等问题，研究的核心在于如何协调各个子分区状态估计结果，以期获得与集中式估计相似的计算精度。

目前，国内外已有许多专家学者对分布式状态估计进行研究，并取得了大量成果。按照各分区状态估计协调阶段的不同，可以分为迭代收敛后协调的分布式状态估计和迭代期间协调的分布式状态估计两种。

(1)迭代收敛后协调的分布式状态估计最大的优点是能够将残差污染控制在本区域内，各区域的收敛性互不影响，因此，该种分布式状态估计研究比较多，其核心在于如何协调各个子分区状态估计的估计结果，以期获得与集中式估计相近的计算精度，各个分区的协调机制对算法的计算量与计算结果影响较大，一般各区域需要反复协调，但是该种分布式状态估计对边界节点的注入量测处理较为粗糙。

(2)迭代期间协调的分布式状态估计在各个区域每次迭代后，均要进行一次统一协调，其优点是反复协调可以使得解更趋近最优解，缺点是协调侧与子分区或子分区之间需要频繁地交换数据，无法使残差污染控制在本区域内，容易造成残差污染的传播与扩大，各区域间的收敛性容易相互影响。

主配网一体化实时分析和计算对数据断面的基本要求是全网状态量的一致性，主配独立进行状态估计的最大问题是得到的全网状态量不一致，比如主配网互相的等值注入的估计结果不匹配，因此有必要对主配网进行一体化状态估计，由于主配网的网络特性差异较大、量测覆盖率相差较大，主网只需单相分析、配网需要三相分析，如果将主配网作为一个整体进行状态估计，在数学上收敛性无法保证，况且主网的状态估计精度也会因配网影响大大降低，所以将主配网的实时量测整合到一个系统中，进行一体化的分布式状态估计尤为必要。

发明内容

有鉴于此，本发明提供的一种基于两阶段法的主配网一体化分布式状态估计方法，该方法降低了计算量，有效且准确地解决了由于配网网络规模庞大造成的主配网一体化状态估计维数灾的问题，而且能够得到主配网一致的、精确的实时状态，为主配网一体化的分析与计算提供基础数据，进而提高了主配网一体化的运行可靠性。

本发明的目的是采用下述技术方案实现的：

一种基于两阶段法的主配网一体化分布式状态估计方法，其改进之处在于，包括：

步骤1.建立主配网一体化状态估计模型；

步骤2.在所述主配网一体化状态估计模型中添加零阻抗支路及虚拟母线，将所述主配网一体化状态估计模型的配网模型解耦为若干子系统；

步骤3.建立解耦前所述主配网的状态量与解耦后所述子系统中状态量的线性关系；

步骤4.第一阶段状态估计：对解耦后子系统以分布式的方式独立进行状态估计；

步骤5.第二阶段状态估计：以解耦前后主配网状态量的线性关系为量测方程，建立状态估计模型进行状态估计，获取主配网的一体化状态。

优选的，所述步骤1中，所述主配网一体化状态估计模型包括：主网模型和配网模型；

其中，所述主网模型为单相模型，由主网高压母线、变压器和主网低压母线组成，所述主网高压母线的出线依次连接所述变压器和主网低压母线；

所述配网模型为三相模型，由馈线组成，所述馈线为所述主网低压母线的出线。

优选的，按下式确定所述主配网一体化状态估计模型的量测方程Z：

Z＝H(X_I,X_B)+e (1)

式(1)中，Z为主配网全网的量测向量，H为量测函数，e为量测误差向量，X_I为除边界节点之外的节点状态向量，X_B为边界节点的状态量。

优选的，所述步骤2包括：

将主配网一体化状态估计模型中配网模型的馈线与低压母线之间添加零阻抗支路，并在所述馈线侧添加所述馈线对应的虚拟母线，由所述馈线及所述馈线对应的虚拟母线构成所述子系统；

其中，所述零阻抗支路为所述子系统之间的联络线，按下式确定解耦后所述子系统的内部状态量与解耦前主配网的关系：

$X_I^a=X_I---\left(2\right)$

式(2)中，为解耦后各子系统内部的状态向量；X_I为解耦前主配网的内部变量；

优选的，所述步骤3包括：选择主网中非边界母线作为解耦前主配网全系统及解耦后主网的相角参考点，选择虚拟母线作为解耦后所述子系统的相角参考点，令解耦前主配网节点数为N，边界母线数为M，总馈线数为K，并按下式确定解耦前所述主配网的状态量Y与解耦后所述子系统中状态量X的线性关系：

$Y=\left[\begin{array}{c}{Y}_I\\ Y_B\end{array}\right]=\left[\begin{array}{cc}I& B_I\\ 0& B_B\end{array}\right]\left[\begin{array}{c}{X}_I\\ X_B\end{array}\right]=CX---\left(3\right)$

式(3)中，C为常数；X_B为解耦前主配网边界母线的2M维状态列向量，Y_I为解耦后子系统内部状态量2(N-M)-1维度列向量，N为解耦前主配网节点数，M为解耦前主配网边界母线数，Y_B为解耦后子系统边界状态量组成的4M维列向量，I为单位矩阵，B_I为除对应边界母线相角元素为1外均为零的矩阵，B_B为除对应边界母线幅值的位置为1外均为零的矩阵。

优选的，所述步骤4包括：

确定解耦后主配网的量侧方程z_k，公式为：

z_k＝f_k(y_i,k,y_b,k)+e_k (4)

式(4)中，z_k为子系统k的量测向量，f_k为子系统k的量测函数，e_k为子系统k的量测误差，y_i,k为主配网解耦子系统k的内部各节点变量，i为节点编号；y_b,k为主配网解耦子系统k的边界节点变量；

其中，基于加权最小二乘状态估计方程，式(4)的求解公式为：

$G_{F,k}^{\left(j\right)}{\mathrm{\Δy}}_k^{\left(j\right)}={\[F_k^{\left(j\right)}\]}^T{W}_k\[z_j-f(y_{i,k},y_{b,k})\]---\left(5\right)$

式(5)中，是子系统k的雅可比矩阵；j是迭代次数；子系统k的信息矩阵；为子系统k第j次迭代时状态量的修正量；W_k为加权矩阵；y_i,k为子系统k的内部节点状态向量；y_b,k为子系统k的边界节点状态向量；z_j为第j次迭代时的量测向量。

优选的，所述步骤5包括：

考虑所述子系统状态量的估计误差，按下式确定解耦前所述主配网的状态量Y与解耦后所述子系统中状态量X的线性关系：

Y＝CX+E (6)

式(6)中，E为子系统状态量的误差向量。

进一步的，当不考虑约束时，基于加权最小二乘算法，获取解耦后所述子系统中状态量X，公式为：

$\left(C^T{W}_{Y}C\right)\hat{x}=C^T{W}_Y\hat{y}---\left(7\right)$

式(7)中，是X的估计值，W_Y为Y的量测权矩阵，为Y的估计值，W_Y为的Y量测权矩阵；

其中，Y量测权矩阵W_Y的公式为：

$W_Y={\mathrm{cov}}^{-1}\left(\tilde{y}\right)=G_F---\left(8\right)$

式(8)中，解耦后子系统状态量的估计值，G_F为根据计算出的信息矩阵；

将G_F表示为内部变量与外部变量的形式，公式为：

$G_F=\left[\begin{array}{cc}{G}_{ii}& G_{ib}\\ G_{ib}^T& G_{bb}\end{array}\right]---\left(9\right)$

式(9)中，G_ii、G_ib、及G_bb均为G_F中各元素；

将式(7)进一步表示为：

$\left[\begin{array}{cc}I& \\ B_i^T& B_b^T\end{array}\right]\left[\begin{array}{cc}{G}_{ii}& G_{ib}\\ G_{ib}^T& G_{bb}\end{array}\right]\left[\begin{array}{cc}I& B_i\\ & B_b\end{array}\right]\left[\begin{array}{c}{\hat{x}}_i\\ {\hat{x}}_b\end{array}\right]=\left[\begin{array}{cc}I& \\ B_b^T& B_b^T\end{array}\right]\left[\begin{array}{cc}{G}_{ii}& G_{ib}\\ G_{ib}^T& G_{bb}\end{array}\right]\left[\begin{array}{c}{\hat{y}}_i\\ {\hat{y}}_b\end{array}\right]---\left(10\right)$

式(10)中，B_i与B_b为C中的元素,I为单位阵，为第二阶段状态估计对各子系统内部状态量的估计值，为第二阶段状态估计对各子系统边界状态量的估计值，为第一阶段状态估计对各子系统内部状态量的估计值，为第一阶段状态估计对各子系统边界状态量的估计值。

进一步的，当考虑零注入约束时，求解零注入约束，公式为：

$\left\{\begin{array}{c}\min f\left(x\right)={(Y-CX)}^T{W}_Y(Y-CX)\\ \begin{array}{cc}s.t.& p\left(x\right)=0\end{array}\end{array}\right.---\left(11\right)$

式(11)中，p(x)为边界母线的零注入约束；

添加拉格朗日乘子处理式(11)，则有：

$\min L\left(Y\right)={(\tilde{y}-Cx)}^T{W}_Y(\tilde{y}-Cx)+{\mathrm{\λ}}^{T}p\left(x\right)---\left(12\right)$

对式(12)中的Y与λ分别求极值得，并对极值条件在处展开，则有：

$\left[\begin{array}{cc}{C}^T{W}_{Y}C& P^T\\ P& 0\end{array}\right]\left[\begin{array}{c}\mathrm{\Δ}x\\ \mathrm{\λ}\end{array}\right]=\left[\begin{array}{c}0\\ -p\left(\hat{x}\right)\end{array}\right]---\left(13\right)$

式(13)中，x为待估计的状态量，λ为拉格朗日乘子，P为等式约束的雅可比矩阵；

消去式(13)中的Δx，求解λ：

$\mathrm{\λ}={\[P{\left(C^T{W}_{Y}C\right)}^{-1}{P}^T\]}^{-1}p\left(\hat{x}\right)---\left(14\right)$

进而获取Δx：

Δx＝-(C^TW_YC)^-1P^Tλ (15)。

进一步的，当边界母线存在等值的馈线时，注入量测为：

z_p＝p(x)+e_p (16)

式(16)中，z_p为原系统边界母线的等值注入量测，p(x)为量测函数，e_p为量测误差；

构造拉格朗日函数，公式为：

$\min L\left(Y\right)={(\tilde{y}-Cx)}^T{W}_Y(\tilde{y}-Cx)+\mathrm{\λ}{(z_p-p(x\left)\right)}^T{W}_Y(z_p-p(x\left)\right)---\left(17\right)$

式(17)中，L(Y)为目标函数，x为待估计的状态量，λ为拉格朗日乘子；

在x处进行泰勒公式展开，公式为：

$\left[\begin{array}{cc}{C}^T{W}_{Y}C& P^T\\ P& -W_P^{-1}\end{array}\right]\left[\begin{array}{c}\mathrm{\Δ}x\\ \mathrm{\λ}\end{array}\right]=\left[\begin{array}{c}0\\ z_p-p\left(\hat{x}\right)\end{array}\right]---\left(18\right)$

式(18)中，Δx为状态量的修正量，P为等式约束的雅可比矩阵，W_P为边界量测的量测权矩阵，为边界量测的量测函数；

消除Δx，得到λ的求解公式：

$\[P{\left(C^T{W}_{Y}C\right)}^{-1}{P}^T+W_p^{-1}\]\mathrm{\λ}=p\left(\hat{x}\right)-z_p---\left(19\right)$

根据求得的λ计算Δx：

Δx＝-(C^TW_YC)^-1P^Tλ (20)

从而得到满足式的解为：

$x=\hat{x}+\mathrm{\Δ}x---\left(21\right).$

从上述的技术方案可以看出，本发明提供了一种基于两阶段法的主配网一体化分布式状态估计方法，建立主配网一体化状态估计模型；在主配网一体化状态估计模型中添加零阻抗支路及虚拟母线；根据零阻抗支路，将主配网解耦为多个子系统；建立主配网的状态量与各子系统中状态量的线性关系；进行第一阶段估计及第二阶段估计，得到主配网的一体化状态。本发明提出的方法降低了计算量，有效且准确地解决了由于配网网络规模庞大造成的主配网一体化状态估计维数灾的问题，而且能够得到主配网一致的、精确的实时状态，为主配网一体化的分析与计算提供基础数据，进而提高了主配网一体化的运行可靠性。

本发明的有益效果：

1、本发明所提供的技术方案中，降低了计算量，有效且准确地解决了由于配网网络规模庞大造成的主配网一体化状态估计维数灾的问题。

2、本发明所提供的技术方案，主网、配网各馈线可以分别采用单相模型、三相模型；能够得到主配网一致的、精确的实时状态，为主配网一体化的分析与计算提供基础数据。

3、本发明所提供的技术方案，坏数据处理在各子系统内部处理完成；提高了主配网一体化的运行可靠性。

4、本发明提供的技术方案，应用广泛，具有显著的社会效益和经济效益。

附图说明

图1是本发明一种基于两阶段法的主配网一体化分布式状态估计方法的流程图；

图2是本发明实施例中解耦前主配网一体化状态估计模型示意图；

图3是本发明实施例中解耦后主配网一体化状态估计模型示意图；

图4是本发明实施例中第二阶段状态估计流程图。

具体实施方式

下面结合附图对本发明的具体实施方式作详细说明。

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。

术语及定义：

1)状态估计：状态估计也称滤波，它是利用实时量测系统的冗余度来提高数据精度，排除随机干扰所引起的错误信息，估计或预报系统的运行状态。电力系统状态估计分为静态和动态状态估计。静态状态估计目前应用较为成熟，以最小二乘法等为主，静态状态估计的局限性是没有考虑系统的动态。

2)分布式状态估计：分布式状态估计利用状态估计算法本身可分区域计算的特性，在各分区进行本区域的状态估计，再通过与其他区域的信息交换，从而使所有分区的状态估计结果达到或接近整体计算的效果，相对于集中式状态估计，不存在因增益矩阵的条件数过大、数值稳定性差等问题导致的整体状态估计失败等问题。

3)主配网：本专利中指的是目前我国电网调度运行管辖权的规定，以220kV、110kV、35kV变电站为界，上述变电站及电压等级高于或等于上述变电站的电网，称为主网，10kV或6.3kV变电站与线路称为配网。

4)零阻抗支路：零阻抗支路是指支路阻抗、对地导纳均为零的支路，其特性为支路两端节点电压幅值与相角均相等的支路，在电力系统零阻抗支路是真实存在，比如开关中的母联开关、同一电厂的新建机组与原有机组高压母线之间的联络线，由于电气距离只有几十米，电气特性上也相当于零阻抗支路。

本发明提供的一种基于两阶段法的主配网一体化分布式状态估计方法，如图1所示，包括：

步骤1.建立主配网一体化状态估计模型；

步骤2.在所述主配网一体化状态估计模型中添加零阻抗支路及虚拟母线，将所述主配网一体化状态估计模型的配网模型解耦为若干子系统；

步骤3.建立解耦前所述主配网的状态量与解耦后所述子系统中状态量的线性关系；

步骤4.第一阶段状态估计：对解耦后子系统以分布式的方式独立进行状态估计；

步骤5.第二阶段状态估计：以解耦前后主配网状态量的线性关系为量测方程，建立状态估计模型进行状态估计，获取主配网的一体化状态。

具体的，如图2所示，所述步骤1中，所述主配网一体化状态估计模型包括：主网模型和配网模型；

其中，所述主网模型为单相模型，由主网高压母线、变压器和主网低压母线组成，所述主网高压母线的出线依次连接所述变压器和主网低压母线；

所述配网模型为三相模型，由馈线组成，所述馈线为所述主网低压母线的出线。

按下式确定所述主配网一体化状态估计模型的量测方程Z：

Z＝H(X_I,X_B)+e (1)

式(1)中，Z为主配网全网的量测向量，H为量测函数，e为量测误差向量，X_I为除边界节点之外的节点状态向量，X_B为边界节点的状态量。

其中，

在图2中的配网各馈线与主网母线之间添加零阻抗支路，并为每条馈线添加虚拟母线，如图3所示，所述步骤2包括：

将主配网一体化状态估计模型中配网模型的馈线与低压母线之间添加零阻抗支路，并在所述馈线侧添加所述馈线对应的虚拟母线，由所述馈线及所述馈线对应的虚拟母线构成所述子系统，其中，所述零阻抗支路为所述子系统之间的联络线；

令添加零阻抗后B1、B2母线的状态量分别为X_B，1与X_B，2，母线B3、B4的状量分别为与则与原状态量的关系如下：

$\left[\begin{array}{c}{X}_{B,1}^a\\ X_{B,2}^a\\ X_{B,3}^a\\ X_{B,4}^a\end{array}\right]=\left[\begin{array}{cc}1& 0\\ 0& 1\\ 1& 0\\ 0& 1\end{array}\right]\left[\begin{array}{c}{X}_{B,1}\\ X_{B,2}\end{array}\right]$

通过增加零阻抗支路与对主网边界母线增加新状态量，实现了主网与配网各馈线之间的解耦，则按下式确定解耦后所述子系统的内部状态量与解耦前主配网的关系：

$X_I^a=X_I---\left(2\right)$

式(2)中，为解耦后各子系统内部的状态向量；X_I为解耦前主配网的内部变量；

通过添加零阻抗支路，增广了边界母线的状态量，将原系统解耦为主网、馈线1与馈线2三个子系统，选择主网中非边界母线为解耦前主配网全系统的相角参考点，并且也作为解耦后主网的相角参考点，即对于主网解耦前后节点相角参考点不变，解耦后选择母线B3、B4分别为馈线1与馈线2子系统的相角参考点，则解耦后子系统状态量与原系统状态量的关系为：

$\left[\begin{array}{c}{U}_{M,I}\\ {\mathrm{\δ}}_{M,I}\\ U_{D,I}^1\\ {\mathrm{\δ}}_{D,I}^1\\ U_{D,I}^2\\ {\mathrm{\δ}}_{D,I}^2\\ U_{M,B}^1\\ {\mathrm{\δ}}_{M,B}^1\\ U_{M,B}^2\\ {\mathrm{\δ}}_{M,B}^2\\ U_{D,B}^1\\ {\mathrm{\δ}}_{D,B}^1\\ U_{D,B}^2\\ {\mathrm{\δ}}_{D,B}^2\end{array}\right]=\left[\begin{array}{cccccccccc}{I}_M& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& I_M& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& I_{F1}& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& I_{F1}& 0& 0& 0& I{1}_{F2}& 0& 0\\ 0& 0& 0& 0& I_{F2}& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& I_{F2}& 0& 0& 0& I{1}_{F2}\\ 0& 0& 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}{V}_{M,I}\\ {\mathrm{\θ}}_{M,I}\\ V_{D,I}^1\\ {\mathrm{\θ}}_{D,I}^1\\ V_{D,I}^2\\ {\mathrm{\θ}}_{D,I}^2\\ V_B^1\\ {\mathrm{\θ}}_B^1\\ V_B^2\\ {\mathrm{\θ}}_B^2\end{array}\right]$

基于上式解耦后子系统状态量与原系统状态量的关系，选择主网中非边界母线作为解耦前主配网全系统及解耦后主网的相角参考点，选择虚拟母线作为解耦后所述子系统的相角参考点，令解耦前主配网节点数为N，边界母线数为M，总馈线数为K，并按下式确定解耦前所述主配网的状态量Y与解耦后所述子系统中状态量X的线性关系：

$Y=\left[\begin{array}{c}{Y}_I\\ Y_B\end{array}\right]=\left[\begin{array}{cc}I& B_I\\ 0& B_B\end{array}\right]\left[\begin{array}{c}{X}_I\\ X_B\end{array}\right]=CX---\left(3\right)$

式(3)中，C为常数；X_B为解耦前主配网边界母线的2M维状态列向量，Y_I为解耦后子系统内部状态量2(N-M)-1维度列向量，N为解耦前主配网节点数，M为解耦前主配网边界母线数，Y_B为解耦后子系统边界状态量组成的4M维列向量，I为单位矩阵，B_I为除对应边界母线相角元素为1外均为零的矩阵，B_B为除对应边界母线幅值的位置为1外均为零的矩阵。

所述步骤4包括：

确定解耦后主配网的量侧方程z_k，公式为：

z_k＝f_k(y_i,k,y_b,k)+e_k (4)

式(4)中，z_k为子系统k的量测向量，f_k为子系统k的量测函数，e_k为子系统k的量测误差，y_i,k为主配网解耦子系统k的内部各节点变量，i为节点编号；y_b,k为主配网解耦子系统k的边界节点变量；

其中，基于加权最小二乘状态估计方程，式(4)的求解公式为：

$G_{F,k}^{\left(j\right)}{\mathrm{\Δy}}_k^{\left(j\right)}={\[F_k^{\left(j\right)}\]}^T{W}_k\[z_j-f(y_{i,k},y_{b,k})\]---\left(5\right)$

式(5)中，是子系统k的雅可比矩阵；j是迭代次数；子系统k的信息矩阵；为子系统k第j次迭代时状态量的修正量；W_k为加权矩阵；y_i,k为子系统k的内部节点状态向量；y_b,k为子系统k的边界节点状态向量；z_j为第j次迭代时的量测向量。

所述步骤5，如图4所示，包括：

考虑所述子系统状态量的估计误差，按下式确定解耦前所述主配网的状态量Y与解耦后所述子系统中状态量X的线性关系：

Y＝CX+E (6)

式(6)中，E为子系统状态量的误差向量。

当不考虑约束时，基于加权最小二乘算法，获取解耦后所述子系统中状态量X，公式为：

$\left(C^T{W}_{Y}C\right)\hat{x}=C^T{W}_Y\hat{y}---\left(7\right)$

式(7)中，是X的估计值，W_Y为Y的量测权矩阵，为Y的估计值，W_Y为的Y量测权矩阵；

其中，Y量测权矩阵W_Y的公式为：

$W_Y={\mathrm{cov}}^{-1}\left(\tilde{y}\right)=G_F---\left(8\right)$

式(8)中，解耦后子系统状态量的估计值，G_F为根据计算出的信息矩阵；

将G_F表示为内部变量与外部变量的形式，公式为：

$G_F=\left[\begin{array}{cc}{G}_{ii}& G_{ib}\\ G_{ib}^T& G_{bb}\end{array}\right]---\left(9\right)$

式(9)中，G_ii、G_ib、及G_bb均为G_F中各元素；

将式(7)进一步表示为：

$\left[\begin{array}{cc}I& \\ B_i^T& B_b^T\end{array}\right]\left[\begin{array}{cc}{G}_{ii}& G_{ib}\\ G_{ib}^T& G_{bb}\end{array}\right]\left[\begin{array}{cc}I& B_i\\ & B_b\end{array}\right]\left[\begin{array}{c}{\hat{x}}_i\\ {\hat{x}}_b\end{array}\right]=\left[\begin{array}{cc}I& \\ B_b^T& B_b^T\end{array}\right]\left[\begin{array}{cc}{G}_{ii}& G_{ib}\\ G_{ib}^T& G_{bb}\end{array}\right]\left[\begin{array}{c}{\hat{y}}_i\\ {\hat{y}}_b\end{array}\right]---\left(10\right)$

式(10)中，B_i与B_b为C中的元素,I为单位阵，为第二阶段状态估计对各子系统内部状态量的估计值，为第二阶段状态估计对各子系统边界状态量的估计值，为第一阶段状态估计对各子系统内部状态量的估计值，为第一阶段状态估计对各子系统边界状态量的估计值。

其中，式(10)的方程维数庞大，求解工作量很大，可以采用分步求解的方式，先求解后求解 的求解公式如下：

$\left(B_b^T\right(G_{bb}-G_{ib}^T{G}_{ii}^{-1}{G}_{ib}\left)\right)B_b^T{\hat{x}}_b=B_b^T(G_{bb}-G_{ib}^T{G}_{ii}^{-1}{G}_{ib}){\hat{y}}_b$

根据上式求解得到后，将带入式(10)后，得：

${\hat{x}}_i={\mathrm{\Δx}}_i+{\hat{y}}_i-B_i{\hat{x}}_b$

其中Δx_i由下式求得：

$G_{ii}{\mathrm{\Δx}}_i=G_{ib}({\hat{y}}_b-B_b{\hat{x}}_b)=\hat{Y}$

由于G_ii是对角分块的形式，则对于Δx_i的求解可以分块进行，即：

$G_{ii,k}{\mathrm{\Δx}}_{i,k}={\hat{Y}}_k$

下标k表示第k个子系统。

当主网与配网的边界当没有对馈线进行等值时，考虑零注入约束，求解零注入约束，公式为：

$\left\{\begin{array}{c}\min f\left(x\right)={(Y-CX)}^T{W}_Y(Y-CX)\\ \begin{array}{cc}s.t.& p\left(x\right)=0\end{array}\end{array}\right.---\left(11\right)$

式(11)中，p(x)为边界母线的零注入约束；

添加拉格朗日乘子处理式(11)，则有：

$\min L\left(Y\right)={(\tilde{y}-Cx)}^T{W}_Y(\tilde{y}-Cx)+{\mathrm{\λ}}^{T}p\left(x\right)---\left(12\right)$

对式(12)中的Y与λ分别求极值得，并对极值条件在处展开，则有：

$\left[\begin{array}{cc}{C}^T{W}_{Y}C& P^T\\ P& 0\end{array}\right]\left[\begin{array}{c}\mathrm{\Δ}x\\ \mathrm{\λ}\end{array}\right]=\left[\begin{array}{c}0\\ -p\left(\hat{x}\right)\end{array}\right]---\left(13\right)$

式(13)中，x为待估计的状态量，λ为拉格朗日乘子，P为等式约束的雅可比矩阵；

消去式(13)中的Δx，求解λ：

$\mathrm{\λ}={\[P{\left(C^T{W}_{Y}C\right)}^{-1}{P}^T\]}^{-1}p\left(\hat{x}\right)---\left(14\right)$

进而获取Δx：

Δx＝-(C^TW_YC)^-1P^Tλ (15)。

当边界母线存在等值的馈线时，考虑馈线等值注入，注入量测为：

z_p＝p(x)+e_p (16)

式(16)中，z_p为原系统边界母线的等值注入量测，p(x)为量测函数，e_p为量测误差；

构造拉格朗日函数，公式为：

$\min L\left(Y\right)={(\tilde{y}-Cx)}^T{W}_Y(\tilde{y}-Cx)+\mathrm{\λ}{(z_p-p(x\left)\right)}^T{W}_Y(z_p-p(x\left)\right)---\left(17\right)$

式(17)中，L(Y)为目标函数，x为待估计的状态量，λ为拉格朗日乘子；

在x处进行泰勒公式展开，公式为：

$\left[\begin{array}{cc}{C}^T{W}_{Y}C& P^T\\ P& -W_P^{-1}\end{array}\right]\left[\begin{array}{c}\mathrm{\Δ}x\\ \mathrm{\λ}\end{array}\right]=\left[\begin{array}{c}0\\ z_p-p\left(\hat{x}\right)\end{array}\right]---\left(18\right)$

式(18)中，Δx为状态量的修正量，P为等式约束的雅可比矩阵，W_P为边界量测的量测权矩阵，为边界量测的量测函数；

消除Δx，得到λ的求解公式：

$\[P{\left(C^T{W}_{Y}C\right)}^{-1}{P}^T+W_p^{-1}\]\mathrm{\λ}=p\left(\hat{x}\right)-z_p---\left(19\right)$

根据求得的λ计算Δx：

Δx＝-(C^TW_YC)^-1P^Tλ (20)

从而得到满足式的解为：

$x=\hat{x}+\mathrm{\Δ}x---\left(21\right).$

最后应当说明的是：以上实施例仅用以说明本发明的技术方案而非对其限制，尽管参照上述实施例对本发明进行了详细的说明，所属领域的普通技术人员应当理解：依然可以对本发明的具体实施方式进行修改或者等同替换，而未脱离本发明精神和范围的任何修改或者等同替换，其均应涵盖在本发明的权利要求保护范围之内。

Claims (10)

1.一种基于两阶段法的主配网一体化分布式状态估计方法，其特征在于，所述方法包括：

步骤1.建立主配网一体化状态估计模型；

步骤2.在所述主配网一体化状态估计模型中添加零阻抗支路及虚拟母线，将所述主配网一体化状态估计模型的配网模型解耦为若干子系统；

步骤3.建立解耦前所述主配网的状态量与解耦后所述子系统中状态量的线性关系；

步骤4.第一阶段状态估计：对解耦后子系统以分布式的方式独立进行状态估计；

步骤5.第二阶段状态估计：以解耦前后主配网状态量的线性关系为量测方程，建立状态估计模型进行状态估计，获取主配网的一体化状态。

2.如权利要求1所述的方法，其特征在于，所述步骤1中，所述主配网一体化状态估计模型包括：主网模型和配网模型；

其中，所述主网模型为单相模型，由主网高压母线、变压器和主网低压母线组成，所述主网高压母线的出线依次连接所述变压器和主网低压母线；

所述配网模型为三相模型，由馈线组成，所述馈线为所述主网低压母线的出线。

3.如权利要求2所述的方法，其特征在于，按下式确定所述主配网一体化状态估计模型的量测方程Z：

Z＝H(X_I,X_B)+e (1)

式(1)中，Z为主配网全网的量测向量，H为量测函数，e为量测误差向量，X_I为除边界节点之外的节点状态向量，X_B为边界节点的状态量。

4.如权利要求1所述的方法，其特征在于，所述步骤2包括：

将主配网一体化状态估计模型中配网模型的馈线与低压母线之间添加零阻抗支路，并在所述馈线侧添加所述馈线对应的虚拟母线，由所述馈线及所述馈线对应的虚拟母线构成所述子系统；

其中，所述零阻抗支路为所述子系统之间的联络线，按下式确定解耦后所述子系统的内部状态量与解耦前主配网的关系：

式(2)中，为解耦后各子系统内部的状态向量；X_I为解耦前主配网的内部变量。

5.如权利要求1所述的方法，其特征在于，所述步骤3包括：选择主网中非边界母线作为解耦前主配网全系统及解耦后主网的相角参考点，选择虚拟母线作为解耦后所述子系统的相角参考点，令解耦前主配网节点数为N，边界母线数为M，总馈线数为K，并按下式确定解耦前所述主配网的状态量Y与解耦后所述子系统中状态量X的线性关系：

式(3)中，C为常数；X_B为解耦前主配网边界母线的2M维状态列向量，Y_I为解耦后子系统内部状态量2(N-M)-1维度列向量，N为解耦前主配网节点数，M为解耦前主配网边界母线数，Y_B为解耦后子系统边界状态量组成的4M维列向量，I为单位矩阵，B_I为除对应边界母线相角元素为1外均为零的矩阵，B_B为除对应边界母线幅值的位置为1外均为零的矩阵。

6.如权利要求1所述的方法，其特征在于，所述步骤4包括：

确定解耦后主配网的量侧方程z_k，公式为：

z_k＝f_k(y_i,k,y_b,k)+e_k (4)

式(4)中，z_k为子系统k的量测向量，f_k为子系统k的量测函数，e_k为子系统k的量测误差，y_i,k为主配网解耦子系统k的内部各节点变量，i为节点编号；y_b,k为主配网解耦子系统k的边界节点变量；

其中，基于加权最小二乘状态估计方程，式(4)的求解公式为：

式(5)中，是子系统k的雅可比矩阵；j是迭代次数；子系统k的信息矩阵；为子系统k第j次迭代时状态量的修正量；W_k为加权矩阵；y_i,k为子系统k的内部节点状态向量；y_b,k为子系统k的边界节点状态向量；z_j为第j次迭代时的量测向量。

7.如权利要求1所述的方法，其特征在于，所述步骤5包括：

考虑所述子系统状态量的估计误差，按下式确定解耦前所述主配网的状态量Y与解耦后所述子系统中状态量X的线性关系：

Y＝CX+E (6)

式(6)中，E为子系统状态量的误差向量。

8.如权利要求7所述的方法，其特征在于，当不考虑约束时，基于加权最小二乘算法，获取解耦后所述子系统中状态量X，公式为：

式(7)中，是X的估计值，W_Y为Y的量测权矩阵，为Y的估计值，W_Y为的Y量测权矩阵；

其中，Y量测权矩阵W_Y的公式为：

式(8)中，解耦后子系统状态量的估计值，G_F为根据计算出的信息矩阵；

将G_F表示为内部变量与外部变量的形式，公式为：

式(9)中，G_ii、G_ib、及G_bb均为G_F中各元素；

将式(7)进一步表示为：

式(10)中，B_i与B_b为C中的元素,I为单位阵，为第二阶段状态估计对各子系统内部状态量的估计值，为第二阶段状态估计对各子系统边界状态量的估计值，为第一阶段状态估计对各子系统内部状态量的估计值，为第一阶段状态估计对各子系统边界状态量的估计值。

9.如权利要求8所述的方法，其特征在于，当主网与配网的边界没有对馈线进行等值时考虑零注入约束，求解零注入约束，公式为：

式(11)中，p(x)为边界母线的零注入约束；

添加拉格朗日乘子处理式(11)，则有：

对式(12)中的Y与λ分别求极值得，并对极值条件在处展开，则有：

式(13)中，x为待估计的状态量，λ为拉格朗日乘子，P为等式约束的雅可比矩阵；

消去式(13)中的Δx，求解λ：

进而获取Δx：

Δx＝-(C^TW_YC)^-1P^Tλ (15)。

10.如权利要求9所述的方法，其特征在于，当边界母线存在等值的馈线时，考虑馈线等值注入，注入量测为：

z_p＝p(x)+e_p (16)

式(16)中，z_p为原系统边界母线的等值注入量测，p(x)为量测函数，e_p为量测误差；

构造拉格朗日函数，公式为：

式(17)中，L(Y)为目标函数，x为待估计的状态量，λ为拉格朗日乘子；

在x处进行泰勒公式展开，公式为：

式(18)中，Δx为状态量的修正量，P为等式约束的雅可比矩阵，W_P为边界量测的量测权矩阵，为边界量测的量测函数；

消除Δx，得到λ的求解公式：

根据求得的λ计算Δx：

Δx＝-(C^TW_YC)^-1P^Tλ (20)

从而得到满足式的解为：