CN104184144A - 一种用于多电压等级电网模型的抗差状态估计方法 - Google Patents

Abstract

本发明涉及一种用于多电压等级电网模型的抗差状态估计方法，所述方法包括读取电网模型和SCADA量测；进行网络拓扑分析，将物理电网模型转化成母线-支路计算模型；建立有功加权雅可比矩阵R_a ^-1/2B_a和无功计算加权雅可比矩阵R_r ^-1/2B_r，计算有功系数矩阵A和无功系数矩阵B；利用正交变换获得有功系数矩阵和无功系数矩阵的上三角矩阵L和L₁，进行状态估计迭代求解；计算量测残差v_i，选择量测类型基准值S_i，计算遥测估计值误差计算量测权因子，判断是否有量测处于降权区或淘汰区，没有则状态估计计算结束；有则根据量测权因子修正量测权重矩阵R_a ^-1和R_r ^-1，重新生成加权雅可比矩阵R_a ^-1/2B_a和R_r ^-1/2B_r，计算A和B，返回步骤(4)。

Description

一种用于多电压等级电网模型的抗差状态估计方法

技术领域

本发明涉及一种抗差状态估计方法，具体讲涉及一种用于多电压等级电网模型的抗差状态估计方法。

背景技术

随着智能电网的建设和发展，电网规模越来越大，电网建模更加精细化，对电力系统在线分析计算，尤其是状态估计计算带来了新的挑战，如何从计算精度和计算速度两个方面对状态估计算法进行提升，快速准确的计算出大电网全网实时运行状态，为其他各类分析应用提供准确可靠的数据源，已经成为大电网状态估计计算所面临的新问题。因此，具有更强的不良数据辨识能力、更高计算精度和计算速度的抗差状态估计模型与算法已经成为电力系统研究中的一个热点方向。

现有的M估计通过等价权原理化为加权最小二乘状态估计的形式，引入权函数对量测权重进行灵活调整，继承了最小二乘的求解形式，在迭代计算过程中，根据残差动态调节量测权重，具有一定的抗差效果。但是，随着我国大规模可再生能源的集约化开发，特高压交、直流输电技术的快速发展，以及智能电网一体化运行的深入开展，各级调度中心对在线分析和仿真计算的精度提出了更高的要求，近年来已有部分省级调度中心试验将省、地电网模型拼接，并基于拼接后的大电网模型进行在线分析和调度员培训模拟。省、地电网模型拼接后电网模型电压等级跨度非常大，不同电压等级量测量的数值及其满刻度值往往也相差数倍甚至数十倍。抗差最小二乘估计中常用的权函数有Huber权函数、Hampel权函数、Turkey权函数、丹麦权函数、IGG权函数等。现有的基于IGG权函数法的抗差状态估计采用量测残差绝对值作为判断是否进行降权处理的依据，对不同电压等级的具有不同特点的量测没有进行区分，其对低电压等级设备量测或小容量设备量测的坏数据辨识能力有限。

另外，在抗差状态估计数学模型进行求解迭代过程中，需要根据权重调整情况重新建立加权雅可比矩阵及正交分解，因此，基于等价权原理的抗差最小二乘状态估计的计算效率较传统加权最小二乘方法有所下降，限制了该算法在在线分析领域中的实用化进程。

综上所述，要将抗差最小二乘状态估计算法更好的应用到实际电力系统中去，还需要对抗差状态估计算法进行深入研究。

发明内容

针对现有技术的不足，本发明提供一种用于多电压等级电网模型的抗差状态估计方法，在IGG分段权函数中引入量测类型基准值，实现对不同电压等级量测的区分，对处于降权区的量测权因子计算公式进行修正，整体提高抗差状态估计对不同电压等级量测坏数据的抗差性能。

本发明的目的是采用下述技术方案实现的：

一种用于多电压等级电网模型的抗差状态估计方法，其改进之处在于，所述方法包括

(1)读取电网模型和SCADA量测；

(2)进行网络拓扑分析，将物理电网模型转化成母线-支路计算模型；

(3)建立有功加权雅可比矩阵R_a ^-1/2B_a和无功计算加权雅可比矩阵R_r ^-1/2B_r，计算有功系数矩阵A和无功系数矩阵B；

(4)利用正交变换获得有功系数矩阵和无功系数矩阵的上三角矩阵L和L₁，进行状态估计迭代求解；

(5)计算量测残差v_i，选择量测类型基准值S_i，计算遥测估计值误差

(6)计算量测权因子，判断是否有量测处于降权区或淘汰区，没有则状态估计计算结束；有则根据量测权因子修正量测权重矩阵R_a ^-1和R_r ^-1，重新生成加权雅可比矩阵R_a ^-1/2B_a和R_r ^-1/2B_r，计算A和B，返回步骤(4)。

优选的，所述步骤(1)包括从调度控制系统模型库获取电网模型，从SCADA应用获取实时遥信遥测数据。

优选的，所述步骤(2)包括根据电网模型电气联接关系和开关/刀闸运行状态，进行厂站母线拓扑分析，形成闭合开关/刀闸相连的结点集合，分配计算母线编号；根据投运支路联接的计算母线集合，进行电气岛划分，判断电气岛带电状态。

优选的，所述步骤(3)包括

