CN103413053A - 一种基于内点法的电力系统抗差状态估计方法 - Google Patents

Abstract

本发明涉及一种基于内点法的电力系统抗差状态估计方法，属于电力系统运行和控制技术领域。该方法包括：获取量测数据、拓扑分析以及计算电网参数；建立指数函数抗差状态估计模型；引入残差变量，对抗差状态估计模型进行等价转换，规范成内点法容易求解的形式；利用内点法对等价转换后的状态估计模型进行求解，并对海森矩阵进行数值近似。本发明通过对状态估计模型的等价转换，降低了电力系统抗差状态估计的难度；通过对海森矩阵进行近似，大大提高了电力系统抗差状态估计的计算效率；状态估计结果严格满足零注入等式约束；计算方法具有很强的鲁棒性。

Description

一种基于内点法的电力系统抗差状态估计方法

技术领域

本发明涉及一种基于内点法的电力系统抗差状态估计方法，属于电力系统运行和控制技术领域。

背景技术

电力系统的状态估计就是以数据采集与监视控制系统采集的实时量测为数据输入源，剔除其中的不良数据之后，结合电网模型，按照特定的状态估计模型，对电网的状态量进行最优估算。状态估计在能量管理系统中处于极为重要的地位，是后续一系列高级应用的重要基石。

传统的电力系统状态估计需要进行不良数据的辨识以及剔除，才能进行准确的估计。但现在希望能在进行状态估计的过程中自动降低不良数据所导致的估计偏差，因此可以采用抗差状态估计的模型。

在求解的算法方面，传统的牛顿拉夫逊方法存在收敛域小、对初值要求苛刻、无法处理不等式约束等缺点，在现场应用中，在重载的情况下，牛顿法有可能出现即便平启动也无法收敛的情况，同时也无法考虑节点负荷上下限等等的不等式约束，因此必须寻求更优的解法。而现代内点法是公认的优秀的算法。

发明内容

本发明的目的是提出一种基于内点法的电力系统抗差状态估计方法，将现代内点法应用到抗差状态估计中，以解决实际工程中状态估计收敛域窄的问题。

本发明提出一种基于内点法的电力系统抗差状态估计方法，包括以下步骤：

(1)从电力系统的数据采集中心实时获取量测数据，量测数据包括节点电压幅值、节点有功功率、节点无功功率、支路有功功率和支路无功功率，并存储量测数据，根据量测数据进行拓扑分析，得到电力系统结构图，根据量测数据计算电力系统导纳矩阵Y；

(2)根据步骤(1)中量测数据、电力系统结构图和电力系统导纳矩阵Y，建立初始电力系统抗差状态估计模型如下：

$\underset{s}{\min }J\left(s\right)=-\underset{i=1}{\overset{n_{\mathrm{meas}}}{\mathrm{\Σ}}}\exp (-\frac{{(m_i\left(s\right)-m_i^{\mathrm{meas}})}^2}{2{\mathrm{\σ}}^2})$

s.t.  c(s)=O

其中，s为电力系统状态量，s=[θ，V]^T，m_i(s)为电力系统第i个量测方程的估计值，为上述量测数据中的第i个量测值，σ是核估计法的窗宽，窗宽σ的取值范围为：O<σ<1，c(s)=O表示零注入等式约束；

(3)对上述初始电力系统抗差状态估计模型进行等价转换，得到电力系统抗差状态估计模型如下：

$\min f(s,t)=-\underset{i}{\mathrm{\&Sigma;}}{e}^{-\frac{{t_i}^2}{2{\mathrm{\&sigma;}}^2}}$

s.t.c(s，t)=O

$g_i(s,t)=m_i\left(s\right)-m_i^{\mathrm{meas}}-t_i=0,i=\mathrm{1,2},...,n_{\mathrm{meas}}$

其中，m_i(s)为电力系统第i个量测方程的估计值，t为残差变量，残差变量中的第i个分量为：

i=1，2，...，n_meas，n_meas为上述量测数据的个数，s为电力系统状态量，s=[θ，V]^T，夕是节点电压的相角，V是节点电压的幅值，T是矩阵转置；

(4)采用原-对偶内点法，对上述电力系统抗差状态估计模型进行求解，具体实现过程如下：

(4-1)将上述电力系统抗差状态估计模型改写为以下标准形式：

min  f(x)

s.f  h(x)=O

$\underset{\&OverBar;}{g}\&le;g\left(x\right)\&le;\stackrel{\&OverBar;}{g}$

其中，x是电力系统状态量s和残差变量t合并后的复合向量，x=[s，t]^T，函数h(x)包含量测方程m_i(s)的约束以及零注入约束，g(x)为电力系统节点负荷与复合向量x之间的函数关系，为用户设定的电力系统负荷上限，g为用户设定的电力系统负荷下限；

(4-2)将上述标准形式中的不等式约束简化为等式约束，得到简化后的电力系统抗差状态估计模型：

min  f(x)

s.t  h(x)=O

