一种基于DTW的改进K-means的日负荷曲线聚类方法与流程

本发明属于电力系统分析与控制技术领域，特别涉及一种基于dtw的改进k-means的日负荷曲线聚类方法。

背景技术：

日负荷曲线聚类是配用电大数据挖掘的基础，对负荷预测、电网规划和需求侧响应均有一定的指导意义。随着智能电网的不断推进，电力系统信息化程度不断提高，用电信息采集系统、配网gis系统、配网自动化系统等逐渐完善，配用电数据呈现出数据量大、类型多、增长快等大数据特征。如何采取有效的数据挖掘技术，在大数据背景下对不同类型的海量用户进行精细化划分，从而挖掘出不同类型负荷间的内在联系及对应的用电行为、用电特性等信息，无疑对电网公司及电力用户具有重要的意义。

传统的日负荷曲线聚类方法通常以日负荷曲线各采样时刻点的功率值经极大值归一化后，采用k-means、模糊c均值等算法以欧式距离作为相似性判据对日负荷曲线进行聚类。该类方法存在如下两个弊端：1)作为时间序列的负荷曲线来说，曲线之间的相似性易受气温气候、收入、电价政策等许多因素的影响，不能简单通过单纯的欧式距离得到充分反映；2)随着负荷数据规模的不断增长，该类方法在计算效率上面临着巨大的挑战。

并且现在传统聚类方法面临一个重要问题：相似度衡量准确度。传统方法，通常只是单纯选取欧氏距离衡量曲线间相似度，而忽视了曲线形态变化对曲线相似度的影响。因此如何选择合适的相似度衡量指标可以在很大程度上提升日负荷曲线聚类结果的准确性及效率。

技术实现要素：

本发明所解决的技术问题是，针对现有日负荷曲线聚类方法中存在的问题，提出一种基于dtw的改进k-means的日负荷曲线聚类方法，结合dtw与欧氏距离综合考虑日负荷曲线的分布特性与动态特性，并以这两种特征参数分别进行各自特征聚类中心的计算，最后采用加权使每条曲线与自身的分布特性聚类中心和动态特性聚类中心的相似度衡量，以此为依据采用k-means算法进行电力系统中日负荷曲线的聚类。

本发明采取的技术方案为：

一种基于dtw的改进k-means的日负荷曲线聚类方法，包括以下步骤：

步骤1)对日负荷功率曲线数据组中的异常数据进行识别与修正，获得原始数据矩阵，然后对修正后的日负荷功率曲线数据进行标幺化处理，得到标幺化数据矩阵及相应动态数据矩阵，并确定初始聚类数目、迭代次数和最大聚类数目；

步骤2)基于标幺化数据矩阵，以日负荷曲线数据组中数据对象的密度参数集合为依据，将密度参数集合中的密度最大值作为初始分布特性聚类中心曲线，进而得出当前聚类数目下的动态特性聚类中心曲线；

步骤3)以当前聚类数目开始进行循环计算：采用欧氏距离与dtw，分别计算每条日负荷曲线关于分布特性聚类中心曲线及动态特性聚类中心曲线的相似度，再以相似度数值为根据，进行数据组的归类，然后对分布特性与动态特性聚类中心曲线进行更新，并将分布特性与动态特性相似度值的加权值作为k-means算法中曲线相似度的衡量依据；

步骤4)结合步骤3)所得更新后的聚类中心曲线，计算日负荷曲线的综合聚类离散度，若聚类中心曲线或聚类离散度的误差值已满足要求，则停止聚类中心曲线的更新，否则返回步骤3)继续对聚类中心曲线进行更新，直至满足要求为止；然后得到当前循环对应聚类数目下的聚类中心曲线及对应dbi的值，即完成本次循环，再对聚类数目加一后检查是否达到最大聚类数目，没有达到则返回步骤2)，计算下一个指定聚类数目的聚类中心曲线及对应dbi值，达到则执行步骤5)，其中dbi是类内距离之和与类外距离的比值；

5)结合迭代完成后所得不同聚类数目下的聚类结果，基于dbi确定最佳聚类数目与最终聚类结果，其中dbi越小则聚类结果越好，并以选取dbi数值最小时所对应的聚类数目为最佳聚类数目。

