CN101707373B - 基于自动微分的电力系统状态估计方法 - Google Patents

Abstract

本发明公布了一种基于自动微分的电力系统状态估计方法，将AD技术分别应用于线性和非线性状态估计中。首先，为AD分配内存，声明活跃变量。然后，在迭代前将雅克比矩阵和海森矩阵中不变元素的位置和数值存到一个链表中。接着迭代开始，运用AD工具计算雅克比矩阵和海森矩阵的可变元素，同时读取链表中相应矩阵的不变元素，以此获得状态估计所需的雅克比矩阵和海森矩阵。本发明由于用AD替代了传统的手工编写微分代码计算雅克比矩阵和海森矩阵，有效避免了截断误差，提高了计算效率，而且本发明便于在已有的状态估计软件上实现。

Description

基于自动微分的电力系统状态估计方法

技术领域

发明涉及一种基于自动微分的电力系统状态估计方法，属于电力系统运行和控制技术领域。

背景技术

能量管理系统(energy manage system，EMS)是以计算机为基础的现代电力调度自动化系统，其带来的最根本的改变就是由传统经验型调度上升到分析型调度，从而提高了电力系统运行的安全性和经济性。随着电网规模的不断扩大，电力系统的结构和运行方式日趋复杂，电力系统调度中心的自动化水平也需要逐步由低级向高级发展，因此EMS得到了广泛的应用。状态估计软件是EMS的核心软件，它为其他高级应用软件提供了一个可靠而完整的电力系统实时数据库，被誉为应用软件的“心脏”，因此状态估计是电力系统运行、控制和安全评估等方面的基础。

状态估计中量测函数是状态矢量的非线性矢量函数，过去往往采取将非线性函数泰勒级数展开并取至一次项，即化为线性模型来处理，从而忽略了模型的非线性程度。本发明采用双一次项展开方法，保留了泰勒级数的二阶项，提出了非线性最小二乘法。

状态估计中需要计算量测函数的雅克比(Jacobian)矩阵和海森(Hessian)矩阵。而在常规状态估计中，形成Jacobian矩阵和Hessian矩阵最直接的方法是手工对量测函数的代码进行微分，得到微分代码。手工编写微分代码工作过于繁琐且容易出错，自动微分(automatic differentiation，AD)技术的出现克服了这个缺点，它是计算机数值计算和分析领域内的一项完全崭新的技术。AD将微分定义为代数运算，与其他微分方法(如数值微分、符号微分)相比，它可以自动计算函数的任意阶导数，而且避免了截断误差。

发明内容

本发明所要解决的技术问题是针对现有技术存在的缺陷提供一种基于自动微分的电力系统状态估计方法。

本发明为实现上述目的，采用如下技术方案：

本发明为基于自动微分的电力系统状态估计方法，其特征在于包括以下步骤：

(1)获取电力系统的网络参数，包括：输电线路的支路号、首端节点和末端节点编号、串联电阻、串联电抗、并联电导、并联电纳、变压器变比和阻抗；

(2)初始化，包括：对状态量设置初值、节点次序优化、形成节点导纳矩阵、设置门槛值，分配内存、声明活跃变量；

(3)输入遥测数据z，包括电压幅值、发电机有功功率、发电机无功功率、负荷有功功率、负荷无功功率、线路首端有功功率、线路首端无功功率、线路末端有功功率以及线路末端无功功率；

(4)将雅克比矩阵和海森矩阵中不变元素的位置和数值存到一个链表中；

(5)恢复迭代计数器迭代次数k＝1；

(6)由现有的状态量x^(k)，应用自动微分计算雅克比矩阵和海森矩阵中的可变元素，同时读取步骤(4)所述链表中相应矩阵的不变元素，以此获得所需的雅克比矩阵和海森矩阵；

(7)求取状态修正量

选取

并修正状态量得到

$\mathrm{\Δ}{\hat{x}}_1^{\left(k\right)}=[H^T\left({\hat{x}}^{\left(k\right)}\right)R^{-1}H\left({\hat{x}}^{\left(k\right)}\right){]}^{-1}\×H^T\left({\hat{x}}^{\left(k\right)}\right)R^{-1}[z-h\left({\hat{x}}^{\left(k\right)}\right)]$

