CN102623993A - 一种分布式电力系统状态估计方法 - Google Patents

Abstract

本发明公开了一种分布式电力系统状态估计方法。本发明建立了考虑DG特性的分布式电力系统状态估计的模型，不仅可以得到电网中各节点的实时准确的运行状态，还可以得到DG自身的实时准确的运行状态；其次，将保留非线性法引入状态估计中，避免了截断误差，可以使得状态估计具有更高的计算精度，以及提高算法的迭代效率；最后，由于保留非线性法状态估计中需要计算量测函数的Jacobian矩阵和Hessian矩阵，而手工计算数量巨大的微分函数和编写微分代码的工作量繁重，工作过于繁琐且容易出错，所以本发明还通过使用自动微分技术替代传统的手工编写微分代码计算Jacobian矩阵和Hessian矩阵，减少了手写代码量，提高了程序的开发效率。

Description

一种分布式电力系统状态估计方法

技术领域

本发明涉及一种分布式电力系统状态估计方法，属于电力系统运行和控制技术领域。

背景技术

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

状态估计是EMS的重要组成部分，其结果直接影响电网调度的智能化分析与决策，它是远动装置与数据库之间的重要一环。在这里，要对电网中DG自身的运行状态进行监视，是保证分布式电力系统安全稳定运行基础，然而要获取它们准确的运行状态，那就必须要求对状态估计的功能进行完善，同时也要求状态估计根据电网的发展而不断的提高。现有分布式电力系统状态估计方法是将DG看做负荷处理，未考虑DG的类型、模型及自身状态变量，即只把DG看做一般的负荷节点，只会估计出节点的电压幅值U和相角θ，而无法得到DG自身的运行状态，也无法区分DG的类型等特征。

发明内容

本发明所要解决的技术问题在于克服现有技术不足，提供一种分布式电力系统状态估计方法，该方法充分考虑DG自身运行特点，建立了考虑DG特性的分布式电力系统状态估计模型，不仅可以得到电网中各节点的实时准确的运行状态，还可以得到DG的实时准确的运行状态。

本发明具体采用以下技术方案解决上述技术问题。

一种分布式电力系统状态估计方法，通过对以下分布式电力系统状态估计模型求解，估计出分布式电力系统状态量的值：

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

式中，z为实时量测量，包括：电压幅值、发电机有功功率、发电机无功功率、风机有功功率、风机无功功率、光伏有功功率、光伏无功功率、负荷有功功率、负荷无功功率、线路首端有功功率、线路首端无功功率、线路末端有功功率以及线路末端无功功率；状态量x＝[θ，U，s，M，α]^T，其中θ和U分别为电网中各节点的相角、电压幅值，s为风机发电模型中的滑差，M和α分别为光伏发电模型中的逆变器幅值调制比、相位调制角；R为以实时量测量的方差为对角元素的量测误差对角阵；h(x)为量测方程，其中，风机注入功率量测方程如下：

$P_w=-\frac{{\mathrm{sR}}_r{U}^2}{s^2{(X_s+X_r)}^2+R_r^2}$

$Q_w=-\frac{U^2}{X_m}-\frac{s^2{U}^2(X_s+X_r)}{s^2{(X_s+X_r)}^2+R_r^2}$

式中，P_w和Q_w分别为风机发电的有功功率和无功功率；s为滑差；X_s为定子绕组的电抗；R_r和X_r分别为转子绕组电阻和电抗；X_m为励磁绕组的电抗；

光伏注入功率量测方程如下：

式中，P_pv和Q_pv分别为光伏发电的有功功率和无功功率；V_cell为光伏电池的电压；N_s为单个光伏模块中光伏电池串联的数量；N_ss为串联光伏模块的数量；式中，P_pv和Q_pv分别为光伏发电的有功功率和无功功率；V_cell为光伏电池的电压；N_s为单个光伏模块中光伏电池串联的数量；N_ss为串联光伏模块的数量；M和α分别为光伏发电模型中的逆变器幅值调制比、相位调制角；z₁₂、

z₂₃、

z₁₃、

分别为光伏发电模型等值电路中的AC Part的Y型电路转换成的△型等值电路的等值阻抗幅值和相角。

在对状态估计模型求解时，传统的基于加权最小二乘算法的状态估计是线性化处理，会造成截断误差，为解决这个问题，本发明进一步采用保留非线性法对所述分布式电力系统状态估计模型进行求解，具体如下：

步骤A、设定估计迭代次数初始值k＝1；

步骤B、根据所述状态量x^(k)，计算所需的Jacobian矩阵；其中含有状态量s，M，α的分块扩展Jacobian矩阵为

也即为：

$H\left(x^{\left(k\right)}\right)={\left[\begin{array}{ccccc}\frac{\∂U}{\∂\mathrm{\θ}}& \frac{\∂U}{\∂U}& \frac{\∂U}{\∂s}& \frac{\∂U}{\∂M}& \frac{\∂U}{\∂\mathrm{\α}}\\ \frac{\∂P}{\∂\mathrm{\θ}}& \frac{\∂P}{\∂U}& \frac{\∂P}{\∂s}& \frac{\∂P}{\∂M}& \frac{\∂P}{\∂\mathrm{\α}}\\ \frac{\∂Q}{\∂\mathrm{\θ}}& \frac{\∂Q}{\∂U}& \frac{\∂Q}{\∂s}& \frac{\∂Q}{\∂M}& \frac{\∂Q}{\∂\mathrm{\α}}\\ \frac{\∂P_w}{\∂\mathrm{\θ}}& \frac{\∂P_w}{\∂U}& \frac{\∂P_w}{\∂s}& \frac{\∂P_w}{\∂M}& \frac{\∂P_w}{\∂\mathrm{\α}}\\ \frac{\∂Q_w}{\∂\mathrm{\θ}}& \frac{\∂Q_w}{\∂U}& \frac{\∂Q_w}{\∂s}& \frac{\∂Q_w}{\∂M}& \frac{\∂Q_w}{\∂\mathrm{\α}}\\ \frac{\∂P_{\mathrm{pv}}}{\∂\mathrm{\θ}}& \frac{\∂P_{\mathrm{pv}}}{\∂U}& \frac{\∂P_{\mathrm{pv}}}{\∂s}& \frac{\∂P_{\mathrm{pv}}}{\∂M}& \frac{\∂P_{\mathrm{pv}}}{\∂\mathrm{\α}}\\ \frac{\∂Q_{\mathrm{pv}}}{\∂\mathrm{\θ}}& \frac{\∂Q_{\mathrm{pv}}}{\∂U}& \frac{\∂Q_{\mathrm{pv}}}{\∂s}& \frac{\∂Q_{\mathrm{pv}}}{\∂M}& \frac{\∂Q_{\mathrm{pv}}}{\∂\mathrm{\α}}\end{array}\right]}_{x=x^{\left(k\right)}}$

