电力系统状态的双曲余弦型最大指数平方抗差估计方法

摘要

本发明公开了属于电力系统调度自动化领域的电力系统状态的双曲余弦型最大指数平方抗差估计方法。该方法是基于抗差性好、计算效率高的双曲余弦型最大指数平方抗差状态估计；包括如下步骤：步骤A.提取电力系统节点注入有功和无功、支路有功和无功功率，以及节点电压幅值参数；以此建立双曲余弦型最大指数平方抗差状态估计模型；步骤B.利用原-对偶内点算法，对所述双曲余弦型最大指数平方抗差状态估计模型求解。本发明有效抑制包括一致性不良数据在内的多个不良数据，显示了良好的抗差性，并有很高的计算效率，得到广泛应用。

说明书

技术领域

本发明属于电力系统调度自动化领域。特别涉及一种电力系统状态的双曲余 弦型最大指数平方抗差估计方法。

背景技术

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

目前，在国内外应用最为广泛的状态估计是加权最小二乘法(Weightedleast squares,WLS)。WLS模型简洁，求解容易，但是其抗差性很差。为了增强抗差性， 一般有两种方法。第一种是在WLS估计之后加入不良数据辨识环节，例如最大正 则化残差检验法(LNR)或估计辨识方法等；另一种是采用抗差状态估计方法。目 前，国内外学者已经提出的抗差状态估计方法(Robuststateestimation)包括 加权最小绝对值估计(Weightedleastabsolutevalue,WLAV)、非二次准则法 (QL、QC等)以合格率最大为目标的状态估计(Maximumnormalmeasurementrate, MNMR)等。但是这些抗差状态估计方法的估计性能仍有待提高。

发明内容

本发明的目的是提出一种电力系统状态的双曲余弦型最大指数平方抗差估 计方法，其特征在于，该方法是基于抗差性好、计算效率高的双曲余弦型最大指 数平方抗差状态估计；包括如下步骤：

步骤A.提取电力系统节点注入有功和无功、支路有功和无功功率，以及节点 电压幅值参数；以此建立双曲余弦型最大指数平方抗差状态估计模型；

步骤B.利用原-对偶内点算法，对所述双曲余弦型最大指数平方抗差状态估 计模型求解。

所述步骤A建立双曲余弦型最大指数平方抗差状态估计模型为： 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包括：

步骤B1：令x为平启动状态变量；选择r⁽⁰⁾＝λ⁽⁰⁾＝π⁽⁰⁾＝0；设置收敛判据ε＝10^-6， 置迭代计数器k＝0；

步骤B2：求解修正方程，得到[dx^Tdr^Tdλ^Tdπ^T]；

步骤B3，修正变量$\left(\begin{array}{c}{x}^{(k+1)}\\ r^{(k+1)}\\ {\lambda }^{(k+1)}\\ {\pi }^{(k+1)}\end{array}\right)=\left(\begin{array}{c}{x}^{\left(k\right)}\\ r^{\left(k\right)}\\ {\lambda }^{\left(k\right)}\\ {\pi }^{\left(k\right)}\end{array}\right)+\left(\begin{array}{c}dx\\ dr\\ d\lambda \\ d\pi \end{array}\right);$

步骤B4：判断是否收敛，若max(dx)＜ε,则转步骤B6，否则进入步骤B5；

步骤B5：令迭代计数器k＝k+1，进入步骤B2；

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

所述步骤B2包括：

步骤B21：形成量测方程以及零注入功率约束对应的雅克比矩阵及 形成量测方程以及零注入功率约束对应的海森矩阵▽²h(x)及▽²g(x)；

其中h(x)为状态矢量到量测矢量的映射，即为量测估计值；z＝h(x)+ε为量测 方程，z和x如步骤A所示，ε∈R_m为量测误差向量；h(x)的具体表达式介绍如下： 对于节点i的电压幅值量测，v_i＝v_i；对于注入有功量测， j∈N_i意思为与节点i相联的所有节点；G_ij,B_ij分别为 节点导纳矩阵第i行第j列的实部和虚步，v_i和v_j分别为节点i和节点j的电压 幅值；θ_ij为节点i和节点j的相角差；对于注入无功量测， 对于支路ij的有功量测， 其中_gsi为支路ij的首端对地电导，其中g_ij和b_ij分别为支路ij的串联电导和串联电纳；对于支路ij的无功量测， $Q_{ij}=-v_i^2(b_{si}+b_{ij})+v_i{v}_j{b}_{ij}{\mathrm{cos\theta }}_{ij}-v_i{v}_j{g}_{ij}{\mathrm{sin\theta }}_{ij},$其中b_si为支路ij的首端对地电纳。

