基于GPU并行计算的电力系统PQ分解状态估计方法与流程

本发明涉及电力系统状态估计的技术领域，更具体地，涉及一种基于gpu并行计算的电力系统pq分解状态估计方法。

背景技术：

电力系统状态估计的功能是根据电力系统的各种量测信息，估计出电力系统当前的运行状态。电力系统状态估计是电网调度系统中各大高级应用的基础，更是维持电力系统稳定、确保电力系统安全高效运行的基础，因此，电力系统状态估计的计算响应速度，决定了电网调度系统能否及时地反映电网的真实状态，也影响着其他电网自动控制系统的控制效果，对整个电力系统的稳定运行至关重要。

目前，电力调度系统常采用最小二乘法以及pq分解法对电力系统进行状态估计。在以上两种方法的计算过程中，存在大量相互独立的计算过程以及矩阵运算，无法保证电力调度系统具有及时反映电网真实状态的能力，为提高电力调度系统反映电网真实状态的实时性，在状态估计的计算过程中进行并行化处理十分有必要，公开号为cn1070696a的中国专利在2017年8月公开了一种电力系统状态估计的并行计算方法，整合建立状态估计并行模型与并行分解两个过程，克服了现有采用优化算法进行状态估计时的复杂性，此方法方便结合相关并行平台实现并行化处理，可以提高状态估计计算响应速度，但该专利提出的并行计算方法涉及两个过程，过程也比较复杂，而且并没有统一封装整合改进。

技术实现要素：

为克服传统电力系统状态估计的方法具有无法保证电力调度系统及时反映电网真实状态的缺陷，而目前公开的电力系统状态估计并行计算的方法涉及的过程比较复杂，且存在未统一封装整合改进的弊端；本发明提出一种基于gpu并行计算的电力系统pq分解状态估计方法，提高电力系统状态估计的速度，保证电力调度系统及时反映电网的真实运行状态，确保电力系统安全高效运行。

为了达到上述技术效果，本发明的技术方案如下：

本发明提出一种基于gpu并行计算的电力系统pq分解状态估计方法，至少包括：

s1.采集电力系统的量测参数和电力系统的网络参数；

s2.根据电力系统的网络参数生成节点导纳矩阵及支路导纳列表；

s3.以电力系统的节点电压幅值及电压相角作为状态变量向量，对状态变量向量初始化，将节点电压幅值初始化为1.0pu，节点电压相角初始化为0rad；

s4.根据电力系统的量测参数和电力系统的网络参数，分别计算有功雅各比矩阵hp及无功雅各比矩阵hq；

s5.计算有功残差方程矩阵ap和无功残差方程矩阵aq，对有功残差方程矩阵ap和无功残差方程矩阵aq分别进行llt分解得到有功下三角矩阵apl及无功下三角矩阵aql；

s6.将有功雅各比矩阵hp、无功雅各比矩阵hq、有功下三角矩阵apl、无功下三角矩阵aql及支路导纳列表复制到gpu的主内存中；

s7.基于gpu并行计算，在gpu中进行有功迭代，判断有功迭代收敛条件是否满足，若是，执行步骤s8，否则，更新节点电压相角，执行步骤s8；

s8.在gpu中进行无功迭代，判断无功迭代收敛条件是否满足，若是，执行步骤s9；否则，更新节点电压幅值，执行步骤s9；

s9.判断有功迭代及无功迭代是否均收敛，若是，输出状态估计结果；否则，执行步骤s10；

s10.判断迭代次数kk是否大于最大迭代次数kth，若是，状态估计迭代不收敛；否则，返回执行步骤s7。

在此，电力系统的量测参数通过数据采集系统采集,电力系统的网络参数通过能量管理系统(ems)采集，步骤s2所述的根据电力系统的网络参数生成节点导纳矩阵按传统计算方法生成；另外，对有功残差方程矩阵ap和无功残差方程矩阵aq分别进行llt分解是比较成熟的传统技术；在整个状态估计中，如果在有功迭代收敛条件及无功迭代收敛条件均不满足时，将当前的迭代次数与设置的最大迭代次数相比，若当前的迭代次数已经超过了最大迭代次数，但状态估计过程仍不满足有功迭代收敛条件及无功迭代收敛条件，则表示状态估计过程不收敛；在有功迭代收敛条件及无功迭代收敛条件均满足时，输出的状态估计结果为电力系统的节点电压幅值及电压相角数据，即状态变量向量的当前值。

