一种考虑零注入量测等式约束的电力系统状态估计方法

摘要

本发明涉及一种考虑零注入量测等式约束的电力系统状态估计方法，属于电力系统调度自动化与电网仿真技术领域。该方法包括：将零注入量测作为等式约束加入到状态估计的优化模型中，形成一个含等式约束的优化模型。对该模型将零注入节点的电压幅值和相角作为基变量，而将非零注入节点的电压幅值和相角作为非基变量，参与寻优，采用既约梯度法求解，零注入节点电压根据非零注入节点电压确定。本发明方法可使零注入约束方程严格满足，可以解决目前电力系统状态估计中零注入节点功率失配量较大导致估计结果不满足潮流方程的问题。

说明书

技术领域

本发明涉及一种考虑零注入量测等式约束的电力系统状态估计方法，属于电力系统调度自动化与电网仿真技术领域。

背景技术

电力系统状态估计是能量管理系统(EMS)的核心基础功能之一，EMS的几乎所有在线和离线应用都是以状态估计结果为基础进行的。实际电力系统中存在大量注入功率为0的零注入节点。在目前调度中心常用的最小二乘法状态估计(WLS)程序中，通常的处理方法是加入量测值为0的节点注入功率伪量测，并对这些伪量测设置大权重。如果零注入量测权重设置得过小，最终的估计结果中零注入节点将存在较大的功率失配量，从而导致估计结果不符合潮流方程；反之，如果零注入量测权重过大，最小二乘法中的信息矩阵条件数将恶化，可能导致状态估计发散。

此外，为了解决最小二乘法状态估计结果可能受不良数据影响的问题，我们提出了新的指数型抗差估计模型，具体方法已申请中国发明专利(申请号：200910082501.8)。这类新型的估计模型也存在类似的零注入量测问题。

数学上处理等式约束的常规算法是拉格朗日乘子法，即将等式约束优化问题：

$\underset{x}{\min }J\left(x\right)$

s.t.c(x)＝0

加入拉格朗日乘子，转化为无约束优化问题：

$\underset{x,\lambda }{\min }J\left(x\right)+{\lambda }^{T}c\left(x\right)$

在实际的电力系统状态估计问题中，几乎所有220kV及以上电压等级节点，以及三卷变压器中心节点都属于零注入节点。对于省级电网，零注入节点的比例往往在50％左右。因此，传统的拉格朗日乘子法会大幅地增加了优化变量数目，造成计算量急剧上升。这也是拉格朗日乘子法长期没有在电力系统状态估计领域实际应用的主要原因。

发明内容

本发明的目的是提出一种考虑零注入量测等式约束的电力系统状态估计方法，在电力系统调度中心，将零注入节点的注入功率为0作为等式约束加入到优化方程中，并采用既约梯度法对此含等式约束的优化模型进行求解。

包括以下步骤：

(1)将电力系统中零注入节点的注入功率作为等式约束，建立一个含等式约束的电力系统状态估计模型：

min J(x)

s.t c(x)＝0

其中J(x)为状态估计的目标函数，该目标函数为最小二乘目标函数或可微非最小二乘目标函数，c(x)为电力系统零注入伪量测的量测方程，其中的分量c_i(x)为：

$\left(\begin{array}{c}{V}_i\underset{j\in i}{\Sigma }{V}_j(G_{\mathrm{ij}}\cos {\theta }_{\mathrm{ij}}+B_{\mathrm{ij}}{\sin \theta }_{\mathrm{ij}})=0\\ V_i\underset{j\in i}{\Sigma }{V}_j(G_{\mathrm{ij}}{\sin \theta }_{\mathrm{ij}}-B_{\mathrm{ij}}{\cos \theta }_{\mathrm{ij}})=0\end{array}\right)$

其中V_i是节点i的电压幅值，V_j是节点j的电压幅值，θ_ij是节点i和j之间的相角差值，G_ij和B_ij分别是节点导纳矩阵中的元素；

(2)对上述步骤(1)中建立的含等式约束的电力系统状态估计模型，采用既约梯度法进行求解，求解步骤如下：

(2-1)将状态变量x划分为基变量x_B和非基变量x_N，其中x_B为零注入节点的电压和相角，x_N为电力系统中非零注入节点的电压和相角；

(2-2)设定状态变量x的初值x⁽⁰⁾和收敛精度ε，设置迭代次数计数器k，k＝0；

(2-3)计算在当前点x^(k)处目标函数J(x)对非基变量x_N的全梯度如下：

