CN113595072A - 基于gpu并行加速共轭梯度算法的配网抗差状态估计方法 - Google Patents

Abstract

本发明公开了一种基于GPU并行加速共轭梯度算法的配网抗差状态估计方法，包括：构建配电网系统的量测方程、基于指数加权最小二乘模型、正常测点模型以及考虑量测数据不确定性的区间分析状态估计模型；求解所述基于指数加权最小二乘模型，将求解结果代入量测方程并计算量测残差，根据所述正常测点模型，对于异常测点采用量测计算值替代异常测点的实际量测值；求解所述区间分析状态估计模型，得到状态变量的不确定值区间；其中在指数加权最小二乘模型求解过程中，利用GPU并行加速的方法进行求解。本方法解决配电网状态估计中大量低精度数据造成估计精度不高甚至估计结果严重偏离实际的问题，同时减少状态估计的计算时间。

Description

基于GPU并行加速共轭梯度算法的配网抗差状态估计方法

技术领域

本发明涉及配电网控制领域，具体涉及一种基于GPU并行加速共轭梯度算法的配网抗差状态估计方法。

背景技术

配电网是指从输电网或地区发电厂接受电能，通过配电设施就地分配或按电压逐级分配给各类用户的电力网。配电网具有电压等级多，网络结构复杂，设备类型多样，使得配电网的运行状态多变。并且随着社会生产和人民生活对电能的需求和要求与日俱增，配电网发展迅速，小风电、小水电、光伏发电等分布式能源装机容量不断提高，对配电网的稳定运行提出了更高的要求。

电力系统状态估计是电力系统调度中心的能量管理系统(EMS)的核心功能之一，其功能是根据电力系统的各种量测信息，估计出电力系统当前的运行状态。现代电网的安全经济运行依赖于能量管理系统(EMS)，而能量管理系统的众多功能又可分成针对电网实时变化进行分析的在线应用和针对典型潮流断面进行分析的离线应用两大部分。电力系统状态估计可以说是大部分在线应用的高级软件的基础。如果电力系统状态估计结果不准确，后续的任何分析计算将不可能得到准确的结果。

配电网状态估计是利用冗余的量测量估计配电网系统的运行状态，配电网结构复杂，可再生能源并网、电动汽车入网等增加了配电网的不确定性，使得配电网运行和控制日益复杂，并且配电网支路繁多，只有少部分支路装有实时量测装置，大量数据需要通过伪量测来补充，伪量测数据通常精度不高，对传统最小二乘法状态估计的准确度和精度影响较大。在输电网中，为了获取准确的量测值，通常在重要节点或支路安装PMU同步向量量测单元，但是PMU价格昂贵目前尚不能在配电网中大量安装。为了在众多低精度数据进行状态估计，往往要求量测数据的个数远远大于状态变量个数，即要求数据冗余度不能太低，当量测数据较大时会造成求解时间过长的问题。

发明内容

本发明的目的是提供一种基于GPU并行加速共轭梯度算法的配网抗差状态估计方法，目的是辨识出不良数据并给出考虑量测数据不确定性的状态估计。

为了实现上述任务，本发明采用以下技术方案：

一种基于GPU并行加速共轭梯度算法的配网抗差状态估计方法，包括：

构建配电网系统的量测方程，包括状态变量和量测量；

利用量测量与量测计算量的量测残差，构建基于指数加权最小二乘模型；将测点分为正常测点和异常测点两类，同类型的测点采用不同准则赋予区间值，利用线性规划求解状态变量的边界信息，构建正常测点模型；

利用区间值分析所述正常测点模型求解结果的不确定性，构建考虑量测数据不确定性的区间分析状态估计模型；

求解所述基于指数加权最小二乘模型，将求解结果代入量测方程并计算量测残差，根据所述正常测点模型，对于异常测点采用量测计算值替代异常测点的实际量测值；求解所述区间分析状态估计模型，得到状态变量的不确定值区间；其中在指数加权最小二乘模型求解过程中，利用GPU并行加速的方法进行量测函数生成、量测雅克比矩阵生成以及信息矩阵的并行求解。