式中下标w和pv分别表示风电场节点对应的功率和光伏电场节点对应的功率，

的维数与电网中风电场的个数相同，

的维数与电网中光伏电场的个数相同；

步骤C、根据所述状态量x^(k)，计算所需的Hessian矩阵；其中含有状态量s，M，α的分块扩展Hessian矩阵为

(i，j＝1，2，…，n)，也即为：

$W\left(x^{\left(k\right)}\right){=\left[\begin{array}{ccccc}\frac{{\∂}^{2}U}{\∂\mathrm{\θ}\∂x}& \frac{{\∂}^{2}U}{\∂U\∂x}& \frac{{\∂}^{2}U}{\∂s\∂x}& \frac{{\∂}^{2}U}{\∂M\∂x}& \frac{{\∂}^{2}U}{\∂\mathrm{\α}\∂x}\\ \frac{{\∂}^{2}P}{\∂\mathrm{\θ}\∂x}& \frac{{\∂}^{2}P}{\∂U\∂x}& \frac{{\∂}^{2}P}{\∂s\∂x}& \frac{{\∂}^{2}P}{\∂M\∂x}& \frac{{\∂}^{2}P}{\∂\mathrm{\α}\∂x}\\ \frac{{\∂}^{2}Q}{\∂\mathrm{\θ}\∂x}& \frac{{\∂}^{2}Q}{\∂U\∂x}& \frac{{\∂}^{2}Q}{\∂s\∂x}& \frac{{\∂}^{2}Q}{\∂M\∂x}& \frac{{\∂}^{2}Q}{\∂\mathrm{\α}\∂x}\\ \frac{{\∂}^2{P}_w}{\∂\mathrm{\θ}\∂x}& \frac{{\∂}^2{P}_w}{\∂U\∂x}& \frac{{\∂}^2{P}_w}{\∂s\∂x}& \frac{{\∂}^2{P}_w}{\∂M\∂x}& \frac{{\∂}^2{P}_w}{\∂\mathrm{\α}\∂x}\\ \frac{{\∂}^2{Q}_w}{\∂\mathrm{\θ}\∂x}& \frac{{\∂}^2{Q}_w}{\∂U\∂x}& \frac{{\∂}^2{Q}_w}{\∂s\∂x}& \frac{{\∂}^2{Q}_w}{\∂M\∂x}& \frac{{\∂}^2{Q}_w}{\∂\mathrm{\α}\∂x}\\ \frac{{\∂}^2{P}_{\mathrm{pv}}}{\∂\mathrm{\θ}\∂x}& \frac{{\∂}^2{P}_{\mathrm{pv}}}{\∂U\∂x}& \frac{{\∂}^2{P}_{\mathrm{pv}}}{\∂s\∂x}& \frac{{\∂}^2{P}_{\mathrm{pv}}}{\∂M\∂x}& \frac{{\∂}^2{P}_{\mathrm{pv}}}{\∂\mathrm{\α}\∂x}\\ \frac{{\∂}^2{Q}_{\mathrm{pv}}}{\∂\mathrm{\θ}\∂x}& \frac{{\∂}^2{Q}_{\mathrm{pv}}}{\∂U\∂x}& \frac{{\∂}^2{Q}_{\mathrm{pv}}}{\∂s\∂x}& \frac{{\∂}^2{Q}_{\mathrm{pv}}}{\∂M\∂x}& \frac{{\∂}^2{Q}_{\mathrm{pv}}}{\∂\mathrm{\α}\∂x}\end{array}\right]}_{x=x^{\left(k\right)}}$

步骤D、解下述方程，求得状态修正量Δx^(k)，选取并修正状态量得到x^(k+1)；

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

$\mathrm{\Δ}{x}_{\left(2\right)}^{\left(k\right)}={\left[H^T\left(x^{\left(k\right)}\right)R^{-1}H\left(x^{\left(k\right)}\right)\right]}^{-1}{\left[W\left(x^{\left(k\right)}\right)\mathrm{\Δ}{x}_{\left(1\right)}^{\left(k\right)}\right]}^T{R}^{-1}\left[\right[z-h\left(x^{\left(k\right)}\right)]-H\left(x^{\left(k\right)}\right)\mathrm{\Δ}{x}_{\left(1\right)}^{\left(k\right)}];$

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

x^(k+1)＝x^(k)+Δx^(k).

步骤E、判断

是否小于预设的收敛标准ε_x，如果是，则结束计算；否则，返回步骤B进行第k+1次估计。

更进一步地，所述Jacobian矩阵和Hessian矩阵中的可变元素利用自动微分(Automatic Diferentiation，AD)技术自动计算获得。从而避免了人工编写微分代码计算Jacobian矩阵和Hessian矩阵所存在的繁琐且效率低的问题。

相比现有技术，本发明方法具有以下有益效果：