g(x)＝0为零注入功率约束；

步骤B22：引入拉格朗日函数

$L\equiv \underset{i=1}{\overset{m}{\Sigma }}{w}_i\exp (-\frac{r_i^2}{2{\sigma }_0^2})[\exp \left(\frac{r_i}{{\sigma }_1}\right)+\exp (-\frac{r_i}{{\sigma }_1})]-{\lambda }^{T}g\left(x\right)-{\pi }^T(z-h(x)-r)---\left(4\right)$

式中：λ∈R^c及π∈R^m为拉格朗日乘子矢量；计算

$L_{r_i}=w_i\exp (-\frac{r_i^2}{2{\sigma }_0^2})(-\frac{r_i}{{\sigma }_0^2})[\exp \left(\frac{r_i}{{\sigma }_1}\right)+\exp (-\frac{r_i}{{\sigma }_1})]+\frac{w_i}{{\sigma }_1}\exp (-\frac{r_i^2}{2{\sigma }_0^2})[\exp \left(\frac{r_i}{{\sigma }_1}\right)-\exp (-\frac{r_i}{{\sigma }_1})]+{\pi }_i,$

式中，W为R^m×m的对角阵，其对角元素为

$\begin{array}{c}{W}_{ii}=-\frac{w_i}{{\sigma }_0^2}\left\{\exp (-\frac{r_i^2}{2{\sigma }_0^2})(-\frac{r_i}{{\sigma }_0^2})r_i\right[\exp \left(\frac{r_i}{{\sigma }_1}\right)+\exp (-\frac{r_i}{{\sigma }_1})]+\exp (-\frac{r_i^2}{2{\sigma }_0^2})[\exp \left(\frac{r_i}{{\sigma }_1}\right)+\exp (-\frac{r_i}{{\sigma }_1})\left]\right\}\\ +(-\frac{w_i}{{\sigma }_0^2})\exp (-\frac{r_i^2}{2{\sigma }_0^2})\frac{r_i}{{\sigma }_1}[\exp \left(\frac{r_i}{{\sigma }_1}\right)-\exp (-\frac{r_i}{{\sigma }_1})]\\ +\frac{w_i}{{\sigma }_1}\left\{\exp (-\frac{r_i^2}{2{\sigma }_0^2})(-\frac{r_i}{{\sigma }_0^2})\right[\exp \left(\frac{r_i}{{\sigma }_1}\right)-\exp (-\frac{r_i}{{\sigma }_1})]+\exp (-\frac{r_i^2}{2{\sigma }_0^2})\frac{1}{{\sigma }_1}[\exp \left(\frac{r_i}{{\sigma }_1}\right)+\exp (-\frac{r_i}{{\sigma }_1})\left]\right\}\end{array};$

步骤B23：求解方程$\left(\begin{array}{cccc}{▿}^{2}g\left(x\right)\lambda -{▿}^{2}h\left(x\right)\pi & 0& G^T& -H^T\\ G& 0& 0& 0\\ H& I& 0& 0\\ 0& W& 0& I\end{array}\right)\left(\begin{array}{c}dx\\ dr\\ d\lambda \\ d\pi \end{array}\right)=\left(\begin{array}{c}-L_x\\ -L_{\lambda }\\ -L_{\pi }\\ -L_r\end{array}\right)$得到 [dx^Tdr^Tdλ^Tdπ^T]。

本发明的有益效果对电力系统状态采用双曲余弦型最大指数平方抗差状态 估计，在估计过程中可有效抑制包括一致性不良数据在内的多个不良数据，显示 了良好的抗差性，并具有很高的计算效率，非常适宜于实际工程应用。

具体实施方式