进一步地，所述构建配电网系统的量测方程，包括状态变量和量测量；

其中，状态变量包括电压幅值和相角，量测量包括节点注入有功功率和无功功率、支路测量有功功率和无功功率、电压幅值。

进一步地，所述利用量测量与量测计算量的量测残差，构建基于指数加权最小二乘模型，表示为：

J(X)＝[Z-h(X)]^TW[Z-h(X)]

其中，X为状态变量，Z为量测量，h(X)为量测方程向量，Z-h(X)为量测残差矢量；W为指数型权函数对角阵。

进一步地，所述指数型权函数对角阵W中，对角元素计算公式如下：

式中，z_i为量测量Z的第i个元素；h_i(X)为量测方程向量h(X)中的第i个量测方程；R_i为权矩阵R的第i个元素；(KR)_ii为矩阵(KR)的主对角线元素；σ为标准化残差的尺度参数，其初值的计算公式为：

σ²＝(Z-h(X₀))^TR(Z-h(X₀))/(m-n)

其中，X₀为X的初值；m为量测量Z的维数，n为状态变量的维数；K为残差灵敏度矩阵，计算公式为：

K＝E-H(H^TRH)^-1H^TR

其中

E为m阶单位阵。

进一步地，所述将测点分为正常测点和异常测点两类，同类型的测点采用不同准则赋予区间值，利用线性规划求解状态变量的边界信息，构建正常测点模型，包括：

利用扩展不确定度对基于指数加权最小二乘求出的状态估计对应测点进行分类，每个确定的状态变量X代入量测函数得到的量测计算值h(X)与实际测量值Z具有差值；有如下定义：

P(|z_i-h_i(X)|≤U_i)＝p

式中，U_i为量测点i在置信概率p下的扩展不确定度值，为预设值，P表示概率；

在给定置信区间的条件下，满足|z_i-h_i(X)|≤U_i定义为正常测点，不满足定义为异常测点。

进一步地，所述利用区间值分析所述正常测点模型求解结果的不确定性，构建考虑量测数据不确定性的区间分析状态估计模型，包括：

式中，x_i为X的第i个元素；x _i为x_i的下限，

为x_i的上限；Ω₁为异常测点集合，Ω₂为正常测点集合，h_i(X)为第i个量测方程，h_i(X₁)表示量测计算值。

进一步地，所述求解所述基于指数加权最小二乘模型，包括：

对基于指数加权最小二乘模型的目标函数求极小值应满足

为了使求解方便，忽略权函数对状态变量的偏导这一项，则有：

为了便于求解上述非线性方程组，将h(X)在X₀处一阶泰勒展开并略去余项，有h(X)＝h(X₀)+H(X₀)ΔX，代入式2H^TW[Z-h(X)]＝0，并令ΔZ＝Z-h(X₀)，得到线性方程组：

H^TWHΔX＝H^TWΔZ

将H^TWH进行LU数值分解，利用前推回代计算X的增量ΔX；令X＝X+ΔX，计算量测雅可比矩阵H、权矩阵W和量测矩阵h(X)，代入上述线性方程组，再次求解；如此循环，直到收敛。

进一步地，求解所述区间分析状态估计模型，包括：

在X^(k)处对h(X)进行泰勒展开，忽略高次项后可得h(X)表达式：

h(X)≈h(X^(k))+HΔX

式中，X^(k)为方程迭代区间近似解的中间值，即

利用量测函数的一阶泰勒展开公式将区间分析状态估计模型的修正方程转化为如下两个线性优化问题：

ΔZ _i分别为ΔZ的上下限；

Δx _i分别为x_i增量Δx_i的上下限；

上述两个线性优化问题采用线性规划方法进行求解，经过线性优化求解后可得

进一步地，利用GPU并行加速的方法进行量测函数生成，包括：