所述的一种基于dtw的改进k-means的日负荷曲线聚类方法，所述步骤1)中，对日负荷功率曲线中的异常数据进行识别，是计算日负荷功率曲线中各采样时刻点的负荷功率的变化率并与预设阈值比较，若超出阈值范围则视为异常数据；对异常数据进行进行修正，是首先判断负荷曲线的数据缺失量和异常量是否小于10％，若是，则先将异常量置为0即置为缺失量，再以一元三点抛物线插值算法对缺失量进行插值拟合，否则直接删除该条负荷曲线。

所述的一种基于dtw的改进k-means的日负荷曲线聚类方法，述步骤1)中，对修正后的日负荷功率曲线数据进行标幺化处理，得到动态数据矩阵的方法为：

记pk＝[pk1,...,pki,...,pkm]∈r^1×m为修正后第k条日负荷功率曲线的m点原始有功功率矩阵，k＝1,2,3,…,n，n为日负荷功率曲线总条数，pki为第k条日负荷功率曲线的第i点原始有功功率，i＝1,2,…,m，m为采样点个数；则p＝[p1,...,pk,...,pn]^t∈r^n×m为n条日负荷功率曲线的m点原始有功功率矩阵，其中t为矩阵转置符号；

取日负荷功率曲线的功率最大值pk.max＝max{pk1,pk2,...,pki,...,pkm}为基准值，根据下式对原始数据样本进行标幺化处理，得到功率点pki的标幺值p'ki：

p'ki＝pki/pk.max

进而得到归一化的日负荷功率曲线有功功率标幺值矩阵p'k：

p'k＝[p'k1,p'k2,...,p'ki,...,p'km]∈r^1×m，则标幺化后的日负荷功率曲线矩阵a＝[p1',...,pk',...,pn']^t∈r^n×m；

然后由上式获取日负荷曲线有功功率动态数据矩阵：

f'k＝[fk1,fk2,...,fki,...,fk(m-1)]∈r^1×(m-1)；

其中，fki为第k条日负荷功率曲线的第i点的日负荷曲线斜率数据，δt为相邻有功功率采样点的时间间隔。

所述的一种基于dtw的改进k-means的日负荷曲线聚类方法，所述的初始聚类数目lmin为2、迭代次数为1000，并预设最大聚类数目n为日负荷功率曲线数据组中的总曲线数目。

所述的一种基于dtw的改进k-means的日负荷曲线聚类方法，所述的步骤2)包括以下过程：

基于标幺化数据矩阵a，得到初始分布特性聚类中心：

首先计算任意两负荷曲线px＝(px1,px2,...,pxn)，py＝(py1,py2,...,pyn)之间的欧氏距离d：

再计算得到日负荷曲线之间的平均距离dist：

其中为n条日负荷曲线之间的组合数；

然后计算以每条日负荷曲线为中心基于平均距离的密度参数den，

其中当dist-d(pi,p)＜0时，u(dist-d(pi,p))＝0；当dist-d(pi,p)≥0时，u(dist-d(pi,p))＝1；

进而得到日负荷曲线的密度参数集合d，根据当前聚类数目l选取前l个密度最大者，且每次选取一个密度最大者之后，将该曲线和与该曲线欧式距离小于dist的负荷曲线从标幺化数据矩阵a中删除，从而得出日负荷曲线的初始分布特性聚类中心曲线集合cl＝(c1,c2,...,cl)；

通过下式得到动态特性聚类中心曲线集合cl'＝(c1',c2',...,cl')：

ci,j'＝ci,j+1-ci,j

式中，ci,j为负荷曲线的第i类分布特性聚类中心ci的第j维数据，ci,j'为负荷曲线的第i类动态特性聚类中心ci'的第j维数据。

所述的一种基于dtw的改进k-means的日负荷曲线聚类方法，所述的步骤3)中，采用欧氏距离与dtw，来计算每条日负荷曲线与分布特性聚类中心曲线以及动态特性聚类中心曲线的相似度的过程为：

先计算每条日负荷曲线与分布特性聚类中心曲线相似度d0和动态特性聚类中心曲线的相似度d1：

d0(i,j)＝d(pi,cj)

d1(i,j)＝dtw(fi',cj')