当自动微分运用于线性状态估计时，

取为零；

当自动微分运用于非线性状态估计时，如下：

${\mathrm{\Δ}\hat{x}}_2^{\left(k\right)}=[H^T\left({\hat{x}}^{\left(k\right)}\right)R^{-1}H\left({\hat{x}}^{\left(k\right)}\right){]}^{-1}\×[W\left({\hat{x}}^{\left(k\right)}\right)\mathrm{\Δ}{\hat{x}}_1^{\left(k\right)}{]}^T\left[R^{-1}\right[z-h\left({\hat{x}}^{\left(k\right)}\right)]-R^{-1}H\left({\hat{x}}^{\left(k\right)}\right)\mathrm{\Δ}{\hat{x}}_1^{\left(k\right)}]$

$\mathrm{\Δ}{\hat{x}}^{\left(k\right)}=\mathrm{\Δ}{\hat{x}}_1^{\left(k\right)}+\mathrm{\Δ}{\hat{x}}_2^{\left(k\right)}$

${\hat{x}}^{(k+1)}={\hat{x}}^{\left(k\right)}+{\mathrm{\Δ}\hat{x}}^{\left(k\right)}$

其中，x为电力系统的n维状态变量；

是状态量的估计值，i＝1，2，…n；z为量测值矢量即遥测数据；h(x)是x的非线性函数；H(x)为量测函数的雅克比矩阵；R^-1为量测矢量的加权阵；T为转秩符号；(k)表示迭代序号。(8)当

小于收敛标准，结束状态估计，否则返回步骤(6)进行第k+1次状态估计。

基本加权最小二乘(WLS)状态估计方法是电力系统状态估计的最基本方法，该方法模型简单，收敛性能好，估计质量高，是目前应用最为广泛的方法之一。这种方法通过围绕着状态估计值将量测函数进行线性化而把WLS方法应用于非线性系统，量测函数的线性化处理决定了WLS方法的两点不足：一是线性化带来的误差；二是手工编写微分代码计算雅克比矩阵的繁琐易错。

本发明提出的基于自动微分的电力系统状态估计方法，在WLS方法的基础上推导建立了非线性最小二乘的模型，并且通过AD替代传统的手工编写微分代码来计算雅克比矩阵和海森矩阵，有效避免了截断误差，提高了计算效率。

附图说明

图1：本发明方法流程图。

图2：本发明采用的元件等值电路图，其中：图(a)是线路∏形等值电路图，图(b)是变压器∏形等值电路图。

图3：本发明提出的基于AD技术的状态估计方法所应用的三个算例系统，其中：图(a)是IEEE-14节点系统，图(b)是IEEE-30节点系统，图(c)是IEEE-57节点系统。

具体实施方式

下面结合附图对发明的技术方案进行详细说明：

1969年美国麻省理工学院的许怀丕(F.C.Schweppe)等人提出了电力系统状态估计的最基本算法——基本加权最小二乘(WLS)状态估计算法，其基本思想是以量测量和量测估计值之差的平方和最小为目标准则的估计方法。这种方法通过围绕着状态估计值将量测函数进行线性化而把WLS方法应用于非线性系统，量测函数的线性化处理忽略了模型的非线性程度。本发明在WLS方法的基础上推导建立了非线性最小二乘的数学模型。

如图1所示，在给定网络接线、支路参数和量测系统的条件下，非线性量测方程可表示为：

z＝h(x)+v

式中，z为量测值矢量即遥测数据，绝大多数是通过遥测得到的，也有一小部分是人工设定的；h(x)为由基尔霍夫等基本电路定律所建立的量测函数；x为n维系统状态变量；v为量测随机误差，假设其服从均值为零、方差为σ²的正态分布。

在电力系统状态估计中，量测量配置的类型要比常规潮流多，不仅包括了各节点的注入功率量测P_i、Q_i，还可以包括支路的功率量测P_ij、Q_ij、P_ji、Q_ji以及节点的电压量测V_i。以节点注入有功功率量测P_i、支路始端有功功率量测P_ij以及节点电压量测V_i为例，量测方程如下式所示：

节点注入有功功率：