量测函数的计算可分为P_i、Q_i、P_ij、Q_ij四部分，每个部分利用一个内核函数完成相应的计算，一个线程负责一个量测功率，其中，P_i、Q_i表示节点注入有功功率和无功功率，P_ij、Q_ij表示支路量测有功功率和无功功率。

进一步地，利用GPU并行加速的方法进行量测雅可比矩阵生成，包括：

将获取的配电网节点支路数据进行SCR格式存储；所述数据包括支路阻抗、有功功率和无功功率；

量测雅可比矩阵的第一类子矩阵的稀疏性与节点导纳矩阵相同，第二类子矩阵的稀疏性与支路参数矩阵相同，找到各自非零元的对应关系，此过程在CPU中进行；

根据已经计算好的量测函数，在GPU中利用8个内核函数分别计算8个子矩阵中非零元素；

将非零元根据对应关系填入，得到完整的量测雅可比矩阵。

进一步地，利用GPU并行加速的方法进行信息矩阵的并行求解，包括：

对于求解过程中的H^TWH和H^TRH，首先计算H^TH，再与对角阵相乘，降低重复计算；计算H^TH时，若令H＝(H₁ ^T,H₂ ^T,…,H_m ^T)^T，则有

各行向量计算相互独立，适用于GPU加速；

在求非主对角元素时，利用CSR存储中的Indices向量分成的m段，将相邻两端中对应列号相同的Data元素相乘后相加便可得到。

一种终端设备，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，处理器执行计算机程序时实现前述基于GPU并行加速共轭梯度算法的配网抗差状态估计方法的步骤。

一种计算机可读存储介质，计算机可读存储介质存储有计算机程序，计算机程序被处理器执行时实现前述基于GPU并行加速共轭梯度算法的配网抗差状态估计方法的步骤。

与现有技术相比，本发明具有以下技术特点：

本发明方法能够得到比潮流计算可信度更高、精度更高的低压配电网运行状态，这种基于GPU并行加速的低压配电网区间分析状态估计方法减少了存储空间和状态估计的运行时间，同时降低不良数据对状态估计的影响，为低压配电网的运行状态提供一个可信度高的合理区间，给智能电网调度运行提供必要的技术支撑，解决配电网状态估计中大量低精度数据造成估计精度不高甚至估计结果严重偏离实际的问题，同时减少状态估计的计算时间，为配电网的控制与运行提供分析模型。

附图说明

图1为变压器等值电路示意图；

图2为CSR存储格式；

图3为本发明方法的流程示意图。

具体实施方式

参见图1，本发明提供了一种基于GPU并行加速共轭梯度算法的配网抗差状态估计方法，包括以下步骤：

S1低压配电网量测方程

配电网系统量测量有节点注入有功功率、节点注入无功功率、支路有功功率、支路无功功率、电压幅值等，这些量测量之间都满足基尔霍夫定律、欧姆定律或功率方程，可以选择不同的量测量作为状态估计的量测量，状态变量可以选取功率、节点电压、电流等。本发明选取电压幅值V和相角θ作为状态变量X，节点注入有功功率P_i和无功功率Q_i、支路量测有功功率P_ij和无功功率Q_ij、电压幅值V_i作为量测量Z，关系如下：

Z＝[P_i,Q_i,P_ij,Q_ij,V_i]

X＝[V,θ]

在上述状态变量X中，V表示电压幅值矩阵，V_i为矩阵里的元素，表示i节点的电压幅值；θ为电压相角矩阵，θ_i为矩阵里的元素，表示i节点的电压相角；而θ_ij＝θ_i-θ_j。

S1.1节点注入量测方程

节点注入功率方程包括节点注入有功功率方程和节点注入无功功率方程，具体如下：

式中，V_i和V_j分别为节点i和节点j的电压幅值；i,j＝1,2,...,n，n为节点个数；P_i和Q_i分别节点i的有功功率和无功功率；G_ij为系统节点导纳矩阵中节点i和节点j对应导纳实数部分；B_ij为系统节点导纳矩阵中节点i和节点j对应导纳虚数部分；θ_ij为节点i和节点j之间的相角差。

