说明书一种电力系统自适应核密度抗差状态估计方法
技术领域
本发明涉及电力系统状态估计领域，尤其是涉及一种电力系统自适应核密度抗差状态估计方法。
背景技术
电力系统状态估计作为EMS中最核心的功能之一，其性能的优劣直接影响着所有基于它而运行的高级应用软件运行及计算的可靠性和可信性。当前国家电网正全面实施状态估计精度的行业对标，以使得各省级电网的状态估计软件运行的可靠性及精度提高到一个新的高度。相应地，在此背景下，提高不良数据辨识的鲁棒性、快速性和精确性，不仅可以使调度人员能够准确、全面地掌握电力系统的实际运行状态，大大提高电力系统安全与经济运行水平；同时也能增强电力系统运行数据库的可靠性，为安全分析和运行计划等提供准确的参考数据，对提升电力公司的生产、管理效率以及使得电力公司在状态估计的行业对标中名列前茅，均具有重要的现实意义。
对于不同的状态估计模型，其求解方法有主要有牛顿迭代法、拉格朗日乘子法、内点法和智能算法。其中牛顿迭代法主要应用在基于加权最小二乘法的状态估计模型中，通过构造量测雅可比矩阵和信息矩阵相关的线性方程组实现状态变量的迭代求解。后三种方法主要应用于基于非二次准则而建立的状态估计模型中，其中拉格朗日乘子法通常应用于求解含连续型目标函数的状态估计模型，通常是为了避免因权重值相差悬殊导致方程数值不稳定而引入拉格朗日乘子将零注入功率等虚拟量测作为等式约束进行处理，然后采用牛顿迭代法或者修正牛顿迭代法进行求解；内点法和智能算法主要用于含非连续型目标函数的状态估计模型中。这些方法为实现不同状态估计模型的快速求解提供了很好的基础。
可以看出，为了保证解的有效性，不同的状态估计模型需要有相应的求解方法与之相适应，而目前只有基于一维核密度的抗差状态估计模型及其求解方法，并且 其未给出系统的核带宽确定策略。
发明内容
本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种收敛性好、抗差能力好的电力系统自适应核密度抗差状态估计方法，可用于求解基于自适应核密度抗差状态估计模型、核带宽确定策略及模型，保证自适应核密度抗差状态估计方法在实际应用中能够有效得到推广。
本发明的目的可以通过以下技术方案来实现：
一种电力系统自适应核密度抗差状态估计方法，包括以下步骤：
1)读取电力系统网络参数、量测数据及其标准差；
2)根据所述网络参数将电力系统划分为M个可观测岛；
3)依次对各观测岛进行状态估计：
a、建立自适应核密度抗差状态估计数学模型：
maxxJ(x)=Σi=1nωiσ~~iexp(-[hi(x)-zi]22σ~~i2)s.t.c(x)=0]]>
其中，i表示量测序号，i＝1,2,…,Im，Im为量测总数，J(x)为目标函数，ωi为第i个量测的权重，zi为第i个量测的量测值，hi(x)为第i个量测的估计值，ri＝hi(x)-zi表示第i个量测的量测残差，x为电力系统状态变量，表示节点电压幅值、相角向量，c(x)＝0为零注入节点功率约束方程，n表示量测数量，为第i个量测所对应的自适应核函数带宽；
b、采用牛顿迭代法对自适应核密度抗差状态估计数学模型进行求解，获得各观测岛的状态估计结果。
所述权重ωi的计算公式为：
ωi=α+exp(-si2s2),s=1nΣi=1nsi2]]>
其中，α为正常数，si为第i个量测标准差，s为全部量测标准差的几何平均值。
所述第i个量测所对应的自适应核函数带宽通过以下公式获得：
σ~~i=max(σ~~θj,σi,nor)zi=Pjmeamax(σ~~Vj,σi,nor)zi=Qjmea∪zi=Vjmeamax(min(σ~~θj,σ~~θk),σi,nor)zi=Pjkmea∪zi=Pkjmeamax(min(σ~~Vj,σ~~Vk),σi,nor)zi=Qjkmea∪zi=Qkjmea]]>
其中，为节点j注入功率量测或者支路j-k有功量测的状态变量分量综合带宽，为节点j注入功率量测或电压量测或支路j-k无功量测的状态变量分量综合带宽，分别为支路k-j的有功量测和无功量测的状态变量分量综合带宽，σi,nor为与量测zi的估计误差的方差相关的带宽，σi,nor＝βsi，si为第i个量测标准差，β为一常数，3≤β≤5。
所述综合带宽通过以下公式获得：
σθj~=max(σθj,ini~,σθj,opt~)k≤2max(σθj,cre~,σθj,opt~)k>2]]>
σVj~=max(σVj,ini~,σVj,opt~)k≤2max(σVj,cre~,σVj,opt~)k>2]]>
其中，k为牛顿迭代法的迭代次数，分别为相角、电压初始带宽，分别为相角、电压修正带宽，分别为相角、电压近似最优带宽。
所述相角、电压初始带宽具体为：
σθj,ini~=maxi=1,2,...,nj|riθjk|-2lnδgdθj,σVj,ini~=maxi=1,2,...,nj|riVjk|-2lnδgdVj]]>
其中，为与节点j的相角强相关的量测组中的第i个量测在第k次迭代中的残差，为与节点j的电压强相关的量测组中的第i个量测在第k次迭代中的残差，为正常量测的核函数阈值。
所述相角、电压修正带宽具体为：
σθj,cre~=rmidθjk-2lnδgdθj,σVj,cre~=rmidVjk-2lnδgdVj]]>
其中，为与节点j的相角强相关的nj个量测的量测残差中位数，为与节点j的电压强相关的nj个量测的量测残差中位数，为正常量测的核函数阈值。
所述相角近似最优带宽通过以下公式计算：
4Aσθj5~+2Bσθj3~-C=0]]>
其中，A=14(nj)2(Σi=1nj1ci2)2[∫u2ψ(u)du]2∫f′′2(θj)dθj]]>
B=1njΣi=1nj1ci3(Σi=1nj1ci-1)∫u2ψ(u)du∫f(θj)f′′(θj)dθj]]>
C=f(θj)∫ψ2(θj)dθj(nj)2Σi=1nj1ci]]>
f(θj)≈12πsθjexp(-(θj-θ‾j)22sθj2)]]>
θ‾j≈1njΣi=1njzi-bici,sθj=1njΣi=1njsi2ci2]]>
式中，nj为与节点j的相角强相关的量测个数，ψ(·)为高斯核函数，u表示任意变量，θj为状态变量的第j个分量，为相角估计值的均值，为θj的方差，n为节点j的邻接节点总数，参数ci,bi的具体表达式为：
ci≈1zi=VjmeaBjjzi=Pjmea∪zi=Qjmeabjkzi=Pjkmea∪zi=Qjkmea-bjkzi=Pkjmea∪zi=Qkjmea]]>
bi≈0zi=Vjmea∪zi=Pjmea∪zi=Pjkmea∪zi=Pkjmea-Σk=1,k≠jnBjkzi=Qjmea-bjkzi=Qjkmeabjkzi=Qkjmea]]>
式中，Bjj、Bjk为节点j的自导纳及节点j与节点k之间的互导纳，bjk为j-k支路的线路电纳值，为节点j的电压幅值量测；
同理获得电压近似最优带宽
所述步骤b中，采用牛顿迭代法对自适应核密度抗差状态估计数学模型进行求解，获得观测岛m的状态估计结果，m＝1,…,M，具体为：
b1)初始化，置牛顿迭代次数k＝0，电压和相角采用平值启动；
b2)对观测岛m的量测zi，根据其相关节点的电压和相角计算该量测的计算值hi(xk)，xk为状态变量的第k次迭代值，xk＝(Vm1,Vm2,…,VmNm,θm1,θm2…,θmNm)，Nm为观测岛m内的节点总数；
b3)建立自适应核密度抗差状态估计数学模型的拉格朗日方程式：
maxx,λL(x,λ)=Σi=1nωiσ~~ifi(x)+Σj=1pλjcj(x)]]>
其中，fi(x)=exp(-ri22σ~~i2),]]>λi为拉格朗日乘子；
b4)获取拉格朗日方程式的极值条件：
∂L(x,λ)∂x=Σi=1nωiσ~~i∂fi(x)∂x+Σj=1p∂cj(x)∂xλj=HTF(x)[z-h(x)]+CTλ=0∂L(x,λ)∂λ=c(x)=0F(x)=diag{Fii(x)}=diag{ωi,fi(x)/σ~~i3}]]>
其中C为p×d维零注入雅可比矩阵，H为量测雅可比矩阵；
b5)获得拉格朗日函数L(x,λ)的海森矩阵中的各项：
∂2L(x,λ)∂x2=D,∂2L(x,λ)∂λ2=0]]>
∂2L(x,λ)∂x∂λ=CT,∂2L(x,λ)∂λ∂x=C]]>
其中，D=HTF(x)(I-diag{ri2σ~~i2})H,]]>
则极值条件相应的线性方程组为：
DCTC0Δx(k)Δλ(k)=HTF(x)(z-h(x))+CTλc(x);]]>
b6)求解线性方程组获得状态变量与拉格朗日乘子的增量Δx(k)、Δλ(k)，并通过以下公式修正状态变量x(k)、拉格朗日乘子λ(k)：
x(k+1)=xk+Δx(k)λ(k+1)=λ(k)+Δλ(k)]]>
b7)判断状态变量与拉格朗日乘子的增量Δx(k)、Δλ(k)是否满足以下条件：
max{|Δx(k)|}＜ε且max{|Δλ(k)|}＜ε
ε为收敛判据，若是，则以修正后的状态变量作为当前观测岛的状态估计结果输出，若否，则k＝k+1，返回步骤b2)，进行下一次迭代。
与现有技术相比，本发明具有以下优点：
(1)本发明采用了自适应核密度带宽的方法进行电力系统状态估计，对于正常量测具有与WLS算法类似的收敛性好、渐进无偏的特点；对于不良数据，在迭代过程中通过自适应地降低其核密度带宽而等效降低其权重直至为零，相应减小或消除其影响而具有良好的抗差能力；对于可疑量测，在迭代中自适应核带宽逐渐减小，部分量测的等效权值逐步增大而被辨识为好数据，部分量测的等效权值减小而逐步被辨识为不良数据。
(2)本发明实现了可疑数据的平滑逐步辨识，避免了“非好即坏”的判定而提高了对于不良数据的辨识能力，具有很好的工程应用前景。
附图说明
图1为本发明的流程示意图；
图2为实施例中9节点系统及量测配置图；
图中，■表示有功、无功量测，●表示电压量测，1-9表示9个节点，G1、G2、 G3表示发电机；
图3为9节点系统在迭代过程部分量测的核带宽变化过程；
图4为9节点系统在迭代过程部分量测的核函数值变化过程；
图5为9节点系统计算之后各节点电压幅值误差曲线对比图；
图6为9节点系统计算之后各节点电压相角误差曲线对比图。
具体实施方式
下面结合附图和具体实施例对本发明进行详细说明。本实施例以本发明技术方案为前提进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。
如图1所示，本实施例提供一种电力系统自适应核密度抗差状态估计方法，包括以下步骤：
1)读取电力系统网络参数、量测数据及其标准差。
2)根据所述网络参数将电力系统划分为M个可观测岛，对于第m个可观测岛，其节点总数为Nm，量测总数为Im；置m＝1。
3)初始化：判断m是否大于M，若是，则完成所有可观测岛的状态估计；若否，则对第m个可观测岛，置牛顿迭代次数k＝0，电压和相角采用平值启动，即Vmj＝1，θmj＝0(j＝1，2，…，Nm)。
4)计算近似最优核带宽：对于节点j(j＝1,2,…,Nm)，对该节点的相角强相关量测，根据下式(1)的近似最优核带宽方程求解出该节点的相角强相关量测近似最优带宽
4Aσθj5~+2Bσθj3~-C=0---(1)]]>
式中，A=14(nj)2(Σi=1nj1ci2)2[∫u2ψ(u)du]2∫f′′2(θj)dθj]]>
B=1njΣi=1nj1ci3(Σi=1nj1ci-1)∫u2ψ(u)du∫f(θj)f′′(θj)dθj]]>
C=f(θj)∫ψ2(θj)dθj(nj)2Σi=1nj1ci]]>
f(θj)≈12πsθjexp(-(θj-θ‾j)22sθj2)]]>
θ‾j≈1njΣi=1njzi-bici,sθj=1njΣi=1njsi2ci2]]>
式中，nj为与节点j的相角强相关的量测个数，ψ(·)为高斯核函数，u表示任意变量，θj为状态变量的第j个分量，为相角估计值的均值，为θj的方差，n为节点j的邻接节点总数，参数ci,bi的具体表达式为：
ci≈1zi=VjmeaBjjzi=Pjmea∪zi=Qjmeabjkzi=Pjkmea∪zi=Qjkmea-bjkzi=Pkjmea∪zi=Qkjmea]]>
bi≈0zi=Vjmea∪zi=Pjmea∪zi=Pjkmea∪zi=Pkjmea-Σk=1,k≠jnBjkzi=Qjmea-bjkzi=Qjkmeabjkzi=Qkjmea]]>
式中，Bjj、Bjk为节点j的自导纳及节点j与节点k之间的互导纳，bjk为j-k支路的线路电纳值，为节点j的电压幅值量测。
节点j的电压近似最优核带宽的计算过程与上述类似。
5)计算量测值与其计算值的偏差：对该岛内的量测zi，根据其相关节点的电压和相角计算该量测的计算值hi(xk)，得到当前牛顿法迭代次数k下的残差rik＝zi-hi(xk)，其中状态变量的第k次迭代值xk＝(Vm1,Vm2,…,VmNm,θm1,θm2…,θmNm),i＝1,2,…,Im。
6)计算初始带宽和修正带宽：对于节点j(j＝1,2,…,Nm)，若k<3，则根据该节点的相角强相关量测残差和电压强相关量测残差，由式(2)的初始带宽计算公式得到相应的初始带宽和若k≥3，根据式(3)得到相应的修正带宽 和
σθj,ini~=maxi=1,2,...,nj|riθjk|-2lnδgdθj,σVj,ini~=maxi=1,2,...,nj|riVjk|-2lnδgdVj---(2)]]>
其中，为与节点j的相角强相关的量测组中的第i个量测在第k次迭代中的残差，为与节点j的电压强相关的量测组中的第i个量测在第k次迭代中的残差，为正常量测的核函数阈值。
σθj,cre~=rmidθjk-2lnδgdθj,σVj,cre~=rmidVjk-2lnδgdVj---(3)]]>
其中，为与节点j的相角强相关的nj个量测的量测残差中位数，为与节点j的电压强相关的nj个量测的量测残差中位数，为正常量测的核函数阈值。
7)计算综合带宽：根据迭代次数k和式(4)确定节点j(j＝1,2,…,Nm)的状态变量分量综合带宽和
σ~~θj=max(σθj,ini~,σθj,opt~)k≤2max(σθj,cre~σθj,opt~)k>2σ~~Vj=max(σVj,ini~,σV,opt~)k≤2max(σVj,cre~,σVj,opt~)k>2---(4)]]>
式中，(或)称为初始带宽；(或)为修正带宽；(或)为近似最优带宽。
8)计算量测核带宽及其权重：对该岛内的第i个量测zi(i＝1,2,…,Im)，根据该量测的类型及其相关节点的状态变量分量综合带宽，由式(5)确定该量测的带宽由式(7)计算其权重ωi。
σ~~i=max(σ~~θj,σi,nor)zi=Pjmeamax(σ~~Vj,σi,nor)zi=Qjmea∪zi=Vjmeamax(min(σ~~θj,σ~~θk),σi,nor)zi=Pjkmea∪zi=Pkjmeamax(min(σ~~Vj,σ~~Vk),σi,nor)zi=Qjkmea∪zi=Qkjmea---(5)]]>
式中，为节点j注入功率量测或者支路j-k有功量测的状态变量分量综合带宽；为节点j注入功率量测或电压量测或支路j-k无功量测的状态变量分量综合带宽；与分别为支路j-k的有功量测和无功量测的状态变量分量综合带宽；σi,nor为与量测zi的估计误差的方差相关的带宽，由式(6)确定，其中3≤β≤5。
σi,nor＝βsi (6)
量测权重ωi公式：
ωi=α+exp(-si2s2),s=1nΣi=1nsi2---(7)]]>
其中α为一很小的正数，si为第i个量测标准差；s为全部量测标准差的几何平均值。
9)计算量测雅可比矩阵：对于该岛内的每个量测zi，根据相应量测计算方程对电压和相角的偏导数，计算得到量测雅可比矩阵H。
10)根据自适应核密度抗差状态估计数学模型式(8)，采用牛顿法进行求解，建立拉格朗日方程式(9)：
maxxJ(x)=Σi=1nωiσ~~iexp(-[hi(x)-zi]22σ~~i2)s.t.c(x)=0---(8)]]>
其中，i表示量测序号，ωi为第i个量测的权重；zi为第i个量测的量测值；hi(x)为第i个量测的估计值；ri＝hi(x)-zi表示第i个量测的量测残差；x为系统状态变量，表示节点电压幅值、相角向量；c(x)＝0为零注入节点功率约束方程；n表示量测数量；为第i个量测所对应的自适应核函数带宽。
maxx,λL(x,λ)=Σi=1nωiσ~~ifi(x)+Σj=1pλjcj(x)---(9)]]>
其中，λi为拉格朗日乘子，其它变量及符号的含义与式(8)相同。
式(9)的极值条件为：
∂L(x,λ)∂x=Σi=1nωiσ~~i∂fi(x)∂x+Σj=1p∂cj(x)∂xλj=HTF(x)[z-h(x)]+CTλ=0∂L(x,λ)∂λ=c(x)=0F(x)=diag{Fii(x)}=diag{ωi,fi(x)/σ~~i3}---(10)]]>
其中C为p×d维零注入雅可比矩阵。
采用牛顿法求解式(10)，则拉格朗日函数L(x,λ)的海森矩阵中的各项为：
∂2L(x,λ)∂x2=D,∂2L(x,L)∂λ2=0∂2L(x,λ)∂x∂λ=CT,∂2L(x,λ)∂λ∂x=C---(11)]]>
其中D=HTF(x)(I-diag{ri2σ~~i2})H.]]>
求解式(10)的相应线性方程组为：
DCTC0Δx(k)Δλ(k)=HTF(x)(z-h(x))+CTλc(x)---(12)]]>
11)求解迭代方程：计算式(12)中的方程系数并求解该线性方程组获得Δx(k)、Δλ(k)，按式(13)修正状态变量x(k)、拉格朗日乘子λ(k)。
x(k+1)=xk+Δx(k)λ(k+1)=λ(k)+Δλ(k)---(13)]]>
12)校验状态变量与拉格朗日乘子的增量Δx(k)、Δλ(k)，若max{|Δx(k)|}＜ε且max{|Δλ(k)|}＜ε，则岛m计算结束，输出岛m的状态估计结果，m＝m+1，转至步骤3)进行下一个可观测岛的状态估计；否则k＝k+1,则转至步骤5)。其中ε为收敛判据，本实施例中，ε＝0.0001。
为了验证上述电力系统自适应核密度抗差状态估计方法(以KD为缩写)的有效性，本实施例以9节点系统为例，对本发明方法(KD)、指数型目标函数法(MES)、 平方-最小绝对值法(QLAV)、多段法(IGG)、最小绝对值(WLAV)、平方-常数法(QC)、合格率最大法(MNMR)等各种状态估计方法进行对比。
9节点系统的网络结构、量测配置如图2所示，量测信息见表1，网络参数见表2。节点5的注入有功、支路5-4的传输有功以及支路6-9、9-3的传输无功设置为不良数据，则两个有功不良数据为强相关不良数据，两个无功不良数据为强相关不良数据。
表1 9节点系统量测
表2 9节点系统网络参数
1)不良数据核带宽、核函数及其目标函数值在迭代过程的变化情况
图3和图4为利用本发明方法在9节点系统中节点5相关的部分量测及不良数据在迭代过程其核带宽和核函数值的变化过程，表3列出了相应的目标函数值变化情况。从图3和图4可以看出，随着迭代的进行，不良数据的核函数值会变小，其中由于Q6-9的核带宽逐渐增大而残差也逐渐增大，在迭代后期其核函数值迅速下降，其中Q6-9由于与真值偏离过大，因此其核函数值降为零以减小其对状态估计的影响；虽然Q5和Q5-4在迭代过程中核带宽略大，但由于其残差较小，使得其核函数值较大而不被误辨识为不良数据；由表3可知，虽然不良数据P5、Q5-4和Q9-3的核函数值和目标值略大，但其对整体目标函数的贡献相对值很小，所有不良数据的目标函数之和仅占全部量测目标函数总和的0.07％，因此其对状态估计结果影响不大，也体现了算法的抗差能力。
表3量测的目标函数值变化情况
2)与其它算法的估计结果比较
计算结果示于表4中，表4中KD表示本发明方法，Vmax表示最大电压幅值误差，Vavg表示平均电压幅值误差，θmax表示相角最大误差，θavg表示相角平均误差，ITE表示迭代次数，WBD表示残差污染，FBD表示残差淹没。
表4不同抗差算法9节点系统计算结果对比
由表4可以看出，全部算法均出现了残差污染(本发明方法KD只有支路3-9的无功误辨识为不良数据，为所有算法中最少)，而KD的残差淹没个数为零，即成功辨识了全部不良数据，也为所有算法中最少；故KD具有最强的抗差能力，即获得最高的计算精度。从迭代次数看，MES、IGG由于需要两阶段计算而导致迭代次数多于其它方法；KD与QLAV迭代次数为最少；WLAV由于目标函数在原点附近不可微导致迭代次数相对于KD及QLAV有所增加；QC由于拒绝域过宽可观测性被破坏导致估计失败；MNMR由于残差位于其近似方型的评价函数的底部时，评价函数值近似相同，即对残差较小的量测不具备滤波能力而导致了严重的污染；且该方法的目标函数在残差接近零或残差过大时变化过于缓慢，不得不采用收敛性好但迭代次数多且计算速度慢的BFGS算法或内点法，因此，其迭代次数明显多于其他方法。
图5、图6给出了KD、MES及QLAV算法对该系统计算之后节点电压幅值误差和节点电压相角误差曲线对比图。图5、图6中的横坐标表示9节点系统的节点 编号，纵坐标表示不同计算结果与真值间的偏差。
从图5、图6及表4可以看出，MES算法由于全系统仅有1个带宽，不能适应各量测残差不同的情况，未能成功辨识全部不良数据，有两个不良数据没有辨识出来，同时误辨识一个量测，故误差是最大的；而QLAV算法由于门槛固定，也未能成功辨识全部不良数据，有一个不良数据没有辨识出来，4个正常量测出现了误辨识，误差相对于MES要小。KD由于根据不同的量测残差自动选定不同的合理带宽而具有适应于不同量测的能力，实现了所有不良数据的成功辨识，由于该算例的不良数据均为强相关不良数据，因此，表明本发明方法具有辨识强相关不良数据的能力。故本发明方法具有最好的计算精度。
以上所述的具体实施例仅为说明本发明的实现效果，并不用以限制本发明。凡在本发明所提出的方法的基本思路和框架之内所作的任何非实质性的修改、转换和改进，均应包含在本发明的保护范围之内。