首先，建立了考虑DG特性的分布式电力系统状态估计的模型，充分考虑了DG自身的运行特点，不再是简单的将其看成负荷来处理，不仅可以得到电网中各节点的实时准确的运行状态，还可以得到DG的实时准确的运行状态；其次，将保留非线性法引入状态估计中，避免了截断误差，可以使得状态估计具有更高的计算精度，以及提高算法的迭代效率；最后，由于保留非线性法状态估计中需要计算量测函数的Jacobian矩阵和Hessian矩阵，而手工计算数量巨大的微分函数和编写微分代码的工作量繁重，工作过于繁琐且容易出错，所以本发明还通过使用AD技术替代传统的手工编写微分代码计算Jacobian矩阵和Hessian矩阵，减少了手写代码量，提高了程序的开发效率。

附图说明

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

图2：本发明采用的风机发电模型等值电路图；

图3：本发明采用的光伏发电模型等值电路图，其中：图3(a)是光伏发电模型等值电路图，图3(b)是将图3(a)中AC Part的Y型电路转换成△型的等值电路图；

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

具体实施方式

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

电力系统的实时运行和控制需要了解系统的真实运行工况，由于测量和传输等方面的原因，得到的生数据难免存在误差，状态估计能在一定程度上提高数据的精度。自1969年美国麻省理工学院的许怀丕(F.C.Schweppe)等人提出了电力系统状态估计的最基本算法——基本加权最小二乘(Weighted Least Squares WLS)状态估计算法以来，加权最小二乘法成为电力系统状态估计中应用最多的算法。其基本思想是以量测量和量测估计值之差的平方和最小为目标准则的估计方法。本发明提出了一种考虑DG特性的新的分布式电力系统状态估计方法。该方法是建立在非线性WLS状态估计的模型基础上，考虑DG自身的运行特性，运用AD计算非线性WLS状态估计方法中的Jacobian矩阵和Hessian矩阵。而非线性WLS状态估计方法就是对量测函数采用双一次项展开来获取至泰勒级数的二阶项，建立非线性状态估计的数学模型。

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

z＝h(x)+v

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

在分布式电力系统状态估计中，量测量配置的类型要比常规潮流多，不仅包括了各节点的注入功率量测P_i、Q_i、P_wi、Q_wi、P_pvi、Q_pvi，还可以包括支路的功率量测P_ij、Q_ij、P_ji、Q_ji以及节点的电压幅值量测U_i。它们的量测函数分别如下所示：

节点i电压幅值(以下都是以极坐标表示)：

U_i＝U_i

节点注入有功功率和无功功率：

$P_i=U_i\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j(G_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}}+B_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}})$

$Q_i=U_i\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{U}_j(G_{\mathrm{ij}}\sin {\mathrm{\θ}}_{\mathrm{ij}}-B_{\mathrm{ij}}\cos {\mathrm{\θ}}_{\mathrm{ij}})$

式中θ_ij＝θ_i-θ_j，P_i和Q_i分别为节点i有功注入功率和无功注入功率，其方向规定：流入节点i为正，流出节点i为负。

其中Y_ij为节点导纳矩阵中对应节点i和j之间的元素，当i＝j时为自导纳，当i≠j时为互导纳，其通式为Y_ij＝G_ij+jB_ij，G_ij为电导，B_ij为电纳。

对于图4(a)所示的非变压器支路，其线路潮流量测方程如下所示：

非变压器支路i-j上始端有功功率和无功功率：

$P_{\mathrm{ij}}=U_i^{2}g-U_i{U}_{j}g\cos {\mathrm{\θ}}_{\mathrm{ij}}-U_i{U}_{j}b\sin {\mathrm{\θ}}_{\mathrm{ij}}$

$Q_{\mathrm{ij}}=-U_i^2(b+y_c)-U_i{U}_{j}g\sin {\mathrm{\θ}}_{\mathrm{ij}}+U_i{U}_{j}b\cos {\mathrm{\θ}}_{\mathrm{ij}}$

式中P_ij和Q_ij分别为线路i-j上始端有功功率和无功功率，其方向规定：由i流向j为正，由j流向i为负。

非变压器支路i-j上终端有功功率和无功功率：

$P_{\mathrm{ji}}=U_j^{2}g-U_i{U}_j\mathrm{gcso}{\mathrm{\θ}}_{\mathrm{ij}}+U_i{U}_{j}b\sin {\mathrm{\θ}}_{\mathrm{ij}}$

$Q_{\mathrm{ji}}=-U_j^2(b+y_c)+U_i{U}_{j}g\sin {\mathrm{\θ}}_{\mathrm{ij}}+U_i{U}_{j}b\cos {\mathrm{\θ}}_{\mathrm{ij}}$

式中P_ji和Q_ji分别为线路i-j上终端有功功率和无功功率，其方向规定：由j流向i为正，由i流向j为负。

其中y_ij为线路的导纳值，有

g为线路导纳，b为线路电纳，R为线路电阻，X为线路电抗，y_c为线路对地导纳值。

对于如图4(b)所示变压器支路，其线路潮流量测方程如下所示：

变压器支路i-j上始端有功功率和无功功率：

$P_{\mathrm{ij}}=-\frac{1}{K}{U}_i{U}_j{b}_T\sin {\mathrm{\θ}}_{\mathrm{ij}}$

$Q_{\mathrm{ij}}=-\frac{1}{K^2}{U}_i^2{b}_T+\frac{1}{K}{U}_i{U}_j{b}_T\cos {\mathrm{\θ}}_{\mathrm{ij}}$

其方向规定：由i流向j为正，由j流向i为负。

变压器支路i-j上终端有功功率和无功功率：

$P_{\mathrm{ji}}=\frac{1}{K}{U}_i{U}_j{b}_T\sin {\mathrm{\θ}}_{\mathrm{ij}}$