其中d0(i,j)表示第i条日负荷曲线pi关于第j个分布特性聚类中心cj的相似度，d1(i,j)表示第i条日负荷动态曲线fi'关于第j个动态特性聚类中心cj'的相似度；

然后通过加权得到日负荷曲线的综合相似度dall，其中dall越小表示相似度越高，最后将日负荷曲线归类于其与所有聚类中心曲线中dall最小的类别：

dall(i,j)＝a*d0(i,j)+b*d1(i,j)

式中dall(i,j)表示第i条日负荷曲线pi与第j类聚类中心曲线的综合相似度；a与b分别是日负荷曲线分布特性相似度与动态特性相似度的权重；

将每条曲线根据归于不同类别之后，对数据组分布特性与动态特性聚类中心曲线进行更新：

按照下式对分布特性聚类中心曲线进行更新；

式中cj(i+1)表示第i+1次迭代所得的分布特性聚类中心曲线；nj表示属于第j类聚类中心的日负荷曲线条数；pn表示属于第j类聚类中心的日负荷曲线；

按照下式对动态特性聚类中心曲线进行更新；

ci,j'＝ci+1,j-ci,j

式中，ci,j'表示第j类动态特性聚类中心曲线的第i点；ci+1,j表示第j类分布特性聚类中心曲线的第i+1点；ci,j表示第j类分布特性聚类中心曲线的第i点；欧氏距离衡量相似度的方法如所示；

其中基于dtw衡量相似度的方法为：

设有两负荷曲线px＝(px1,px2,...,pxn)，py＝(py1,py2,...,pym)，则有：

a)构建n*m的距离矩阵dn*m，其中元素d(x,y)如下式所示；

b)将矩阵d中每一组相邻元素组成的集合称为弯曲路径，记为k＝(k1,k2,k3,..ks..,kg)，其中g为路径中元素的总个数，元素ks为路径上第s个点的坐标，即ks＝(x,y)；

c)上述的路径k有多条，然后通过dtw找到1条最优弯曲路径，使得序列px和py的弯曲总代价最小，即：

式中，d(ks)为弯曲路径的累积距离，通过动态规划方法来构造一个累积代价矩阵g来求解上式，即gxy＝d(x,y)+min(g(x-1,y-1),g(x,y-1),g(x-1,y))，其中，x＝1,2,...,n；y＝1,2,...m；g(0,0)＝0；g(x,0)＝g(o,y)＝+∞，可知，时间序列px和py的动态时间弯曲距离为dtw(px,py)＝g(n,m)。

所述的一种基于dtw的改进k-means的日负荷曲线聚类方法，所述的步骤4)中，综合聚类离散度通过以下公式计算：

式中，jn为第n次迭代所得聚类结果的综合聚类离散度；ci(n)为第n次迭代所得的第i个聚类中心；pn为属于聚类中心曲线ci(n)这一类的日负荷曲线；l为聚类数目；

若jn+1-jn＜ω，则聚类停止，其中ω为预设的收敛极限值。

所述的一种基于dtw的改进k-means的日负荷曲线聚类方法，所述的步骤5)中以基于dbi指标确定最佳聚类数与最终聚类结果的的方法如下：

式中，ri用来衡量第i类与第j类的相似度；l为聚类数目；

其中，si用来度量第i个类中数据点的分散程度，通过下式计算：

其中，xl为第i类中第l个数据点；ci为第i类的中心；t为第i类中数据点的个数；q取1时si为各点到中心的距离的均值，q取2时si为各点到中心距离的标准差，均用于衡量类内分散程度；

其中mij为第i类中心与第j类中心的距离；ci为第i类的中心；adi为第i类的中心点ci的第d个属性的值；p取1时表示1-范数，p取2时表示2-范数，即两个类中心的欧式距离。

本发明的技术效果在于，综合考虑了日负荷曲线的分布特征与动态特征，在大数据背景下可以很大程度上提升日负荷曲线聚类的效率及质量。聚类结果与工程实际相符，能够为电网公司分析用户用电行为，制定合理的用电计划提供有力的支撑。具有良好的应用前景。

附图说明

图1为本方法流程示意图。

图2为dtw路径示意图。

图3为改进k-means算法流程图。

图4为基于dbi指标确定最佳聚类数及聚类结果流程图。