本发明提出一种电力系统状态的双曲余弦型最大指数平方抗差状态估计方 法，下面结合实施例详细描述本发明。

本发明是基于抗差性好、计算效率高的双曲余弦型最大指数平方抗差状态估 计(Hyperboliccosinemaximumexponentialsquarestateestimation, COSH-MES)方法，包括下列步骤：

步骤A：提取电力系统节点注入有功和无功、支路有功和无功功率，以及节 点电压幅值参数；以此建立双曲余弦型最大指数平方抗差状态估计模型；具体地， 本发明提出的COSH-MES的模型如下所示

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：利用原-对偶内点算法，对所述双曲余弦型最大指数平方抗差状态估 计模型求解。

(1)COSH-MES模型的求解方法

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

引入拉格朗日函数

$L\equiv \underset{i=1}{\overset{m}{\Sigma }}{w}_i\exp (-\frac{r_i^2}{2{\sigma }_0^2})[\exp \left(\frac{r_i}{{\sigma }_1}\right)+\exp (-\frac{r_i}{{\sigma }_1})]-{\lambda }^{T}g\left(x\right)-{\pi }^T(z-h(x)-r)---\left(4\right)$

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

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

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

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

$L_{\pi }\equiv \partial L/\partial \pi =z-h\left(x\right)-r=0---\left(7\right)$

$L_{r_i}\equiv \partial L/\partial r_i=w_i\exp (-\frac{r_i^2}{2{\sigma }_0^2})(-\frac{r_i}{{\sigma }_0^2})[\exp \left(\frac{r_i}{{\sigma }_1}\right)+\exp (-\frac{r_i}{{\sigma }_1})]+\frac{w_i}{{\sigma }_1}\exp (-\frac{r_i^2}{2{\sigma }_0^2})[\exp \left(\frac{r_i}{{\sigma }_1}\right)-\exp (-\frac{r_i}{{\sigma }_1})]+{\pi }_i=0---\left(8\right)$

式中：$H=\partial h\left(x\right)/\partial x,G=\partial g\left(x\right)/\partial x.$

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

[▽²g(x)λ-▽²h(x)π]dx+G^Tdλ-H^Tdπ＝-L_x(9)

Gdx＝-L_λ(10)

-Hdx-dr＝-L_π(11)

$\left(\begin{array}{c}-\frac{w_i}{{\sigma }_0^2}\left\{\exp (-\frac{r_i^2}{2{\sigma }_0^2})(-\frac{r_i}{{\sigma }_0^2})r_i\right[\exp \left(\frac{r_i}{{\sigma }_1}\right)+\exp (-\frac{r_i}{{\sigma }_1})]+\exp (-\frac{r_i^2}{2{\sigma }_0^2})[\exp \left(\frac{r_i}{{\sigma }_1}\right)+\exp (-\frac{r_i}{{\sigma }_1})\left]\right\}\\ +(-\frac{w_i}{{\sigma }_0^2})\exp (-\frac{r_i^2}{2{\sigma }_0^2})\frac{r_i}{{\sigma }_1}[\exp \left(\frac{r_i}{{\sigma }_1}\right)-\exp (-\frac{r_i}{{\sigma }_1})]\\ +\frac{w_i}{{\sigma }_1}\left\{\exp (-\frac{r_i^2}{2{\sigma }_0^2})(-\frac{r_i}{{\sigma }_0^2})\right[\exp \left(\frac{r_i}{{\sigma }_1}\right)-\exp (-\frac{r_i}{{\sigma }_1})]+\exp (-\frac{r_i^2}{2{\sigma }_0^2})\frac{1}{{\sigma }_1}[\exp \left(\frac{r_i}{{\sigma }_1}\right)+\exp (-\frac{r_i}{{\sigma }_1})\left]\right\}\end{array}\right){\mathrm{dr}}_i+{\mathrm{d\pi }}_i=-L_{r_i}---\left(12\right)$

式(12)可表示为矩阵形式，为

Wdr+dπ＝-L_r(13)

式中：W为R^m×m的对角阵，其对角元素为