$Q_{\mathrm{ji}}=-b_T{U}_j^2+\frac{1}{K}{U}_i{U}_j{b}_T\cos {\mathrm{\θ}}_{\mathrm{ij}}$

其方向规定：由j流向i为正，由i流向j为负。

上述变压器支路量测方程中，K为变压器非标准变比。j为标准侧，变比为1；i为非标准侧，变比为K；b_T为变压器标准侧(j侧)的电纳，有

其中X_T为变压器标准侧之电抗。

本发明考虑了DG自身的运行特点。根据图2所示的风机发电模型等值电路，得出风机发电的有功功率和无功功率如下：

$P_w=-\frac{{\mathrm{sR}}_r{U}^2}{s^2{(X_s+X_r)}^2+R_r^2}$

$Q_w=-\frac{U^2}{X_m}-\frac{s^2{U}^2(X_s+X_r)}{s^2{(X_s+X_r)}^2+R_r^2}$

其中s为滑差；X_s为定子绕组的电抗；R_r和X_r分别为转子绕组电阻和电抗；X_m为励磁绕组的电抗；

在图3(a)所示的光伏发电模型等值电路中，把光伏发电模型等值电路中的AC Part的Y型电路转换成如图3(b)所示的△型等值电路，根据光伏发电模型等值电路可以得出光伏发电的有功功率和无功功率如下：

其中在光伏发电模型等值电路的PV arrays中的PV cell满足下面的I-V方程：

$I_{\mathrm{cell}}=I_L-I_O[\exp \left(\frac{V_{\mathrm{cell}}+I_{\mathrm{cell}}{R}_s}{a}\right)-1]-\frac{V_{\mathrm{cell}}+I_{\mathrm{cell}}{R}_s}{R_{\mathrm{sh}}}$

I_cell和V_cell需满足以下条件方程，使光伏发电工作在最大功率处：

$\frac{d\left(V_{\mathrm{cell}}{I}_{\mathrm{cell}}\right)}{{\mathrm{dV}}_{\mathrm{cell}}}{|}_{\mathrm{mpp}}=I_{\mathrm{cell},\mathrm{mpp}}+V_{\mathrm{cell},\mathrm{mpp}}\frac{\mathrm{dI}}{\mathrm{dV}}{|}_{\mathrm{cell},\mathrm{mpp}}=I_{\mathrm{cell},\mathrm{mpp}}+V_{\mathrm{cell},\mathrm{mpp}}\frac{-I_o{R}_{\mathrm{sh}}\exp \left(\frac{V_{\mathrm{cell},\mathrm{mpp}}+I_{\mathrm{cell},\mathrm{mpp}}{R}_s}{a}\right)-a}{{\mathrm{aR}}_{\mathrm{sh}}+L_o{R}_s{R}_{\mathrm{sh}}\exp \left(\frac{V_{\mathrm{cell},\mathrm{mpp}}+I_{\mathrm{cell},\mathrm{mpp}}{R}_s}{a}\right)+a{R}_s}=0$

由I-V方程和条件方程可以求出I_cell和V_cell，接着求出DC Part中的U_PV-arrays＝N_sN_ssV_cell和I_PV-arrays＝N_ppI_cell，以及Inverter Part中的 ${\stackrel{\·}{U}}_{\mathrm{IP}}=U_{\mathrm{IP}}\∠\mathrm{\α}=\frac{\sqrt{2}}{4}{\mathrm{MU}}_{\mathrm{PV}-\mathrm{arrays}}\∠\mathrm{\α}=\frac{\sqrt{2}}{4}{\mathrm{MN}}_s{N}_{\mathrm{ss}}{V}_{\mathrm{cell}}\∠\mathrm{\α};$ 最后得出 $U_{\mathrm{IP}}=\frac{\sqrt{2}}{4}{\mathrm{MN}}_s{N}_{\mathrm{ss}}{V}_{\mathrm{cell}};$ 其中I_L为光电流；I_O为二极管反向饱和电流；R_s为串联电阻；R_sh为并联电阻；a为校正理想因子；I_cell为光伏电池的电流；V_cell为光伏电池的电压；N_s为单个光伏模块中光伏电池串联的数量；N_ss为串联光伏模块的数量；N_pp为并联光伏模块的数量；U_PV-arrays为光伏阵列的电压；I_PV-arrays为光伏阵列的电流；M为逆变器的幅值调制比；α为逆变器的相位调制角；R为逆变器功率消耗等值电阻；L_f为滤波器的电感；C_f为滤波器的电容；R_T、X_T、G_T和B_T分别为升压变压器的电阻、电抗、电导和电纳。

分别将风机发电模型中的滑差s、光伏发电模型中的逆变器幅值调制比M和相位调制角α作为增广状态量引入修正方程，进行状态估计；状态估计状态量扩展为x＝[θ，U，s，M，α]^T，状态修正量扩展为Δx＝[Δθ，ΔU，Δs，ΔM，Δα]^T。这样即得到本发明的分布式电力系统状态估计模型。

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

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

达到最小时的x的值。其中，R是以

为对角元素的量测误差对角阵，在状态估计中取其逆矩阵为量测矢量的加权阵。

根据目标函数J(x)＝[z-h(x)]^TR^-1[z-h(x)]，由极值函数可以得出

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

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

为了求取

首先分别将H(x)和h(x)进行泰勒级数展开并取至一次项，然后迭代求解。假定状态量初值为x⁽⁰⁾，将

和

在x⁽⁰⁾处泰勒级数展开取至一次项，可得：