(3.1)根据节点相关支路数进行节点排序优化，采用支路电抗形成导纳矩阵，建立有功量测的雅可比矩阵B_a；形成导纳矩阵，建立无功量测的雅可比矩阵B_r；

(3.2)获取系统量测初始权重设置，建立有功加权雅可比矩阵R_a ^-1/2B_a和无功计算加权雅可比矩阵R_r ^-1/2B_r；其中，R_a ^-1和R_r ^-1分别为有功量测权重矩阵和无功量测权重矩阵；

(3.3)根据公式 $A=V_0^4\left[{(-B_a)}^T{R}_a^{-1}\left({-B}_a\right)\right]$ 和 $B=V_0^2\left[{(-B_r)}^T{R}_r^{-1}\left({-B}_r\right)\right]$ 分别计算有功系数矩阵A和计算无功系数矩阵B，其中，V₀取系统参考节点电压。

优选的，所述步骤(4)包括

(4.1)利用正交变换分别获得A和B的因子表L和L₁，置迭代计数器l＝1；

(4.2)按当前状态量电压幅值V^(l)和电压相角θ^(l)计算有功量测h_a，根据公式计算a^(l)；其中，a^(l)为有功自由矢量，z_a为有功量测矢量，h_a为用状态量V^(l)和θ^(l)表达的有功量测；

(4.3)利用A的因子表L由公式AΔθ^(l)＝a^(l)解出电压相角修正矢量Δθ^(l)，判断|Δθ^(l)|_max是否达到收敛精度；若不收敛，则计算θ^(l+1)＝θ^(l)+Δθ^(l)，若收敛，则有功收敛标志位1；

(4.4)按当前状态量电压幅值V^(l)和电压相角θ^(l)计算无功量测h_r，根据公式计算b^(l)；其中，b^(l)为无功自由矢量，z_r为无功量测矢量，h_r为用状态量V^(l)和θ^(l)表达的无功量测；

(4.5)利用B的因子表L₁由公式BΔV^(l)＝b^(l)解出电压幅值修正矢量ΔV^(l)，判断|ΔV^(l)|_max是否达到收敛精度；若不收敛，则计算V^(l+1)＝V^(l)+ΔV^(l)，若收敛，则无功收敛标志位1；

(4.6)判断有功和无功收敛标志位是否都为1，为1则状态估计迭代收敛，否则迭代计数器l＝l+1；判断是否达到最大迭代次数，达到最大迭代次数，则状态估计不收敛，计算结束；未达到最大迭代次数则返回步骤(4.2)。

进一步地，采用Givens正交三角化的列超前并行消去算法求解L和L₁，采用4进程并行模式，所述步骤(4.1)包括

(4.1.1)记加权雅可比矩阵R_a ^-1/2B_a为B_W，Givens变换将B_W的某一元素b_il零化，可以看成矩阵B_W的两行左乘一个二阶平面旋转阵，即 $\left(\begin{array}{cc}c& s\\ -s& c\end{array}\right)\left(\begin{array}{c}{\mathrm{\β}}_j\\ {\mathrm{\β}}_i\end{array}\right)=\left(\begin{array}{c}{\mathrm{\β}}_j^{\′}\\ {\mathrm{\β}}_i^{\′}\end{array}\right),$ 其中β_i，β_j为B_W矩阵的i行和j行，i,j∈{1,…,m}，i≠j，c＝sinθ_ij，s＝cosθ_ij，θ_ij表示旋转的角度；由β_i'中的b_il’＝0，可确定c和s；β_j'的元素和β_i'的其余元素为：

$b_{\mathrm{jl}}^{\′}={(b_{\mathrm{jl}}^2+b_{\mathrm{il}}^2)}^{\frac{1}{2}},c=b_{\mathrm{jl}}/b_{\mathrm{jl}}^{\′},s=b_{\mathrm{il}}/b_{\mathrm{jl}}^{\′}$

b'_jt＝cb_jt+sb_it,(t＝1,…n)

b'_it＝-sb_jt+cb_it,(t＝1,…,n,t≠l)；

(4.1.2)B_W矩阵分块存储，B_W矩阵储存在进程0中，将B_W矩阵按行分解为4块，序号为0、1、2、3，将第1、2、3块分别发送给进程1、2和3，进程0保存并处理第0块B_W矩阵数据；

(4.1.3)进程0从B_W矩阵取一行数据，经消去操作后，将此行数据发送给进程1，进程1完成相应的列消去操作后，再将此行数组发送给进程2，以此类推；进程0将数据发送后可再次从B_W矩阵取下一行数据，进行消去操作，并重复上述过程，以实现并行计算；

(4.1.4)进程3完成列消去操作后，将此行数组保存到L矩阵当中；

(4.1.5)全部消去完成后，进程3将最终生成的L矩阵发送给进程0。

优选的，所述步骤(5)包括

(5.1)根据公式ν_i＝z_i-h_i(V,θ)计算各量测残差ν_i；

(5.2)根据量测电压等级和量测设备类型选择量测类型基准值S_i，根据S_i计算量测估计误差

优选的，所述步骤(6)包括