优选地，步骤s1所述电力系统的量测参数包括：节点电压u、节点注入有功功率pinj_m、节点注入无功功率qinj_m、支路有功功率pft_m、ptf_m支路无功功率qft_m、qtf_m及电流幅值；通过量测参数形成有功测值向量zp_m、无功量测值向量zq_m，以及有功量测权矩阵rp和无功量测权矩阵rq；采集的无效或实际未经传送的量测参数设置为0，权系数也为0。

所述电力系统的网络参数包括：节点数nodenum、支路数量lnnum、支路首节点编号f、支路末节点编号t、支路电阻r、支路电抗x、支路非标准变比k、对地导纳b、节点并联电导gsh、节点并联电纳bsh，无并联元件节点的节点并联电导gsh及节点并联电纳bsh为0。

在此，采集的网络参数里的节点电压u信息包括节点电压的幅值和电压相角；另外，为适应后续在gpu中进行的并行计算，考虑cpu内存以128个浮点数或整数为单位对参数进行读取，因此步骤s1采集的参数数据均作对齐处理：将量测量向量的长度扩展为128的整数倍，节点数nodenum及支路数量lnnum也扩展为128的整数倍。

步骤s2所述的支路导纳列表由一侧至另一侧的导纳存放结构依次为：支路导纳列：ys[1]～ys[lnnum]、支路f侧对地导纳列：yf[1]～yf[lnnum]及支路t侧对地导纳列：yt[1]～yt[lnnum]；

支路导纳列中任意一个支路k导纳ys(k)的求解公式为：

其中，ys(k)表示该支路导纳，r(k)表示该支路电阻，x(k)表示该支路电抗；k(k)表示该支路非标准变比；gs(k)为该支路电导；bs(k)为该支路电纳；j为复数计算中虚部标识；

支路f侧对地导纳列中第k个支路f侧对地导纳yf(k)的求解公式为：

其中，b(k)表示该线路对地电容，gf(k)为该支路f侧对地电导；bf(k)为该支路f侧对地电纳；

支路t侧对地导纳列中第k个支路t侧对地导纳yt(k)的求解公式为：

其中，gt(k)为该支路t侧对地电导；bt(k)为该支路t侧对地电纳。

在此，为了提升gpu获取参数数据的速度，防止gpu各线程之间的访问冲突，保证gpu并行计算的速度，将支路导纳列表按照由一侧至另一侧的导纳存放结构依次存放，实际不存在的支路的导纳参数为0。

优选地，步骤s3所述的状态变量向量的长度为：nodenum*2。

优选地，步骤s4所述的有功雅各比矩阵hp为nodenum+lnnum*2行nodenum列的稀疏矩阵，有功雅各比矩阵hp中的元素计算公式为：

θft＝θf(k)-θt(k)

其中，pinj(i)为点i的注入有功功率；bij表示节点i与节点j之间的支路电导，bii表示节点i的自电导、θi为节点i的电压相角、θj为节点j的电压相角；uf(k)为支路k起始节点的电压幅值；ut(k)为支路k末节点的电压幅值；θf(k)支路k起始节点的电压相角，θt(k)为支路k末节点的电压相角。f(k)为支路k起始节点编号、t(k)为支路k末节点编号,θft表示节点i与节点j的电压相角差。

无功雅各比矩阵hq为nodenum*2+lnnum*2行nodenum列的稀疏矩阵，无功雅各比矩阵hq中的元素计算公式为：

θft＝θf(k)-θt(k)

其中，qinj(i)为节点i的注入无功功率，bij表示节点i与节点j之间的支路电导，bii表示节点i的自电导、ui为节点i的电压幅值、uj为节点j的电压幅值；uf(k)为支路k起始节点的电压幅值；ut(k)为支路k末节点的电压幅值；θf(k)支路k起始节点的电压相角，θt(k)为支路k末节点的电压相角。f(k)为支路k起始节点编号、t(k)为支路k末节点编号,θft表示节点i与节点j的电压相角差。

优选地，步骤s5所述的有功残差方程矩阵ap的计算公式为：

其中，为有功雅各比矩阵hp的转置，为有功量测权矩阵rp的逆矩阵；无功残差方程矩阵aq的计算公式为：

其中，为无功雅各比矩阵hq的转置，为无功功率量测权矩阵rq的逆矩阵。