$g_N^{\left(k\right)}=\frac{\mathrm{dJ}\left(x\right)}{{\mathrm{dx}}_N^{\left(k\right)}}=\frac{\partial J\left(x\right)}{{\partial x}_N^{\left(k\right)}}+\frac{\partial J\left(x\right)}{{\partial x}_B^{\left(k\right)}}\frac{{\partial x}_B^{\left(k\right)}}{{\partial x}_N^{\left(k\right)}}$

$=\frac{\partial J\left(x\right)}{{\partial x}_N^{\left(k\right)}}-{\left[{\left(\frac{\partial c\left(x\right)}{{\partial x}_B^{\left(k\right)}}\right)}^{-1}\frac{\partial c\left(x\right)}{{\partial x}_N^{\left(k\right)}}\right]}^T\frac{\partial J\left(x\right)}{{\partial x}_B^{\left(k\right)}}$

(2-4)设定若则状态估计收敛，退出计算，若则用共轭梯度方向取代梯度方向如下：

$d^{\left(k\right)}=-g_N^{\left(k\right)}+{\beta }_k{d}^{(k-1)},k\ge 1$

$d^{\left(k\right)}={-g}_N^{\left(0\right)},k=0$

其中系数β_k采用PRP(Polak，Ribiere，Polyak)共轭梯度法求取：

(2-5)根据上述步骤(2-4)确定的寻优方向进行一维搜索，求出最优步长则k＝k+1，重复步骤(2-3)-(2-5)。

本发明提出的考虑零注入量测等式约束的电力系统状态估计方法，其优点是：将零注入量测作为等式约束加入到状态估计的优化模型中，并采用既约梯度法进行求解，仅将非零注入节点的节点电压x_N作为寻优变量，零注入节点节点电压x_B根据非零注入节点的节点电压x_N确定。方法与状态估计模型的目标函数无关。因此本发明具有以下优点：

1、在本发明方法中，零注入节点的注入功率作为等式约束处理，其估计值严格为0，状态估计结果严格满足潮流方程。

2、本发明方法中，寻优变量仅为非零注入节点电压x_N，寻优空间维度小，计算效率高。

3、本发明方法中计算全梯度时，只涉及等式约束函数c(x)对基变量x_B求导的雅可比矩阵求逆，该雅可比矩阵稀疏，采用因子分解法计算，计算快速。

4、本发明方法通用性强，适用于最小二乘估计和其他目标函数可微的非最小二乘估计。

5、本发明方法采用梯度类算法，收敛可靠，迭代不会发散。在极端情况下，可在估计未收敛时中止状态估计计算，得到的估计结果仍严格满足潮流方程。

具体实施方式

本发明提出的考虑零注入量测等式约束的电力系统状态估计方法，包括以下步骤：

1、将零注入节点的节点注入功率为0看成等式约束，并建立一个含等式约束的状态估计模型如(1)：

min J(x)            (1)

s.t c(x)＝0

(1)中J(x)为状态估计的目标函数，c(x)是没有挂接负荷和发电机的节点的有功和无功零注入伪量测的量测方程。其分量C_i(x)可表示成：

$\left(\begin{array}{c}{V}_i\underset{j\in i}{\Sigma }{V}_j(G_{\mathrm{ij}}\cos {\theta }_{\mathrm{ij}}+B_{\mathrm{ij}}{\sin \theta }_{\mathrm{ij}})=0\\ V_i\underset{j\in i}{\Sigma }{V}_j(G_{\mathrm{ij}}{\sin \theta }_{\mathrm{ij}}-B_{\mathrm{ij}}{\cos \theta }_{\mathrm{ij}})=0\end{array}\right)---\left(2\right)$

其中V_i是节点i的电压幅值，V_j是节点j的电压幅值，θ_ij是节点i和节点j间的相角差值，G_ij和B_ij分别是节点导纳矩阵中的元素。

状态估计模型(1)对于目标函数的形式没有特殊的要求，如果采用常规的最小二乘法状态估计，则状态估计模型可写成：

s.t c(x)＝0

或最大负指数平方抗差估计模型可写成：(具体方法已申请中国专利(申请号为200910082501.8))

$\max J\left(x\right)=\underset{i}{\Sigma }\exp (-\frac{4}{R_{\mathrm{ii}}}{(z_i-h_i\left(x\right))}^2)---\left(4\right)$

s.t c(x)＝0

在上述(3)和(4)式中，z_i是量测值，包括线路或变压器的有功功率P_ij和无功功率Q_ij、节点的电压幅值V_i、发电机和负荷的有功功率P_i和无功功率Q_i等，R_ii是实时量测的方差，x是电力系统的状态变量，包括所有节点的电压幅值V和相角θ，h_i(x)是系统量测方程，对于常见的量测类型，其定义为：