$\begin{array}{c}{W}_{ii}=-\frac{w_i}{{\sigma }_0^2}\left\{\exp (-\frac{r_i^2}{2{\sigma }_0^2})(-\frac{r_i}{{\sigma }_0^2})r_i\right[\exp \left(\frac{r_i}{{\sigma }_1}\right)+\exp (-\frac{r_i}{{\sigma }_1})]+\exp (-\frac{r_i^2}{2{\sigma }_0^2})[\exp \left(\frac{r_i}{{\sigma }_1}\right)+\exp (-\frac{r_i}{{\sigma }_1})\left]\right\}\\ +(-\frac{w_i}{{\sigma }_0^2})\exp (-\frac{r_i^2}{2{\sigma }_0^2})\frac{r_i}{{\sigma }_1}[\exp \left(\frac{r_i}{{\sigma }_1}\right)-\exp (-\frac{r_i}{{\sigma }_1})]\\ +\frac{w_i}{{\sigma }_1}\left\{\exp (-\frac{r_i^2}{2{\sigma }_0^2})(-\frac{r_i}{{\sigma }_0^2})\right[\exp \left(\frac{r_i}{{\sigma }_1}\right)-\exp (-\frac{r_i}{{\sigma }_1})]+\exp (-\frac{r_i^2}{2{\sigma }_0^2})\frac{1}{{\sigma }_1}[\exp \left(\frac{r_i}{{\sigma }_1}\right)+\exp (-\frac{r_i}{{\sigma }_1})\left]\right\}\end{array};$$\begin{array}{cc}d\pi ={\left(\begin{array}{cccc}{\mathrm{d\pi }}_1^T& {\mathrm{d\pi }}_2^T& ...& {\mathrm{d\pi }}_m^T\end{array}\right)}^T,& L_r={\left(\begin{array}{cccc}{L}_{r_1}^T& L_{r_1}^T& ...& L_{r_m}^T\end{array}\right)}^T\end{array}.$

根据式(9)、(10)、(11)及(13)，可得修正方程为

$\left(\begin{array}{cccc}{▿}^{2}g\left(x\right)\lambda -{▿}^{2}h\left(x\right)\pi & 0& G^T& -H^T\\ G& 0& 0& 0\\ H& I& 0& 0\\ 0& W& 0& I\end{array}\right)\left(\begin{array}{c}dx\\ dr\\ d\lambda \\ d\pi \end{array}\right)=\left(\begin{array}{c}-L_x\\ -L_{\lambda }\\ -L_{\pi }\\ -L_r\end{array}\right)---\left(14\right)$

求解式(14)可得[dx^Tdr^Tdλ^Tdπ^T]，则迭代即可持续进行。

(2)COSH-MES模型的求解步骤

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

步骤B1:进行初始化，令x为平启动状态变量；选择r⁽⁰⁾＝λ⁽⁰⁾＝π⁽⁰⁾＝0；设置 收敛判据ε＝10^-6，置迭代计数器k＝0。

具体地，令x⁽⁰⁾∈Rⁿ代表由所有节点电压幅值和相角组成的的平启动状态变量 (参考节点相角除外)；选择r⁽⁰⁾＝λ⁽⁰⁾＝π⁽⁰⁾＝0，其中λ∈R^c及π∈R^m为拉格朗日乘子矢 量，m为量测量的个数，而c为零注入功率约束的个数；置迭代计数器k＝0。

步骤B2：求解式(14)修正方程，以完成对变量的修正，得到[dx^Tdr^Tdλ^Tdπ^T]。

步骤B3：修正变量为：$\left(\begin{array}{c}{x}^{(k+1)}\\ r^{(k+1)}\\ {\lambda }^{(k+1)}\\ {\pi }^{(k+1)}\end{array}\right)=\left(\begin{array}{c}{x}^{\left(k\right)}\\ r^{\left(k\right)}\\ {\lambda }^{\left(k\right)}\\ {\pi }^{\left(k\right)}\end{array}\right)+\left(\begin{array}{c}dx\\ dr\\ d\lambda \\ d\pi \end{array}\right);$

步骤B4：判断是否收敛，若max(dx)＜ε,则转步骤B6，否则进入步骤B5；

步骤B5：令迭代计数器k＝k+1，进入步骤B2；以及

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

实施例

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

1.抗差性能的比较

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

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

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

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

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

2.计算效率的比较

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

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

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

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