在此，由于有功量测权矩阵rp与无功量测权矩阵rq均为对角阵，而且在采集电力系统量测参数时，采集的无效或实际未经传送的量测参数设置为0，因此，有功量测权矩阵rp与无功量测权矩阵rq中有将部分对角元设为0，可能存在无逆矩阵的情况，对有功量测权矩阵rp与无功量测权矩阵rq中非零元素取倒数来等效求逆。

优选地，步骤s6所述的有功雅各比矩阵hp、无功雅各比矩阵hq、有功下三角矩阵apl、无功下三角矩阵aql复制到gpu的主内存时采用列优先稀疏矩阵压缩存储格式；支路导纳列表复制到gpu的主内存中采用列优先的存储格式，当有功雅各比矩阵hp、无功雅各比矩阵hq、有功下三角矩阵apl、无功下三角矩阵aql及支路导纳列表中任意一参数列的行数不等于gpu线程数的整数倍时，插入全为0的参数空行补全，从而保证触发gpu对复制到主内存的参数的合并访问。

优选地，步骤s7所述的有功迭代的过程为：

s701.计算有功功率量测估计值；

首先计算支路有功量测估计值，支路有功功率量测估计值的计算公式为：

θft＝θf(k)-θt(k)

其中，pft_est(k)为该支路由f侧到t侧支路有功功率量测估计值；ptf_est(k)为该支路由t侧到f侧支路有功功率量测估计值；uf(k)表示该支路f侧的节点电压幅值；ut(k)表示该支路t侧的节点电压幅值；θf(k)支路k起始节点的电压相角，θt(k)为支路k末节点的电压相角,θft表示节点i与节点j的电压相角差。

s702.遍历所有节点，对于任意点i，设置有功注入量测估计值变量pinj_est(i)，并初始化为0。

s703.遍历所有支路，对于任意支路k，将支路两侧有功量测估计值叠加到量测节点的注入有功功率量测估计值。例如，对于第k个支路

pinj_est(f(k))＝pinj_est(f(k))+pft_est(k)

pinj_est(t(k))＝pinj_est(t(k))+ptf_est(k)

其中，pft_est(f)为该支路f侧节点的节点注入有功量测估计值，ptf_est(t)为该支路t侧节点的节点注入有功量测估计值。

s704.遍历所有节点，将节点并联导纳有功量测估计值叠加到节点注入有功量测估计值，对于任意节点i：

其中，gsh(i)为节点i并联电导

s705.将pinj_est、pft_est、ptf_est依次排列形成有功量测估计值向量zp_est；利用有功量测估计值向量zp_est与有功测值向量zp_m求取有功量测估计误差向量：

zp_err＝zp_m-zp_est

其中，zp_err表示有功量测估计误差向量；

s706.根据有功残差方程矩阵ap及有功量测估计误差向量zp_err，求取状态变量θ的变化量δθ，状态变量θ的变化量δθ满足：

优选地，步骤s7所述的无功迭代的过程为：

s711.计算无功功率量测估计值；

首先计算支路无功量测估计值，支路无功功率量测估计值的计算公式为：

θft＝θf(k)-θt(k)

其中，qft_est(f)为该支路f侧节点的节点注入无功量测估计值，qtf_est(t)为该支路t侧节点的节点注入无功量测估计值。uf(k)表示该支路f侧的节点电压幅值；ut(k)表示该支路t侧的节点电压幅值；θf(k)支路k起始节点的电压相角，θt(k)为支路k末节点的电压相角,θft表示节点i与节点j的电压相角差。

s712.遍历所有节点，对于任意点i，设置有功注入量测估计值变量qinj_est(i)，并初始化为0。

s713.遍历所有支路，对于任意支路k，将支路两侧无功量测估计值叠加到量测节点的注入无功功率量测估计值。例如，对于第k个支路

qinj_est(f(k))＝qinj_est(f(k))+qft_est(k)

qinj_est(t(k))＝qinj_est(t(k))+qtf_est(k)

其中，qft_est(f)为该支路f侧节点的节点注入无功量测估计值，qtf_est(t)为该支路t侧节点的节点注入无功量测估计值。

s714.遍历所有节点，将节点并联导纳无功量测估计值叠加到节点注入无功量测估计值，对于任意节点i：

其中，bsh(i)为节点i并联电纳。