$\left\{\begin{array}{c}H\left(\hat{x}\right)=H\left(x^{\left(0\right)}\right)+W\left(x^{\left(0\right)}\right)\mathrm{\Δ}{\hat{x}}^{\left(0\right)}\\ h\left(\hat{x}\right)=h\left(x^{\left(0\right)}\right)+H\left(x^{\left(0\right)}\right)\mathrm{\Δ}{\hat{x}}^{\left(0\right)}\end{array}\right.---(**)$

式中， $\mathrm{\Δ}{\hat{x}}^{\left(0\right)}=\hat{x}-x^{\left(0\right)}.$

H(x⁽⁰⁾)和W(x⁽⁰⁾)分别是非线性函数方程h(x)的Jacobian矩阵和Hessian矩阵，H(x⁽⁰⁾)是一个二维矩阵，W(x⁽⁰⁾)是一个三维矩阵，其元素分别是

$H\left(x^{\left(0\right)}\right)=\frac{\∂h\left(x\right)}{\∂x}{|}_{x=x^{\left(0\right)}};$

$W\left(x^{\left(0\right)}\right)=\frac{{\∂}^{2}h\left(x\right)}{\∂x_i\∂x_j}{|}_{x=x^{\left(0\right)}},$ (i，j＝1，2，…，n).

在算法程序中H(x⁽⁰⁾)和W(x⁽⁰⁾)是运用AD技术自动求得的。

将方程(**)代入方程(*)中，可得：

$\left[H^T\left(x^{\left(0\right)}\right)R^{-1}H\left(x^{\left(0\right)}\right)\right]\mathrm{\Δ}{\hat{x}}^{\left(0\right)}=H^T\left(x^{\left(0\right)}\right)R^{-1}\mathrm{\Δz}+{\left[W\left(x^{\left(0\right)}\right)\mathrm{\Δ}{\hat{x}}^{\left(0\right)}\right]}^T{R}^{-1}[\mathrm{\Δz}-H\left(x^{\left(0\right)}\right)\mathrm{\Δ}{\hat{x}}^{\left(0\right)}]---(***)$

式中，Δz＝z-h(x⁽⁰⁾)。

方程(***)等号右边

仍然存在，无法直接求解，根据WLS状态估计的结论，这里将等式右边的

用[H^T(x⁽⁰⁾)R^-1H(x⁽⁰⁾)]^-1H^T(x⁽⁰⁾)R^-1[z-h(x⁽⁰⁾)]代替，由此可以得

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

$\mathrm{\Δ}{\hat{x}}_{\left(2\right)}^{\left(0\right)}={\left[H^T\left(x^{\left(0\right)}\right)R^{-1}H\left(x^{\left(0\right)}\right)\right]}^{-1}{\left[W\left(x^{\left(0\right)}\right)\mathrm{\Δ}{\hat{x}}_{\left(1\right)}^{\left(0\right)}\right]}^T{R}^{-1}\left[\right[z-h\left(x^{\left(0\right)}\right)]-H\left(x^{\left(0\right)}\right)\mathrm{\Δ}{\hat{x}}_{\left(1\right)}^{\left(0\right)}];$

$\mathrm{\Δ}{\hat{x}}^{\left(0\right)}=\mathrm{\Δ}{\hat{x}}_{\left(1\right)}^{\left(0\right)}+\mathrm{\Δ}{\hat{x}}_{\left(2\right)}^{\left(0\right)}.$

因此可以求出

$\hat{x}=x^{\left(0\right)}+\mathrm{\Δ}{\hat{x}}^{\left(0\right)}.$

由此可见，若以(k)表示迭代次数，则非线性WLS状态估计的迭代计算公式如下：

$\left\{\begin{array}{c}\mathrm{\Δ}{\hat{x}}_{\left(1\right)}^{\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{\Δ}{x}_{\left(2\right)}^{\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}}_{\left(1\right)}^{\left(k\right)}\right]}^T{R}^{-1}\left[\right[z-h\left({\hat{x}}^{\left(k\right)}\right)]-H\left({\hat{x}}^{\left(k\right)}\right)\mathrm{\Δ}{\hat{x}}_{\left(1\right)}^{\left(k\right)}];\\ \mathrm{\Δ}{\hat{x}}^{\left(k\right)}={\mathrm{\Δ}\hat{x}}_{\left(1\right)}^{\left(k\right)}+\mathrm{\Δ}{\hat{x}}_{\left(2\right)}^{\left(k\right)};\\ {\hat{x}}^{(k+1)}={\hat{x}}^{\left(k\right)}+\mathrm{\Δ}{\hat{x}}^{\left(k\right)}.\end{array}\right.$

其中

和

分别是非线性函数方程h(x)的Jacobian矩阵和Hessian矩阵，

是一个二维矩阵，

是一个三维矩阵，其元素分别是

$H\left({\hat{x}}^{\left(k\right)}\right)=\frac{\∂h\left(x\right)}{\∂x}{|}_{x={\hat{x}}^{\left(k\right)}};$

$W\left({\hat{x}}^{\left(k\right)}\right)=\frac{{\∂}^{2}h\left(x\right)}{{\∂x}_i\∂x_j}{|}_{x={\hat{x}}^{\left(k\right)}},$ (i，j＝1，2，…，n).

在算法程序中

和

是运用AD技术自动求得的，不需要手工推导微分公式和编写微分代码。Jacobian矩阵和Hessian矩阵中有些元素的值随着状态量变化而变化，但也有很大一部分元素的值均为常数或零，其在迭代过程中不随状态量的变化而变化。而AD工具在生成Jacobian矩阵和Hessian矩阵时，处理可变元素和不变元素的过程是相同的，因此其对Jacobian矩阵和Hessian矩阵中的不变元素应用链式法则求导是在重复运算，降低了程序的计算效率。针对这个问题，本发明提出的处理方法是：在迭代前，将Jacobian矩阵中不变元素各自的位置和数值存储在程序中开辟的链表Link1中，而将Hessian矩阵中不变元素各自的位置和数值存储在程序中开辟的链表Link2中，在生成Jacobian矩阵和Hessian矩阵时直接将其从链表Link1和Link2中分别读出放入对应的位置，这样就可以避免对Jacobian矩阵和Hessian矩阵中不变元素的重复计算，并且只需要对Jacobian矩阵和Hessian矩阵中可变元素，在迭代中运用AD技术自动求出，最终自动获取Jacobian矩阵和Hessian矩阵，提高程序的计算效率。