$P_i=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}(G_{\mathrm{ij}}{e}_i{e}_j-B_{\mathrm{ij}}{e}_i{f}_j+B_{\mathrm{ij}}{e}_j{f}_i+G_{\mathrm{ij}}{f}_i{f}_j)$

线路i-j上始端有功功率：

$P_{\mathrm{ij}}=(e_i^2+f_i^2)g-(e_i{e}_j+f_i{f}_j)g+(e_i{f}_j-e_j{f}_i)b$

节点i电压幅值：

$V_i^2=e_i^2+f_i^2$

上述量测方程中，e_i，e_j为节点i，j的电压实部；f_i，f_j为节点i，j的电压虚部；G_ij、B_ij为导纳矩阵的实部和虚部；g、b为线路的电导和电纳。

给定量测矢量z以后，状态估计问题就是求使目标函数

J(x)＝[z-h(x)]^TR^-1[z-h(x)]

达到最小时的x的值。其中，R是以σ_i ²为对角元素的量测误差方差阵，在状态估计中取其逆矩阵为量测矢量的加权阵。

由于h(x)是x的非线性函数，所以无法直接计算得出x的估计值

为了求取

以往常用的方法是对h(x)进行线性化处理，即进行泰勒级数展开并忽略二次以上非线性项，由此得到线性状态估计的迭代公式：

$\mathrm{\Δ}{\hat{x}}^{\left(k\right)}=[H^T\left({\hat{x}}^{\left(k\right)}\right)R^{-1}H\left({\hat{x}}^{\left(k\right)}\right){]}^{-1}\×H^T\left({\hat{x}}^{\left(k\right)}\right)R^{-1}[z-h\left({\hat{x}}^{\left(k\right)}\right)]$

${\hat{x}}^{(k+1)}={\hat{x}}^{\left(k\right)}+\mathrm{\Δ}{\hat{x}}^{\left(k\right)}$

式中(k)表示迭代序号， $H=\frac{\∂h\left(x\right)}{\∂x}.$

不同于线性状态估计，本发明对量测函数采用双一次项展开的方法取至泰勒级数的二阶项，建立非线性状态估计的数学模型。

由极值函数有

$H^T\left(\hat{x}\right)R^{-1}[z-h\left(\hat{x}\right)]=0$

令x₀是x的某一近似值，将

和在x₀附近泰勒级数展开并取至一次项，得

$H\left(\hat{x}\right)=H\left(x_0\right)+W\left(x_0\mathrm{\Δx}\right)$

$h\left(\hat{x}\right)=h\left(x_0\right)+H\left(x_0\right)\mathrm{\Δx}$

其中， $H=\frac{\∂h\left(x\right)}{\∂x},$ $W=\frac{{\∂}^{2}h\left(x\right)}{{\∂x}_i{\∂x}_j},(i,j=\mathrm{1,2},\·\·\·n).$

将上式代入可得：

[H^TR^-1H]Δx＝H^TR^-1Δz+[WΔx]^T×R^-1[Δz-HΔx]

其中，Δz＝z-h(x₀)。

上式等号右边Δx仍然存在，无法直接求解，这里将等式右边的Δx用[H^TR^-1H]^-1H^TR^-1Δz代替，由此得到：

${\mathrm{\Δ}\hat{x}}_1^{\left(k\right)}=[H^T\left({\hat{x}}^{\left(k\right)}\right)R^{-1}H\left({\hat{x}}^{\left(k\right)}\right){]}^{-1}\×H^T\left({\hat{x}}^{\left(k\right)}\right)R^{-1}[z-h\left({\hat{x}}^{\left(k\right)}\right)]$

${\mathrm{\Δ}\hat{x}}_2^{\left(k\right)}=[H^T\left({\hat{x}}^{\left(k\right)}\right)R^{-1}H\left({\hat{x}}^{\left(k\right)}\right){]}^{-1}\×[W\left({\hat{x}}^{\left(k\right)}\right)\mathrm{\Δ}{\hat{x}}_1^{\left(k\right)}{]}^T\left[R^{-1}\right[z-h\left({\hat{x}}^{\left(k\right)}\right)]-R^{-1}H\left({\hat{x}}^{\left(k\right)}\right)\mathrm{\Δ}{\hat{x}}_1^{\left(k\right)}]$