g(x)-1-g＝0

$g\left(x\right)+u-\stackrel{\&OverBar;}{g}=0$

(1，u)≥O

其中，(1，u)∈R^r为将不等式约束转换为等式约束的松弛向量，R为实数集，r为函数g(x)的维数；

(4-3)根据上述简化后的电力系统抗差状态估计模型，构建一个增广拉格朗日函数如下：

$L(x,l,u,y,z,w,\tilde{z},\tilde{w})\&equiv;f\left(x\right)-y^{T}h\left(x\right)-z^T(g\left(x\right)-l-\underset{\&OverBar;}{g})$

$-w^T(g\left(x\right)+u-\stackrel{\&OverBar;}{g})-{\tilde{z}}^{T}l{-\tilde{w}}^{T}u$

其中，x为电力系统状态量s与残差变量t合并后的复合向量，对应于原-对偶内点法计算中的原始变量，(l，u)∈R^r为松弛向量，r为函数g(x)的维数，y，z，w，

分别为原-对偶内点法计算中的对偶变量；

(4-4)设置迭代初始值，将原始变量x中，电力系统状态量s中的电压幅值初值设为1，电压相角初值设为0，将原始变量x中的残差变量t的初始值设为0，将对偶变量的初始值设为0，设置迭代次数的初始值k=0，设置最大迭代次数K_max，K_max的取值范围为：10～100，对于电力系统中的零注入母线，设定残差变量的第i个分量

(4-5)对迭代次数k进行判断，若k＜K_max，则进行步骤(4-6)，若k≥K_max，则结束计算，并输出计算不收敛；

(4-6)根据下式计算简化后的电力系统抗差状态估计模型的互补间隙C_Gap：

$C_{\mathrm{Gap}}=\underset{i=1}{\overset{r}{\mathrm{\&Sigma;}}}(l_i{z}_i-u_i{w}_i)$

其中，l和u为松弛变量，z和w为对偶变量，r为函数g(x)的维数；

设定一个电力系统抗差状态估计计算精度ε，一般可设为10^-6，将互补间隙C_Gap与计算精度ε进行比较，若C_Gap＜ε，则输出x，包括电力系统状态量和残差变量，并结束计算，若C_Gap≥ε，则进行步骤(4-7)；

(4-7)根据下式计算扰动因子μ：

$\mathrm{\&mu;}\&equiv;\mathrm{\&tau;}\frac{C_{\mathrm{Gap}}}{2r}$

其中τ为中心参数，取值范围为τ∈(0，1]，r为函数g(x)的维数，C_Gap为上述互补间隙；

(4-8)按照以下卡罗需-库恩-塔克方程组，计算上述增广拉格朗日函数对原-对偶内点法中各原始变量和对偶变量的一阶偏导数：

$({\&dtri;}_x^{2}f\left(x\right)-{\&dtri;}_x^{2}h\left(x\right)y-{\&dtri;}_x^{2}g\left(x\right)(z+w))\mathrm{\&Delta;x}$

$-{\&dtri;}_{x}h\left(x\right)\mathrm{\&Delta;y}-{\&dtri;}_{x}g\left(x\right)(\mathrm{\&Delta;z}+\mathrm{\&Delta;w})=-{\&dtri;}_x{L}_0$

${\&dtri;}_{x}h{\left(x\right)}^T\mathrm{\&Delta;x}=-{\&dtri;}_y{L}_0$

${\&dtri;}_{x}g{\left(x\right)}^T\mathrm{\&Delta;x}-\mathrm{\&Delta;l}=-{\&dtri;}_z{L}_0$

${\&dtri;}_{x}g{\left(x\right)}^T\mathrm{\&Delta;x}+\mathrm{\&Delta;u}=-{\&dtri;}_w{L}_0$

$\mathrm{Z\&Delta;l}+\mathrm{L\&Delta;z}=-{\&dtri;}_l^{\mathrm{\&mu;}}{L}_0$

$\mathrm{W\&Delta;u}+\mathrm{U\&Delta;w}=-{\&dtri;}_u^{\mathrm{\&mu;}}{L}_0$

其中，

分别是卡罗需-库恩-塔克方程组对应的扰动方程的残差，

和分别是步骤(4-1)电力系统抗差状态估计模型标准形式中的f(x)、h(x)和g(x)的海森矩阵；

(4-9)根据以下修正方程，计算得到在第k次迭代时原始变量x的修正量Δx、对偶变量y、z、w的修正量Δy、Δz、Δw和松弛变量1、u的修正量Δl、Δu：