s715.将qinj_est、qft_est、qtf_est、u依次排列形成无功量测估计值向量zq_est；利用有功量测估计值向量zq_est与有功测值向量zq_m求取有功量测估计误差向量：

zp_err＝zp_m-zp_est

其中，zp_err表示有功量测估计误差向量；

s716.根据有功残差方程矩阵aq及有功量测估计误差向量zq_err，求取状态变量u的变化量δu，状态变量u的变化量δu满足：

在此，无论是有功迭代还是无功迭代的过程中，在计算有功或无功功率量测估计值时，由于gpu各线程以单指令多数据流的形式实现并行计算，每个节点连接的支路数量不一样，用传统按节点计算节点注入有功或无功功率量测估计值的方法会导致各个线程之间的等待，严重影响了并行度，所以首先计算支路有功、无功量测估计值，以及所有节点的并联导纳有功、无功量测估计值，对支路功率求和来得到节点注入功率，提高计算的并行度；另外，在状态变量u的变化量δu或状态变量θ的变化量δθ的求解过程中，状态估计是基于pq分解法的，这时有功残差方程矩阵ap及无功残差方程矩阵aq均为常矩阵，在迭代过程中不发生变化，因此利用llt分解分别得到的下三角阵apl、aql前推回代，进一步得到δθ、δu，以简化求解过程。

优选地，有功迭代收敛条件为：

max(δθ)<θth

其中,δθ表示状态变量θ的变化量；θth表示有功迭代收敛阈值；当有功迭代收敛条件不满足时，更新节点电压相角，节点电压相角的更新公式为：

θ^(kk+1)＝θ^(kk)+δθ

其中，θ^(kk)表示第kk次迭代时的节点电压相角状态变量；θ^(kk+1)表示第kk+1次迭代时的节点电压相角状态变量；无功迭代收敛条件为：

max(δu)<uth

其中,δu表示状态变量节点电压幅值u的变化量；uth表示无功迭代收敛阈值；当无功迭代收敛条件不满足时，更新节点电压幅值，节点电压幅值的更新公式为：

u(^kk+1)＝u(^kk)+δu

其中，u^(kk)表示第kk次迭代时的节点电压幅值状态变量；u^(kk+1)表示第kk+1次迭代时的节点电压幅值状态变量。

与现有技术相比，本发明技术方案的有益效果是：

本发明提出一种基于gpu并行计算的电力系统pq分解状态估计方法，基于传统电力系统pq分解法的状态估计，将状态估计过程中涉及到的参数统一复制封装到gpu的主内存中，触发gpu对复制到主内存的参数的合并访问，然后基于gpu进行并行计算，克服传统电力系统状态估计的方法具有无法保证电力调度系统及时反映电网真实状态的缺陷，提高电力系统状态估计的速度，保证电力调度系统及时反映电网的真实运行状态，确保电力系统安全高效运行。

附图说明

图1表示本发明提出的基于gpu并行计算的电力系统pq分解状态估计方法的流程示意图；

图2表示本发明实施例中提出的支路导纳列表封装结构图。

具体实施方式

附图仅用于示例性说明，不能理解为对本专利的限制；

为了更好地说明本实施例，附图某些部位会有省略、放大或缩小，并不代表实际尺寸；

对于本领域技术人员来说，附图中某些公知内容说明可能省略是可以理解的。

下面结合附图和实施例对本发明的技术方案做进一步的说明。

实施例1

图1为本发明提出的基于gpu并行计算的电力系统pq分解状态估计方法的流程示意图，参见图1，其中，图示中hp表示有功雅各比矩阵；hq表示无功雅各比矩阵；ap表示有功残差方程矩阵；aq表示无功残差方程矩阵；apl表示有功下三角矩阵；aql表示无功下三角矩阵，所述方法包括：

s1.采集电力系统的量测参数和电力系统的网络参数；在具体实施时，电力系统的量测参数通过数据采集系统采集,电力系统的网络参数通过能量管理系统(ems)采集；

电力系统的量测参数包括：节点电压u、节点注入有功功率pinj_m、节点注入无功功率qinj_m、支路两侧有功功率pft_m、ptf_m支路两侧无功功率qft_m、qtf_m及电流幅值，采集的无效或实际未经传送的量测参数设置为0，权系数也为0；通过量测参数形成有功测值向量zp_m、无功量测值向量zq_m，以及有功量测权矩阵rp和无功量测权矩阵rq；采集的网络参数里的节点电压u信息包括节点电压的幅值和电压相角；