$\mathrm{\Δ}{\hat{x}}^{\left(k\right)}=\mathrm{\Δ}{\hat{x}}_1^{\left(k\right)}+\mathrm{\Δ}{\hat{x}}_2^{\left(k\right)}$

${\hat{x}}^{(k+1)}={\hat{x}}^{\left(k\right)}+{\mathrm{\Δ}\hat{x}}^{\left(k\right)}$

式中(k)表示迭代序号。

在状态估计中，需要求取量测函数的雅克比矩阵和海森矩阵，以节点注入有功功率P_i、线路i-j上的始端有功功率P_ij以及节点电压V_i为例，量测矢量的雅克比矩阵和海森矩阵对应的元素为：

雅克比矩阵：

节点注入功率

$\left\{\begin{array}{c}\frac{{\∂P}_i}{{\∂e}_i}=\underset{\underset{j\≠i}{j=1}}{\overset{n}{\mathrm{\Σ}}}(G_{\mathrm{ij}}{e}_j-B_{\mathrm{ij}}{f}_j)+{2G}_{\mathrm{ii}}{e}_i\\ \frac{{\∂P}_i}{{\∂f}_i}=\underset{\underset{j\≠i}{j=1}}{\overset{n}{\mathrm{\Σ}}}(B_{\mathrm{ij}}{e}_j+G_{\mathrm{ij}}{f}_j)+{2G}_{\mathrm{ii}}{f}_i\end{array}\right.$

$\left\{\begin{array}{c}\frac{{\∂P}_i}{{\∂e}_j}=G_{\mathrm{ij}}{e}_i+B_{\mathrm{ij}}{f}_i\\ \frac{{\∂P}_i}{{\∂f}_j}=G_{\mathrm{ij}}{f}_i-B_{\mathrm{ij}}{e}_i\end{array}\right.$

线路i-j上的始端功率

$\left\{\begin{array}{c}\frac{{\∂P}_{\mathrm{ij}}}{{\∂e}_i}=2{\mathrm{ge}}_i-{\mathrm{ge}}_j+{\mathrm{bf}}_j\\ \frac{{\∂P}_{\mathrm{ij}}}{{\∂f}_i}={2\mathrm{gf}}_i-{\mathrm{gf}}_j-{\mathrm{be}}_j\end{array}\right.$

$\left\{\begin{array}{c}\frac{{\∂P}_{\mathrm{ij}}}{{\∂e}_j}=-{\mathrm{ge}}_i-{\mathrm{bf}}_i\\ \frac{{\∂P}_{\mathrm{ij}}}{{\∂f}_j}=-{\mathrm{gf}}_i+{\mathrm{be}}_i\end{array}\right.$

节点电压

$\left\{\begin{array}{c}\frac{{\∂V}_i^2}{{\∂e}_i}={2e}_i\\ \frac{\∂V_i^2}{{\∂f}_i}={2f}_i\end{array}\right.$

由此可见，在直角坐标下雅克比矩阵的元素均是可变元素。

海森矩阵：

节点注入功率

$\frac{{\∂}^2{P}_i}{{\∂e}_i^2}=\frac{{\∂}^2{P}_i}{{\∂f}_i^2}={2G}_{\mathrm{ii}}$

$\frac{{\∂}^2{P}_i}{{\∂e}_i{\∂e}_j}=\frac{{\∂}^2{P}_i}{{\∂f}_i{\∂f}_j}=G_{\mathrm{ij}}$

$\frac{{\∂}^2{P}_i}{{\∂e}_i{\∂f}_j}=-\frac{{\∂}^2{P}_i}{{\∂f}_i{\∂e}_j}=-B_{\mathrm{ij}}$

线路i-j上的始端功率

$\frac{{\∂}^2{P}_{\mathrm{ij}}}{{\∂e}_i^2}=\frac{{\∂}^2{P}_{\mathrm{ij}}}{{\∂f}_i^2}=2g$

$\frac{{\∂}^2{P}_{\mathrm{ij}}}{{\∂e}_i{\∂e}_j}=\frac{{\∂}^2{P}_{\mathrm{ij}}}{{\∂f}_i{\∂f}_j}=-g$

$\frac{{\∂}^2{P}_{\mathrm{ij}}}{{\∂e}_i{\∂f}_j}=-\frac{{\∂}^2{P}_{\mathrm{ij}}}{{\∂f}_i{\∂e}_j}=b$

节点电压

$\frac{{\∂}^2{V}_i^2}{{\∂e}_i^2}=\frac{{\∂}^2{V}_i^2}{{\∂f}_i^2}=2$

由此可见，在直角坐标系下，海森矩阵的所有元素均为常数，在迭代中保持不变。

在状态估计中，形成雅克比矩阵和海森矩阵最直接的方法是手工对量测函数进行微分，推导出导数的解析表达式，然后再编制出计算机程序，这种方法可以称为“手工方法”。虽然这种方法很容易得到效率最高的代码，在状态估计中经常被使用，但是手工编写微分代码的工作太繁琐且容易出错，尤其是对于复杂的大型问题。

自动微分(AD)技术是计算机数值计算和分析领域内的一项完全崭新的技术。1961年A.Robinson引入“无穷大量”和“无穷小量”，建立了“非标准分析”的理论，提出在实数的延拓空间中，微分可以定义为代数运算。这种思想同计算机技术相结合就形成了今天的“自动微分”技术。

AD技术依赖于以下事实，无论函数多么复杂，总是由一系列基本元函数(如sin，cos，exp等)和四则运算符(+，-，×，÷)经过有限次组合而成。根据微分计算的链式法则，不需要推导导数计算公式，就可以以完全机械的方式获得函数的微分，而且不含截断误差。

目前存在两种自动微分的实现方法：源代码变换方法和运算符重载方法。源代码变换是显示的重写源代码以得到导数代码的方法；运算符重载是基于C++等高级语言中操作符重载的机理，通过重载每个操作以达到计算导数和根据链式法则传递导数的目的，代表软件有ADOL-C等。

ADOL-C是由Dresden技术大学科学计算学院开发的自动微分系统，使用运算符重载的方法对C++程序进行自动微分，能够以正模式和逆模式计算任意阶导数。该软件可以免费使用。

本发明从http://www.coin-or.org/projects/ADOL-C.xml获得ADOL-C软件，运用该软件自动求取雅克比矩阵和海森矩阵。通过上面的分析，发现在直角坐标系中海森矩阵的所有元素均为常数，在迭代过程中不随自变量的变化而变化。而ADOL-C软件在生成相应矩阵时，处理可变元素和不变元素的过程是相同的，因此其对海森矩阵的元素应用链式法则求导是在重复运算，降低了程序的性能。

针对这个问题，本发明提出的方法是这样处理的：对于不变元素，在迭代前先将其各自的位置和数值存储在一个链表中；对于可变元素，在迭代中运用ADOL-C软件自动求导。

如图2所示，图(a)是线路∏形等值电路图，节点i和节点j之间串接导纳g+j′b，节点i、j的输出端分别串接一个接地电纳j′y_c后接地。

图(b)是变压器∏形等值电路图，节点i和节点j之间串接

节点i串接一个

后接地，j的输出端串接一个

后接地。j′表示虚部。

下面介绍本发明的两个实施例：

实施例一：

本发明采用图3所示的IEEE-14节点、IEEE-30节点、IEEE-57节点的标准算例，将AD技术运用于线性状态估计中，并与手工编写微分代码的线性状态估计方法进行了对比，仿真结果如下表所示：

表1 AD技术应用于线性状态估计中的仿真结果

其中， $J=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\left[\frac{z_i-h_i\left(\hat{x}\right)}{{\mathrm{\σ}}_i}\right]}^2,$ $S_M={\left[\frac{1}{m}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\left(\frac{z_i-s_i}{{\mathrm{\σ}}_i}\right)}^2\right]}^{\frac{1}{2}},$ $S_B={\left\{\frac{1}{m}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\left[\frac{h_i\left(\hat{x}\right)-s_i}{{\mathrm{\σ}}_i}\right]}^2\right\}}^{\frac{1}{2}}$

式中，J为目标函数值；S_M为量测误差统计值；S_B为估计误差统计值；m为量测量数目；s_i为量测量i的真值。

实施例二：

在非线性状态估计中，应用手工编写微分代码的方法计算雅克比矩阵和海森矩阵，复杂程度明显增大。本发明采用图3所示的IEEE-14节点、IEEE-30节点、IEEE-57节点的标准算例，将AD技术运用于非线性状态估计中，仿真结果如下表所示：

表2 AD技术应用于非线性状态估计中的仿真结果

本发明方法具有如下优势，一方面，由表1所示的仿真结果可以看出，在线性状态估计中，相对于采用手工方法，采用AD技术的方法计算速度明显加快，而且避免了截断误差，计算精度也得到了提高，估计结果更加符合电力系统的实际状态；另一方面，由表2所示的仿真结果可以看出，将AD技术运用于非线性状态估计中，大大降低了应用的难度，提高了计算速度，尤其在大系统中，速度提高了一倍多。由此可见，AD技术在电力系统中的应用是可行的，优越的。

Claims (1)

1.一种基于自动微分的电力系统状态估计方法，其特征在于包括以下步骤：

(1)获取电力系统的网络参数，包括：输电线路的支路号、首端节点和末端节点编号、串联电阻、串联电抗、并联电导、并联电纳、变压器变比和阻抗；

(2)初始化，包括：对状态量设置初值、节点次序优化、形成节点导纳矩阵、设置门槛值、分配内存、声明活跃变量；

(3)输入遥测数据z，包括电压幅值、发电机有功功率、发电机无功功率、负荷有功功率、负荷无功功率、线路首端有功功率、线路首端无功功率、线路末端有功功率以及线路末端无功功率；

(4)将雅克比矩阵和海森矩阵中不变元素的位置和数值存到一个链表中；

(5)恢复迭代计数器迭代次数k＝1；

(6)由现有的状态量x^(k)，应用自动微分计算雅克比矩阵和海森矩阵中的可变元素，同时读取步骤(4)所述链表中相应矩阵的不变元素，以此获得所需的雅克比矩阵和海森矩阵；

(7)求取状态修正量选取

并修正状态量得到

$\mathrm{\Δ}{\hat{x}}_1^{\left(k\right)}={\left[H^T\left({\hat{x}}^{\left(k\right)}\right)R^{-1}H\left({\hat{x}}^{\left(k\right)}\right)\right]}^{-1}\×H^T\left({\hat{x}}^{\left(k\right)}\right)R^{-1}[z-h\left({\hat{x}}^{\left(k\right)}\right)],$

当自动微分运用于线性状态估计时，

取为零；

当自动微分运用于非线性状态估计时，

如下：

$\mathrm{\Δ}{\hat{x}}_2^{\left(k\right)}={\left[H^T\left({\hat{x}}^{\left(k\right)}\right)R^{-1}H\left({\hat{x}}^{\left(k\right)}\right)\right]}^{-1}\×{\left[W\left({\hat{x}}^{\left(k\right)}\right)\mathrm{\Δ}{\hat{x}}_1^{\left(k\right)}\right]}^T\left[R^{-1}\right[z-h\left({\hat{x}}^{\left(k\right)}\right)]-R^{-1}H\left({\hat{x}}^{\left(k\right)}\right)\mathrm{\Δ}{\hat{x}}_1^{\left(k\right)}]$

$\mathrm{\Δ}{\hat{x}}^{\left(k\right)}=\mathrm{\Δ}{\hat{x}}_1^{\left(k\right)}+{\mathrm{\Δ}\hat{x}}_2^{\left(k\right)},$

${\hat{x}}^{(k+1)}={\hat{x}}^{\left(k\right)}+\mathrm{\Δ}{\hat{x}}^{\left(k\right)},$

其中，x为电力系统的n维状态变量；

是状态量的估计值，i＝1，2，Λn，n为电力系统节点数；z为量测值矢量即遥测数据；h(x)是x的非线性函数，h(x)为量测函数；H(x)为量测函数的雅克比矩阵；W(x)为量测函数的海森矩阵；R^-1为量测矢量的加权阵；T为转秩符号；(k)表示迭代序号；

(8)当小于设定的收敛标准，则结束状态估计，否则返回步骤(6)进行第k+1次状态估计。

Images