线路的实时量测方程为：

$\left(\begin{array}{c}{P}_{\mathrm{ij}}^s=V_i^2{g}_{\mathrm{ij}}-V_i{V}_j(g_{\mathrm{ij}}{\cos \theta }_{\mathrm{ij}}+b_{\mathrm{ij}}{\sin \theta }_{\mathrm{ij}})\\ Q_{\mathrm{ij}}^s=-V_i^2(b_{\mathrm{ij}}+y_c)-V_i{V}_j(g_{\mathrm{ij}}{\sin \theta }_{\mathrm{ij}}-b_{\mathrm{ij}}{\cos \theta }_{\mathrm{ij}})\end{array}\right)---\left(5\right)$

变压器的实时量测方程为：

$\left(\begin{array}{c}{P}_{\mathrm{ij}}^s=-\frac{1}{k}{V}_i{V}_j{b}_{\mathrm{ij}}{\sin \theta }_{\mathrm{ij}}\\ Q_{\mathrm{ij}}^s={-V}_i^2{b}_{\mathrm{ij}}+\frac{1}{k}{V}_i{V}_j{b_{\mathrm{ij}}\cos \theta }_{\mathrm{ij}}\end{array}\right)---\left(6\right)$

上式中，P_ij^s和Q_ij^s分别是线路或变压器支路的有功功率和无功功率估计值，V_i和V_j分别是节点i和节点j的电压幅值，θ_ij是节点i和节点j之间的相角差值；对(5)式的线路，g_ij，b_ij和y_c分别是线路电导、电纳和充电电纳；对(6)式的变压器支路，b_ij和k分别是变压器支路的电纳和变比，规定j侧为非标准变比侧。

节点i的电压实时量测方程：

$V_i^s=V_i---\left(7\right)$

V_i^s为节点i的电压估计值。

节点i的注入实时量测方程：

$\left(\begin{array}{c}{P}_i^s=V_i\underset{j\in i}{\Sigma }{V}_j(G_{\mathrm{ij}}\cos {\theta }_{\mathrm{ij}}+B_{\mathrm{ij}}{\sin \theta }_{\mathrm{ij}})\\ Q_i^s=V_i\underset{j\in i}{\Sigma }{V}_j(G_{\mathrm{ij}}\sin {\theta }_{\mathrm{ij}}+B_{\mathrm{ij}}{\cos \theta }_{\mathrm{ij}})\end{array}\right)---\left(8\right)$

上式中，P_i、Q_i是节点i的有功注入功率和无功注入功率估计值。

2、对步骤1中建立的含等式约束的状态估计模型，采用既约梯度法进行求解。根据等式约束，将系统变量x划分为基变量x_B和非基变量x_N。基变量x_B通过等式约束由非基变量x_N确定，即式(1)可写成以下的无约束优化模型：

min J(x_N，x_B(x_N))            (9)

其中x_B(x_N)表示基变量x_B的取值是非基变量x_N的函数，(8)式的优化模型中，仅将非基变量x_N作为优化变量，基变量x_B的取值由非基变量x_N确定。这样不但保证了优化结果满足等式约束c(x)＝0，而且大幅减小了寻优空间维数，提高了计算效率。

在电力系统状态估计模型(1)中，x_B取为零注入节点的节点电压和相角，x_N取为非零注入节点的节点电压和相角。则可使用梯度法，对(8)式的无优化约束问题进行求解。具体步骤如下：

(2-1)给状态变量x设初值给定收敛精度ε，设置迭代次数计数器k，k＝0。

(2-2)计算在当前点x^(k)处目标函数J(x)对非基变量x_N的全梯度：

$g_N^{\left(k\right)}=\frac{\mathrm{dJ}\left(x\right)}{{\mathrm{dx}}_N^{\left(k\right)}}=\frac{\partial J\left(x\right)}{{\partial x}_N^{\left(k\right)}}+\frac{\partial J\left(x\right)}{{\partial x}_B^{\left(k\right)}}\frac{{\partial x}_B^{\left(k\right)}}{{\partial x}_N^{\left(k\right)}}---\left(10\right)$

$=\frac{\partial J\left(x\right)}{{\partial x}_N^{\left(k\right)}}-{\left[{\left(\frac{\partial c\left(x\right)}{{\partial x}_B^{\left(k\right)}}\right)}^{-1}\frac{\partial c\left(x\right)}{{\partial x}_N^{\left(k\right)}}\right]}^T\frac{\partial J\left(x\right)}{{\partial x}_B^{\left(k\right)}}$