S1.2普通非变压器支路量测方程

对于普通非变压器线路，线路l_ij上i侧的有功功率和无功功率为：

P_ij＝-V_i ²G_ij+V_iV_j(G_ijcosθ_ij+B_ijsinθ_ij)

Q_ij＝V_i ²B_ij+V_iV_j(G_ijsinθ_ij-B_ijcosθ_ij)

S1.3变压器支路量测方程

对于变压器线路，变压器支路的等值电路与普通线路不同，变压器的主要参数是电抗和变比，等值电路如附图1所示，其中k是变压器的变比，i代表非标准侧，j代表标准侧，b_T表示标压器的电纳。

变压器支路i侧的有功功率和无功功率为：

变压器支路j侧的有功功率和无功功率为：

S2基于指数加权最小二乘模型的构建

利用基于指数加权二乘求出一组不受粗差影响的状态估计。

量测残差指量测量Z与量测计算向量h(X)的差值，采用指数加权最小二乘法状态估计，不仅能够在量测差值较小时保持最小二乘的简洁性，还能抑制残差很大时对状态估计的影响。基于指数加权最小二乘目标模型：

J(X)＝[Z-h(X)]^TW[Z-h(X)]

式中，X为n维状态变量；Z为m维量测量；h(X)为量测方程向量；Z-h(X)为量测残差矢量；W为m阶指数型权函数对角阵，对角元素计算公式如下：

式中，z_i为量测量Z的第i个元素；h_i(X)为量测方程向量h(X)的第i个元素，即第i个量测方程；R_i为权矩阵R的第i个元素；(KR)_ii为矩阵(KR)的主对角线元素；σ为标准化残差的尺度参数，其初值的计算公式为：

σ²＝(Z-h(X₀))^TR(Z-h(X₀))/(m-n)

其中，X₀为X的初值；K为残差灵敏度矩阵，计算公式为：

K＝E-H(H^TRH)^-1H^TR

其中

E为m阶单位阵。

S3正常测点模型的构建

根据正常测点的定义，将测点分为正常测点和异常测点两类，不同类型的测点采用不同准则赋予区间值，利用线性规划求解状态变量的边界信息。所述测点对应于量测量Z中每个元素数值的位置。

利用扩展不确定度对基于指数加权最小二乘求出的状态估计对应测点进行分类。每个确定的状态变量X代入量测函数h得到的量测计算值h(X)与实际测量值Z具有差值。由于量测装置本身的精度的局限，每个量测点的差值不能直接比较。借用扩展不确定度的定义，有如下定义：

P(|z_i-h_i(X)|≤U_i)＝p

式中，U_i为量测点i在置信概率p下的扩展不确定度值，为预设值。

在给定置信区间的条件下，满足|z_i-h_i(X)|≤U_i定义为正常测点，不满足定义为异常测点。

指数加权最小二乘模型，虽然考虑了粗差对状态估计结果的影响，但是无法排除所有不良数据对状态估计的影响。根据指数加权最小二乘模型的求解得到的状态变量X₁，对每个测点判断。正常测点则有-U_i≤z_i-h_i(X₁)≤U_i，对于异常测点，量测数据z_i对状态估计的准确度影响较大，选用量测计算值h_i(X₁)代替量测值。

S4考虑量测数据不确定性的区间分析状态估计模型的构建

运用正常测点模型将测点分成两类，由于无法穷举每种状态变量，所以异常测点只是没有理由拒绝定义为异常测点。考虑实际为正常测点却有可能被定义为异常测点对结果的影响，利用区间值分析上述模型求解结果的不确定性，具体如下：

根据上述给定置信概率下的正常测点模型，正常测点的不确定性可以描述为：

(z_i-h_i(X))∈[-U_i,U_i]

异常测点选用计算值h_i(X₁)代替量测值，则其不确定性可以描述为：

(h_i(X₁)-h_i(X))∈[-U_i,U_i]