本发明方法的详细流程如图1所示，按照以下各步骤：

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

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

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

(4)构造DG注入功率量测方程：

风机注入功率量测方程如下：

$P_w=-\frac{{\mathrm{sR}}_r{U}^2}{s^2{(X_s+X_r)}^2+R_r^2}$

$Q_w=-\frac{U^2}{X_m}-\frac{s^2{U}^2(X_s+X_r)}{s^2{(X_s+X_r)}^2+R_r^2}$

光伏注入功率量测方程如下：

(5)将Jacobian矩阵中不变元素的位置和数值存到程序中开辟的一个链表Link1中；

(6)将Hessian矩阵中不变元素的位置和数值存到程序中开辟的另一个链表Link2中；

(7)恢复迭代计数器：k＝1；

(8)根据本发明所采用的扩展状态量x^(k)，运用AD技术计算Jacobian矩阵中的可变元素，同时读取链表Link1中相应矩阵的不变元素，获得所需的Jacobian矩阵；其中含有状态量s，M，α的分块扩展Jacobian矩阵为

也即为：

$H\left(x^{\left(k\right)}\right)={\left[\begin{array}{ccccc}\frac{\∂U}{\∂\mathrm{\θ}}& \frac{\∂U}{\∂U}& \frac{\∂U}{\∂s}& \frac{\∂U}{\∂M}& \frac{\∂U}{\∂\mathrm{\α}}\\ \frac{\∂P}{\∂\mathrm{\θ}}& \frac{\∂P}{\∂U}& \frac{\∂P}{\∂s}& \frac{\∂P}{\∂M}& \frac{\∂P}{\∂\mathrm{\α}}\\ \frac{\∂Q}{\∂\mathrm{\θ}}& \frac{\∂Q}{\∂U}& \frac{\∂Q}{\∂s}& \frac{\∂Q}{\∂M}& \frac{\∂Q}{\∂\mathrm{\α}}\\ \frac{\∂P_w}{\∂\mathrm{\θ}}& \frac{\∂P_w}{\∂U}& \frac{\∂P_w}{\∂s}& \frac{\∂P_w}{\∂M}& \frac{\∂P_w}{\∂\mathrm{\α}}\\ \frac{\∂Q_w}{\∂\mathrm{\θ}}& \frac{\∂Q_w}{\∂U}& \frac{\∂Q_w}{\∂s}& \frac{\∂Q_w}{\∂M}& \frac{\∂Q_w}{\∂\mathrm{\α}}\\ \frac{\∂P_{\mathrm{pv}}}{\∂\mathrm{\θ}}& \frac{\∂P_{\mathrm{pv}}}{\∂U}& \frac{\∂P_{\mathrm{pv}}}{\∂s}& \frac{\∂P_{\mathrm{pv}}}{\∂M}& \frac{\∂P_{\mathrm{pv}}}{\∂\mathrm{\α}}\\ \frac{\∂Q_{\mathrm{pv}}}{\∂\mathrm{\θ}}& \frac{\∂Q_{\mathrm{pv}}}{\∂U}& \frac{\∂Q_{\mathrm{pv}}}{\∂s}& \frac{\∂Q_{\mathrm{pv}}}{\∂M}& \frac{\∂Q_{\mathrm{pv}}}{\∂\mathrm{\α}}\end{array}\right]}_{x=x^{\left(k\right)}}$

式中下标w和pv分别表示风电场节点对应的功率和光伏电场节点对应的功率，

的维数与电网中风电场的个数相同，

的维数与电网中光伏电场的个数相同；

(9)根据本发明所采用的扩展状态量x^(k)，运用AD技术计算Hessian矩阵中的可变元素，同时读取链表Link2中相应矩阵的不变元素，获得所需的Hessian矩阵；其中含有状态量s，M，α的分块扩展Hessian矩阵为(i，j＝1，2，…，n)，也即为：

$W\left(x^{\left(k\right)}\right){=\left[\begin{array}{ccccc}\frac{{\∂}^{2}U}{\∂\mathrm{\θ}\∂x}& \frac{{\∂}^{2}U}{\∂U\∂x}& \frac{{\∂}^{2}U}{\∂s\∂x}& \frac{{\∂}^{2}U}{\∂M\∂x}& \frac{{\∂}^{2}U}{\∂\mathrm{\α}\∂x}\\ \frac{{\∂}^{2}P}{\∂\mathrm{\θ}\∂x}& \frac{{\∂}^{2}P}{\∂U\∂x}& \frac{{\∂}^{2}P}{\∂s\∂x}& \frac{{\∂}^{2}P}{\∂M\∂x}& \frac{{\∂}^{2}P}{\∂\mathrm{\α}\∂x}\\ \frac{{\∂}^{2}Q}{\∂\mathrm{\θ}\∂x}& \frac{{\∂}^{2}Q}{\∂U\∂x}& \frac{{\∂}^{2}Q}{\∂s\∂x}& \frac{{\∂}^{2}Q}{\∂M\∂x}& \frac{{\∂}^{2}Q}{\∂\mathrm{\α}\∂x}\\ \frac{{\∂}^2{P}_w}{\∂\mathrm{\θ}\∂x}& \frac{{\∂}^2{P}_w}{\∂U\∂x}& \frac{{\∂}^2{P}_w}{\∂s\∂x}& \frac{{\∂}^2{P}_w}{\∂M\∂x}& \frac{{\∂}^2{P}_w}{\∂\mathrm{\α}\∂x}\\ \frac{{\∂}^2{Q}_w}{\∂\mathrm{\θ}\∂x}& \frac{{\∂}^2{Q}_w}{\∂U\∂x}& \frac{{\∂}^2{Q}_w}{\∂s\∂x}& \frac{{\∂}^2{Q}_w}{\∂M\∂x}& \frac{{\∂}^2{Q}_w}{\∂\mathrm{\α}\∂x}\\ \frac{{\∂}^2{P}_{\mathrm{pv}}}{\∂\mathrm{\θ}\∂x}& \frac{{\∂}^2{P}_{\mathrm{pv}}}{\∂U\∂x}& \frac{{\∂}^2{P}_{\mathrm{pv}}}{\∂s\∂x}& \frac{{\∂}^2{P}_{\mathrm{pv}}}{\∂M\∂x}& \frac{{\∂}^2{P}_{\mathrm{pv}}}{\∂\mathrm{\α}\∂x}\\ \frac{{\∂}^2{Q}_{\mathrm{pv}}}{\∂\mathrm{\θ}\∂x}& \frac{{\∂}^2{Q}_{\mathrm{pv}}}{\∂U\∂x}& \frac{{\∂}^2{Q}_{\mathrm{pv}}}{\∂s\∂x}& \frac{{\∂}^2{Q}_{\mathrm{pv}}}{\∂M\∂x}& \frac{{\∂}^2{Q}_{\mathrm{pv}}}{\∂\mathrm{\α}\∂x}\end{array}\right]}_{x=x^{\left(k\right)}}$