具体实施方式

下面结合附图对本发明作进一步说明：

本实施例的总体思路框图如图1所示，包括以下步骤：

1)对日负荷功率曲线数据进行预处理，获得标幺化数据矩阵a∈r^n×m及其动态数据矩阵a'∈r^n×(m-1)其中n为日负荷曲线条数，m为维数，并确定初始聚类数目为lmin为2、迭代次数为1000，并预设n为日负荷功率曲线数据组中的总曲线数目。算法以lmin为初始聚类数目开始进行循环算法，每次循环完成之后聚类数目加一，直至聚类数目达到lmax；

2)结合步骤1)所得标幺化数据矩阵a，以日负荷曲线数据组的密度函数为依据选取当前聚类数目下的初始分布特性聚类中心cl＝(c1,c2,...,cl)与动态特性聚类中心cl'＝(c1',c2',...,cl')；

3)将每条日负荷曲线分别采用欧氏距离与dtw衡量其关于分布特性聚类中心曲线与动态特性聚类中心曲线的相似度，根据此相似度参数进行数据组分布特性与动态特性聚类中心曲线的更新，并将此两相似度的加权值作为k-means算法中曲线相似度的衡量依据；

4)结合步骤3)所得更新后的聚类中心曲线，计算日负荷曲线的综合聚类离散度，若聚类中心曲线或聚类离散度的误差值已满足要求，则停止聚类中心曲线的更新，否则返回步骤3)继续对聚类中心曲线进行更新，直至满足要求为止；然后得到当前循环对应聚类数目下的聚类中心曲线及对应dbi的值，即完成本次循环，再对聚类数目加一后检查是否达到最大聚类数目，没有达到则返回步骤2)，计算下一个指定聚类数目的聚类中心曲线及对应dbi值，达到则执行步骤5)，其中dbi是类内距离之和与类外距离的比值；

5)结合步骤4)所得不同聚类数目下的聚类结果，基于dbi指标确定最佳聚类数与最终聚类结果。

其中步骤1)包括以下步骤：

1-1)对日负荷功率曲线中的异常数据进行识别与修正；

1-2)对修正后的日负荷功率曲线数据进行标幺化处理，并得到相应动态数据矩阵a'；

对于以上步骤进行相关解释如下：

所述步骤1-1)中异常数据的识别方法具体为：

记pk＝[pk，1,pk,2…,pk,m]^t为某条负荷曲线在各采样时刻点的功率值，以公式(1)对异常数据进行识别。

式中：δk,i为负荷曲线在第i点的负荷功率变化率，当其超过预设的阀值ε后视为异常数据，不失一般性ε可取0.5～0.8。

所述步骤1-1)中异常数据的修正方法具体为：

若某条负荷曲线的数据缺失量和异常量达到10％或以上时，认定该曲线无效直接删除该条负荷曲线。

若某条负荷曲线的数据缺失量和异常量低于10％时，将异常量置为0，再以一元三点抛物线插值算法对缺失量进行插值拟合。一元三点抛物线插值算法的原理为：

设n个节点xi(i＝0,1,…,n-1)的函数值为yi＝f(xi)，有x0＜x1＜…＜xn-1，对应函数值y0＜y1＜…＜yn-1。为计算指定的插值点t的近似函数值z＝f(t)，选择最靠近t的3个节点：xk-1、xk、xk+1(xk＜t＜xk+1)，然后根据抛物线插值公式(2)计算z的值，即

式中，当|xk-t|＜|t-xk+1|时，m＝k-1；当|xk-t|＞|t-xk+1|时，m＝k。

若插值点t不在包含n个节点的区间内，则只选取区间某一端的2个节点来进行线性插值。

所述步骤1-2)中对修正后的日负荷功率曲线数据进行标幺化处理的方法具体为：

记pk＝[pk1,...,pki,...,pkm]∈r^1×m为修正后第k条日负荷功率曲线的m点原始有功功率矩阵，k＝1,2,3,…,n，n为日负荷功率曲线总条数，pki为第k条日负荷功率曲线的第i点原始有功功率，i＝1,2,…,m，m为采样点个数，一般为48；则p＝[p1,...,pk,...,pn]^t∈r^n×m为n条日负荷功率曲线的m点原始有功功率矩阵；