即正常、异常测点量测值区间分别为[z_i-U_i,z_i+U_i]、[h_i(X₁)-U_i,h_i(X₁)+U_i]。根据量测值区间得到状态变量的不确定值区间，如下：

式中，x_i为X的第i个元素；x _i为x_i的下限，

为x_i的上限；Ω₁为异常测点集合，Ω₂为正常测点集合。

合理的区间模型可以提高低压配电网区间状态估计的可信性与精确性和降低状态估计结果区间保守性。区间状态估计得到电压幅值和相角的范围，表明配电网电压一定在所求的范围内，以便对配电网调整。

S5模型的求解

基于指数加权最小二乘区间分析状态估计模型的求解，主要分成三部分，第一步，利用牛顿-拉夫逊法将确定性非线性方程组转换为线性方程组，采用LU分解前推回代法求解线性方程组；第二步，将第一步得到的状态变量代入量测方程，根据定义对测点分类，不同类型测点根据不同的准则赋予区间值；第三步，采用类似于牛顿-拉夫逊法求解确定性非线性方程组的处理，将含非线性区间约束的超定方程组求解问题转换成相应的线性修正方程组，运用线性规划迭代求解。

S5.1基于指数加权最小二乘模型的求解

对基于指数加权最小二乘模型的目标函数求极小值应满足

为了使求解方便，忽略权函数对状态变量的偏导这一项，则有：

为了便于求解上述非线性方程组，将h(X)在X₀处一阶泰勒展开并略去余项，有h(X)＝h(X₀)+H(X₀)ΔX，代入式2H^TW[Z-h(X)]＝0，并令ΔZ＝Z-h(X₀)，得到线性方程组：

H^TWHΔX＝H^TWΔZ

将H^TWH进行LU数值分解，利用前推回代计算X的增量ΔX。令X＝X+ΔX，计算量测雅可比矩阵H、权矩阵W和量测矩阵h(X)，代入上述线性方程组，再次求解。如此循环，直到收敛。

S5.2正常测点模型的实施

将基于指数加权最小二乘模型的求解结果代入量测方程并求解残差ΔZ＝Z-h(X)，根据正常测点模型，用量测计算值h_i(X)替代异常测点的实际测量值z_i。

S5.3考虑量测数据不确定性的区间分析状态估计模型的求解

对于一个含n节点主动配电网，需要求解2n个带非线性约束的优化问题，建立2m个方程，由于低压配电网状态估计要求远m大于n，不利于非线性方程组求解。类似于牛顿-拉夫逊求解确定性非线性方程组，在X^(k)处对h(X)进行泰勒展开，忽略高次项后可得h(X)表达式：

h(X)≈h(X^(k))+HΔX

式中，X^(k)为方程迭代区间近似解的中间值，即

X⁽⁰⁾采用基于指数加权二乘模型得到的状态变量X。

利用量测函数的一阶泰勒展开公式将区间分析模型的修正方程转化为如下两个线性优化问题：

ΔZ _i分别为ΔZ的上下限。

上述两个线性优化问题采用简单线性规划方法进行求解，经过2n次线性优化求解后可得

S6 GPU加速

随着低压配电网规模的扩大和愈加复杂的计算，传统的在线计算不能满足实时性的要求。近年来GPU(Graphics Processing Unit)价格的下降使得利用GPU加速状态估计成为可能。本方案中状态估计GPU加速部分主要在量测函数生成、量测雅可比矩阵生成、信息矩阵求解过程。

S6.1稀疏矩阵存储

对量测雅可比矩阵是一个非对称的稀疏矩阵，而信息矩阵是一个对称的系数矩阵，在并行计算前通过压缩只存储非零元素节省内存，降低稀疏矩阵在计算过程中的内存消耗。目前最为常用的压缩格式主要有坐标压缩法(Coordinate Storage Scheme，COO)、对角压缩法(Diagonal Storage Scheme，DIA)和行格式压缩法(Compressed Row StorageScheme，CSR)等，三种格式中，COO格式简单但消耗内存最大，DIA格式适用于对角元素为非零元素的稀疏矩阵，对于无规律的稀疏矩阵，CSR较为有效。综合考虑，选用CSR作为量测雅可比矩阵和信息矩阵的存储方法。