注意这里的Hessian矩阵W是一个三维矩阵，其中x＝[θ，U，s，M，α]^T。

(10)解下述方程，求得状态修正量Δx^(k)，选取

并修正状态量得到x^(k+1)；

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

$\mathrm{\Δ}{x}_{\left(2\right)}^{\left(k\right)}={\left[H^T\left(x^{\left(k\right)}\right)R^{-1}H\left(x^{\left(k\right)}\right)\right]}^{-1}{\left[W\left(x^{\left(k\right)}\right)\mathrm{\Δ}{x}_{\left(1\right)}^{\left(k\right)}\right]}^T{R}^{-1}\left[\right[z-h\left(x^{\left(k\right)}\right)]-H\left(x^{\left(k\right)}\right)\mathrm{\Δ}{x}_{\left(1\right)}^{\left(k\right)}];$

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

x^(k+1)＝x^(k)+Δx^(k).

(11)判断是否小于收敛标准ε_x(一般ε_x取10^-5)，如果是，结束计算，否则返回步骤(8)进行第k+1次估计。

Claims (3)

1.一种分布式电力系统状态估计方法，其特征在于，通过对以下分布式电力系统状态估计模型求解，估计出分布式电力系统状态量的值：

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

式中，z为实时量测量，包括：电压幅值、发电机有功功率、发电机无功功率、风机有功功率、风机无功功率、光伏有功功率、光伏无功功率、负荷有功功率、负荷无功功率、线路首端有功功率、线路首端无功功率、线路末端有功功率以及线路末端无功功率；状态量x＝[θ，U，s，M，α]^T，其中θ和U分别为电网中各节点的相角、电压幅值，s为风机发电模型中的滑差，M和α分别为光伏发电模型中的逆变器幅值调制比、相位调制角；R为以实时量测量的方差为对角元素的量测误差对角阵；h(x)为量测方程，其中，风机注入功率量测方程如下：

$P_w=-\frac{{\mathrm{sR}}_r{U}^2}{s^2{(X_s+X_r)}^2+R_r^2}$

$Q_w=-\frac{U^2}{X_m}-\frac{s^2{U}^2(X_s+X_r)}{s^2{(X_s+X_r)}^2+R_r^2}$

式中，P_w和Q_w分别为风机发电的有功功率和无功功率；s为滑差；X_s为定子绕组的电抗；R_r和X_r分别为转子绕组电阻和电抗；X_m为励磁绕组的电抗；

光伏注入功率量测方程如下：

式中，P_pv和Q_pv分别为光伏发电的有功功率和无功功率；V_cell为光伏电池的电压；N_s为单个光伏模块中光伏电池串联的数量；N_ss为串联光伏模块的数量；M和α分别为光伏发电模型中的逆变器幅值调制比、相位调制角；z₁₂、

z₂₃、

z₁₃、

分别为光伏发电模型等值电路中的AC Part的Y型电路转换成的△型等值电路的等值阻抗幅值和相角。

2.如权利要求1所述分布式电力系统状态估计方法，其特征在于，采用保留非线性法对所述分布式电力系统状态估计模型进行求解，具体如下：

步骤A、设定估计迭代次数初始值k＝1；

步骤B、根据所述状态量x^(k)，计算所需的Jacobian矩阵；其中含有状态量s，M，α的分块扩展Jacobian矩阵为也即为：

