估计不确定性（共6篇）
估计不确定性 篇1
0 引言
随着电力系统规模的日益庞大，系统的不确定性也逐渐增加[1]。不确定性的存在使得电力系统的规划和运行控制都面临着各种风险，在电力系统的安全稳定运行方面引起的弊端尤为显著。当前特高压交直流输电技术的迅速发展促进了全国电网的互联[2]，不确定性因素影响的范围更加广泛，产生的危害也更加严重。电力系统的不确定性主要表现在系统潮流的不确定性、发电机组及线路等电力元件运行状态的不确定性、风电等可再生能源出力的不确定性等。其中，由于电力元件运行状态以及可再生能源出力的不确定性会直接或间接地反映在系统潮流中，故对系统潮流的不确定性进行分析成为当前电力系统不确定性研究的基础与重点。
针对系统潮流的不确定性，概率性的求解方法可较好地适应这一需求，当前这类方法主要可以归为模拟法、理论分析法和近似法3类。以基于简单随机采样的蒙特卡洛模拟（Monte Carlo simulation,MCS）法[3]为代表的模拟方法计算精度最高，但在计算时间方面消耗巨大[4]，因此在系统规模较大的情况下该方法几乎不适用；但作为一种算法准确度的验证标准，蒙特卡洛模拟法的作用已得到公认。之后有学者提出改进的基于拉丁超立方采样的蒙特卡洛法[5,6]，在一定程度上降低了采样次数，但并未从根本上解决模拟法因保证计算精度而需要进行大量采样的问题。另一方面，以采用快速傅里叶变换[7]和半不变量方法[8,9]为代表的理论计算方法，凭借其较好的理论基础，在计算时可获得较高的效率，但会引入较大的误差，且半不变量法难以处理注入变量非正态分布尤其是离散度很大时的不确定性问题。
相比之下，在诸多近似法中，最早由Rosenblueth[10]提出的点估计法（point estimate method,PEM），由于其对注入量的分布类型没有限制且计算效率上较模拟法有极大提高，故成为工程应用上一种较实用的方法。在电力系统领域，点估计法也逐渐被用来研究与随机变量相关联的不确定潮流问题[11,12,13]。在这些应用中，点估计法基本仍以Hong[14]提出的两点或三点形式为主，但随着当前可再生能源的多点接入，随机注入量分布的非正态性凸显，该传统点估计法的计算精度会因此受到影响。由于点估计法的精度与取点数关系密切，提高点估计法精度的最根本办法是增加估计点的数目，文献[15]指出Rosenblueth的点估计理论在本质上与Gauss-Hermite的数值积分一致，这也为多点估计方法的应用提供了理论基础。文献[16]在文献[17]的研究基础上，提出一种多点估计法用于不确定性潮流计算，有效地提高了计算精度，但其在估计点选取时采用了Rosenblatt变换，该变换要求较高，需要知道输入随机变量完整的联合概率密度函数。而当随机变量的边际概率分布函数和相关系数矩阵已知时，可采用Nataf变换[18,19]来替代Rosenblatt变换进行概率建模。Nataf变换作为一种空间转换方法，在电力系统领域已有一些研究，文献[20]采用Nataf变换来获取不同风电场之间具有相关性的风速分布样本，进而分析风电场相关性对最优潮流的影响，文献[21]则将Nataf变换与拉丁超立方采样结合来研究不确定性潮流问题，而在其他领域，Nataf变换与点估计法的结合[22]也用来研究边坡稳定可靠度的相关问题[23]。
本文研究了Nataf变换以及Gauss-Hermite积分与点估计之间的关系，以Gauss-Hermite数值积分为基础，将Nataf变换的空间转换过程与点估计法结合，从而提出了一种改进的多点估计法。同时本文考虑了风力发电机、常规发电机出力的波动性，以及服从不同类型分布的负荷随机性，使得在电网中存在大量非正态分布和离散分布的随机变量的情景下，计算精度较传统点估计法得以有效提高。
1 基于Nataf变换的多点估计算法研究思路
1.1 Nataf变换的基本思想
Nataf变换建立的是一个实现空间转换的数学模型，可以完成从原始变量空间到独立标准正态空间的转换。设服从任意分布类型的相关非正态随机向量：
式中：m为随机变量的个数；xi(i=1,2，…，m）的概率密度函数和累积分布函数已知。
标准正态随机向量Y=[y1,y2，…，ym]可通过式（2）表征的等概率转换原则[18,19]引入：
式中：i=1,2，…，m;Fi(xi）为xi的累积分布函数；Ψ（·）和Ψ-1（·）分别为标准正态累积分布函数和逆累积分布函数。
根据隐函数的求导法则，易推导出关于X的联合概率密度函数：
式中：f(xi）为xi的概率密度函数；Y为对应于期望值为0、标准差为1及相关系数为ρ0的m维标准正态分布。
设ρ为输入随机向量X的相关系数矩阵，ρ0为标准正态随机向量Y的相关系数矩阵，则根据相关系数的定义及式（3）和式（4），可得到相关系数矩阵中各元素的计算公式：
式中：ρ0ij为ρ0中的元素；xi和xj由式（2）可分别转化为Fi-1（Ψ（yi））和Fj-1（Ψ（yj））。
由式（5）所建立的ρ与ρ0的映射关系可知，在获得ρ和输入随机变量xi的边缘概率密度函数后，可通过式（5）求解非线性方程来确定ρ0。又由于ρ0是一对称矩阵，通过式（6）可对其进行Choleskey分解：
式中：L0为ρ0经Choleskey分解后得到的下三角阵。
利用L0并根据式（7），可将相关的标准正态随机向量Y转换为独立的标准正态随机向量V:
同样由Nataf变换的过程，根据式（8）和式（9）也可以得到其对应的逆变换。
Nataf的正变换和逆变换所完成的空间转换过程如图1所示。而本文所采用的多点估计法正是通过Nataf逆变换的过程来获取在原始变量空间的估计点。由式（8）和式（9）可将Nataf逆变换的过程写成：
式中：N-1（·）表示Nataf逆变换。
1.2 基于Nataf变换的单变量多点估计法
对于一个单变量响应函数Z=Z(x），首先由Gauss-Hermite积分的要求，在独立标准正态空间上进行多点的选取。由于此空间上的积分点和权重系数可通过查表的方式获取，故不需要进行单独计算。以3点、5点和7点为例，表1列出了其对应的积分点和权重系数。表中：vj为第j个积分点；Pj为第j个积分节点所对应的权重系数。
根据式（10）可将响应函数Z表示为：
响应函数Z统计矩的积分满足式（12）至式（14)[22]，由Gauss-Hermite积分公式可将式（12）至式（14）转化成式（15）至式（17）来求解响应函数Z的统计信息：
式中：μZ，σZ，αkZ分别为Z(x）的期望值、标准差和k阶中心矩；n为Gauss-Hermite积分中所选取的节点数。
1.3 基于Nataf变换的多变量多点估计法
式（15）至式（17）针对的是单变量下的点估计算法，而在实际的工程应用中，常遇到多变量的响应问题，即响应函数Z是多个输入随机变量X=[x1,x2，…，xm]的函数，故针对多变量的多点估计法更具实用价值。多变量函数所对应的点估计法可通过单变量函数的点估计法进行扩展。在该情形下，若采用直接计算，每个变量取n个估计点，则m个随机变量共构成了nm种组合，进而响应函数Z将进行nm次计算，计算量将会随着随机变量数目的增长呈现指数增长趋势。针对该问题，本文采用文献[17]提出的逼近方法来进行简化计算。
式中：ZΣ为含多变量的响应函数；Zμ=Z（μ1，μ2，…，μm）为响应函数在所有随机变量均取为期望值时的响应值；Zi为第i个随机变量xi的函数，而除xi外其余随机变量均取期望值，即Zi=Z（μ1，μ2，…，xi，…，μm），此时Zi可视为仅关于xi的单变量响应函数，根据式（11）可将Zi表示为
从而由式（15）至式（17）可以得到对应的关于Zi的期望值、标准差和k阶中心矩。
在获取Zi的统计信息后，即可进一步根据式（20）至式（22）将单变量响应函数的点估计方法扩展到多变量响应函数上：
式中：为多变量函数Z的统计矩。该逼近方法估计点数为nm，对应计算量为nm次，但其中有n次是在各变量均取期望值时进行计算，故实际计算可降至（n-1)m+1次。
2 基于Nataf变换的不确定性潮流多点估计算法实现
2.1 不确定性潮流建模
设某电网有s个节点、b条支路，并且其中有p个PQ节点、q个PV节点、1个平衡节点。采用极坐标形式时，对应的系统潮流方程为：
式中：Pi和Qi分别为节点i处的注入有功功率和注入无功功率；θij为节点i与节点j之间的相角差；Ui为节点i的电压幅值；Gij和Bij分别为导纳矩阵各元素的实部与虚部；Pij和Qij分别为连接节点i与节点j线路的有功功率和无功功率；tij为支路的非标准变比；Yij0为支路对地导纳。
对式（23）和式（24）进行简化，可得到式（25）和式（26）的形式：
式中：X=[P1,P2，…，Ps-1,Q1,Q2，…，Qp]为节点注入向量，由节点注入有功功率Pi和节点注入无功功率Qi构成；Z为任一支路的有功或无功功率；g（·）和h（·）分别为节点功率和支路功率关于节点电压的表达式。
根据式（25）可以得到节点电压U与节点注入量X的关系：
进而根据式（26）得到：
即支路潮流可表征为关于多变量xi的响应函数，式中m表示注入随机变量的个数，可根据实际情况确定。如在含有p个PQ节点、q个PV节点的网络中，当注入有功和无功功率均作为随机变量考虑时，m=2p+q。
2.2 不确定性潮流的多点估计法实现步骤
由式（28）可知，支路潮流表征为多变量xi的响应函数后，即可根据多变量下的多点估计算法来进行计算。算法步骤如下。
步骤1：获取m维注入随机变量X的信息。
步骤2：选择多点估计法中每个变量的估计点数n，对应于3点、5点及7点估计法时n分别为3,5,7。
步骤3：根据表1中列出的Gauss-Hermite积分点和权重系数，得到在独立标准正态空间V～上的多个估计点vi,j和Pj，其中i=1,2，…，m;j=1,2，…，n。
步骤4：由Nataf逆变换得到vi,j在原始变量空间X～上的映射xi,j。
步骤5：在每个估计点xi,j处进行确定性潮流计算，并根据式（15）至式（17）得到单个随机变量下的μZi和σZi。
步骤6：根据式（20）至式（22）的逼近方法得到多变量响应函数下的期望值μZΣ和σZΣ。
3 算例分析
3.1 计算网络随机性参数设定
为验证基于Nataf变换的多点估计法的准确性，本文采用MATLAB语言，在Windows环境下对IEEE 14节点系统算例进行了测试。该系统拓扑结构如图2所示，共有20条线路，其中节点1和节点2上接有发电机。节点1为平衡节点，节点10和节点14上接入风电机组，并假定各节点注入量相互独立。
该系统中常规发电机出力假定为二项分布，随机特性见附录A表A1。单台风电机组的额定出力设为1.5MW，节点10和节点14上各接有10台风电机组，风电机组参数设定见文献[24]，风速统计模型则采用双参数Weibull分布，其中尺度参数c取为10m/s，形状参数k取为2。负荷节点方面除节点9服从附录A表A2所示的离散分布外，其余负荷节点均服从附录A表A3所示的正态分布。
3.2 计算标准选取
本文以蒙特卡洛模拟法的计算结果作为比较标准，分别将传统的3点估计法和改进的多点估计法进行比较分析。为直观比较各方法的计算精度，根据式（29）定义了相对误差指标：
式中：RMCS和RPEM分别为蒙特卡洛模拟法和各点估计法的计算结果。
为得到精确的比较标准，首先采用蒙特卡洛模拟法对该系统分别进行5 000,10 000,50 000,60 000,70 000次仿真。附录B中列出了蒙特卡洛模拟法各仿真次数下部分具有代表性的计算值。
通过对蒙特卡洛模拟法的灵敏度分析可知，当进行60 000次仿真时其精度已合适，继续增加仿真次数后其计算结果的精度没有太大提高，故选取仿真60 000次的蒙特卡洛模拟法计算结果作为传统3点估计法及改进多点估计法的评价标准。
3.3 计算结果分析
表2和表3分别列出了本文提出的基于Nataf变换和Gauss-Hermite积分技术的3点估计法（N3PEM）、5点估计法（N5PEM）、7点估计法（N7PEM），蒙特卡洛模拟法和传统的3点估计法（3PEM）这5种方法下，部分支路潮流有功和无功功率的计算值。由于各方法下节点电压幅值的差异极小，故在表4中只列出部分节点电压相角的计算值。
根据计算结果进行分析可知，在电压幅值和相角以及支路有功和无功功率的期望值求取上，4种点估计法的计算结果基本一致，且与蒙特卡洛模拟法的计算标准相比也非常接近。但在标准差的求取上，4种点估计法下节点电压幅值和相角，以及支路有功和无功功率的标准差仍有较大差异。为方便比较各方法的计算精度，表5中列出了几种方法下各标准差计算值与蒙特卡洛模拟法计算标准的最大相对误差。
从表5中的计算结果可以清晰地看到，本文所提出的基于Nataf变换和Gauss-Hermite积分技术的点估计法计算精度较传统的点估计法更高，且随着估计点数的增多，各指标的计算值精度整体上有所提升，特别是在支路有功和无功功率标准差的计算值上，多点估计法与蒙特卡洛模拟法的最大相对误差控制得很好。
另一方面，蒙特卡洛模拟法（60 000次）的计算时间为778.10s，而改进的3点估计法、5点估计法及7点估计法的计算时间分别为0.66,0.96,1.65s，与蒙特卡洛模拟法相比，计算效率优势明显，且与传统3点估计法的计算时间0.23s相比，改进多点估计法增加的计算时间也完全可以接受。
4 结语
本文提出一种改进的多点估计法来研究不确定性潮流问题。首先建立了不确定性潮流的计算模型，将潮流的不确定性问题转化为一个求解多变量响应函数的问题；然后根据Gauss-Hermite积分技术选择对应的在标准正态空间分布的估计点和权重，并通过Nataf逆变换实现空间转换得到原始空间的估计点；最后根据转换后的估计点和权重得到支路有功功率和无功功率的估计信息。
含多点接入可再生能源的算例分析结果表明：与传统的点估计法相比，本文算法可获得更高的精度，特别是在系统中存在大量非正态分布和离散分布的随机变量情况下，计算结果的精度会随着选取估计点数的增加而提高。因此，本文提出的基于Nataf变换和Gauss-Hermite积分的多点估计法对于考虑负荷、发电机出力波动，以及可再生能源的多点接入等因素后的不确定性潮流研究具有一定的实用价值。本文中注入随机变量是在假定为相互独立的前提下考虑的，而在实际电网运行中，电源与电源以及电源与负荷之间存在一定的相关性，故下一步将针对这一问题进行深入研究。
附录见本刊网络版（http：//www.aeps-info.com/aeps/ch/index.aspx）。
估计不确定性 篇2
基于非线性约束条件研究了一类不能求解析解的微分方程系统模型的`参数估计方法.通过非线性系统模型和观测模型的线性化,用线性模型理论推导了非线性约束条件下的参数估计方法,设计了全局最优解的高斯牛顿迭代算法.某卫星轨道确定算例表明:该参数估计法的卫星轨道确定精度较传统最小二乘法提高约50%.?
作 者：刘靖 潘晓刚 LIU Jing PAN Xiao-gang 作者单位：刘靖,LIU Jing(湖南涉外经济学院,湖南,长沙,410205)
潘晓刚,PAN Xiao-gang(国防科学技术大学,理学院,湖南,长沙,410073)
估计不确定性 篇3
1 聚类有效性分析与聚类个数估计
1.1 聚类有效性分析
现在常用的估计聚类个数的有效性评价方法大概有3种:有效性指标、检测稳定聚类结构的稳定性方法和系统演化方法。其中,有效性指标在实际中应用广泛,包括Silhouette指标、Davies-Bouldin指标、Calinski-Harabasz指标、加权的类间类内相似度比率和Homogeneity-Separation指标等。
1.1.1 平均Silhouette(Sil)指标。
设a(t)为聚类Cj中的样本t与类内所有其他样本的平均不相似度或距离,d(t,Ci)为样本t到另一个类Ci的所有样本的平均不相似度或距离,则b(t)=min{d(t,Ci)},i=1,2,……,k;i≠j[1]。Sil指标计算每个样本与同一聚类中样本的不相似度以及与其他聚类中样本的不相似度,其每个样本t的计算公式如下:
一般以一个数据集的所有样本的平均Sil值来评价聚类结果的质量,Sil指标越大,表示聚类质量越好,其最大值对应的类数作为最优的聚类个数。
1.1.2 Davies-Bouldin(DBI)指数。
DBI指数是基于样本的类内散度与各聚类中心间距的测度,进行类数估计时其最小值对应的类数作为最优的聚类个数。设DWi表示聚类Ci的所有样本到其聚类中心的平均距离,DCij表示聚类Ci和聚类Cj中心之间的距离,则DB指标如下:
1.1.3 Calinski-Harabasz(CH)指标。
CH指标(伪F统计量)是基于全部样本的类内离差矩阵与类间离差矩阵的测度,其最大值对应的类数作为最优的聚类个数,则CH指标如下:
其中,SSB和SSw分别表示类间离差矩阵和类内离差矩阵,k为聚类数目。类间离差矩阵SSB定义如下:
其中,mi,m分别表示第i个类的中心和类内平均,||mi-m||定义了向量之间的欧式距离。类内离差矩阵SSw的定义如下(也称为Hartigan指标):
其中,x,Ci分别表示数据点和第i个类。
从定义可以看出:好的聚类结果,是希望SSB越大越好和SSw越小越好。因此,CH(k)的值越大,选择的k值越接近最优。
1.2 聚类数目对有效性的影响
对于盲聚类分析,聚类数目和聚类初始中心对于结果的影响很明显,因为目前采用的聚类算法,类似于K-均值聚类,都对类数和初始聚类中心敏感。
不失一般性,对同一组数据集合,分别采用目标类数为3~5进行K-均值聚类,仿真的对比结果如图1所示。
为了对比不同聚类数目下的聚类结果的性能指标,对上述实验数据,分别采用2.1节中的指标进行评估,结果如图2所示。从图2可以看出,前3种指标显示,选择类数为4最接近最优的聚类目标。但是采用Hartigan指标,其结果不然。
上述实验中,不同的指标所侧重的测度是不一致的,从而适用的应用场景也是有区别的。另外,从相邻类数下指标的走势可以看出,K-均值算法对于类数敏感性显而易见。因此,在聚类之前,能够确定有效的聚类数目是非常必要和重要的。
2 基于局部密度估计的聚类个数确定
给定点集S,对于每一个样本点i,为其定义2个属性参量:局部密度ρi和距离高密度点集的距离δi。样本点i的局部目的ρi定义如下:
其中,dij为样本点i和样本点j之间的“距离”度量,并且满足三角不等式;dc是常数参数,定义为截止距离,其物理意义是控制样本点的影响最远距离。式(6)中函数R(x)的定义如式(7):
其中,ρi表示距离i的距离小于dc的样本点数量。i的参数δi计算如下:
其中,max(ρ)为点集S的最大局部密度。从式(8)可以看出,对于非最大局部密度样本点,δi定义为样本点i距离其他任意更高局部密度样本点的最小距离;而对于最大局部密度样本点,δi定义为其他点到此点的最大距离。
通过对点集S中的样本经行二维空间的处理,以ρi和δi两个指标对点集降维处理。其中,ρi用来限制点集在多维测度空间内的局部密度,而δi用来度量局部密度较大的多维点集在测度空间内的相对聚类。距离其他高局部密度较“远”的点,称为在决策图上的“孤点”。“孤点”的个数即是聚类的估计类数。在局部密度的估计中,截止距离dc的选择非常重要,其影响了局部密度测度空间的边界。对于dc的选择准则,还需要进一步的理论分析,但是多次实验证明,dc的选择需要保证局部密度的最小取值不小于点集总数的1%~2%。
3 实验结果与分析
对于2.2中所采样的同样的数据数据集合,使用2中介绍的基于局部密度估计的方法,计算点集中每个样本的ρi和δi两个指标,绘出决策图如图3所示。
从图3可以看出,在ρi和δi两个指标都很“突出”的区域有4个明显的“孤”点。所以,估计算法确定出的聚类数目为k=4。对于2.2节的结论,两者的结论是一致的。通过对多组数据的实验,该文所提出的算法估计出的聚类类数和Sil指标、DBI指标和CH指标的一致率分别是99.4%、99.1%和98.8%。
4 结论
从K-均值算法对于聚类个数和初始中心敏感的缺陷出发,基于聚类中心具有高密度且距高局部密度聚类中心具有较远距离的特点,提出一种基于局部密度估计的聚类个数的估计方法。经过仿真实验,验证了该算法具有良好的有效性和鲁棒性。
摘要：随着人工智能和数据挖掘技术的兴起,聚类分析已被广泛应用于通信、文本数据统计、生物信息学和图像处理中。对于非监督聚类分析,聚类的分类数目是决定聚类质量的关键因素。通常聚类个数事先无法确定,随即选择的初始聚类中心容易使聚类结果不稳定。针对此,基于聚类中心具有高局部密度且距高局部密度聚类中心距离较远的特点,提出一种基于局部密度估计的聚类个数的估计方法。经过仿真实验,验证了该算法具有良好的有效性和鲁棒性。
估计不确定性 篇4
现有状态估计[1,2,3,4]和不良数据辨识[5,6,7]方法主要建立在测点残差概念基础上,但是,测点残差并不反映测点真值任何信息,对于测量结果的可靠程度无法进行定量描述。这给进一步研究状态估计问题带来了较大困难,导致一些基础性问题,如真值未知情况下状态估计结果评价问题,都不能得到解决。
本文基于文献[8]提出的正常测点概念,构造了测点评价函数,用于评价测点是正常测点还是异常测点;然后以测点正常率最大为目标函数,提出了新的状态估计模型及其求解算法,并对该状态估计方法的修正模型进行了进一步讨论。
1 测点评价函数
设x为状态估计结果,对测点i,称di为测点i在估计x下的相对偏离:
式中:zi为测点i的量测值;hi(·)为测点i的量测函数;Ui为与置信概率p对应的测点i的扩展不确定度。
由正常测点定义可知:测点i为异常测点时,|di|>1;测点i为正常测点时,|di|≤1。
根据上述思想,可以建立如下测点评价函数:
定义1:称
显然,测点正常率最高的估计,其状态正常率评价指标最小;反之,状态正常率评价指标最小的估计,其测点正常率最高。因而,可将求测点正常率最高的估计问题转化为求正常率评价指标最小的估计问题。
g(di)的缺点是:在±1处有间断点,并且不是处处连续可导,不便于实际应用。为此,建立g(di)的近似函数f(di),显然f(di)应具有以下特性:
式中:λ为大于0的常量。
λ越小,式(3)与式(2)在形状上越相符,但在实际应用中,若λ过小,相当于函数f(di)在|di|=1处有跃变,对以
式(3)的意义为:当|di|≤1时,认为测点i为正常测点,f(di)接近于0;当|di|≥1+λ时,认为测点i为异常测点,f(di)接近于1;当1≤|di|≤1+λ时,认为测点i为可疑测点。
本文中,f(di)取以下形式:
式中:
k一般取大于2的自然数。k取不同值时,f(±1)和f(±(1+λ))的取值见表1。由表1可知,k越大,f(±1)的值越接近于0,f(±(1+λ))的值越接近于1,式(3)与式(2)在形式上越相符。但k过大会给以
λ=5,k=3时,函数f(di)的示例图形见图1。由图1可见,f(di)没有间断点,处处连续可导,且f(di)与g(di)的特性接近,即当测点i正常时,f(di)接近于0,反之f(di)接近于1。因此,本文采用式(4)定义的f(di)作为测点评价函数。
测点评价函数也可采用其他形式,只要保证处处连续可导,且当测点i为正常测点时,测点评价函数值接近于0,反之接近于1即可。
2 方法介绍
2.1 数学模型
设某一系统状态下,系统中异常测点数目为q,则
由此,得到新的状态估计模型:
由于实际运行中,真实系统状态必然满足潮流约束和其他实际物理约束(如发电机功率上下界约束)。考虑这些约束后,得到改进的状态估计模型:
式中:g(x)=0代表潮流约束;l(x)≤0代表实际物理约束。
2.2 求解算法
式(6)中,将di代入目标函数中,即得:
式(7)是一最优潮流问题,其求解方法较多。对于大规模最优潮流问题,用现代内点法[9]求解较好,其收敛性能得到保证,计算速度快,且当问题规模增大时,计算时间增加不多[10,11]。
本文基于开源现代内点法优化计算软件包[11,12]成功求解了以上问题,收敛性好,计算速度快,适合于在线工程应用。详细求解过程与文献[10,11]类似,此处不赘述。
3 方法特点
与以往状态估计方法相比,本文方法特点如下:
1)以测点正常率最大为目标,因而所求得的状态估计结果的测点正常率高。
2)可在求得状态估计解的同时自动辨识系统中的不良数据。本文方法以测点正常率最大,即以测点异常率最小为目标。可以认为,在本文方法所求得的估计下异常测点即为系统中的不良数据,因而,该方法在求得状态估计解的同时,也自动对系统中的不良数据进行了辨识。
3)估计准确性不易受不良数据影响,具有很强的抗差性。以往状态估计方法容易受不良数据影响,且估计值偏离量测值越远,受影响越大。本文方法中,当测点为异常测点时,无论估计值偏离量测值多远,其反映在目标函数中大小都是1,因而估计准确性不易受不良数据影响,具有很强的抗差性。
4)本文方法还具有以下特点:①所求得状态估计结果为潮流解,且满足各种物理约束,更加接近实际;②计算中无需进行不良数据校验、权重因子设置,调试和维护极为简单;③求解算法为现代内点法,收敛性好,计算速度快。
4 进一步讨论
本文系列文章(共3篇)第1篇[8]中指出,测点正常率越高,状态估计结果越合理。基于此,本文以测点正常率最大为目标,提出了相应状态估计模型和求解算法。可以想见,本文方法求得结果的测点正常率较高,与真实状态较为接近。但是,本文方法所求得结果是否可以进一步改进?是否可以求得更好的解?
本文系列文章第1篇[8]中指出:当2个估计的测点正常率相同时可以降低置信概率,然后在新的置信概率下比较2个估计的测点正常率,测点正常率较高的估计其结果更为合理。
因此,要得到更合理的状态估计结果,可以采用如下方法:设求解模型(6)所得估计为
式中:
模型(8)与模型(6)相比,主要有2点不同:
1)增添了约束:dj≤1,∀j∈A。这样可以保证
2)测点评价函数由式(5)改为式(9)。即对异常测点评价不变,对正常测点更希望其估计值分布在量值区间
同时,模型(8)与模型(6)形式类似,可用同样方法求解,并且还可将模型(6)的解作为模型(8)的初始可行解,从而加快问题求解。
通过以上调整,按本文系列文章第1篇[8]中提出的状态估计结果评价方法,可以认为模型(8)得到的估计结果较模型(6)所得估计结果更为合理。
模型(8)中,去掉约束dj≤1,∀j∈A,得到以下模型:
由于去掉了约束dj≤1,∀j∈A,模型(10)可以直接求解,无需在求解模型(6)后再求解。那么,该模型所得解与模型(6)、模型(8)所得解相比,其合理性如何呢?结论如下:
1)模型(10)所得估计结果的测点正常率应不大于模型(6)所得结果的测点正常率。模型(10)与模型(6)相比,仅在于测点评价函数不同,模型(10)评价函数采用了式(9),而模型(6)采用了式(5)。模型(10)所得估计结果的正常率应不大于模型(6)所得结果的测点正常率,因为以式(5)为目标函数,即希望更多测点分布在量值区间[zi-Ui,zi+Ui]内,而以式(9)为目标函数,即希望有更多测点分布在量值区间
2)模型(10)所得估计结果的测点正常率应不大于模型(8)所得结果的测点正常率。模型(10)与模型(8)相比,少了一个约束,即可行区域变大,因而模型(10)应能求得更优目标函数值,但如前所述,这并不能使更多测点估计值分布在[zi-Ui,zi+Ui]内。而模型(8)由于加上了约束dj≤1,∀j∈A,其测点正常率不会低于模型(6)的测点正常率,因而也就不会低于模型(10)的测点正常率。
综上所述,当系统测点中不存在不良测点或不良测点数目较少时,3种模型求得的解相差不大;而当系统中不良测点较多时,模型(8)求得的解的合理性略高。
5 算例分析
本文所提出的模型中用到了参数k和λ。参数k和λ取值对算法的收敛性和状态估计结果合理性有多大影响呢?本文对此进行了研究,详细结果参见5.1节和5.2节。同时还对第4节中讨论的模型(6)、模型(8)、模型(10)的估计结果进行了对比分析,详细结果参见5.3节。
5.1 参数λ取值试验研究
为检验λ取不同值时对算法的收敛性及估计结果合理性的影响,采用3份生数据对IEEE 118节点系统进行了试验。第1份生数据中不含不良量测(量测值与真值之差绝对值大于3倍标准差的量测);第2份数据中含有32个不良量测,占量测总数的3.036%;第3份数据中含有52个不良量测,占量测总数的4.934%。试验中,测点扩展不确定度取该测点测量误差3倍标准差,各次试验使用初值相同,计算模型均采用模型(6),k均取2。
不含不良数据情况下试验结果见表2,需要说明的是,
5.2 参数k取值试验研究
为检验k取不同值时对算法的收敛性及估计结果合理性影响,采用与5.1节相同的3份生数据对IEEE 118 节点系统进行了试验。试验中,测点扩展不确定度取该测点测量误差3倍标准差,各次试验使用初值相同,计算模型均采用模型(6),λ均取4。
不含不良数据情况下试验结果见表3,含不良数据情况下试验结果见附录A表A3和表A4。从试验结果中可以看出:当系统不含不良数据或不良数据较少(小于3%)时,迭代次数少,解的合理性好,且迭代次数和合理性指标变化不大;当系统中不良数据较多(大于5%)时,解的合理性有所下降,进一步,若k取值较大(大于5),迭代次数有所增加。
5.3 模型(6),(8),(10)的比较
为检验模型(6),(8),(10)的收敛性和估计结果的合理性,采用与5.1节相同的3份生数据对IEEE 118节点系统进行了试验。试验中,k均取2,λ均取4,模型(8)和模型(10)中
不含不良数据情况下试验结果见表4,含不良数据情况下试验结果见附录A表A5和表A6。从试验结果可以看出:在不含不良量测、含较少不良量测和含较多不良量测3种情况下,3个模型的各项指标均相差不大。因此,本文系列文章第3篇(关于算法比较和实际应用)中均采用模型(6)进行计算。
与其他状态估计准则,包括加权最小二乘准则、加权最小绝对值准则、非二次准则等的比较参见本文系列文章第3篇。
6 结语
基于正常测点和异常测点概念,构造了相应测点评价函数。当测点为正常测点时,测点评价函数取0,而当测点为异常测点时,测点评价函数取1。
进一步建立了以测点正常率最大为目标函数的状态估计模型,并提出了相应求解算法。该方法具有以下特点:①所求得状态估计结果的测点正常率较高;②可自动对不良数据点进行辨识;③估计准确性不易受不良数据影响,具有很强的抗差性;④所求得状态估计结果为潮流解,且满足各种物理约束,更加接近实际;⑤计算中无需进行不良数据检验、权重因子设置,调试和维护极为简单;⑥求解算法为现代内点法,收敛性好,计算速度快。
本文还给出了所提出模型的修正模型,但算例分析表明,各模型所得结果的合理性指标基本相同。
与以往状态估计方法以及不良数据辨识方法相比,本文方法优劣如何?本文系列文章第3篇将对此进行详细的对比研究。
参考文献
[1]SCHWEPPE F C,WILDES J,ROMD B.Power systemstatic-state esti mation:PartⅠ,Ⅱ,Ⅲ.IEEE Trans on PowerApparatus and Systems,1970,89(1):120-125;125-130;130-135.
[2]BALDICK R,CLEMENTS K A,PI NJO-LAZIGAL Z,et al.I mplementing non-quadratic objective functions for stateesti mation and bad data rejection.IEEE Trans on PowerSystems,1997,12(1):376-382.
[3]CELIK M K,ABUR A.Arobust WLAVstate esti mator usingtransformations.IEEE Trans on Power Systems,1992,7(1):106-113.
[4]SI NGH H,ALVARADO F L,LI U W-H E.Constrained LAVstate esti mation using penalty functions.IEEE Trans on PowerSystems,1997,12(1):383-388.
[5]ZHANG B M,LO K L.A recursive measurement erroresti mationidentification method for bad data analysis in powersystemstate esti mation.IEEE Trans on Power Systems,1991,6(1):191-197.
[6]ZHANG B M,WANG S Y,XI ANG N D.Alinear recursivebad data identification method with real-ti me application topower systemstate esti mation.IEEE Trans on Power Systems,1992,7(3):1378-1385.
[7]KOTIGUA W W,VIDYASGAR M.Bad data rejectionproperties of weighted least absolute value techniques applied tostatic state esti mation.IEEE Trans on Power Apparatus andSystems,1982,101(4):844-851.
[8]何光宇,董树锋.基于测量不确定度的电力系统状态估计:(一)结果评价.电力系统自动化,2009,33(19):21-24,35.HE Guangyu,DONG Shufeng.Power system static-stateesti mation based on uncertainty of measurement:Part oneresult evaluation.Automation of Electric Power Systems,2009,33(19):21-24,35.
[9]WEI H,SASAKI H,KUBOKAWAJ,et al.An interior pointnonlinear programming for opti mal power flow problems with anovel data structure.IEEE Trans on Power Systems,1998,13(3):870-877.
[10]WEI H,SASAKI H,KUBOKAWA J,et al.Large scalehydrothermal opti mal power flow problems based on interiorpoint nonlinear programming.IEEE Trans on Power Systems,2000,15(1):396-403.
[11]W CHTER A,BIEGLER L T.On the i mplementation of aninterior-point filter line-search algorithm for large-scalenonlinear programming.Mathematical Programming,2006,106(1):25-57.
估计不确定性 篇5
经典测量数据处理的主要对象大多是来源同类观测量的单一观测数据,数据的利用也仅仅以单一精度观测量为依据。随着测量技术的发展,测量数据处理与利用对象已从过去的单一来源、单一精度的类型观测量发展为多类型、不同精度的多源观测量[1]。通过多源数据融合可以获取更全面、更有效的信息,产生比单一信息源更精确、更可靠的状态判断[2]。随着大地测量数据的积累,融合多源大地观测数据逐渐成为研究的热点之一。
为了融合多源观测数据,获取最优的融合估值,目前测量数据融合大多需要精确确定各类观测量的先验方差,需合理确定各类数据权重。若观测量的权存在误差或观测量间的权比失调,融合结果便不是最优解。为此,有学者提出一种基于方差分量估计的方法来解决多源观测数据融合中不同数据的权重问题[3]。该方法最大的优点是避免了烦琐且不可靠的多源观测量的先验方差估计,但是未能考虑多源观测数据含有粗差的情况。观测值中含有的粗差,会导致融合结果偏离真值。本文在基于Helmert方差分量估计定权的基础上,将抗差方差分量估计引入到多源观测数据融合中,以确定各类融合数据的权值,获得最优的融合结果。
1 基于最小二乘多源观测数据的融合模型
对于r类多源观测数据观测信息L1,L2,…,Lr,其权阵分别为P1,P2,…,Pr,误差方差分别为:
式(1)中X^为t×1维待估参数向量,A1,A2,…,Ar分别为L1,L2,…,Lr的设计矩阵,V1,V2,…,Vr分别为L1,L2,…,Lr的残差向量,L1,L2,…,Lr的维数分别为n1,n2,…,nr。
各类观测量之间随机独立,即:
根据最小二乘原理,参数求解需满足VTPV=min的要求,转置后并将式(1)带入消去V1,V2,…,Vr,求得:
将式(3)中每个式子相加后,可以求得残数解:
式(4)即为基于最小二乘多源观测信息融合模型,通过模型不难发现各类不同精度的观测信息进行融合,不同的权比将会导致不同的融合结果。
2 基于IGGIII权函数抗差Helmert方差估计
多源数据处理过程中通常采用方差—协方差分量估计的方法以确定不同种类观测量的权重;若各类观测量之间相互独立,则采用方差分量估计的方法确定各类观测量的权重[4,5]。若含有粗差的各类观测量之间相互独立,则采用抗差方差分量估计的方法确定各类观测值的权重。
同样,对于可列公式(1)所示误差方差的r类(或r种精度)的观测量,各类观测量之间随机独立,则r类观测量的方差分量估计公式为[6]:
式(5)中:
如果第i类观测量Li含有粗差,则参数估值及其方差分量估值会被扭曲。为抵制异常观测量对X^,σ2oi的影响,建立如下目标函数[7]:
式(7)中,ρ为一个凸的连续函数;(Pj)i,(Vj)i分别为Lj的权元素和残差元素。
基于等价权原理,有:
式(8)中,为Li的等价权矩阵,,wi为降权因子。选择合适等价权函数可确定相应的等价权矩阵。采用IGGIII方案[8]:
式(9)中,为标准残差;k0和k1为常量,,vi0表示通过法方程得到的v的第一次估值,li表示误差方差的常数项。
则有r类(或r种精度)观测量的抗差方差分量估计公式为[9]:
其中:
求出σ201,σ202,…,σ20r后,根据公式(12)进行迭代判断,直到:
其中,ε为一给定的正小数,本文给定0.0005。
满足条件后,根据公式(13)分别求出各类观测量的权阵P1,P2,…Pr
3 算例与讨论
以山东某煤矿工作面沉降观测数据为例说明本方法的可行性。该煤矿地表岩移监测线共布设了38个监测点,编号分别为1~38。根据《煤矿测量规程》(2011版)要求,为保证监测线首次采集水准数据的可靠性,一般需要进行至少两次重复测量以确保每个监测点之间高差的准确性,并以此作为岩移沉降监测的基准数据。本次研究首次监测点高程采集共进行了三次重复观测,三次数据采集结果每两次高差中误差分别是2.48mm、2.65mm、2.89mm,精度满足要求。由于三次观测量不等精度,因此可以将其作为三类观测量。以三类数据的加权平均值作为参考值,为考察粗差对结果的影响,分别在第二类高差段为4~5、10~11、12~13、14~15、20~21、23~24的监测数据上分别加上-2.0cm、2.5cm、2.8cm、3.1cm、2.9cm、2.7cm的粗差(见表1),作为观测异常点。
采用五种方案进行计算:(1)无粗差基于加权平均法计算平均值;(2)无粗差基于Helmert方差分量估计法定权融合;(3)有粗差的基于加权平均法计算平均值;(4)有粗差基于Helmert方差分量估计定权融合;(5)有粗差的基于抗差Helmert方差分量估计定权融合。
(1)抗差Helmert方差分量估计的单位权方差迭代结果见表2。
(2)各个方案实验结果的对比情况见表3,图1、2。
方案(1)、(4)、(5)的结果对比如图2所示。
计算结果分析:对监测点观测数据分别进行在无粗差情况下的加权平均法求平均值和基于Helmert方差分量估计法定权融合,解算各点高程。图1(a)所示的两者结果有差异,原因是各类观测数据精度不同;当加入粗差后无论加权平均法求均值还是基于Helmert方差分量估计法定权融合,都无法抵御粗差的影响,而采用基于抗差Helmert方差分量估计确定权重融合的结果则比较显著。这表明:
(1)对同一条监测线分别进行三次重复观测,精度因操作人员而异,即使同一操作人员不同时期观测,其精度也不同。对不同精度的测量数据进行融合,参与融合的各类测量数据权比不合理,对数据结果有影响。
(2)当观测量不包含粗差时,无论是基于加权平均法计算平均值还是基于Helmert方差分量法定权后再进行融合计算,其计算结果是基本一致的(根据《煤矿测量规程》(2011版),同一点高程重复观测不大于4mm,可取平均值。实验中,两种方法处理结果较差绝对值不大于4mm,为方便分析,我们认为结果一致),说明在不含粗差的多源水准观测数据融合中两种方法对结果均没有显著影响;当观测量包含粗差时,无论是基于加权平均法计算平均值,还是基于Helmert方差分量法定权后的融合,都无法抑制粗差对结果的影响,且粗差对结果影响较明显。
(3)当观测量包含粗差时,基于抗差Helmert方差分量法定权后的融合在一定程度上有效地抑制了粗差对结果的影响。
4 结论
对多类观测数据进行融合时,应考虑各类观测数据对融合结果的影响。通过对方案(1)和(2)的研究表明,融合时根据观测量精度对结果的贡献,合理地确定权重。当观测数据中含有粗差或者观测资料受到噪声的污染时,通过对方案(1)、(2)、(3)、(4)和(5)的研究表明,采用抗差方差分量估计的方法确定多源数据的权比可以有效地改善粗差对结果的影响。
摘要：为合理确定多源数据的权重,同时为减轻多源观测数据融合过程中观测资料粗差的影响,将抗差方差分量估计引入多源数据融合中。采用山东某煤矿某工作面岩移的初次监测数据进行抗差方差分量估计,以说明多源观测数据融合权比算法的有效性。结果表明,基于抗差方差分量估计确定多源观测数据权重的方法既可以有效地减弱观测资料中粗差的影响,同时可以较好地解决多源观测数据加权的不合理性,提高了数据融合精度。
估计不确定性 篇6
关键词：综采工作面,稳定率,x2检验
综合机械化采煤工作面由于机械化程度高, 工作连续性强, 效率高, 一般占全矿产量比例很高。从2005年开始, 陕西省境内新开发的矿区综采面一般占到全矿产量的75%以上。综采工作面工序接替的稳定性直接影响全矿产量的稳定性。综采工作面的主要工序有落煤, 运煤, 推移刮板运输机, 移架, 采空区处理等主要工序, 这些工序可以依据工作面作业循环表的编制事先计划好。但是, 由于采煤工作面是一个不断连续移动的工作空间, 其顶底板, 煤壁, 采空区, 地质因素不断变化, 或者是设备故障, 或者是人员原因, 或者是采煤工艺本身设计的问题, 都会引起工序接替不顺畅, 由于上道工序的延误, 挤占下道工序的时间, 使得作业循环表安排的工作时间变化, 影响工作面出煤效率。如何找到一个能从总体上衡量一个综采工作面的稳定性指标就显得很重要。同时要求构造的指标具有现场可操作性强, 衡量程度高, 也便于纵向和横向比较的特点。
1 综采面稳定率
综采面的稳定率可用如下公式表达:
其中:T正常代表工作面正常工作时间;T非正常代表工作面非正常工作时间;24代表综采面一个正常工作循环时间, 即每天的24小时。
工作面正常时间的含义是指按工作循环表安排工艺时间, 而非正常工作时间就是未按工作循环表运行工序的时间。
用稳定率衡量综采面稳定性的优点首先在于统计上的可操作性。由于影响工作面非正常工作因素众多, 无论是哪个因素影响了循环作业的进度, 都会通过工序时间滞后而体现出来, 在现场管理时只需记录未按工作循环的工序时间即可;其次, 稳定率指标也可用于纵向和横向比较。用于纵向比较时, 只要计算出被评价的综采队的以前和现在的稳定率就可以直接做出评价, 用于横向比较时只要比较两个不同综采队的稳定率即可。最后, 由于稳定率是一个没有量纲的数字, 可以直接用来衡量综采队的生产工艺稳定性。
2 稳定率的算例
下列数据为某煤矿综采面连续100天记录的非正常工作时间
2.1 做出频数直方图
对100个数据适宜分为10个组, 组距应为0.03, 为了解决端点的数据归属, 原数据在1.27~1.55之间, 现取为1.265~1.555之间做出频数分布表如下。
2.2 非正常工作时间的总体分布x2检验
2.2.1 x2检验方法
从非正常作业时间频数直方图形状看, 应该为服从正态分布。设H0:T~N (μ, σ2) , 当μ, σ2未知时, 由数据求出μ, σ2的极大似然估计值, 令
将原始数据分组组限代入上式, 求出Z组限, 比如原始数据某一组组限为 (ti, ti+1) , 相应组限为, 再将其分别记做 (ai, ai+1) , i=1…n-1。由所设正态分布计算出Z组限内每一组概率Pi=P (ai<Z<ai+1) 。
如果原假设H0:为真, 那么应该fi/n与pi差异不大, 即
皮尔逊统计量:
其中K为分组数, R为数据数。
在给出显著水平α下, 有
所以, 如果事件一旦发生。则否定H0, 即理论总体不是正态分布, 否则接受原假设, 即非正常工作时间的总体分布为正态分布。
2.2.2计算
由表2
则H0:T~N (1.406, 0.002322)
2.2.3 数据的x2检验
首先由计算出组限;再查表求出npi组合后计算出fi-npi, (fi-npi) 2, (fi-npi) 2/npi具体计算列表如下。
查x2分布表, 最后并组后, K=6, 估计了两个未知参数的值, 所以自由度为6-1-2=3:
当给出置信水平α=0.1时, 则x20.1 (3) =6.251
由于2.5235<6.251
所以在水平a=0.1下接受H0, 认为非正常工作时间分布服从正态分布。
3 非正常作业时间区间估计
由于综采面地质条件, 采煤工艺, 作业人员, 采煤设备, 是一个系统, 而非正常作业时间是此系统非稳定运行的指标, 因此, 统计作业的非正常作业时间是综采面采煤系统非正常时间的一个样本。以上x2检验, 说明综采面采煤系统非正常时间为正态分布, 因此, 可以给出非正常作业时间的区间统计, 即给出一个概率水平α=0.05有
已知σ=0.0482, 查正态分布表C=1.96代入上式即
即事件发生的概率是0.95
上式说明我们有95%的把握得出该综采队非正常工作时间的区间1.112<μ<1.700。
4 结论以及应用
(1) 对稳定率的估计, 经过x2检验, 非正常作业时间分布服从正态分布。因而对于综采工作面工序稳定性可以估计, 为煤矿评价综采工作面工序稳定性提供了客观评价的定量依据。也可以作为综采工作面外包的煤矿对分包综采工程的定量评价依据。
(2) 由于计算量大, 本文计算结果均以非正常工作作业时间的数据进行处理, 没有直接换算成稳定率计算。另外用于实际当中时原始数据统计量大, 开发相关的统计程序软件将是需要解决的问题。
参考文献
[1]张海彬.非能源矿产资源安全开发分类评价研究[D].中南大学, 2014.
[2]胡敬朋.综采工作面人-机-环境系统可靠性研究[D].安徽理工大学, 2005.