CSR存储m维稀疏矩阵A共需3个一维向量完成，第一个向量Data用来存储A中的非零元素，按行可分为m段，第二个向量Indices存储Data中各个元素的列号，第三个向量Ptr用于存储Data中m个段的段首在系数矩阵A中非零元元素的顺序号，示例如附图2所示。

S6.2量测函数生成

与量测雅可比矩阵和信息矩阵不同，量测函数h具有稠密性，且各个量测点的有功功率和无功功率计算相互独立，利于GPU的并行计算。量测函数的计算可分为P_i、Q_i、P_ij、Q_ij四部分，每个部分利用一个内核函数完成相应的计算，一个线程负责一个量测功率。在计算P_ij、和Q_ij时，由于变压器支路与一般支路的计算公式不同，所以在计算时需要加以判断。以节点注入有功功率P_i为例算法如表1所示：

S6.3量测雅可比矩阵生成

鉴于潮流计算的雅可比矩阵分块形式，可将量测雅可比矩阵修改形式修改如下：

其中K、N、M、L与潮流计算类似，K1、N1、M1、L1是根据增加的量测支路得到。在进行并行运算内核函数设计时，可将计算公式较为类似的K、N、M、L归作第一类，K1、N1、M1、L1作为第二类。普通支路和变压器支路的功率方程不同，具体如下：

上式中的前两个式子为矩阵K的具体表达式，后两个式子为矩阵N的具体表达式。

上式中的前两个式子为矩阵M的具体表达式，后两个式子为矩阵L的具体表达式。

上式是针对普通支路，前两个式子为矩阵K1的具体表达式，后两个式子为矩阵N1的具体表达式。

上式是针对变压器支路，前两个式子为矩阵K1的具体表达式，后两个式子为矩阵N1的具体表达式。

上式子是针对普通支路，前两个式子为矩阵M1的具体表达式，后两个式子为矩阵L1的具体表达式。

上式子是针对变压器支路，前两个式子为矩阵M1的具体表达式，后两个式子为矩阵L1的具体表达式。

GPU加速的量测雅可比矩阵并行生成算法过程如下：

1)将获取的配电网节点支路数据进行SCR格式存储；所述数据包括支路阻抗、有功功率和无功功率。

2)量测雅可比矩阵H的第一类子矩阵的稀疏性与节点导纳矩阵相同，第二类子矩阵的稀疏性与支路参数矩阵相同，找到各自非零元的对应关系，此过程在CPU中进行；

3)根据上已经计算好的量测函数h，在GPU中利用8个内核函数分别计算8个子矩阵中非零元素；

4)将非零元根据对应关系填入H，得到完整的量测雅可比矩阵。

S6.4信息矩阵H^TWH的并行求解

信息矩阵H^TWH是一个对称矩阵且具有很高的稀疏性，虽然H^TWH是一个n阶矩阵，但分析计算公式知，H是一个m*n阶的矩阵，计算过程涉及维数m的计算，由于m远远大于n，所以计算量和内存都非常巨大。

信息矩阵的生成过程的计算量集中在残差灵敏度矩阵K和H^TWH的计算，分析可知，权矩阵W的求解只需K的主对角线元素，所以在求解过程中只需要求解m个元素值，无需求出所有的m*m个元素值。在求解过程中出现了H^TWH和H^TRH，由于W和R都是m阶对角矩阵，所以可以首先计算H^TH，再与对角阵相乘，降低重复计算。计算H^TH时，若令H＝(H₁ ^T,H₂ ^T,…,H_m ^T)^T，则有

各行向量计算相互独立，适用于GPU加速。

H^TH为大型稀疏矩阵相乘，主对角线元素可以通过H列向量的非零元素平方和求得，在求非主对角元素时，利用CSR中的Indices向量分成的m段，将相邻两端中对应列号相同的Data元素相乘后相加便可得到。