$H\left(x^{\left(k\right)}\right)={\left[\begin{array}{ccccc}\frac{\∂U}{\∂\mathrm{\θ}}& \frac{\∂U}{\∂U}& \frac{\∂U}{\∂s}& \frac{\∂U}{\∂M}& \frac{\∂U}{\∂\mathrm{\α}}\\ \frac{\∂P}{\∂\mathrm{\θ}}& \frac{\∂P}{\∂U}& \frac{\∂P}{\∂s}& \frac{\∂P}{\∂M}& \frac{\∂P}{\∂\mathrm{\α}}\\ \frac{\∂Q}{\∂\mathrm{\θ}}& \frac{\∂Q}{\∂U}& \frac{\∂Q}{\∂s}& \frac{\∂Q}{\∂M}& \frac{\∂Q}{\∂\mathrm{\α}}\\ \frac{\∂P_w}{\∂\mathrm{\θ}}& \frac{\∂P_w}{\∂U}& \frac{\∂P_w}{\∂s}& \frac{\∂P_w}{\∂M}& \frac{\∂P_w}{\∂\mathrm{\α}}\\ \frac{\∂Q_w}{\∂\mathrm{\θ}}& \frac{\∂Q_w}{\∂U}& \frac{\∂Q_w}{\∂s}& \frac{\∂Q_w}{\∂M}& \frac{\∂Q_w}{\∂\mathrm{\α}}\\ \frac{\∂P_{\mathrm{pv}}}{\∂\mathrm{\θ}}& \frac{\∂P_{\mathrm{pv}}}{\∂U}& \frac{\∂P_{\mathrm{pv}}}{\∂s}& \frac{\∂P_{\mathrm{pv}}}{\∂M}& \frac{\∂P_{\mathrm{pv}}}{\∂\mathrm{\α}}\\ \frac{\∂Q_{\mathrm{pv}}}{\∂\mathrm{\θ}}& \frac{\∂Q_{\mathrm{pv}}}{\∂U}& \frac{\∂Q_{\mathrm{pv}}}{\∂s}& \frac{\∂Q_{\mathrm{pv}}}{\∂M}& \frac{\∂Q_{\mathrm{pv}}}{\∂\mathrm{\α}}\end{array}\right]}_{x=x^{\left(k\right)}}$

式中下标w和pv分别表示风电场节点对应的功率和光伏电场节点对应的功率，

的维数与电网中风电场的个数相同，

的维数与电网中光伏电场的个数相同；

步骤C、根据所述状态量x^(k)，计算所需的Hessian矩阵；其中含有状态量s，M，α的分块扩展Hessian矩阵为

(i，j＝1，2，…，n)，也即为：

$W\left(x^{\left(k\right)}\right){=\left[\begin{array}{ccccc}\frac{{\∂}^{2}U}{\∂\mathrm{\θ}\∂x}& \frac{{\∂}^{2}U}{\∂U\∂x}& \frac{{\∂}^{2}U}{\∂s\∂x}& \frac{{\∂}^{2}U}{\∂M\∂x}& \frac{{\∂}^{2}U}{\∂\mathrm{\α}\∂x}\\ \frac{{\∂}^{2}P}{\∂\mathrm{\θ}\∂x}& \frac{{\∂}^{2}P}{\∂U\∂x}& \frac{{\∂}^{2}P}{\∂s\∂x}& \frac{{\∂}^{2}P}{\∂M\∂x}& \frac{{\∂}^{2}P}{\∂\mathrm{\α}\∂x}\\ \frac{{\∂}^{2}Q}{\∂\mathrm{\θ}\∂x}& \frac{{\∂}^{2}Q}{\∂U\∂x}& \frac{{\∂}^{2}Q}{\∂s\∂x}& \frac{{\∂}^{2}Q}{\∂M\∂x}& \frac{{\∂}^{2}Q}{\∂\mathrm{\α}\∂x}\\ \frac{{\∂}^2{P}_w}{\∂\mathrm{\θ}\∂x}& \frac{{\∂}^2{P}_w}{\∂U\∂x}& \frac{{\∂}^2{P}_w}{\∂s\∂x}& \frac{{\∂}^2{P}_w}{\∂M\∂x}& \frac{{\∂}^2{P}_w}{\∂\mathrm{\α}\∂x}\\ \frac{{\∂}^2{Q}_w}{\∂\mathrm{\θ}\∂x}& \frac{{\∂}^2{Q}_w}{\∂U\∂x}& \frac{{\∂}^2{Q}_w}{\∂s\∂x}& \frac{{\∂}^2{Q}_w}{\∂M\∂x}& \frac{{\∂}^2{Q}_w}{\∂\mathrm{\α}\∂x}\\ \frac{{\∂}^2{P}_{\mathrm{pv}}}{\∂\mathrm{\θ}\∂x}& \frac{{\∂}^2{P}_{\mathrm{pv}}}{\∂U\∂x}& \frac{{\∂}^2{P}_{\mathrm{pv}}}{\∂s\∂x}& \frac{{\∂}^2{P}_{\mathrm{pv}}}{\∂M\∂x}& \frac{{\∂}^2{P}_{\mathrm{pv}}}{\∂\mathrm{\α}\∂x}\\ \frac{{\∂}^2{Q}_{\mathrm{pv}}}{\∂\mathrm{\θ}\∂x}& \frac{{\∂}^2{Q}_{\mathrm{pv}}}{\∂U\∂x}& \frac{{\∂}^2{Q}_{\mathrm{pv}}}{\∂s\∂x}& \frac{{\∂}^2{Q}_{\mathrm{pv}}}{\∂M\∂x}& \frac{{\∂}^2{Q}_{\mathrm{pv}}}{\∂\mathrm{\α}\∂x}\end{array}\right]}_{x=x^{\left(k\right)}}$

步骤D、解下述方程，求得状态修正量Δx^(k)，选取

并修正状态量得到x^(k+1)；

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

$\mathrm{\Δ}{x}_{\left(2\right)}^{\left(k\right)}={\left[H^T\left(x^{\left(k\right)}\right)R^{-1}H\left(x^{\left(k\right)}\right)\right]}^{-1}{\left[W\left(x^{\left(k\right)}\right)\mathrm{\Δ}{x}_{\left(1\right)}^{\left(k\right)}\right]}^T{R}^{-1}\left[\right[z-h\left(x^{\left(k\right)}\right)]-H\left(x^{\left(k\right)}\right)\mathrm{\Δ}{x}_{\left(1\right)}^{\left(k\right)}];$

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

x^(k+1)＝x^(k)+Δx^(k).

步骤E、判断

是否小于预设的收敛标准ε_x，如果是，则结束计算；否则，返回步骤B进行第k+1次估计。

3.如权利要求2所述分布式电力系统状态估计方法，其特征在于，所述Jacobian矩阵和Hessian矩阵中的可变元素利用自动微分技术自动计算获得。

Images