根据改进IGG权函数 $\mathrm{\&omega;}\left(v\right)=\left\{\begin{array}{cc}1& \frac{\left|v_i\right|}{S_i}<k_1\mathrm{\&sigma;}\\ \frac{{\left(k_1\mathrm{\&sigma;}\right)}^2}{{v_i}^2}& k_1\mathrm{\&sigma;}\&le;\frac{\left|v_i\right|}{S_i}\&le;k_2\mathrm{\&sigma;}\\ 0& \frac{\left|v_i\right|}{S_i}>k_2\mathrm{\&sigma;}\end{array}\right.$ 计算所有量测权因子，其中，为保权区，为降权区，为淘汰区；判断是否有量测处于降权区或淘汰区，没有则状态估计计算结束；有则根据量测权因子修正量测权重矩阵R_a ^-1和R_r ^-1，重新生成有功加权雅可比矩阵R_a ^-1/2B_a和无功计算加权雅可比矩阵R_r ^-1/2B_r，计算有功系数矩阵A和计算无功系数矩阵B，返回步骤(4)，直至没有量测处于降权区或淘汰区，则状态估计结束。

进一步地，所述公式 $\mathrm{\&omega;}\left(v\right)=\left\{\begin{array}{cc}1& \frac{\left|v_i\right|}{S_i}<k_1\mathrm{\&sigma;}\\ \frac{{\left(k_1\mathrm{\&sigma;}\right)}^2}{{v_i}^2}& k_1\mathrm{\&sigma;}\&le;\frac{\left|v_i\right|}{S_i}\&le;k_2\mathrm{\&sigma;}\\ 0& \frac{\left|v_i\right|}{S_i}>k_2\mathrm{\&sigma;}\end{array}\right.$ 的变换过程如下：

(9.1)IGG权函数 $\mathrm{\&omega;}\left(v\right)=\left\{\begin{array}{cc}1& \left|v\right|\&le;k_1\mathrm{\&sigma;}\\ \frac{k_1\mathrm{\&sigma;}}{\left|v\right|}& k_1\mathrm{\&sigma;}\&le;\left|v\right|k_2\mathrm{\&sigma;}\\ 0& \left|v\right|>k_2\mathrm{\&sigma;}\end{array}\right.$ 采用量测残差绝对值对权函数进行分区，引入量测类型基准值S_i，计算遥测数据估计值误差v_i/S_i，对权函数进行分区，区分了不同电压等级的量测；其中，k1，k2∈R，0<k1<k2，σ为电力系统遥测量的标准误差，|ν|为量测残差绝对值；

(9.2)IGG权函数 $\mathrm{\&omega;}\left(v\right)=\left\{\begin{array}{cc}1& \left|v\right|\&le;k_1\mathrm{\&sigma;}\\ \frac{k_1\mathrm{\&sigma;}}{\left|v\right|}& k_1\mathrm{\&sigma;}\&le;\left|v\right|k_2\mathrm{\&sigma;}\\ 0& \left|v\right|>k_2\mathrm{\&sigma;}\end{array}\right.$ 通过等价权化作用于加权最小二乘法目标函数时，处于降权区的量测的目标函数分量为J_i(x)＝r_ik₁σ|v|，修正原有权函数中处于降权区的量测权因子为

(3)通过(9.1)和(9.2)得到改进的IGG权函数为 $\mathrm{\&omega;}\left(v\right)=\left\{\begin{array}{cc}1& \frac{\left|v_i\right|}{S_i}<k_1\mathrm{\&sigma;}\\ \frac{{\left(k_1\mathrm{\&sigma;}\right)}^2}{{v_i}^2}& k_1\mathrm{\&sigma;}\&le;\frac{\left|v_i\right|}{S_i}\&le;k_2\mathrm{\&sigma;}\\ 0& \frac{\left|v_i\right|}{S_i}>k_2\mathrm{\&sigma;}\end{array}\right.;$ 其中，为遥测数据估计值误差，遥测数据估计值误差门槛k₁，k₂取值根据实际系统量测情况设置，k₁大于0.03，k₂大于k₁。

与现有技术比，本发明的有益效果为：

(1)使用遥测估计值相对误差对不同电压等级量测进行区分。将量测满刻度值(因各量测满刻度值难以统计，实际系统使用量测类型基准值代替量测满刻度值)引入到抗差状态估计的权函数中，不同电压等级的量测具有相同的IGG权函数形式。但由于各量测满刻度值不同，其分区策略也不同，实现了对不同电压等级量测的区分，在现有抗差状态估计抗差性能的基础上，提高了对低电压等级量测坏数据的辨识能力和抗差能力，整体提高抗差状态估计的计算精度。

(2)充分利用各级调度中心自动化系统的多核CPU配置高性能服务器的优势，使用Givens正交三角化的列超前并行消去算法和矩阵向量并行乘法对状态估计修正方程的求解进行加速，从而显著缩短每次调节权重后加权雅可比矩阵的正交变换时间，提高抗差状态估计的计算速度。

附图说明

图1为本发明提供的一种用于多电压等级电网模型的抗差状态估计方法流程图；

图2为本发明提供的并行Givens变化程序设计。

具体实施方式

下面结合附图对本发明的具体实施方式作进一步的详细说明。

本发明一种用于多电压等级电网模型的抗差状态估计方法，具体为基于等价权原理的抗差最小二乘状态估计保留了加权最小二乘状态估计的求解形式。改进IGG权函数法抗差最小二乘状态估计先利用加权最小二乘状态估计进行求解，根据WLS的结果计算量测估计值误差，确定量测全因子，对量测权重进行修正，重新建立加权雅可比矩阵、进行正交变换和快速分解迭代求解，直至无量测权重调整为止，状态估计计算结束。

1、抗差最小二乘状态估计数学模型

电力系统状态估计的量测方程是：

z＝h(x)+v                (1)

式中：z为量测向量，h(x)为量测方程，v为量测误差向量，x为状态向量。

由于各量测量的精度不同，对不同量测取不同权重R_i，精度高的取权重大些，精度低的取权重小些，当状态量的估计值为最优时，目标函数J(x)最小。

加权最小二乘法(WLS)目标函数为：

$J\left(x\right)=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{r}_i{v_i}^2,$ 其中 $r_i=\frac{1}{{\mathrm{\&sigma;}}_i^2}---\left(2\right)$

式中：r_i为量测i的权重(即方差的倒数)，v_i残差v的第i个分量，m为量测个数，σ为标准差。

式(2)的迭代修正公式如下：

(H^TR^-1H)Δx^(k)＝H^TR^-1v^(k)                   (3)

抗差估计的形式如下：

$J\left(x\right)=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{l}_i\frac{v_i}{{\mathrm{\&sigma;}}_i}\frac{\mathrm{\&Psi;}\left(u_i\right)}{u_i}---\left(4\right)$

式中：v_i为残差v的第i个分量；P_i为先验权；Ψ函数是待选择的，u_i＝v_i/(σ_iv(l_i))；v(l_i)为已知函数，l_i ^T为R^-1/2H矩阵的行矢量；H为m×n阶雅可比矩阵，m，n分别表示量测量及状态量的个数，R^-1/2H＝[l₁,l₂,...,l_m]^T则表示加权雅可比矩阵。将式(4)整理成矩阵形式：

H^TR^-1Qv＝0                  (5)

$Q=\mathrm{diag}\left\{\frac{\mathrm{\&psi;}\left(u_i\right)}{u_i}\right\}---\left(6\right)$

式(6)的最小二乘解如下：

(H^TR^-1Q^(k)H)Δx^(k)＝H^TR^-1Q^(k)v^(k)              (7)

适当地选择Ψ函数可使对x的M估计具有较高的效率及较强的抗差能力。Ψ函数确定后，等价权Q^(k)也就确定了，式(6)便是修正量求解方程。本专利抗差最小二乘状态估计的等价权函数采用IGG权函数的形式。

2、IGG权函数

本方法在利用等加权原理调节量测量权重时选用IGG权函数形式作为抗差最小二乘状态估计的权函数形式。IGG法的权函数分为保权区、降权区和淘汰区。在状态估计过程中，充分利用主体可靠信息，使其保持原权值不变。对可疑量测值依其可信度分别降权处理，对显著异常量测值则取权值为0，相当于删除该量测。

电力系统遥测量的标准误差σ大约为量测正常范围的0.5％～2.0％，误差大于±3σ的量测就可以认为是不良数据，在正常量测条件下，误差大于±3σ的量测值出现的概率仅为0.27％，实际系统中，可根据量测冗余度及全网量测整体精度调整不良数据检测标准。

设定两个界限k₁σ和k₂σ(k₁,k₂∈R，0<k₁<k₂)，当|v|≤k₁σ时，量测处于保权区，权因子为1，其目标函数分量J_i(x)＝r_iv_i ²保留了最小二乘的形式；当k₁σ＜|v|≤k₂σ时，量测处于降权区，需要降低量测的观测权，其目标函数分量为J_i(x)＝r_ik₁σ|v|，小于最小二乘形式J_i(x)＝r_iv_i ²，达到量测降权目的；对于|v|＞k₂σ的情况，这些量测不应作为观测信息，其权因子取0，相当于删除该量测重新进行状态估计。权因子取值形式如下：

$\mathrm{\&omega;}\left(v\right)=\left\{\begin{array}{cc}1& \left|v\right|\&le;k_1\mathrm{\&sigma;}\\ \frac{k_1\mathrm{\&sigma;}}{\left|v\right|}& k_1\mathrm{\&sigma;}\&le;\left|v\right|k_2\mathrm{\&sigma;}\\ 0& \left|v\right|>k_2\mathrm{\&sigma;}\end{array}\right.---\left(8\right)$

当式(8)的权函数作用于加权最小二乘法(WLS)目标函数时，处于降权区的量测的目标函数分量为J_i(x)＝r_ik₁σ|v|，为残差绝对值的一次函数，一定程度上起到降权目的，但其目标函数分量仍然随残差绝对值的增大而增大，说明在此区间内，量测坏数据的污染依然随残差绝对值增大而增大。另外，式(8)所示的IGG权函数未对不同电压等级量测进行区分。对此，本方法使用改进的IGG权函数。

3、改进IGG权函数

(1)由于不同电压等级、不同电气设备的量测量的满刻度值不同，甚至相差几倍至数十倍，其测量精度及误差也就不同。通常，较低电压等级的量测值较小，其坏数据的绝对残差也较小。如公式(8)中所示的利用量测残差绝对值作为分区界限，不能很好的处理不同量测特点的坏数据，尤其对低电压等级设备量测坏数据辨识能力有限。为此，本专利改进IGG权函数法抗差最小二乘状态估计，在权函数中引入量测类型基准值S_i，对不同电压等级量测进行区分；

量测类型基准值通常是指该量测量的满刻度值。为使计算方便，线路功率量测、零阻抗支路量测、串联补偿器功率量测及母线电压量测可根据其电压等级选取相应的基准值；发电机功率量测及变压器功率量测则可使用设备额定容量作为量测基准值。

利用量测类型基准值S_i计算遥测数据估计值误差，根据估计值误差对IGG权函数重新进行分段，增加了不同电压等级不良数据辨识的区分度。

(2)将权函数中处于降权区的量测权因子修改为此方法具有更强的抗差性能。

改进后权函数如下所示。

其中，为遥测数据估计值误差。遥测数据估计值误差门槛k₁，k₂取值根据实际系统量测情况设置，通常情况，k₁大于0.03，k₂大于k₁。

式(9)与式(8)相比，有两点改进：第一、引入量测类型基准值S_i对权函数重新进行分段，虽然不同电压等级的量测具有相同的权函数形式，但其权函数分段策略却不同，实现对不同电压等级量测的区分。对低电压等级量测，其具有较小的量测类型基准值，其处于降权区的量测的残差绝对值较小，符合调度自动化系统对不同电压等级量测精度的要求。第二、处于降权区的量测权因子形式与式(8)不同，利用式(9)的权函数，代入加权最小二乘法(WLS)目标函数，则处于降权区的量测的目标函数分量为常数：J_i(x)＝r_ik₁ ²σ²，即，当量测处于降权区时，其残差污染不再随残差增大而扩大，有效抑制了残差污染。

如图1所示，本发明一种用于多电压等级电网模型的抗差状态估计方法流程如下：

步骤1：读取电网模型和SCADA实时数据。从调度控制系统模型库获取电网模型，从SCADA应用获取实时遥信遥测数据。

步骤2：进行网络拓扑分析，通过拓扑分析将物理电网模型转化成母线-支路计算模型。根据电网模型电气联接关系和开关/刀闸运行状态，进行厂站母线拓扑分析，形成闭合开关/刀闸相连的结点集合(即计算母线)，分配计算母线编号；根据投运支路联接的计算母线集合，进行电气岛划分，判断电气岛带电状态。

步骤3：根据量测量、量测权重及导纳矩阵分别建立有功加权雅可比矩阵R_a ^-1/2B_a和无功计算加权雅可比矩阵R_r ^-1/2B_r，计算有功系数矩阵A和无功系数矩阵B；

步骤3-1：根据节点相关支路数进行节点排序优化，用支路电抗形成导纳矩阵，建立有功量测的雅可比矩阵B_a；形成正常导纳矩阵，建立无功量测的雅可比矩阵B_r；

步骤3-2：获取系统量测初始权重设置，建立有功加权雅可比矩阵R_a ^-1/2B_a和无功计算加权雅可比矩阵R_r ^-1/2B_r，其中R_a ^-1和R_r ^-1分别为有功量测权重矩阵和无功量测权重矩阵；

步骤3-3：根据公式 $A=V_0^4\left[{(-B_a)}^T{R}_a^{-1}\left({-B}_a\right)\right]$ 和 $B=V_0^2\left[{(-B_r)}^T{R}_r^{-1}\left({-B}_r\right)\right]$ 分别计算有功系数矩阵A和计算无功系数矩阵B，其中V₀取系统参考节点电压。

步骤4：利用正交变换获得有功和无功迭代上三角矩阵L和L₁，进行状态估计迭代求解；

步骤4-1：利用Givens正交三角化的列超前并行消去算法分别获得A和B的因子表L和L₁，置迭代计数器l＝1。以4进程并行模式为例说明，所述步骤4-1包括：

步骤4-1-1：记加权雅可比矩阵R_a ^-1/2B_a为B_W，Givens变换将B_W的某一元素b_il零化，可以看成矩阵B_W的两行左乘一个二阶平面旋转阵，即 $\left(\begin{array}{cc}c& s\\ -s& c\end{array}\right)\left(\begin{array}{c}{\mathrm{\&beta;}}_j\\ {\mathrm{\&beta;}}_i\end{array}\right)=\left(\begin{array}{c}{\mathrm{\&beta;}}_j^{\&prime;}\\ {\mathrm{\&beta;}}_i^{\&prime;}\end{array}\right),$ 其中β_i，β_j为B_W矩阵的i行和j行，i,j∈{1,…,m}，i≠j，c＝sinθ_ij，s＝cosθ_ij，θ_ij表示旋转的角度。由β_i'中的b_il’＝0，可以确定c和s。β_j'的元素和β_i'的其余元素为：

$b_{\mathrm{jl}}^{\&prime;}={(b_{\mathrm{jl}}^2+b_{\mathrm{il}}^2)}^{\frac{1}{2}},c=b_{\mathrm{jl}}/b_{\mathrm{jl}}^{\&prime;},s=b_{\mathrm{il}}/b_{\mathrm{jl}}^{\&prime;}---\left(10\right)$

b'_jt＝cb_jt+sb_it,(t＝1,…n)               (11)

b'_it＝-sb_jt+cb_it,(t＝1,…,n,t≠l)        (12)

步骤4-1-2：B_W矩阵分块存储。B_W矩阵储存在进程0中，将B_W矩阵按行分解为4块，序号为0、1、2、3，将第1、2、3块分别发送给进程1、2和3，进程0保存并处理第0块B_W矩阵数据。

步骤4-1-3：进程0从B_W矩阵取一行数据，经消去操作后，将此行数据发送给进程1，进程1完成相应的列消去操作后，再将此行数组发送给进程2，以此类推。进程0将数据发送后可再次从B_W矩阵取下一行数据，进行消去操作，并重复上述过程，以实现并行计算。

步骤4-1-4：进程3完成列消去操作后，将此行数组保存到L矩阵当中。

步骤4-1-5：全部消去完成后，进程3将最终生成的L矩阵发送给进程0。

步骤4-2：按当前状态量电压幅值V^(l)和电压相角θ^(l)计算有功量测h_a，根据公式计算a^(l)。其中：a^(l)为有功自由矢量，z_a为有功量测矢量，h_a为用状态量V^(l)和θ^(l)表达的有功量测；

步骤4-3：利用A的因子表L由公式AΔθ^(l)＝a^(l)解出电压相角修正矢量Δθ^(l)，判断|Δθ^(l)|_max是否达到收敛精度；若不收敛，则计算θ^(l+1)＝θ^(l)+Δθ^(l)，若收敛，则有功收敛标志位1；

步骤4-4：按当前状态量电压幅值V^(l)和电压相角θ^(l)计算无功量测h_r，根据公式计算b^(l)。其中：b^(l)为无功自由矢量，z_r为无功量测矢量，h_r为用状态量V^(l)和θ^(l)表达的无功量测；

步骤4-5：利用B的因子表L₁由公式BΔV^(l)＝b^(l)解出电压幅值修正矢量ΔV^(l)，判断|ΔV^(l)|_max是否达到收敛精度；若不收敛，则计算V^(l+1)＝V^(l)+ΔV^(l)，若收敛，则无功收敛标志位1；

步骤4-6：判断有功和无功收敛标志位是否都为1，若都是1则状态估计迭代收敛，否则迭代计数器l＝l+1；判断是否达到最大迭代次数，若达到最大迭代次数，则状态估计不收敛，计算结束；未达到最大迭代次数则返回步骤4-2。

步骤5：计算遥测估计值误差，根据式(9)计算各量测权因子，若有量测权重调整，则修正量测权重矩阵，重新建立加权雅可比矩阵，重新计算A和B，返回步骤(4)；若无量测权重调整，则状态估计结束。

步骤5-1：根据公式ν_i＝z_i-h_i(V,θ)计算各量测残差ν_i；

步骤5-2：根据量测电压等级和量测设备类型选择量测类型基准值S_i，具体S_i取值可参考2013年8月7日国家电网公司国家电力调度控制中心印发《智能电网调度控制系统实用化验收办法(试行)》，根据S_i计算量测估计误差|ν_i|/S_i；

步骤6：根据改进IGG权函数 $\mathrm{\&omega;}\left(v\right)=\left\{\begin{array}{cc}1& \frac{\left|v_i\right|}{S_i}<k_1\mathrm{\&sigma;}\\ \frac{{\left(k_1\mathrm{\&sigma;}\right)}^2}{{v_i}^2}& k_1\mathrm{\&sigma;}\&le;\frac{\left|v_i\right|}{S_i}\&le;k_2\mathrm{\&sigma;}\\ 0& \frac{\left|v_i\right|}{S_i}>k_2\mathrm{\&sigma;}\end{array}\right.$ 计算所有量测权因子，其中，为保权区，为降权区，为淘汰区；判断是否有量测处于降权区或淘汰区，没有则状态估计计算结束；有则根据量测权因子修正量测权重矩阵R_a ^-1和R_r ^-1，重新生成有功加权雅可比矩阵R_a ^-1/2B_a和无功计算加权雅可比矩阵R_r ^-1/2B_r，计算有功系数矩阵A和计算无功系数矩阵B，返回步骤(4)，直至没有量测处于降权区或淘汰区，则状态估计结束。

最后应当说明的是：以上实施例仅用以说明本发明的技术方案而非对其限制，所属领域的普通技术人员参照上述实施例依然可以对本发明的具体实施方式进行修改或者等同替换，这些未脱离本发明精神和范围的任何修改或者等同替换，均在申请待批的本发明的权利要求保护范围之内。

Claims (9)

1.一种用于多电压等级电网模型的抗差状态估计方法，其特征在于，所述方法包括

(1)读取电网模型和SCADA量测；

(2)进行网络拓扑分析，将物理电网模型转化成母线-支路计算模型；

(3)建立有功加权雅可比矩阵R_a ^-1/2B_a和无功计算加权雅可比矩阵R_r ^-1/2B_r，计算有功系数矩阵A和无功系数矩阵B；

(4)利用正交变换获得有功系数矩阵和无功系数矩阵的上三角矩阵L和L₁，进行状态估计迭代求解；

(5)计算量测残差v_i，选择量测类型基准值S_i，计算遥测估计值误差

(6)计算量测权因子，判断是否有量测处于降权区或淘汰区，没有则状态估计计算结束；有则根据量测权因子修正量测权重矩阵R_a ^-1和R_r ^-1，重新生成加权雅可比矩阵R_a ^-1/2B_a和R_r ^-1/2B_r，计算A和B，返回步骤(4)。

2.如权利要求1所述的一种用于多电压等级电网模型的抗差状态估计方法，其特征在于，所述步骤(1)包括从调度控制系统模型库获取电网模型，从SCADA应用获取实时遥信遥测数据。

3.如权利要求1所述的一种用于多电压等级电网模型的抗差状态估计方法，其特征在于，所述步骤(2)包括根据电网模型电气联接关系和开关/刀闸运行状态，进行厂站母线拓扑分析，形成闭合开关/刀闸相连的结点集合，分配计算母线编号；根据投运支路联接的计算母线集合，进行电气岛划分，判断电气岛带电状态。

4.如权利要求1所述的一种用于多电压等级电网模型的抗差状态估计方法，其特征在于，所述步骤(3)包括

(3.1)根据节点相关支路数进行节点排序优化，采用支路电抗形成导纳矩阵，建立有功量测的雅可比矩阵B_a；形成导纳矩阵，建立无功量测的雅可比矩阵B_r；

(3.2)获取系统量测初始权重设置，建立有功加权雅可比矩阵R_a ^-1/2B_a和无功计算加权雅可比矩阵R_r ^-1/2B_r；其中，R_a ^-1和R_r ^-1分别为有功量测权重矩阵和无功量测权重矩阵；

(3.3)根据公式 $A=V_0^4\left[{(-B_a)}^T{R}_a^{-1}\left({-B}_a\right)\right]$ 和 $B=V_0^2\left[{(-B_r)}^T{R}_r^{-1}\left({-B}_r\right)\right]$ 分别计算有功系数矩阵A和计算无功系数矩阵B，其中，V₀取系统参考节点电压。

5.如权利要求1所述的一种用于多电压等级电网模型的抗差状态估计方法，其特征在于，所述步骤(4)包括

(4.1)利用正交变换分别获得A和B的因子表L和L₁，置迭代计数器l＝1；

(4.2)按当前状态量电压幅值V^(l)和电压相角θ^(l)计算有功量测h_a，根据公式计算a^(l)；其中，a^(l)为有功自由矢量，z_a为有功量测矢量，h_a为用状态量V^(l)和θ^(l)表达的有功量测；

(4.3)利用A的因子表L由公式AΔθ^(l)＝a^(l)解出电压相角修正矢量Δθ^(l)，判断|Δθ^(l)|_max是否达到收敛精度；若不收敛，则计算θ^(l+1)＝θ^(l)+Δθ^(l)，若收敛，则有功收敛标志位1；

(4.4)按当前状态量电压幅值V^(l)和电压相角θ^(l)计算无功量测h_r，根据公式计算b^(l)；其中，b^(l)为无功自由矢量，z_r为无功量测矢量，h_r为用状态量V^(l)和θ^(l)表达的无功量测；

(4.5)利用B的因子表L₁由公式BΔV^(l)＝b^(l)解出电压幅值修正矢量ΔV^(l)，判断|ΔV^(l)|_max是否达到收敛精度；若不收敛，则计算V^(l+1)＝V^(l)+ΔV^(l)，若收敛，则无功收敛标志位1；

(4.6)判断有功和无功收敛标志位是否都为1，为1则状态估计迭代收敛，否则迭代计数器l＝l+1；判断是否达到最大迭代次数，达到最大迭代次数，则状态估计不收敛，计算结束；未达到最大迭代次数则返回步骤(4.2)。

6.如权利要求5所述的一种用于多电压等级电网模型的抗差状态估计方法，其特征在于，采用Givens正交三角化的列超前并行消去算法求解L和L₁，采用4进程并行模式，所述步骤(4.1)包括

(4.1.1)记加权雅可比矩阵R_a ^-1/2B_a为B_W，Givens变换将B_W的某一元素b_il零化，可以看成矩阵B_W的两行左乘一个二阶平面旋转阵，即 $\left(\begin{array}{cc}c& s\\ -s& c\end{array}\right)\left(\begin{array}{c}{\mathrm{\&beta;}}_j\\ {\mathrm{\&beta;}}_i\end{array}\right)=\left(\begin{array}{c}{\mathrm{\&beta;}}_j^{\&prime;}\\ {\mathrm{\&beta;}}_i^{\&prime;}\end{array}\right),$ 其中β_i，β_j为B_W矩阵的i行和j行，i,j∈{1,…,m}，i≠j，c＝sinθ_ij，s＝cosθ_ij，θ_ij表示旋转的角度；由β_i'中的b_il’＝0，可确定c和s；β_j'的元素和β_i'的其余元素为：

$b_{\mathrm{jl}}^{\&prime;}={(b_{\mathrm{jl}}^2+b_{\mathrm{il}}^2)}^{\frac{1}{2}},c=b_{\mathrm{jl}}/b_{\mathrm{jl}}^{\&prime;},s=b_{\mathrm{il}}/b_{\mathrm{jl}}^{\&prime;}$

b'_jt＝cb_jt+sb_it,(t＝1,…n)

b'_it＝-sb_jt+cb_it,(t＝1,…,n,t≠l)；

(4.1.2)B_W矩阵分块存储，B_W矩阵储存在进程0中，将B_W矩阵按行分解为4块，序号为0、1、2、3，将第1、2、3块分别发送给进程1、2和3，进程0保存并处理第0块B_W矩阵数据；

(4.1.3)进程0从B_W矩阵取一行数据，经消去操作后，将此行数据发送给进程1，进程1完成相应的列消去操作后，再将此行数组发送给进程2，以此类推；进程0将数据发送后可再次从B_W矩阵取下一行数据，进行消去操作，并重复上述过程，以实现并行计算；

(4.1.4)进程3完成列消去操作后，将此行数组保存到L矩阵当中；

(4.1.5)全部消去完成后，进程3将最终生成的L矩阵发送给进程0。

7.如权利要求1所述的一种用于多电压等级电网模型的抗差状态估计方法，其特征在于，所述步骤(5)包括

(5.1)根据公式ν_i＝z_i-h_i(V,θ)计算各量测残差ν_i；

(5.2)根据量测电压等级和量测设备类型选择量测类型基准值S_i，根据S_i计算量测估计误差

8.如权利要求1所述的一种用于多电压等级电网模型的抗差状态估计方法，其特征在于，所述步骤(6)包括

根据改进IGG权函数 $\mathrm{\&omega;}\left(v\right)=\left\{\begin{array}{cc}1& \frac{\left|v_i\right|}{S_i}<k_1\mathrm{\&sigma;}\\ \frac{{\left(k_1\mathrm{\&sigma;}\right)}^2}{{v_i}^2}& k_1\mathrm{\&sigma;}\&le;\frac{\left|v_i\right|}{S_i}\&le;k_2\mathrm{\&sigma;}\\ 0& \frac{\left|v_i\right|}{S_i}>k_2\mathrm{\&sigma;}\end{array}\right.$ 计算所有量测权因子，其中，为保权区，为降权区，为淘汰区；判断是否有量测处于降权区或淘汰区，没有则状态估计计算结束；有则根据量测权因子修正量测权重矩阵R_a ^-1和R_r ^-1，重新生成有功加权雅可比矩阵R_a ^-1/2B_a和无功计算加权雅可比矩阵R_r ^-1/2B_r，计算有功系数矩阵A和计算无功系数矩阵B，返回步骤(4)，直至没有量测处于降权区或淘汰区，则状态估计结束。

9.如权利要求8所述的一种用于多电压等级电网模型的抗差状态估计方法，其特征在于，所述公式 $\mathrm{\&omega;}\left(v\right)=\left\{\begin{array}{cc}1& \frac{\left|v_i\right|}{S_i}<k_1\mathrm{\&sigma;}\\ \frac{{\left(k_1\mathrm{\&sigma;}\right)}^2}{{v_i}^2}& k_1\mathrm{\&sigma;}\&le;\frac{\left|v_i\right|}{S_i}\&le;k_2\mathrm{\&sigma;}\\ 0& \frac{\left|v_i\right|}{S_i}>k_2\mathrm{\&sigma;}\end{array}\right.$ 的变换过程如下：

(9.1)IGG权函数 $\mathrm{\&omega;}\left(v\right)=\left\{\begin{array}{cc}1& \left|v\right|\&le;k_1\mathrm{\&sigma;}\\ \frac{k_1\mathrm{\&sigma;}}{\left|v\right|}& k_1\mathrm{\&sigma;}\&le;\left|v\right|k_2\mathrm{\&sigma;}\\ 0& \left|v\right|>k_2\mathrm{\&sigma;}\end{array}\right.$ 采用量测残差绝对值对权函数进行分区，引入量测类型基准值S_i，计算遥测数据估计值误差v_i/S_i，对权函数进行分区，区分了不同电压等级的量测；其中，k1，k2∈R，0<k1<k2，σ为电力系统遥测量的标准误差，|ν|为量测残差绝对值；

(9.2)IGG权函数 $\mathrm{\&omega;}\left(v\right)=\left\{\begin{array}{cc}1& \left|v\right|\&le;k_1\mathrm{\&sigma;}\\ \frac{k_1\mathrm{\&sigma;}}{\left|v\right|}& k_1\mathrm{\&sigma;}\&le;\left|v\right|k_2\mathrm{\&sigma;}\\ 0& \left|v\right|>k_2\mathrm{\&sigma;}\end{array}\right.$ 通过等价权化作用于加权最小二乘法目标函数时，处于降权区的量测的目标函数分量为J_i(x)＝r_ik₁σ|v|，修正原有权函数中处于降权区的量测权因子为

(3)通过(9.1)和(9.2)得到改进的IGG权函数为 $\mathrm{\&omega;}\left(v\right)=\left\{\begin{array}{cc}1& \frac{\left|v_i\right|}{S_i}<k_1\mathrm{\&sigma;}\\ \frac{{\left(k_1\mathrm{\&sigma;}\right)}^2}{{v_i}^2}& k_1\mathrm{\&sigma;}\&le;\frac{\left|v_i\right|}{S_i}\&le;k_2\mathrm{\&sigma;}\\ 0& \frac{\left|v_i\right|}{S_i}>k_2\mathrm{\&sigma;}\end{array}\right.;$ 其中，为遥测数据估计值误差，遥测数据估计值误差门槛k₁，k₂取值根据实际系统量测情况设置，k₁大于0.03，k₂大于k₁。