本申请实施例进一步提供一种终端设备，该终端设备可以为计算机、服务器；包括存储器、处理器以及存储在存储器中并可在处理器上运行的计算机程序，处理器执行计算机程序时实现上述基于GPU并行加速共轭梯度算法的配网抗差状态估计方法的步骤。

计算机程序也可以被分割成一个或多个模块/单元，一个或者多个模块/单元被存储在存储器中，并由处理器执行，以完成本申请。一个或多个模块/单元可以是能够完成特定功能的一系列计算机程序指令段，该指令段用于描述计算机程序在终端设备中的执行过程。

本申请的实施提供了一种计算机可读存储介质，计算机可读存储介质存储有计算机程序，计算机程序被处理器执行时实现上述基于GPU并行加速共轭梯度算法的配网抗差状态估计方法的步骤。

以上实施例仅用以说明本申请的技术方案，而非对其限制；尽管参照前述实施例对本申请进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本申请各实施例技术方案的精神和范围，均应包含在本申请的保护范围之内。

Claims (10)

1.一种基于GPU并行加速共轭梯度算法的配网抗差状态估计方法，其特征在于，包括：

构建配电网系统的量测方程，包括状态变量和量测量；

利用量测量与量测计算量的量测残差，构建基于指数加权最小二乘模型；

将测点分为正常测点和异常测点两类，同类型的测点采用不同准则赋予区间值，利用线性规划求解状态变量的边界信息，构建正常测点模型；

利用区间值分析所述正常测点模型求解结果的不确定性，构建考虑量测数据不确定性的区间分析状态估计模型；

求解所述基于指数加权最小二乘模型，将求解结果代入量测方程并计算量测残差，根据所述正常测点模型，对于异常测点采用量测计算值替代异常测点的实际量测值；求解所述区间分析状态估计模型，得到状态变量的不确定值区间；其中在指数加权最小二乘模型求解过程中，利用GPU并行加速的方法进行量测函数生成、量测雅克比矩阵生成以及信息矩阵的并行求解。

2.根据权利要求1所述的基于GPU并行加速共轭梯度算法的配网抗差状态估计方法，其特征在于，所述构建配电网系统的量测方程，包括状态变量和量测量；

其中，状态变量包括电压幅值和相角，量测量包括节点注入有功功率和无功功率、支路测量有功功率和无功功率、电压幅值。

3.根据权利要求1所述的基于GPU并行加速共轭梯度算法的配网抗差状态估计方法，其特征在于，所述利用量测量与量测计算量的量测残差，构建基于指数加权最小二乘模型，表示为：

J(X)＝[Z-h(X)]^TW[Z-h(X)]

其中，X为状态变量，Z为量测量，h(X)为量测方程向量，Z-h(X)为量测残差矢量；W为指数型权函数对角阵。

4.根据权利要求3所述的基于GPU并行加速共轭梯度算法的配网抗差状态估计方法，其特征在于，所述指数型权函数对角阵W中，对角元素计算公式如下：

式中，z_i为量测量Z的第i个元素；h_i(X)为量测方程向量h(X)中的第i个量测方程；R_i为权矩阵R的第i个元素；(KR)_ii为矩阵(KR)的主对角线元素；σ为标准化残差的尺度参数，其初值的计算公式为：

σ²＝(Z-h(X₀))^TR(Z-h(X₀))/(m-n)

其中，X₀为X的初值；m为量测量Z的维数，n为状态变量的维数；K为残差灵敏度矩阵，计算公式为：

K＝E-H(H^TRH)^-1H^TR

其中

E为m阶单位阵。

5.根据权利要求1所述的基于GPU并行加速共轭梯度算法的配网抗差状态估计方法，其特征在于，所述将测点分为正常测点和异常测点两类，同类型的测点采用不同准则赋予区间值，利用线性规划求解状态变量的边界信息，构建正常测点模型，包括：