$\left\{\begin{array}{c}\mathrm{\&Delta;l}=\&dtri;g{\left(x\right)}^T\mathrm{\&Delta;x}+L_{z0}\\ \mathrm{\&Delta;u}=-(\&dtri;g{\left(x\right)}^T+L_{w0})\end{array}\right.$

$\left\{\begin{array}{c}\mathrm{\&Delta;z}=-L^{-1}Z\&dtri;g{\left(x\right)}^T\mathrm{\&Delta;x}-L^{-1}({\mathrm{ZL}}_{z0}+L_{l0}^{\mathrm{\&mu;}})\\ \mathrm{\&Delta;w}=U^{-1}W\&dtri;g{\left(x\right)}^T\mathrm{\&Delta;x}+U^{-1}({\mathrm{WL}}_{w0}-L_{l0}^{\mathrm{\&mu;}})\end{array}\right.$

其中：

$H(\&CenterDot;)\&equiv;{\&dtri;}_{x}g\left(x\right)(U^{-1}W-L^{-1}Z){\&dtri;}_{x}g{\left(x\right)}^T+$

$(-{\&dtri;}_x^{2}f\left(x\right)+{\&dtri;}_x^{2}h\left(x\right)y+{\&dtri;}_x^{2}g\left(x\right)(z+w))=H_g+H_h$

$J\left(x\right)\&equiv;{\&dtri;}_{x}h{\left(x\right)}^T$

$\mathrm{\&psi;}(g,\mathrm{\&mu;})\&equiv;-{\&dtri;}_{x}f\left(x\right)+{\&dtri;}_{x}h\left(x\right)y-{\&dtri;}_{x}g\left(x\right)((U^{-1}-L^{-1})\mathrm{\&mu;e}$

$+L^{-1}Z(g\left(x\right)-l-\underset{\&OverBar;}{g})-U^{-1}W(g\left(x\right)+u-\stackrel{\&OverBar;}{g}))$

设定一个转换阈值ε_switch，一般设为计算精度ε的2～3倍即可，ε_switch的取值范围为2×10^-6～3×10^-6，对原始-对偶互补间隙C_Gap进行判断，若C_Gap＜ε_switch，则由于x=[s，t]^T，

${\&dtri;}_x^{2}f\left(x\right)={[{\&dtri;}_s^{2}f(s,t),{\&dtri;}_t^{2}f(s,t)]}^T,$ 在计算上述修正方程中的H(·)时，在

中，

各元素采用下式精确表达：

${\&dtri;}_{t_i}^{2}f(s,t)=\frac{e^{-\frac{1}{2}\frac{t_i^2}{{\mathrm{\&sigma;}}^2}}}{{\mathrm{\&sigma;}}^2}(1-\frac{t_i^2}{{\mathrm{\&sigma;}}^2}),i=1,2,...,n_{\mathrm{meas}}$

其中n_meas为上述量测数据的个数，t_i为残差变量的第i个分量，σ为核估计法的窗宽；

若C_Gap≥ε，则各元素采用下式的近似表达：

${\&dtri;}_{t_i}^{2}f(s,t)\&ap;\frac{e^{-\frac{1}{2}\frac{t_i^2}{{\mathrm{\&sigma;}}^2}}}{{\mathrm{\&sigma;}}^2},i=\mathrm{1,2},...,n_{\mathrm{meas}}$

(4-10)根据下式计算在第k次迭代时原始变量的修正步长

和对偶变量的修正步长

${\mathrm{Step}}_{P}0.9995\min \left\{\underset{i}{\min }(\frac{-l_i}{\mathrm{\&Delta;}{l}_i}:{\mathrm{\&Delta;l}}_i<0;\frac{-u_i}{{\mathrm{\&Delta;u}}_i}:{\mathrm{\&Delta;u}}_i<0)\right\}$

${\mathrm{Step}}_{D}0.9995\min \left\{\underset{i}{\min }(\frac{-z_i}{\mathrm{\&Delta;}{z}_i}:{\mathrm{\&Delta;z}}_i<0;\frac{-w_i}{{\mathrm{\&Delta;w}}_i}:{\mathrm{\&Delta;w}}_i<0)\right\}$

(4-11)根据步骤(4-10)的修正步长，更新上述原始变量和对偶变量：

使 $\left[\begin{array}{c}x\\ l\\ u\end{array}\right]=\left[\begin{array}{c}x\\ l\\ u\end{array}\right]+\mathrm{ste}{p}_P\left[\begin{array}{c}\mathrm{\&Delta;x}\\ \mathrm{\&Delta;l}\\ \mathrm{\&Delta;u}\end{array}\right],\left[\begin{array}{c}y\\ z\\ w\end{array}\right]=\left[\begin{array}{c}y\\ z\\ w\end{array}\right]+{\mathrm{step}}_D\left[\begin{array}{c}\mathrm{\&Delta;y}\\ \mathrm{\&Delta;z}\\ \mathrm{\&Delta;w}\end{array}\right],k=k+1$

(4-12)重复步骤(4-5)-步骤(4-12)。

本发明提出的基于内点法的电力系统抗差状态估计方法，其优点是：

1、本发明状态估计方法对目标函数的海森矩阵进行数值近似之后，状态估计的计算效率较高，因而能够基本满足状态估计实时性的需求。

2、本发明状态估计方法中采用了内点法，因此收敛域广，即便在电力系统运行负荷过大或出现故障等异常情况下，也能迅速得到状态估计的最终结果。

3、本发明状态估计方法，计算过程中严格满足零注入等式约束，并能考虑多种复杂的不等式约束，因此能够得到更加符合实际运行状态的电力系统状态估计结果。

具体实施方式

本发明提出一种基于内点法的电力系统抗差状态估计方法，包括以下步骤：

(1)从电力系统的数据采集中心实时获取量测数据，量测数据包括节点电压幅值、节点有功功率、节点无功功率、支路有功功率和支路无功功率，并存储量测数据，根据量测数据进行拓扑分析，得到电力系统结构图，根据量测数据计算电力系统导纳矩阵Y；

(2)根据步骤(1)中量测数据、电力系统结构图和电力系统导纳矩阵Y，建立初始电力系统抗差状态估计模型如下：

$\underset{s}{\min }J\left(s\right)=-\underset{i=1}{\overset{n_{\mathrm{meas}}}{\mathrm{\&Sigma;}}}\exp (-\frac{{(m_i\left(s\right)-m_i^{\mathrm{meas}})}^2}{2{\mathrm{\&sigma;}}^2})$

s.t.  c(s)＝0

其中，s为电力系统状态量，s＝[θ，V]^T，m_i(s)为电力系统第i个量测方程的估计值，

为上述量测数据中的第i个量测值，σ是核估计法的窗宽，窗宽σ的取值范围为：0＜σ＜1，c(s)＝0表示零注入等式约束；

(3)对上述初始电力系统抗差状态估计模型进行等价转换，得到电力系统抗差状态估计模型如下：

$\min f(s,t)=-\underset{i}{\mathrm{\&Sigma;}}{e}^{-\frac{t_i^2}{2{\mathrm{\&sigma;}}^2}}$

s.t.c(s，t)＝0

$g_i(s,t)=m_i\left(s\right)=-m_i^{\mathrm{meas}}-t_i=0,i=\mathrm{1,2},...,n_{\mathrm{meas}}$

其中，m_i(s)为电力系统第i个量测方程的估计值，t为残差变量，残差变量中的第i个分量为：i＝1，2，...，n_meas，n_meas为上述量测数据的个数，s为电力系统状态量，s＝[θ，V]^T，θ是节点电压的相角，V是节点电压的幅值，T是矩阵转置；

(4)采用原-对偶内点法，对上述电力系统抗差状态估计模型进行求解，具体实现过程如下：

(4-1)将上述电力系统抗差状态估计模型改写为以下标准形式：

minf(x)

s.t h(x)＝0

$\underset{\&OverBar;}{g}\&le;g\left(x\right)\&le;\stackrel{\&OverBar;}{g}$

其中，x是电力系统状态量s和残差变量t合并后的复合向量，x＝[s，t]^T，函数h(x)包含量测方程m_i(s)的约束以及零注入约束，g(x)为电力系统节点负荷与复合向量x之间的函数关系，

为用户设定的电力系统负荷上限，g为用户设定的电力系统负荷下限；

(4-2)将上述标准形式中的不等式约束简化为等式约束，得到简化后的电力系统抗差状态估计模型：

min f(x)

s.t h(x)＝0

g(x)-l-g＝0

$g\left(x\right)+u-\stackrel{\&OverBar;}{g}=0$

(l，u)≥0

其中，(l，u)∈R^r为将不等式约束转换为等式约束的松弛向量，R为实数集，r为函数g(x)的维数；

(4-3)根据上述简化后的电力系统抗差状态估计模型，构建一个增广拉格朗日函数如下：

$L(x,l,u,y,z,w,\tilde{z},\tilde{w})\&equiv;f\left(x\right)-y^{T}h\left(x\right)-z^T(g\left(x\right)-l-\underset{\&OverBar;}{g})$

$-w^T(g\left(x\right)+u-\stackrel{\&OverBar;}{g})-{\tilde{z}}^{T}l-{\tilde{w}}^{T}u$

其中，x为电力系统状态量s与残差变量t合并后的复合向量，对应于原-对偶内点法计算中的原始变量，(l，u)∈R^r为松弛向量，r为函数g(x)的维数，y，z，w，

分别为原-对偶内点法计算中的对偶变量；

(4-4)设置迭代初始值，将原始变量x中，电力系统状态量s中的电压幅值初值设为1，电压相角初值设为0，将原始变量x中的残差变量t的初始值设为0，将对偶变量的初始值设为0，设置迭代次数初始值k＝0，设置最大迭代次数K_max，K_max的取值范围为：10～100，对于电力系统中的零注入母线，设定残差变量的第i个分量

(4-5)对迭代次数k进行判断，若k＜K_max，则进行步骤(4-6)，若k≥K_max，则结束计算，并输出计算不收敛；

(4-6)根据下式计算简化后的电力系统抗差状态估计模型的互补间隙C_Gap：

$C_{\mathrm{Gap}}=\underset{i=1}{\overset{r}{\mathrm{\&Sigma;}}}(l_i{z}_i-u_i{w}_i)$

其中，l和u为松弛变量，z和w为对偶变量，r为函数g(x)的维数；

设定一个电力系统抗差状态估计计算精度ε，将互补间隙C_Gap与计算精度ε进行比较，若C_Gap＜ε，则输出x，包括电力系统状态量和残差变量，并结束计算，若C_Gap≥ε，则进行步骤(4-7)；

(4-7)根据下式计算扰动因子μ：

$\mathrm{\&mu;}\&equiv;\mathrm{\&tau;}\frac{C_{\mathrm{Gap}}}{2r}$

其中τ为中心参数，取值范围为τ∈(0，1]，r为函数g(x)的维数，C_Gap为上述互补间隙；

(4-8)按照以下卡罗需-库恩-塔克方程组，计算上述增广拉格朗日函数对原-对偶内点法中各原始变量和对偶变量的一阶偏导数：

$({\&dtri;}_x^{2}f\left(x\right)-{\&dtri;}_x^{2}h\left(x\right)y-{\&dtri;}_x^{2}g\left(x\right)(z+w))\mathrm{\&Delta;x}$

$-{\&dtri;}_{x}h\left(x\right)\mathrm{\&Delta;y}-{\&dtri;}_{x}g\left(x\right)(\mathrm{\&Delta;z}+\mathrm{\&Delta;w})=-{\&dtri;}_x{L}_0$

${\&dtri;}_{x}h{\left(x\right)}^T\mathrm{\&Delta;x}=-{\&dtri;}_y{L}_0$

${\&dtri;}_{x}g{\left(x\right)}^T\mathrm{\&Delta;x}-\mathrm{\&Delta;l}=-{\&dtri;}_z{L}_0$

${\&dtri;}_{x}g{\left(x\right)}^T\mathrm{\&Delta;x}+\mathrm{\&Delta;u}=-{\&dtri;}_w{L}_0$

$\mathrm{Z\&Delta;l}+\mathrm{L\&Delta;z}=-{\&dtri;}_l^{\mathrm{\&mu;}}{L}_0$

$\mathrm{W\&Delta;u}+\mathrm{U\&Delta;w}=-{\&dtri;}_u^{\mathrm{\&mu;}}{L}_0$

其中，

分别是卡罗需-库恩-塔克方程组对应的扰动方程的残差，

和分别是步骤(4-1)电力系统抗差状态估计模型标准形式中的f(x)、h(x)和g(x)的海森矩阵；

(4-9)根据以下修正方程，计算得到在第k次迭代时原始变量x的修正量Δx、对偶变量y、z、w的修正量Δy、Δz、Δw和松弛变量l、u的修正量Δl、Δu：

$\left\{\begin{array}{c}\mathrm{\&Delta;}l=\&dtri;g{\left(x\right)}^T\mathrm{\&Delta;x}+L_{z0}\\ \mathrm{\&Delta;u}=-(\&dtri;g{\left(x\right)}^T+L_{w0})\end{array}\right.$

$\left\{\begin{array}{c}\mathrm{\&Delta;z}=-L^{-1}Z\&dtri;g{\left(x\right)}^T\mathrm{\&Delta;x}-L^{-1}({\mathrm{ZL}}_{z0}+L_{l0}^{\mathrm{\&mu;}})\\ \mathrm{\&Delta;w}=U^{-1}W\&dtri;g{\left(x\right)}^T\mathrm{\&Delta;x}+U^{-1}({\mathrm{WL}}_{w0}-L_{l0}^{\mathrm{\&mu;}})\end{array}\right.$

其中：

$H(\&CenterDot;)\&equiv;{\&dtri;}_{x}g\left(x\right)(U^{-1}W-L^{-1}Z){\&dtri;}_{x}g{\left(x\right)}^T+$

$({-\&dtri;}_x^{2}f\left(x\right)+{\&dtri;}_x^{2}h\left(x\right)y+{\&dtri;}_x^{2}g\left(x\right)(z+w))=H_g+H_h$

$J\left(x\right)\&equiv;{\&dtri;}_{x}h{\left(x\right)}^T$

$\mathrm{\&psi;}(g,\mathrm{\&mu;})\&equiv;{-\&dtri;}_{x}f\left(x\right)+{\&dtri;}_{x}h\left(x\right)y-{\&dtri;}_{x}g\left(x\right)\left((U^{-1}-L^{-1}\right)\mathrm{\&mu;e}$

$+L^{-1}Z(g\left(x\right)-l-\underset{\&OverBar;}{g})-U^{-1}W(g\left(x\right)+u-\stackrel{\&OverBar;}{g}))$

设定一个转换阈值ε_switch，ε_switch的取值范围为2×10^-6～3×10^-6，对原始-对偶互补间隙C_Gap进行判断，若C_Gap＜ε_switch，则由于x＝[s，t]^T，

在计算上述修正方程中的H(·)时，在

中，

各元素采用下式精确表达：

${\&dtri;}_{t_i}^{2}f(s,t)=\frac{e^{-\frac{1}{2}\frac{t_i^2}{{\mathrm{\&sigma;}}^2}}}{{\mathrm{\&sigma;}}^2}(1-\frac{t_i^2}{{\mathrm{\&sigma;}}^2}),i=\mathrm{1,2},...,n_{\mathrm{meas}}$

其中n_meas为上述量测数据的个数，t_i为残差变量的第i个分量，σ为核估计法的窗宽；

若C_Gap≥ε，则

各元素采用下式的近似表达：

${\&dtri;}_{t_i}^{2}f(s,t)\&ap;\frac{e^{-\frac{1}{2}\frac{t_i^2}{{\mathrm{\&sigma;}}^2}}}{{\mathrm{\&sigma;}}^2},i=\mathrm{1,2},...,n_{\mathrm{meas}}$

(4-10)根据下式计算在第k次迭代时原始变量的修正步长

和对偶变量的修正步长

${\mathrm{Step}}_P=0.9995\min \left\{\underset{i}{\min }(\frac{{-l}_i}{{\mathrm{\&Delta;l}}_i}:{\mathrm{\&Delta;l}}_i<0;\frac{{-u}_i}{{\mathrm{\&Delta;u}}_i}:{\mathrm{\&Delta;u}}_i<0)\right\}$

${\mathrm{Step}}_D=0.9995\min \left\{\underset{i}{\min }(\frac{{-z}_i}{{\mathrm{\&Delta;z}}_i}:{\mathrm{\&Delta;z}}_i<0;\frac{{-w}_i}{{\mathrm{\&Delta;w}}_i}:{\mathrm{\&Delta;w}}_i<0)\right\}$

(4-11)根据步骤(4-10)的修正步长，更新上述原始变量和对偶变量：

使 $\left[\begin{array}{c}x\\ l\\ u\end{array}\right]=\left[\begin{array}{c}x\\ l\\ u\end{array}\right]+s{\mathrm{tep}}_P\left[\begin{array}{c}\mathrm{\&Delta;x}\\ \mathrm{\&Delta;l}\\ \mathrm{\&Delta;u}\end{array}\right],$ $\left[\begin{array}{c}y\\ z\\ w\end{array}\right]=\left[\begin{array}{c}y\\ z\\ w\end{array}\right]+{\mathrm{step}}_D\left[\begin{array}{c}\mathrm{\&Delta;y}\\ \mathrm{\&Delta;z}\\ \mathrm{\&Delta;w}\end{array}\right],k=k+1$

(4-12)重复步骤(4-5)-步骤(4-12)。

Claims (1)

1.一种基于内点法的电力系统抗差状态估计方法，其特征在于，该方法包括以下步骤：

(1)从电力系统的数据采集中心实时获取量测数据，量测数据包括节点电压幅值、节点有功功率、节点无功功率、支路有功功率和支路无功功率，并存储量测数据，根据量测数据进行拓扑分析，得到电力系统结构图，根据量测数据计算电力系统导纳矩阵Y；

(2)根据步骤(1)中量测数据、电力系统结构图和电力系统导纳矩阵Y，建立初始电力系统抗差状态估计模型如下：

$\begin{array}{cc}\underset{s}{\min }& J\left(s\right)\end{array}=-\underset{i=1}{\overset{n_{\mathrm{meas}}}{\mathrm{\&Sigma;}}}\exp (-\frac{{(m_i\left(s\right)-m_i^{\mathrm{meas}})}^2}{2{\mathrm{\&sigma;}}^2})$

s.t. c(s)＝0

其中，s为电力系统状态量，s＝[θ，V]^T，m_i(s)为电力系统第i个量测方程的估计值，

为上述量测数据中的第i个量测值，σ是核估计法的窗宽，窗宽σ的取值范围为：0＜σ＜1，c(s)＝0表示零注入等式约束；

(3)对上述初始电力系统抗差状态估计模型进行等价转换，得到电力系统抗差状态估计模型如下：

$\min f(s,t)=-\underset{i}{\mathrm{\&Sigma;}}{e}^{-\frac{{t_i}^2}{{2\mathrm{\&sigma;}}^2}}$

s.t.c(s，t)＝0

$g_i(s,t)=m_i\left(s\right)-m_i^{\mathrm{meas}}-t_i=0,i=\mathrm{1,2},...,n_{\mathrm{meas}}$

其中，m_i(s)为电力系统第i个量测方程的估计值，t为残差变量，残差变量t中的第i个分量为：

i＝1，2，...，n_meas，n_meas为上述量测数据的个数，s为电力系统状态量，s＝[θ，V]^T，θ为节点电压的相角，V是节点电压的幅值，T是矩阵转置；

(4)采用原-对偶内点法，对上述电力系统抗差状态估计模型进行求解，具体实现过程如下：

(4-1)将上述电力系统抗差状态估计模型改写为以下标准形式：

min f(x)

s.t h(x)＝0

$\underset{\&OverBar;}{g}\&le;g\left(x\right)\&le;\stackrel{\&OverBar;}{g}$

其中，x是电力系统状态量s和残差变量t合并后的复合向量，x＝[s，t]^T，函数h(x)包含量测方程m_i(s)的约束以及零注入约束，g(x)为电力系统节点负荷与复合向量x之间的函数关系，

为用户设定的电力系统负荷上限，g为用户设定的电力系统负荷下限；

(4-2)将上述标准形式中的不等式约束简化为等式约束，得到简化后的电力系统抗差状态估计模型：

min f(x)

s.t h(x)＝0

g(x)-1-g＝0

$g\left(x\right)+u-\stackrel{\&OverBar;}{g}=0$

(l，u)≥0

其中，(l，u)∈R^r为将不等式约束转换为等式约束的松弛向量，R为实数集，r为函数g(x)的维数；

(4-3)根据上述简化后的电力系统抗差状态估计模型，构建一个增广拉格朗日函数如下：

$L(x,l,u,y,z,w,\tilde{z}\widetilde{,w})\&equiv;f\left(x\right)-y^{T}h\left(x\right)-z^T(g\left(x\right)-l-\underset{-}{g})$

$-w^T(g\left(x\right)+u-\stackrel{-}{g})-{\tilde{z}}^{T}l-{\tilde{w}}^{T}u$

其中，x为电力系统状态量s与残差变量t合并后的复合向量，对应于原-对偶内点法计算中的原始变量，(l，u)∈R^r为松弛向量，r为函数g(x)的维数，y，z，w，

分别为原-对偶内点法计算中的对偶变量；

(4-4)设置迭代初始值，将原始变量x中，电力系统状态量s中的电压幅值初值设为1，电压相角初值设为0，将原始变量x中的残差变量t的初始值设为0，将对偶变量的初始值设为0，设置迭代次数的初始值k＝0，设置最大迭代次数K_max，K_max的取值范围为：10～100，对于电力系统中的零注入母线，设定残差变量的第i个分量

(4-5)对迭代次数k进行判断，若k＜K_max，则进行步骤(4-6)，若k≥K_max，则结束计算，并输出计算不收敛；

(4-6)根据下式计算上述简化后的电力系统抗差状态估计模型的互补间隙C_Gap：

$C_{\mathrm{Gap}}=\underset{i=1}{\overset{r}{\mathrm{\&Sigma;}}}(l_i{z}_i-u_i{w}_i)$

其中，l和u为松弛变量，z和w为对偶变量，r为函数g(x)的维数；

设定一个电力系统抗差状态估计计算精度ε，将互补间隙C_Gap与计算精度ε进行比较，若C_Gap＜ε，则输出x，包括电力系统状态量和残差变量，并结束计算，若C_Gap≥ε，则进行步骤(4-7)；

(4-7)根据下式计算扰动因子μ：

$\mathrm{\&mu;}\&equiv;\mathrm{\&tau;}\frac{C_{\mathrm{Gap}}}{2r}$

其中τ为中心参数，取值范围为τ∈(0，1]，r为函数g(x)的维数，C_Gap为上述互补间隙；

(4-8)按照以下卡罗需-库恩-塔克方程组，计算上述增广拉格朗日函数对原-对偶内点法中各原始变量和对偶变量的一阶偏导数：

$({\&dtri;}_x^{2}f\left(x\right)-{\&dtri;}_x^{2}h\left(x\right)y-{\&dtri;}_x^{2}g\left(x\right)(z+w))\mathrm{\&Delta;x}$

$-{\&dtri;}_{x}h\left(x\right)\mathrm{\&Delta;y}-{\&dtri;}_{x}g\left(x\right)(\mathrm{\&Delta;z}+\mathrm{\&Delta;w})=-{\&dtri;}_x{L}_0$

${\&dtri;}_{x}h{\left(x\right)}^T\mathrm{\&Delta;x}=-{\&dtri;}_y{L}_0$

${\&dtri;}_{x}g{\left(x\right)}^T\mathrm{\&Delta;x}-\mathrm{\&Delta;l}=-{\&dtri;}_z{L}_0$

${\&dtri;}_{x}g{\left(x\right)}^T+\mathrm{\&Delta;x}+\mathrm{\&Delta;u}=-{\&dtri;}_w{L}_0$

$\mathrm{Z\&Delta;l}+\mathrm{L\&Delta;z}=-{\&dtri;}_l^{\mathrm{\&mu;}}{L}_0$

$\mathrm{W\&Delta;u}+\mathrm{U\&Delta;w}=-{\&dtri;}_u^{\mathrm{\&mu;}}{L}_0$

其中，

分别是卡罗需-库恩-塔克方程组对应的扰动方程的残差，和

分别是步骤(4-1)电力系统抗差状态估计模型标准形式中的f(x)、h(x)和g(x)的海森矩阵；

(4-9)根据以下修正方程，计算得到在第k次迭代时原始变量x的修正量Δx、对偶变量y、z、w的修正量Δy、Δz、Δw和松弛变量l、u的修正量Δl、Δu：

$\left\{\begin{array}{c}\mathrm{\&Delta;l}=\&dtri;g{\left(x\right)}^T\mathrm{\&Delta;x}+L_{z0}\\ \mathrm{\&Delta;u}=-(\&dtri;g{\left(x\right)}^T+L_{w0})\end{array}\right.$

$\left\{\begin{array}{c}\mathrm{\&Delta;z}=-L^{-1}Z\&dtri;g{\left(x\right)}^T\mathrm{\&Delta;x}-L^{-1}({\mathrm{ZL}}_{z0}+L_{l0}^{\mathrm{\&mu;}})\\ \mathrm{\&Delta;w}=U^{-1}W\&dtri;g{\left(x\right)}^T\mathrm{\&Delta;x}+U^{-1}({\mathrm{WL}}_{w0}-L_{l0}^{\mathrm{\&mu;}})\end{array}\right.$

其中：

$H(\&CenterDot;)\&equiv;{\&dtri;}_{x}g\left(x\right)(U^{-1}W-L^{-1}Z){\&dtri;}_{x}g{\left(x\right)}^T+$

$(-{\&dtri;}_x^{2}f\left(x\right)+{\&dtri;}_x^{2}h\left(x\right)y+{\&dtri;}_x^{2}g\left(x\right)(z+w))=H_g+H_h$

$J\left(x\right)\&equiv;{\&dtri;}_{x}h{\left(x\right)}^T$

$\mathrm{\&psi;}(g,\mathrm{\&mu;})\&equiv;-{\&dtri;}_{x}f\left(x\right)+{\&dtri;}_{x}h\left(x\right)y-{\&dtri;}_{x}g\left(x\right)((U^{-1}-L^{-1})\mathrm{\&mu;e}$

$+L^{-1}Z(g\left(x\right)-l-\underset{\&OverBar;}{g})-U^{-1}W(g\left(x\right)+u-\stackrel{\&OverBar;}{g})$

设定一个转换阈值ε_switch，ε_switch的取值范围为2×10^-6～3×10^-6，对原始-对偶互补间隙C_Gap进行判断，若C_Gap＜ε_switch，则由于x＝[s，t]^T，

在计算上述修正方程中的H(·)时，在中，

各元素采用下式精确表达：

${\&dtri;}_{t_i}^{2}f(s,t)=\frac{e^{-\frac{1t_i^2}{2{\mathrm{\&sigma;}}^2}}}{{\mathrm{\&sigma;}}^2}(1-\frac{t_i^2}{{\mathrm{\&sigma;}}^2}),i=\mathrm{1,2},...,n_{\mathrm{meas}}$

其中n_meas为上述量测数据的个数，t_i为残差变量的第i个分量，σ为核估计法的窗宽；

若C_Gap≥ε，则各元素采用下式的近似表达：

${\&dtri;}_{t_i}^{2}f(s,t)\&ap;\frac{e^{-\frac{1}{2}\frac{t_i^2}{{\mathrm{\&sigma;}}^2}}}{{\mathrm{\&sigma;}}^2},i=\mathrm{1,2},...,n_{\mathrm{meas}}$

(4-10)根据下式计算在第k次迭代时原始变量的修正步长

和对偶变量的修正步长

${\mathrm{Step}}_P=0.9995\min \left\{\underset{i}{\min }(\frac{-l_i}{{\mathrm{\&Delta;l}}_i}:{\mathrm{\&Delta;l}}_i<0;\frac{{-u}_i}{{\mathrm{\&Delta;u}}_i}:{\mathrm{\&Delta;u}}_i<0)\right\}$

${\mathrm{Step}}_D=0.9995\min \left\{\underset{i}{\min }(\frac{{-z}_i}{{\mathrm{\&Delta;z}}_i}:{\mathrm{\&Delta;z}}_i<0;\frac{{-w}_i}{{\mathrm{\&Delta;w}}_i}:{\mathrm{\&Delta;w}}_i<0)\right\}$

(4-11)根据步骤(4-10)的修正步长，更新上述原始变量和对偶变量：

使 $\left[\begin{array}{c}x\\ l\\ u\end{array}\right]=\left[\begin{array}{c}x\\ l\\ u\end{array}\right]+{\mathrm{step}}_P\left[\begin{array}{c}\mathrm{\&Delta;x}\\ \mathrm{\&Delta;l}\\ \mathrm{\&Delta;u}\end{array}\right],$ $\left[\begin{array}{c}y\\ z\\ w\end{array}\right]=\left[\begin{array}{c}y\\ z\\ w\end{array}\right]+{\mathrm{step}}_D\left[\begin{array}{c}\mathrm{\&Delta;y}\\ \mathrm{\&Delta;z}\\ \mathrm{\&Delta;w}\end{array}\right],$ k＝k+1

(4-12)重复步骤(4-5)～步骤(4-12)。