所述电力系统的网络参数包括：节点数nodenum、支路数量lnnum、支路首节点编号f、支路末节点编号t、支路电阻r、支路电抗x、支路非标准变比k、对地导纳b、节点并联电导gsh、节点并联电纳bsh，无并联元件节点的节点并联电导gsh及节点并联电纳bsh为0。另外，为适应后续具体实施时在gpu中进行的并行计算，考虑cpu内存以128个浮点数或整数为单位对参数进行读取，因此步骤s1采集的参数数据均作对齐处理：将量测量向量的长度扩展为128的整数倍，节点数nodenum及支路数量lnnum也扩展为128的整数倍。

s2.根据电力系统的网络参数生成节点导纳矩阵及支路导纳列表；在本实施例中，生成节点导纳矩阵按传统计算方法生成；

s3.以电力系统的节点电压幅值及电压相角作为状态变量向量，对状态变量向量初始化，将节点电压幅值初始化为1.0pu，节点电压相角初始化为0rad；

s4.根据电力系统的量测参数和电力系统的网络参数，分别计算有功雅各比矩阵hp及无功雅各比矩阵hq；

s5.计算有功残差方程矩阵ap和无功残差方程矩阵aq，对有功残差方程矩阵ap和无功残差方程矩阵aq分别进行llt分解得到有功下三角矩阵apl及无功下三角矩阵aql；

s6.将有功雅各比矩阵hp、无功雅各比矩阵hq、有功下三角矩阵apl、无功下三角矩阵aql及支路导纳列表复制到gpu的主内存中；

s7.基于gpu并行计算，在gpu中进行有功迭代，判断有功迭代收敛条件是否满足，若是，执行步骤s8，否则，更新节点电压相角，执行步骤s8；

s8.在gpu中进行无功迭代，判断无功迭代收敛条件是否满足，若是，执行步骤s9；否则，更新节点电压幅值，执行步骤s9；

s9.判断有功迭代及无功迭代是否均收敛，若是，输出状态估计结果；否则，执行步骤s10；

s10.判断迭代次数kk是否大于最大迭代次数kth，若是，状态估计迭代不收敛；否则，返回执行步骤s7。

在整个状态估计中，如果在有功迭代收敛条件及无功迭代收敛条件均不满足时，将当前的迭代次数与设置的最大迭代次数相比，若当前的迭代次数已经超过了最大迭代次数，但状态估计过程仍不满足有功迭代收敛条件及无功迭代收敛条件，则表示状态估计过程不收敛；在有功迭代收敛条件及无功迭代收敛条件均满足时，输出的状态估计结果即状态变量向量的当前值，为电力系统的节点电压幅值及电压相角数据。

在本实施例中，支路导纳列表封装结构图如图2所示，为了提升gpu获取参数数据的速度，防止gpu各线程之间的访问冲突，保证gpu并行计算的速度，支路导纳列表由一侧至另一侧的导纳存放结构依次为：支路导纳列：ys[1]～ys[lnnum]、支路f侧对地导纳列：yf[1]～yf[lnnum]及支路t侧对地导纳列：yt[1]～yt[lnnum]；实际不存在的支路的导纳参数为0。

支路导纳列中任意一个支路k导纳ys(k)的求解公式为：

其中，ys(k)表示该支路导纳，r(k)表示该支路电阻，x(k)表示该支路电抗；k(k)表示该支路非标准变比；gs(k)为该支路电导；bs(k)为该支路电纳；j为复数计算中虚部标识；

支路f侧对地导纳列中第k个支路f侧对地导纳yf(k)的求解公式为：

其中，b(k)表示该线路对地电容，gf(k)为该支路f侧对地电导；bf(k)为该支路f侧对地电纳；

支路t侧对地导纳列中第k个支路t侧对地导纳yt(k)的求解公式为：

其中，gt(k)为该支路t侧对地电导；bt(k)为该支路t侧对地电纳。

步骤s3所述的状态变量向量的长度为：nodenum*2。

在本实施例中，步骤s4所述的有功雅各比矩阵hp为(nodenum+lnnum*2)行nodenum列的稀疏矩阵，有功雅各比矩阵hp中的元素计算公式为：

θft＝θf(k)-θt(k)