利用扩展不确定度对基于指数加权最小二乘求出的状态估计对应测点进行分类，每个确定的状态变量X代入量测函数得到的量测计算值h(X)与实际测量值Z具有差值；有如下定义：

P(|z_i-h_i(X)|≤U_i)＝p

式中，U_i为量测点i在置信概率p下的扩展不确定度值，为预设值，P表示概率；

在给定置信区间的条件下，满足|z_i-h_i(X)|≤U_i定义为正常测点，不满足定义为异常测点。

6.根据权利要求1所述的基于GPU并行加速共轭梯度算法的配网抗差状态估计方法，其特征在于，所述利用区间值分析所述正常测点模型求解结果的不确定性，构建考虑量测数据不确定性的区间分析状态估计模型，包括：

式中，x_i为X的第i个元素；x _i为x_i的下限，

为x_i的上限；Ω₁为异常测点集合，Ω₂为正常测点集合，h_i(X)为第i个量测方程，h_i(X₁)表示量测计算值。

7.根据权利要求1所述的基于GPU并行加速共轭梯度算法的配网抗差状态估计方法，其特征在于，所述求解所述基于指数加权最小二乘模型，包括：

对基于指数加权最小二乘模型的目标函数求极小值应满足

为了使求解方便，忽略权函数对状态变量的偏导这一项，则有：

为了便于求解上述非线性方程组，将h(X)在X₀处一阶泰勒展开并略去余项，有h(X)＝h(X₀)+H(X₀)ΔX，代入式2H^TW[Z-h(X)]＝0，并令ΔZ＝Z-h(X₀)，得到线性方程组：

H^TWHΔX＝H^TWΔZ

将H^TWH进行LU数值分解，利用前推回代计算X的增量ΔX；令X＝X+ΔX，计算量测雅可比矩阵H、权矩阵W和量测矩阵h(X)，代入上述线性方程组，再次求解；如此循环，直到收敛。

8.根据权利要求1所述的基于GPU并行加速共轭梯度算法的配网抗差状态估计方法，其特征在于，求解所述区间分析状态估计模型，包括：

在X^(k)处对h(X)进行泰勒展开，忽略高次项后可得h(X)表达式：

h(X)≈h(X^(k))+HΔX

式中，X^(k)为方程迭代区间近似解的中间值，即

利用量测函数的一阶泰勒展开公式将区间分析状态估计模型的修正方程转化为如下两个线性优化问题：

ΔZ _i分别为ΔZ的上下限；

Δx _i分别为x_i增量Δx_i的上下限；

上述两个线性优化问题采用线性规划方法进行求解，经过线性优化求解后可得

9.根据权利要求1所述的基于GPU并行加速共轭梯度算法的配网抗差状态估计方法，其特征在于，利用GPU并行加速的方法进行量测雅可比矩阵生成，包括：

将获取的配电网节点支路数据进行SCR格式存储；所述数据包括支路阻抗、有功功率和无功功率；

量测雅可比矩阵的第一类子矩阵的稀疏性与节点导纳矩阵相同，第二类子矩阵的稀疏性与支路参数矩阵相同，找到各自非零元的对应关系，此过程在CPU中进行；

根据已经计算好的量测函数，在GPU中利用8个内核函数分别计算8个子矩阵中非零元素；

将非零元根据对应关系填入，得到完整的量测雅可比矩阵。

10.根据权利要求1所述的基于GPU并行加速共轭梯度算法的配网抗差状态估计方法，其特征在于，利用GPU并行加速的方法进行信息矩阵的并行求解，包括：

对于求解过程中的H^TWH和H^TRH，首先计算H^TH，再与对角阵相乘，降低重复计算；计算H^TH时，若令H＝(H₁ ^T,H₂ ^T,…,H_m ^T)^T，则有

各行向量计算相互独立，适用于GPU加速；

在求非主对角元素时，利用CSR存储中的Indices向量分成的m段，将相邻两端中对应列号相同的Data元素相乘后相加便可得到。

Images