J(x)对非基变量x_N的梯度由两部分构成，分别是J(x)本身对x_N的梯度，以及由于x_N的变化，通过等式约束c(x)＝0作用下引起基变量x_B变化，从而对目标函数施加的影响。对于电力系统状态估计问题，目标函数J(x)对x_B，x_N的部分梯度由目标函数的形式决定，而等式约束对x_B，x_N的导数为量测方程雅可比矩阵的相应子阵。

(2-3)若则状态估计收敛，退出计算，若则进行步骤(2-4)

(2-4)为改善算法的收敛速度，克服梯度法在接近收敛时会产生震荡现象等缺点，可用共轭梯度方向取代梯度方向。即如下确定寻优方向：

$d^{\left(k\right)}=-g_N^{\left(k\right)}+{\beta }_k{d}^{(k-1)},k\ge 1---\left(11\right)$

$d^{\left(0\right)}={-g}_N^{\left(0\right)},k=0$

其中系数β_k采用PRP(Polak，Ribiere，Polyak)共轭梯度法求取，系数β_k取为：

${\beta }_k=\frac{{g_N^{\left(k\right)}}^T(g_N^{\left(k\right)}-g_N^{(k-1)})}{{g_N^{(k-1)}}^T{g}_N^{(k-1)}}---\left(12\right)$

(2-5)根据(2-4)确定的寻优方向进行一维搜索，求出最优步长。一维搜索的算法没有特殊要求，根据经验，比较有效的一维搜索算法是三次曲线拟合法。一维搜索可能需要非基变量x_N求取沿寻优方向d^(k)前进步长α后的目标函数值J(x_N+αd^(k))，以及在该点处目标函数值对步长的导数

将电力系统的节点划分为零注入节点和非零注入节点后，线性化的网络方程如下：

$\left(\begin{array}{cc}{Y}_{\mathrm{BB}}& Y_{\mathrm{BN}}\\ Y_{\mathrm{NB}}& Y_{\mathrm{NN}}\end{array}\right)\left(\begin{array}{c}{\stackrel{·}{V}}_B\\ {\stackrel{·}{V}}_N\end{array}\right)=\left(\begin{array}{c}{\stackrel{·}{I}}_B\\ {\stackrel{·}{I}}_N\end{array}\right)---\left(13\right)$

其中下标B为零注入节点(基变量)，下标N为非基变量。其中可由当前的寻优方向d^(k)以及步长α确定，为已知量。对于零注入节点，其节点注入电流始终为0。于是我们可以得到：

$Y_{\mathrm{BB}}{\stackrel{·}{V}}_B+Y_{\mathrm{BN}}{\stackrel{·}{V}}_N={\stackrel{·}{I}}_B=0---\left(14\right)$

因此有：

${\stackrel{·}{V}}_B=-Y_{\mathrm{BB}}^{-1}\left(Y_{\mathrm{BN}}{\stackrel{·}{V}}_N\right)---\left(15\right)$

基变量x_B的值可直接通过求解线性方程组确定，无须任何迭代计算。

根据当前的x_B，x_N，非基变量x_N沿寻优方向d^(k)前进步长α后的目标函数值J(x_N+αd^(k))就可以确定了。

而该点处目标函数值对步长的导数为：

$\frac{\mathrm{dJ}(x_N+{\mathrm{\alpha d}}^{\left(k\right)})}{\mathrm{d\alpha }}=\frac{\partial J(x_N+{\mathrm{\alpha d}}^{\left(k\right)})}{\partial (x_N+{\mathrm{\alpha d}}^{\left(k\right)})}\frac{d(x_N+{\mathrm{\alpha d}}^{\left(k\right)})}{\mathrm{d\alpha }}+\frac{\partial J(x_N+{\mathrm{\alpha d}}^{\left(k\right)})}{{\partial x}_B}\frac{{\partial x}_B}{\partial \alpha }---\left(16\right)$

$=\frac{\partial J\left(x\right)}{{\partial x}_B}{d}^{\left(k\right)}-\left[{\left(\frac{\partial h}{{\partial x}_B}\right)}^{-1}\frac{\partial h}{{\partial x}_N}\right]\frac{\partial J\left(x\right)}{{\partial x}_B}{d}^{\left(k\right)}=g(x_N+{\mathrm{\alpha d}}^{\left(k\right)})d^k$

即的数值为，非基变量x_N沿本次寻优方向d^(k)移动距离α后得到的新点的既约梯度和本次寻优方向d^(k)的点乘。

k＝k+1。重复步骤(2-2)-(2-6)。