其中，pinj(i)为点i的注入有功功率；bij表示节点i与节点j之间的支路电导，bii表示节点i的自电导、θi为节点i的电压相角、θj为节点j的电压相角；uf(k)为支路k起始节点的电压幅值；ut(k)为支路k末节点的电压幅值；θf(k)支路k起始节点的电压相角，θt(k)为支路k末节点的电压相角。f(k)为支路k起始节点编号、t(k)为支路k末节点编号,θft表示节点i与节点j的电压相角差。

无功雅各比矩阵hq为(nodenum*2+lnnum*2)行nodenum列的稀疏矩阵，无功雅各比矩阵hq中的元素计算公式为：

θft＝θf(k)-θt(k)

其中，qinj(i)为节点i的注入无功功率，bij表示节点i与节点j之间的支路电导，bii表示节点i的自电导、ui为节点i的电压幅值、uj为节点j的电压幅值；uf(k)为支路k起始节点的电压幅值；ut(k)为支路k末节点的电压幅值；θf(k)支路k起始节点的电压相角，θt(k)为支路k末节点的电压相角。f(k)为支路k起始节点编号、t(k)为支路k末节点编号,θft表示节点i与节点j的电压相角差。

在本实施例中，步骤s5所述的有功残差方程矩阵ap的计算公式为：

其中，为有功雅各比矩阵hp的转置，为有功量测权矩阵rp的逆矩阵；无功残差方程矩阵aq的计算公式为：

其中，为无功雅各比矩阵hq的转置，为无功功量测权矩阵rq的逆矩阵。在具体实施时，由于有功量测权矩阵rp与无功量测权矩阵rq均为对角阵，而且在采集电力系统量测参数时，采集的无效或实际未经传送的量测参数设置为0，因此，有功量测权矩阵rp与无功量测权矩阵rq中有将部分对角元设为0，可能存在无逆矩阵的情况，对有功量测权矩阵rp与无功量测权矩阵rq中非零元素取倒数来等效求逆。

在本实施例中，步骤s6所述的有功雅各比矩阵hp、无功雅各比矩阵hq、有功下三角矩阵apl、无功下三角矩阵aql复制到gpu的主内存时采用列优先稀疏矩阵压缩存储格式；支路导纳列表复制到gpu的主内存中采用列优先的存储格式，当有功雅各比矩阵hp、无功雅各比矩阵hq、有功下三角矩阵apl、无功下三角矩阵aql及支路导纳列表中任意一参数列的行数不等于gpu线程数的整数倍时，插入全为0的参数空行补全，从而保证触发gpu对复制到主内存的参数的合并访问。

在本实施例中，步骤s7所述的有功迭代的过程为：

s701.计算有功功率量测估计值；

首先计算支路有功量测估计值，支路k有功功率量测估计值的计算公式为：

θft＝θf(k)-θt(k)

其中，pft_est(k)为该支路由f侧到t侧支路有功功率量测估计值；ptf_est(k)为该支路由t侧到f侧支路有功功率量测估计值；uf(k)表示该支路f侧的节点电压幅值；ut(k)表示该支路t侧的节点电压幅值；θf(k)支路k起始节点的电压相角，θt(k)为支路k末节点的电压相角,θft表示节点i与节点j的电压相角差。

s702.遍历所有节点，对于任意点i，设置有功注入量测估计值变量pinj_est(i)，并初始化为0。

s703.遍历所有支路，对于任意支路k，将支路两侧有功量测估计值叠加到量测节点的注入有功功率量测估计值。例如，对于第k个支路

pinj_est(f(k))＝pinj_est(f(k))+pft_est(k)

pinj_est(t(k))＝pinj_est(t(k))+ptf_est(k)

其中，pft_est(f)为该支路f侧节点的节点注入有功量测估计值，ptf_est(t)为该支路t侧节点的节点注入有功量测估计值。

s704.遍历所有节点，将节点并联导纳有功量测估计值叠加到节点注入有功量测估计值，对于任意节点i：

其中，gsh(i)为节点i并联电导

s705.将pinj_est、pft_est、ptf_est依次排列形成有功量测估计值向量zp_est；利用有功量测估计值向量zp_est与有功测值向量zp_m求取有功量测估计误差向量：

zp_err＝zp_m-zp_est

其中,zp_err表示有功量测估计误差向量；