取日负荷功率曲线的功率最大值pk.max＝max{pk1,pk2,...,pki,...,pkm}为基准值，根据式(3)对原始数据样本进行标幺化处理，

p'ki＝pki/pk·max(3)

得到归一化的日负荷功率曲线有功功率标幺值矩阵p'k＝[p'k1,p'k2,...,p'ki,...,p'km]∈r^1×m，并令该矩阵为a∈r^n×m。

通过式(4)，获取日负荷曲线有功功率动态数据矩阵f'k＝[fk1,fk2,...,fki,...,fk(m-1)]∈r^1×(m-1)。

式(4)中，fki为第k条日负荷功率曲线的第i点的日负荷曲线斜率数据，δt为相邻有功功率采样点的时间间隔。

2)结合步骤1)所得标幺化数据矩阵a，以日负荷曲线数据组的密度函数为依据选取初始分布特性聚类中心cl＝(c1,c2,...,cl)与动态特性聚类中心cl'＝(c1',c2',...,cl')；

2-1)依据数据对象的密度参数集合d，把d中密度最大值作为初始分布特性聚类中心；

2-2)根据初始分布特性聚类中心相应得出动态特性聚类中心

对以上步骤进行相关解释如下：

所述步骤2-1)中由标幺化数据矩阵a中每条日负荷曲线的密度参数集合得出聚类中心的方法具体为：

根据式(4)计算任意两负荷曲线px＝(px1,px2,...,pxn)，py＝(py1,py2,...,pyn)之间的欧氏距离d，再通过式(5)得到日负荷曲线之间的平均距离dist，根据式(6)计算以每条日负荷曲线为中心基于平均距离的密度参数den，得到日负荷曲线的密度参数集合d，根据当前聚类数目l选取l个密度最大者，且每次选取一个密度最大者之后，将该曲线和与该曲线欧式距离小于dist的负荷曲线从数据库a中删除，从而得出日负荷曲线初始分布特性聚类中心曲线集合cl＝(c1,c2,...,cl)。

式(5)中，为n条日负荷曲线之间的组合数。

式(6)中，当x＜0时，u(x)＝0；当x≥0式，u(x)＝1

所述步骤2-2)中由日负荷曲线初始分布特性聚类中心得出其动态特性聚类中心曲线集合cl'＝(c1',c2',...,cl')的方法具体为：

ci,j'＝ci,j+1-ci,j(7)

式中，ci,j为负荷曲线的第i类分布特性聚类中心ci的第j维数据，ci,j'为负荷曲线的第i类动态特性聚类中心ci'的第j维数据。

3)将每条日负荷曲线分别采用欧氏距离与dtw衡量其关于分布特性聚类中心曲线与动态特性聚类中心曲线的相似度，以此相似度参数进行数据组分布特性与动态特性聚类中心曲线的更新，并将此两特征参数的加权值作为k-means算法中曲线相似度的衡量依据；

所述步骤3)中采用欧氏距离与dtw衡量日负荷曲线分布特性聚类中心曲线与动态特性聚类中心曲线的相似度的方法步骤为：

3-1)先根据式(8)、式(9)计算每条日负荷曲线关于分布特性聚类中心曲线相似度d0和动态特性聚类中心曲线的相似度d1，并由此加权得到日负荷曲线的加权相似度dall，且dall越小表示相似度越高，最后将日负荷曲线归类于其关于所有聚类中心曲线中dall最小的那个类别；

d0(i,j)＝d(pi,cj)(8)

式中d0(i,j)表示第i条日负荷曲线pi与第j个分布特性聚类中心cj的相似度。

d1(i,j)＝dtw(fi',cj')(9)

式中d1(i,j)表示第i条日负荷动态曲线fi'与第j个动态特性聚类中心cj'的相似度。

dall(i,j)＝a*d0(i,j)+b*d1(i,j)(10)

式中dall(i,j)表示第i条日负荷曲线pi与第j类聚类中心曲线的综合相似度；a与b分别是日负荷曲线分布特性相似度与动态特性相似度的权重，经实验验证，通常分别取0.55、0.45聚类效果最佳。

3-2)按照式(11)与式(7)对聚类中心进行更新；

式中cj(i+1)表示第i+1次迭代所得的第j类聚类中心曲线；nj表示属于第j类聚类中心的日负荷曲线条数；pn表示属于第j类聚类中心的日负荷曲线。

所述步骤3-2)中dtw衡量相似度的方法为：

动态时间弯曲(dtw)运用动态规划思想调整时间序列不同时间点对应元素之间的关系来获取一条最优弯曲路径，使沿该路径时间序列间的距离最小，其能很好度量时间序列之间的关系。dtw路径示意图如图2所示。设有两负荷曲线px＝(px1,px2,...,pxn)，py＝(py1,py2,...,pym)，dtw算法步骤如下：

a)构建n*m的距离矩阵dn*m，其中元素d(x,y)如(12)式所示；

b)将矩阵d中每一组相邻元素组成的集合称为弯曲路径，记为k＝(k1,k2,k3,..ks..,kg)，其中g为路径中元素(图2中实心方块)的总个数，元素ks为路径上第s个点的坐标，即ks＝(x,y)；

c)上述的路径k有多条，dtw的目的在于找到1条最优弯曲路径，使得序列px和py的弯曲总代价最小，即：

式中，d(ks)为弯曲路径的累积距离，为了求解上式，通过动态规划方法来构造一个累积代价矩阵g，即gxy＝d(x,y)+min(g(x-1,y-1),g(x,y-1),g(x-1,y))其中，x＝1,2,...,n；y＝1,2,...m；g(0,0)＝0；g(x,0)＝g(0,y)＝+∞，可知，时间序列px和py的动态时间弯曲距离为dtw(px,py)＝g(n,m)。

4)结合步骤3)所得新的聚类中心，计算日负荷曲线的综合聚类离散度，若算法聚类中心不再发生变化或综合聚类离散度已满足要求，则算法跳出循环，得到当前指定聚类数目下的聚类结果；

所述步骤4)如图3所示，其中综合聚类离散度获得的方法如式(14)所示：

式中，jn为第n次迭代所得聚类结果的综合聚类离散度；ci(n)为第n次迭代所得的第i个聚类中心；pn为属于聚类中心曲线ci(n)这一类的日负荷曲线；l为聚类数目；

当jn+1-jn＜ω(ω为收敛极限值)时，则聚类停止，认为此时聚类划分合理。

5)结合步骤4)所得不同聚类数目下的聚类结果，基于dbi指标确定最佳聚类数与最终聚类结果；

dbi指标是类内距离之和与类外距离的比值。dbi指标越小，代表聚类效果越好。同时，该指标也同时用于选择算法的聚类数目，即选取dbi指标数值最小时所对应的聚类数目l。如图4所示，所述步骤5)中以基于dbi指标确定最佳聚类数与最终聚类结果的的方法如下：

式中，ri用来衡量第i类与第j类的相似度；l为聚类数目。

其中，si用来度量第i个类中数据点的分散程度，如下式所示：

其中，xl为第i类中第l个数据点；ci为第i类的中心；t为第i类中数据点的个数；q取1时si为各点到中心的距离的均值，q取2时si为各点到中心距离的标准差，两者皆可用于衡量类内分散程度。

其中mij为第i类中心与第j类中心的距离；ci为第i类的中心；adi为第i类的中心点ci的第d个属性的值；p取1时表示1-范数，p取2时表示2-范数(即两个类中心的欧式距离)。

本实施例首先从sql2012数据库里随机选取1000条日负荷曲线，并在对所选负荷曲线进行筛选之后进行归一化处理，且计算出处理之后的日负荷曲线数据的形态变化序列，采用密度函数分别选取初始分布特性聚类中心与动态特性聚类中心，并初始化聚类数目。随后，分别采用欧氏距离与动态时间序列(dtw)综合衡量负荷曲线的分布特性及衡量负荷曲线的动态特性，并分别以所得分布特性数值与动态特性数值对各自特性聚类中心进行更新，并以每条日负荷曲线此两特性的加权值(dall)作为负荷曲线与聚类中心曲线的相似度衡量手段。最后本算法以聚类中心曲线与聚类离散度的误差值，作为算法是否得出最终聚类结果的判断依据，并基于dbi指标确定最佳聚类数及聚类中心曲线(即分布特性聚类中心曲线)。