s706.根据有功残差方程矩阵ap及有功量测估计误差向量zp_err，求取状态变量θ的变化量δθ，状态变量θ的变化量δθ满足：

将pi_est、pij_est、pji_est依次排列形成有功量测估计值向量zp_est；利用有功量测估计值向量zp_est与有功测值向量zp_m求取有功量测估计误差向量：

zp_err＝zp_m-zp_est

其中,zp_err表示有功量测估计误差向量；

s706.根据有功残差方程矩阵ap及有功量测估计误差向量zp_err，求取状态变量θ的变化量δθ，状态变量θ的变化量δθ满足：

优选地，步骤s7所述的无功迭代的过程为：

s711.计算无功功率量测估计值；

首先计算支路无功量测估计值，支路无功功率量测估计值的计算公式为：

θft＝θf(k)-θt(k)

其中，qft_est(f)为该支路f侧节点的节点注入无功量测估计值，qtf_est(t)为该支路t侧节点的节点注入无功量测估计值。uf(k)表示该支路f侧的节点电压幅值；ut(k)表示该支路t侧的节点电压幅值；θf(k)支路k起始节点的电压相角，θt(k)为支路k末节点的电压相角,θft表示节点i与节点j的电压相角差。

s712.遍历所有节点，对于任意点i，设置有功注入量测估计值变量qinj_est(i)，并初始化为0。

s713.遍历所有支路，对于任意支路k，将支路两侧无功量测估计值叠加到量测节点的注入无功功率量测估计值。例如，对于第k个支路

qinj_est(f(k))＝qinj_est(f(k))+qft_est(k)

qinj_est(t(k))＝qinj_est(t(k))+qtf_est(k)

其中，qft_est(f)为该支路f侧节点的节点注入无功量测估计值，qtf_est(t)为该支路t侧节点的节点注入无功量测估计值。

s714.遍历所有节点，将节点并联导纳无功量测估计值叠加到节点注入无功量测估计值，对于任意节点i：

其中，bsh(i)为节点i并联电纳。

s715.将qi_est、qij_est、qji_est依次排列形成无功量测估计值向量zq_est；利用无功量测估计值向量zq_est与无功测值向量zq_m求取无功量测估计误差向量：

zq_err＝zq_m-zq_est

其中,zq_err表示无功量测估计误差向量；

s713.根据无功残差方程矩阵aq及无功量测估计误差向量zq_err，求取状态变量u的变化量δu，状态变量u的变化量δu满足：

无论是有功迭代还是无功迭代的过程中，在计算有功或无功功率量测估计值时，由于gpu各线程以单指令多数据流的形式实现并行计算，每个节点连接的支路数量不一样，用传统按节点计算节点注入有功或无功功率量测估计值的方法会导致各个线程之间的等待，严重影响了并行度，所以首先计算支路有功、无功量测估计值，以及所有节点的并联导纳有功、无功量测估计值，对支路功率求和来得到节点注入功率，提高计算的并行度；另外，在状态变量u的变化量δu或状态变量θ的变化量δθ的求解过程中，状态估计是基于pq分解法的，这时有功残差方程矩阵ap及无功残差方程矩阵aq均为常矩阵，在迭代过程中不发生变化，因此利用llt分解分别得到的下三角阵apl、aql前推回代，进一步得到δθ、δu，以简化求解过程。

有功迭代收敛条件为：

max(δθ)<θth

其中,δθ表示状态变量θ的变化量；θth表示有功迭代收敛阈值；当有功迭代收敛条件不满足时，更新节点电压相角，节点电压相角的更新公式为：

θ^(kk+1)＝θ^(kk)+δθ

其中，θ^(kk)表示第kk次迭代时的节点电压相角状态变量；θ^(kk+1)表示第kk+1次迭代时的节点电压相角状态变量；无功迭代收敛条件为：

max(δu)<uth

其中,δu表示状态变量节点电压幅值u的变化量；uth表示无功迭代收敛阈值；当无功迭代收敛条件不满足时，更新节点电压幅值，节点电压幅值的更新公式为：

u(^kk+1)＝u(^kk)+δu

其中，u^(kk)表示第kk次迭代时的节点电压幅值状态变量；u^(kk+1)表示第kk+1次迭代时的节点电压幅值状态变量。

附图中描述的关系仅用于示例性说明，不能理解为对本专利的限制；

显然，本发明的上述实施例仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。