CN112149052A - 一种基于plr-dtw的日负荷曲线聚类方法 - Google Patents

Abstract

本发明公开了一种基于PLR‑DTW的日负荷曲线聚类方法。首先，利用PLR算法对被聚类日负荷曲线数据组进行降维处理，使每条曲线根据其曲线特性自适应地降至某一维度；然后，采用密度函数从原始数据库中选取初始聚类中心曲线，并对所选聚类中心曲线进行PLR降维处理；之后采用DTW对降维后的数据组与聚类中心曲线进行相似度对比分析；最后利用k‑means算法对数据组进行划分聚类，并以聚类离散度的误差值作为算法是否得出聚类结果的判断依据，最后基于DBI指标确定最佳聚类数及相应的聚类中心曲线。本方法有效解决了降维之后维度不同的曲线间的相似度衡量问题，改善了传统聚类计算效果与效率，且聚类结果与实际工程相符，具有一定的工程价值。

Description

一种基于PLR-DTW的日负荷曲线聚类方法

技术领域

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

背景技术

用户日负荷曲线可充分体现不同类型用户的负荷特性，因此对用户日负荷曲线进行聚类是对电力系统内负荷进行建模的前提。在当今电力系统中，一个可信赖的负荷模型是进行诸如负荷建模等研究的基础。随着智能电网建设进一步加深，电力系统信息化程度进一步提高，如何从记录着海量用户负荷特性的大数据平台中提取有价值的用户日负荷曲线，并进行识别、聚类，是目前电网大数据平台建设过程中所面临的重要问题，也是解决负荷模型可靠性的关键途径。

传统的日负荷曲线聚类方法通常直接对进行归一化之后的日负荷曲线进行计算处理，并利用诸如K-means、模糊C均值等算法对其进行聚类，且以欧氏距离作为相似度的评判依据。该类方法存在如下两个弊端：1)采用相应采样点之间的欧氏距离作为相似度衡量依据，只考虑了日负荷曲线的分布特性，没有考虑曲线动态特性，所以衡量精确度在极端爬坡情况下偏差较大，且当今负荷曲线的负荷点时间间隔越来越小，使得单纯求取负荷曲线间对应负荷点的欧式距离意义也随之变小；2)当今随着技术进步，日负荷曲线采样点日趋密集，使得日负荷曲线的维度较以往有很大程度提高，这也造成传统算法的运算效率大幅下降。

即现在传统聚类方法面临两个重要问题：1)相似度衡量准确度；2)日负荷曲线降维。具体而言，传统方法通常只是单纯选取欧氏距离衡量曲线间相似度，而忽视了曲线形态变化对曲线相似度的影响，且被聚类日负荷曲线的高维度在直接用于聚类时，严重影响了算法运算效率。因此选择合适的相似度衡量指标并对日负荷曲线进行保持其曲线特征的降维，可以在很大程度上提升日负荷曲线聚类结果的准确性及效率。

发明内容

本发明所解决的技术问题是，针对现有日负荷曲线聚类方法中存在的问题，提出一种基于PLR-DTW的日负荷曲线聚类方法，首先采用PLR对日负荷曲线数据组进行自适应降维处理，得到维度不尽相同的降维数据组，然后采用DTW算法计算降维后的日负荷曲线与聚类中心之间的距离，并以此为聚类算法的相似度衡量依据，最后利用聚类算法k-means对数据组进行划分聚类，以获取电力系统日负荷曲线集的聚类中心曲线。

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

一种基于PLR-DTW的日负荷曲线聚类方法，包括以下步骤：

步骤1)采集日负荷曲线数据以形成日负荷曲线数据组，对日负荷曲线数据进行标幺化处理，得到标幺化数据矩阵，并确定初始聚类数目、迭代次数和最大聚类数目；

步骤2)在标幺化数据矩阵中日负荷曲线的密度参数集合内，选取与当前聚类数目相等且密度最大的这几个值作为聚类中心曲线集，并对标幺化数据矩阵进行PLR降维处理；

步骤3)对聚类中心曲线集进行PLR降维处理，以当前聚类数目，若为第一次循环计算则以初始聚类数目开始进行循环计算，采用DTW算法分别计算降维后的标幺化数据矩阵中每条日负荷曲线关于聚类中心曲线的相似度，再以此相似度数值作为聚类算法k-means的相似度衡量依据，进行降维后的标幺化数据组的归类，并根据归类对聚类中心曲线集进行更新；

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

步骤5)结合迭代完成后所得不同聚类数目下的聚类结果，基于DBI确定最佳聚类数目与最终聚类结果，其中DBI越小则聚类结果越好，即选取DBI数值最小时所对应的聚类数目和聚类中心曲线集作为最佳聚类数目与聚类结果。

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

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

记P_k＝[p_k1,...,p_ki,...,p_km]∈R^1×m为修正后第k条日负荷曲线的m点原始有功功率矩阵，k＝1,2,3,…,N，N为日负荷曲线总条数，p_ki为第k条日负荷曲线的第i点原始有功功率，i＝1,2,…,m，m为采样点个数；则P＝[P₁,...,P_k,...,P_N]^T∈R^N×m为N条日负荷曲线的m点原始有功功率矩阵，其中T为矩阵转置符号；

取日负荷曲线的功率最大值p_k.max＝max{p_k1,p_k2,...,p_ki,...,p_km}为基准值，根据下式对原始数据样本进行标幺化处理，得到功率点p_ki的标幺值p'_ki：

p'_ki＝p_ki/p_k.max

对各功率点计算标幺值，进而得到归一化的日负荷曲线有功功率标幺值矩阵P'_k：

P'_k＝[p'_k1,p'_k2,...,p'_ki,...,p'_km]∈R^1×m，则标幺化后的日负荷曲线矩阵A＝[P₁',...,P_k',...,P_N']^T∈R^N×m。

所述的一种基于PLR-DTW的日负荷曲线聚类方法，所述步骤1)中，最大聚类数目

N为日负荷曲线数据组中的总曲线数目。

所述的一种基于PLR-DTW的日负荷曲线聚类方法，所述的步骤2)中得到聚类中心曲线集包括以下过程：

基于标幺化数据矩阵，首先计算其中任意两负荷曲线P_x＝(p_x1,p_x2,...,p_xn)，P_y＝(p_y1,p_y2,...,p_yn)之间的欧氏距离d：

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

其中

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

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

其中P_i为除P_c外的任意一条日负荷曲线，当dist-d(P_i,P_c)＜0时，u(dist-d(P_i,P_c))＝0；当dist-d(P_i,P_c)≥0时，u(dist-d(P_i,P_c))＝1；

进而得到日负荷曲线的密度参数集合D，根据当前聚类数目L选取前L个密度最大的日负荷曲线，且每选取一个后，将该日负荷曲线和与该日负荷曲线欧式距离小于dist的日负荷曲线从标幺化数据矩阵A中删除，从而得出日负荷曲线的初始聚类中心曲线集合。

所述的一种基于PLR-DTW的日负荷曲线聚类方法，所述的步骤2)和步骤3)中，进行PLR降维处理包括以下过程：

对于待降维的标幺化数据矩阵或聚类中心曲线集中第k条维度为m的标幺化日负荷曲线P'_k＝[p'_k1,...,p'_ki,p'_kj...,p'_km]，首先计算此日负荷曲线的变化序列

即日负荷曲线上每个采样点对应采样时段的变化量也即斜率集，其元素表达式如下：

然后，对于日负荷曲线上的第i时刻对应的采样点，通过计算该采样点关于其左右两侧相邻采样点变化量的差值绝对值，即计算该采样点的左斜率与右斜率差值的绝对值，以组成曲线P'_k的SEEP序列S＝(Δ1,...,Δi,...,Δm)，由于首尾两端元素分别无左斜率与右斜率，所以规定

且中间元素Δi表达式如下：

之后，设定变化量阈值R，并比较SEEP序列中的元素是否大于阈值R，当大于时则对应采样点为特征点；

最后，将通过PLR算法所得的第k条日负荷曲线的所有特征点，按照采样点的先后顺序依次连接，即得到降维后的日负荷曲线h_k＝(p'_k1,...,p'_ki,...,p'_kf)，其中f为此第k条日负荷曲线降维之后的维度；降维后的日负荷曲线的集合为H＝(h₁,...,h_i,...,h_m)。

所述的一种基于PLR-DTW的日负荷曲线聚类方法，所述的步骤3)中，采用DTW来计算降维的标幺化数据矩阵中每条日负荷曲线关于聚类中心曲线的相似度，其过程为：

首先，对所得聚类中心曲线进行PLR降维处理，得到降维聚类中心曲线集合C_L＝(c₁,c₂,...,c_L)；

然后，计算每条降维后的日负荷曲线与聚类中心曲线相似度D₁：

D₁(i,j)＝dtw(h_i,c_j)

其中D₁(i,j)表示，第i条降维后的日负荷曲线h_i关于第j个聚类中心c_j以DTW为衡量依据的相似度值；

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

设有两负荷曲线P_x＝(p_x1,p_x2,...,p_xn)，P_y＝(p_y1,p_y2,...,p_ym)，则有：

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

b)将矩阵D中每一组相邻元素组成的集合称为弯曲路径，记为K＝(K₁,K₂,K₃,..K_s..,K_g)，其中g为路径中元素的总个数，元素K_s为路径上第s个点的坐标，即K_s＝(x,y)；

c)通过DTW算法寻找弯曲路径K中的最优弯曲路径，使得序列P_x和P_y的弯曲总代价最小，即：

式中，D(K_s)为弯曲路径的累积距离，通过动态规划方法来构造一个累积代价矩阵G来求解上式，即G_xy＝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)＝+∞，则时间序列P_x和P_y的动态时间弯曲距离为DTW(P_x,P_y)＝G(n,m)。

同时，弯曲路径还需满足如下的约束：1)所选路径必须从左下角出发，到右上角结束，即K₁＝(1,1)，K_k＝(n,m)；2)每个点必须和相邻的点匹配，即若K_s＝(i,j)则K_s+1＝(a,b)必须满足0≤a-i≤1，0≤b-j≤1，同时为了避免路径在同一水平或垂直方向进行多次连续弯曲而导致过度弯曲的现象，即避免时间序列的一个点对应另一条时间序列过多的点，在已有约束的基础上加上对连续弯曲数的约束，即：

r_x≤r_x-max，r_y≤r_y-max

其中r_x、r_y分别为路径在x轴与y轴上的连续弯曲数；r_x-max，r_y-max分别为在x轴与y轴上所允许的最大连续弯曲数，其值由序列的特征及维数确定。

所述的一种基于PLR-DTW的日负荷曲线聚类方法，所述的步骤3)中，以相似度数值作为聚类算法k-means的相似度衡量依据，进行数据组的归类，并根据归类对聚类中心曲线集进行更新的步骤包括：

将每条标幺化曲线P'_k归类于与其最相似的聚类中心曲线的一类中，之后，对数据组的聚类中心曲线，按照下式进行更新：

式中C_j(I+1)表示第I+1次迭代所得的聚类中心曲线；N_j表示属于第j类聚类中心的日负荷曲线条数；P_n表示属于第j类聚类中心的日负荷曲线。

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

式中，J_N为第N次迭代所得聚类结果的综合聚类离散度；C_i(N)为第N次迭代所得的第i个聚类中心；P_n为属于聚类中心曲线C_i(N)这一类的日负荷曲线；L为聚类数目；

若J_N+1-J_N＜Ω，则认为综合聚类离散度的误差值满足要求，其中Ω为预设的收敛极限值。

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

式中，R_i用来衡量第i类曲线集的紧密程度；L为聚类数目；

其中，R_ij用来衡量第i类与第j类的相似度，S_i用来度量第i个类中数据点的分散程度，通过下式计算：

其中，X_l为第i类中第l个数据点；C_i为第i类的聚类中心曲线；T为第i类中数据点的个数；q取1时S_i为各点到中心的距离的均值，q取2时S_i为各点到中心距离的标准差，均用于衡量类内分散程度；S_j用来度量第j个类中数据点的分散程度，计算方式与S_i相同；

其中M_ij为第i类中心与第j类中心的距离；C_i为第i类的聚类中心曲线；C_j为第j类的聚类中心曲线；a_di为第i类的中心点C_i的第d个属性的值；p取1时表示1-范数，p取2时表示2-范数，即两个类中心的欧式距离。

本发明的技术效果在于，对电力系统中高维度的日负荷曲线集进行了基于PLR的降维处理，且根据DTW算法可以计算不同维度曲线的相似度的特点，利用DTW算法对降维后的数据组进行关于聚类中心曲线的相似度计算，最后聚类算法k-means以此相似度为依据进行聚类划分，并求得最终聚类中心曲线。该发明在大数据背景下可以很大程度上提升日负荷曲线聚类的效率及质量。聚类结果与工程实际相符，能够为电网公司分析用户用电行为，制定合理的用电计划提供有力的支撑。具有良好的应用前景。

附图说明

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

图2为经过PLR降维处理的日负荷曲线示意图

图3为DTW路径示意图。

图4为k-means算法流程图。

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

具体实施方式

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

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

1)对日负荷曲线数据进行预处理和标幺化之后，对其进行PLR降维处理，得到降维数据矩阵A∈R^N×m,其中N为日负荷曲线条数，m为维数，并确定初始聚类数目为L_min为2、迭代次数为1000，并预设

N为日负荷曲线数据组中的总曲线数目。算法以L_min为初始聚类数目开始进行循环算法，每次循环完成之后聚类数目加一，直至聚类数目达到L_max；

2)结合步骤1)所得标幺化数据矩阵A，以日负荷曲线数据组的密度函数为依据选取当前聚类数目下的初始聚类中心C_L＝(c₁,c₂,...,c_L)；

3)对聚类中心进行PLR降维处理，然后对每条降维日负荷曲线，以DTW算法为衡量依据，计算其关于各类聚类中心曲线的相似度，并以此相似度数值为聚类算法K-means的划分依据，对标幺化日负荷曲线进行聚类划分之后，对各类聚类中心曲线进行更新；

4)结合步骤3)所得更新后的聚类中心曲线，计算标幺化日负荷曲线集关于聚类中心的综合聚类离散度，若聚类离散度的误差值未满足要求，则返回步骤3)，对聚类中心曲线进行降维处理之后，继续进行聚类运算，直至所得聚类中心曲线集满足要求为止，否则，停止聚类中心曲线的更新，完成本次循环；然后，计算当前循环对应聚类数目下的聚类中心曲线及对应DBI的值，完成本次循环；最后，对聚类数目加一后，判断聚类数L是否达到最大聚类数目L_max，没有达到，则返回步骤2)，计算下一个指定聚类数目的聚类中心曲线及对应DBI值运算，否则，执行步骤5)；

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

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

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

1-2)对修正后的日负荷曲线数据进行标幺化处理；

1-3)对标幺化处理后的日负荷曲线数据进行PLR降维；

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

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

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

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

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

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

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

设n个节点x_i(i＝0,1,…,n-1)的函数值为y_i＝f(x_i)，有x₀＜x₁＜…＜x_n-1，对应函数值y₀＜y₁＜…＜y_n-1。为计算指定的插值点t的近似函数值z＝f(t)，选择最靠近t的3个节点：x_k-1、x_k、x_k+1(x_k＜t＜x_k+1)，然后根据抛物线插值公式(2)计算z的值，即

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

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

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

记P_k＝[p_k1,...,p_ki,...,p_km]∈R^1×m为修正后第k条日负荷曲线的m点原始有功功率矩阵，k＝1,2,3,…,N，N为日负荷曲线总条数，p_ki为第k条日负荷曲线的第i点原始有功功率，i＝1,2,…,m，m为采样点个数，一般为48；则P＝[P₁,...,P_k,...,P_N]^T∈R^N×m为N条日负荷曲线的m点原始有功功率矩阵；

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

p'_ki＝p_ki/p_k·max (3)

得到归一化的日负荷曲线有功功率标幺值矩阵P'_k＝[p'_k1,p'_k2,...,p'_ki,...,p'_km]∈R^1×m，并令该矩阵为A∈R^N×m。

所述步骤1-3)中对标幺化处理后的日负荷曲线数据进行PLR降维的方法具体为：

采用PLR对标幺化后的数据集中的每条曲线和聚类中心曲线进行自适应降维处理，该过程以采样点的斜率变化率作为参量进行分段线性化降维，降维后的曲线只保留能反映曲线特征的特征点，因此此过程产生的降维数据集的曲线维度不尽相同，步骤如下：

a)对于数据集中第k条维度为m的标幺化日负荷曲线P'_k＝[p'_k1,...,p'_ki,p'_kj...,p'_km]，首先计算此曲线的变化序列

即曲线上每个采样点对应采样时段的变化量(即斜率)集，其元素表达式如下：

b)对于曲线上的第i时刻对应的采样点，通过计算该采样点关于其左右两侧相邻采样点变化量的差值绝对值(即计算该采样点的左斜率与右斜率差值的绝对值)，以组成曲线的SEEP序列S＝(Δ1,...,Δi,...,Δm)，由于首尾两端元素分别无左斜率与右斜率，所以规定

且中间元素Δi表达式如下：

c)设定变化量阈值R，通过比较SEEP序列元素中的元素与阈值R的大小，来判断对应采样点是否为特征点，例如，对于第i时刻的采样点，若Δi＞R，则认为第i时刻对应采样点即为特征点。需要注意的是，阈值R代表着对原始数据集的压缩程度，一般其取值越大，则压缩程度也越大，但对原始数据的过于压缩会导致原始数据的关键信息丢失，所以R的取值至关重要。通常在实际工程试验中，对于标幺化的数据集，该降维程序以0.05为阈值参考值。

d)最后，将通过PLR算法所得的第k条曲线的所有特征点，按照采样点的先后顺序依次连接，即得到降维后的曲线h_k＝(p'_k1,...,p'_ki,...,p'_kf)，其中f为此第k条曲线降维之后的维度,如图2所示；降维后的曲线的集合为H＝(h₁,...,h_i,...,h_m)。

2)结合步骤1)所得标幺化数据矩阵A，以日负荷曲线数据组的密度函数为依据选取初始聚类中心C_L＝(c₁,c₂,...,c_L)；

依据数据对象的密度参数集合D，把D中密度最大值作为初始聚类中心；

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

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

根据式(4)计算任意两负荷曲线P_x＝(p_x1,p_x2,...,p_xn)，P_y＝(p_y1,p_y2,...,p_yn)之间的欧氏距离d，再通过式(5)得到日负荷曲线之间的平均距离dist，根据式(6)计算以每条日负荷曲线为中心基于平均距离的密度参数den，得到日负荷曲线的密度参数集合D，根据当前聚类数目L选取L个密度最大者，且每次选取一个密度最大者之后，将该曲线和与该曲线欧式距离小于dist的负荷曲线从数据库A中删除，从而得出日负荷曲线初始分布特性聚类中心曲线集合C_L＝(c₁,c₂,...,c_L)。

式(5)中，

为n条日负荷曲线之间的组合数。

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

3)对聚类中心进行PLR降维处理，然后对每条降维日负荷曲线，以DTW算法为衡量依据，计算其关于各类聚类中心曲线的相似度，并以此相似度数值为聚类算法K-means的划分依据，对标幺化日负荷曲线进行聚类划分之后，对各类聚类中心曲线进行更新；

所述步骤3)中采用DTW算法衡量降维日负荷曲线与降维聚类中心曲线的相似度方法步骤为：

3-1)先根据式(8)计算每条日负荷曲线关于聚类中心曲线的相似度D₁，D₁越小表示相似度越高，最后将日负荷曲线归类于其关于所有聚类中心曲线中D₁最小的那个类别；

D₁(i,j)＝dtw(h_i,c_j) (7)

其中D₁(i,j)表示，第i条降维后的日负荷曲线h_i关于第j个聚类中心c_j以dtw为衡量依据的相似度；

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

式中C_j(I+1)表示第I+1次迭代所得的第j类聚类中心曲线；N_j表示属于第j类聚类中心的日负荷曲线条数；P_n表示属于第j类聚类中心的日负荷曲线。

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

动态时间弯曲(DTW)运用动态规划思想调整时间序列不同时间点对应元素之间的关系来获取一条最优弯曲路径，使沿该路径时间序列间的距离最小，该算法的最大优势在于其能衡量不同维度间的时间序列的相似度。DTW路径示意图如图3所示。设有两负荷曲线P_x＝(p_x1,p_x2,...,p_xn)，P_y＝(p_y1,p_y2,...,p_ym)，DTW算法步骤如下：

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

b)将矩阵D中每一组相邻元素组成的集合称为弯曲路径，记为K＝(K₁,K₂,K₃,..K_s..,K_g)，其中g为路径中元素的总个数，元素K_s为路径上第s个点的坐标，即K_s＝(x,y)；

c)上述的路径K有多条，需要通过DTW算法寻找1条最优弯曲路径，使得序列P_x和P_y的弯曲总代价最小，即：

式中，D(K_s)为弯曲路径的累积距离，通过动态规划方法来构造一个累积代价矩阵G来求解上式，即G_xy＝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)＝+∞，可知，时间序列P_x和P_y的动态时间弯曲距离为DTW(P_x,P_y)＝G(n,m)。

同时，弯曲路径还需满足如下的约束：1)所选路径必须从左下角出发，到右上角结束，即K₁＝(1,1)，K_k＝(n,m)；2)每个点必须和相邻的点匹配，即若K_s＝(i,j)则K_s+1＝(a,b)必须满足0≤a-i≤1，0≤b-j≤1。并且，为了避免路径在同一水平或垂直方向进行多次连续弯曲而导致过度弯曲的现象(即避免时间序列的一个点对应另一条时间序列过多的点)，在已有约束的基础上加上对连续弯曲数的约束，即：

r_x≤r_x-max，r_y≤r_y-max

其中r_x、r_y分别为路径在x轴与y轴上的连续弯曲数；r_x-max，r_y-max分别为在x轴与y轴上所允许的最大连续弯曲数，其值由序列的特征及维数确定。

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

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

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

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

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

本实例首先从SQL2012数据库里随机选取969条日负荷曲线，在对所选负荷曲线进行筛选之后进行归一化处理，并对归一化的日负荷曲线数据组进行基于PLR的降维处理，之后设定初始化聚类数目L_min；然后，基于密度函数从数据组中选取L条日负荷曲线作为初始聚类中心曲线，并对其进行PLR降维处理；随后，利用DTW算法，以动态时间弯曲距离作为负荷曲线与聚类中心曲线的相似度衡量依据，聚类算法k-means以此相似度值为依据进行日负荷曲线集的划分，并得出相应聚类中心曲线以完成对聚类中心曲线的更新；之后，本算法以聚类中心曲线与聚类离散度的误差值，作为算法是否得出最终聚类结果的判断依据，若未达到允许误差值则继续进行运算，否则则完成此次以L为聚类数目的聚类算法运算；最后，判定此时聚类数目L是否达到规定的最大聚类数目L_max，若未满足则L＝L+1,继续进行聚类算法的运算，否则，跳出循环，以DBI指标最小对应的最佳聚类数及聚类中心曲线作为算法的最终聚类结果。

Claims (10)

1.一种基于PLR-DTW的日负荷曲线聚类方法，其特征在于，包括以下步骤：

步骤1)采集日负荷曲线数据以形成日负荷曲线数据组，对日负荷曲线数据进行标幺化处理，得到标幺化数据矩阵，并确定初始聚类数目、迭代次数和最大聚类数目；

步骤2)在标幺化数据矩阵中日负荷曲线的密度参数集合内，选取与当前聚类数目相等且密度最大的这几个值作为聚类中心曲线集，并对标幺化数据矩阵进行PLR降维处理；

步骤3)对聚类中心曲线集进行PLR降维处理，以当前聚类数目，若为第一次循环计算则以初始聚类数目开始进行循环计算，采用DTW算法分别计算降维后的标幺化数据矩阵中每条日负荷曲线关于聚类中心曲线的相似度，再以此相似度数值作为聚类算法k-means的相似度衡量依据，进行降维后的标幺化数据组的归类，并根据归类对聚类中心曲线集进行更新；

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

步骤5)结合迭代完成后所得不同聚类数目下的聚类结果，基于DBI确定最佳聚类数目与最终聚类结果，其中DBI越小则聚类结果越好，即选取DBI数值最小时所对应的聚类数目和聚类中心曲线集作为最佳聚类数目与聚类结果。

2.根据权利要求1所述的一种基于PLR-DTW的日负荷曲线聚类方法，其特征在于，所述步骤1)中，在进行标幺化处理前，还包括对日负荷曲线数据组中的异常数据进行识别与修正的步骤，其中对日负荷曲线中的异常数据进行识别，即计算日负荷曲线中各采样时刻点的负荷功率的变化率并与预设阈值比较，若超出阈值范围则视为异常数据；对异常数据进行进行修正，是首先判断负荷曲线的数据缺失量和异常量是否小于10％，若是，则先将异常量置为0即置为缺失量，再以一元三点抛物线插值算法对缺失量进行插值拟合，否则直接删除该条负荷曲线。

3.根据权利要求1所述的一种基于PLR-DTW的日负荷曲线聚类方法，其特征在于，所述步骤1)中，对修正后的日负荷曲线数据进行标幺化处理，得到标幺化数据矩阵的方法为：

记P_k＝[p_k1,...,p_ki,...,p_km]∈R^1×m为修正后第k条日负荷曲线的m点原始有功功率矩阵，k＝1,2,3,…,N，N为日负荷曲线总条数，p_ki为第k条日负荷曲线的第i点原始有功功率，i＝1,2,…,m，m为采样点个数；则P＝[P₁,...,P_k,...,P_N]^T∈R^N×m为N条日负荷曲线的m点原始有功功率矩阵，其中T为矩阵转置符号；

取日负荷曲线的功率最大值p_k.max＝max{p_k1,p_k2,...,p_ki,...,p_km}为基准值，根据下式对原始数据样本进行标幺化处理，得到功率点p_ki的标幺值p'_ki：

p'_ki＝p_ki/p_k.max

对各功率点计算标幺值，进而得到归一化的日负荷曲线有功功率标幺值矩阵P'_k：

P'_k＝[p'_k1,p'_k2,...,p'_ki,...,p'_km]∈R^1×m，则标幺化后的日负荷曲线矩阵A＝[P₁',...,P_k',...,P_N']^T∈R^N×m。

4.根据权利要求1所述的一种基于PLR-DTW的日负荷曲线聚类方法，其特征在于，所述步骤1)中，最大聚类数目

N为日负荷曲线数据组中的总曲线数目。

5.根据权利要求1所述的一种基于PLR-DTW的日负荷曲线聚类方法，其特征在于，所述的步骤2)中得到聚类中心曲线集包括以下过程：

基于标幺化数据矩阵，首先计算其中任意两负荷曲线P_x＝(p_x1,p_x2,...,p_xn)，P_y＝(p_y1,p_y2,...,p_yn)之间的欧氏距离d：

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

其中

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

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

其中P_i为除P_c外的任意一条日负荷曲线，当dist-d(P_i,P_c)＜0时，u(dist-d(P_i,P_c))＝0；当dist-d(P_i,P_c)≥0时，u(dist-d(P_i,P_c))＝1；

进而得到日负荷曲线的密度参数集合D，根据当前聚类数目L选取前L个密度最大的日负荷曲线，且每选取一个后，将该日负荷曲线和与该日负荷曲线欧式距离小于dist的日负荷曲线从标幺化数据矩阵A中删除，从而得出日负荷曲线的初始聚类中心曲线集合。

6.根据权利要求1所述的一种基于PLR-DTW的日负荷曲线聚类方法，其特征在于，所述的步骤2)和步骤3)中，进行PLR降维处理包括以下过程：

对于待降维的标幺化数据矩阵或聚类中心曲线集中第k条维度为m的标幺化日负荷曲线P'_k＝[p'_k1,...,p'_ki,p'_kj...,p'_km]，首先计算此日负荷曲线的变化序列

即日负荷曲线上每个采样点对应采样时段的变化量也即斜率集，其元素表达式如下：

然后，对于日负荷曲线上的第i时刻对应的采样点，通过计算该采样点关于其左右两侧相邻采样点变化量的差值绝对值，即计算该采样点的左斜率与右斜率差值的绝对值，以组成曲线P'_k的SEEP序列S＝(Δ1,...,Δi,...,Δm)，由于首尾两端元素分别无左斜率与右斜率，所以规定

且中间元素Δi表达式如下：

之后，设定变化量阈值R，并比较SEEP序列中的元素是否大于阈值R，当大于时则对应采样点为特征点；

最后，将通过PLR算法所得的第k条日负荷曲线的所有特征点，按照采样点的先后顺序依次连接，即得到降维后的日负荷曲线h_k＝(p'_k1,...,p'_ki,...,p'_kf)，其中f为此第k条日负荷曲线降维之后的维度；降维后的日负荷曲线的集合为H＝(h₁,...,h_i,...,h_m)。

7.根据权利要求6所述的一种基于PLR-DTW的日负荷曲线聚类方法，其特征在于，所述的步骤3)中，采用DTW来计算降维的标幺化数据矩阵中每条日负荷曲线关于聚类中心曲线的相似度，其过程为：

首先，对所得聚类中心曲线进行PLR降维处理，得到降维聚类中心曲线集合C_L＝(c₁,c₂,...,c_L)；

然后，计算每条降维后的日负荷曲线与聚类中心曲线相似度D₁：

D₁(i,j)＝dtw(h_i,c_j)

其中D₁(i,j)表示，第i条降维后的日负荷曲线h_i关于第j个聚类中心c_j以DTW为衡量依据的相似度值；

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

设有两负荷曲线P_x＝(p_x1,p_x2,...,p_xn)，P_y＝(p_y1,p_y2,...,p_ym)，则有：

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

b)将矩阵D中每一组相邻元素组成的集合称为弯曲路径，记为K＝(K₁,K₂,K₃,..K_s..,K_g)，其中g为路径中元素的总个数，元素K_s为路径上第s个点的坐标，即K_s＝(x,y)；

c)通过DTW算法寻找弯曲路径K中的最优弯曲路径，使得序列P_x和P_y的弯曲总代价最小，即：

式中，D(K_s)为弯曲路径的累积距离，通过动态规划方法来构造一个累积代价矩阵G来求解上式，即G_xy＝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)＝+∞，则时间序列P_x和P_y的动态时间弯曲距离为DTW(P_x,P_y)＝G(n,m)。

同时，弯曲路径还需满足如下的约束：1)所选路径必须从左下角出发，到右上角结束，即K₁＝(1,1)，K_k＝(n,m)；2)每个点必须和相邻的点匹配，即若K_s＝(i,j)则K_s+1＝(a,b)必须满足0≤a-i≤1，0≤b-j≤1，同时为了避免路径在同一水平或垂直方向进行多次连续弯曲而导致过度弯曲的现象，即避免时间序列的一个点对应另一条时间序列过多的点，在已有约束的基础上加上对连续弯曲数的约束，即：

r_x≤r_x-max，r_y≤r_y-max

其中r_x、r_y分别为路径在x轴与y轴上的连续弯曲数；r_x-max，r_y-max分别为在x轴与y轴上所允许的最大连续弯曲数，其值由序列的特征及维数确定。

8.根据权利要求1所述的一种基于PLR-DTW的日负荷曲线聚类方法，其特征在于，所述的步骤3)中，以相似度数值作为聚类算法k-means的相似度衡量依据，进行数据组的归类，并根据归类对聚类中心曲线集进行更新的步骤包括：

将每条标幺化曲线P'_k归类于与其最相似的聚类中心曲线的一类中，之后，对数据组的聚类中心曲线，按照下式进行更新：

式中C_j(I+1)表示第I+1次迭代所得的聚类中心曲线；N_j表示属于第j类聚类中心的日负荷曲线条数；P_n表示属于第j类聚类中心的日负荷曲线。

9.根据权利要求1所述的一种基于PLR-DTW的日负荷曲线聚类方法，其特征在于，所述的步骤4)中，综合聚类离散度通过以下公式计算：

式中，J_N为第N次迭代所得聚类结果的综合聚类离散度；C_i(N)为第N次迭代所得的第i个聚类中心；P_n为属于聚类中心曲线C_i(N)这一类的日负荷曲线；L为聚类数目；

若J_N+1-J_N＜Ω，则认为综合聚类离散度的误差值满足要求，其中Ω为预设的收敛极限值。

10.根据权利要求1所述的一种基于PLR-DTW的日负荷曲线聚类方法，其特征在于，所述的步骤5)中以基于DBI指标确定最佳聚类数与最终聚类结果的的方法如下：

式中，R_i用来衡量第i类曲线集的紧密程度；L为聚类数目；

其中，R_ij用来衡量第i类与第j类的相似度，S_i用来度量第i个类中数据点的分散程度，通过下式计算：

其中，X_l为第i类中第l个数据点；C_i为第i类的聚类中心曲线；T为第i类中数据点的个数；q取1时S_i为各点到中心的距离的均值，q取2时S_i为各点到中心距离的标准差，均用于衡量类内分散程度；S_j用来度量第j个类中数据点的分散程度，计算方式与S_i相同；

其中M_ij为第i类中心与第j类中心的距离；C_i为第i类的聚类中心曲线；C_j为第j类的聚类中心曲线；a_di为第i类的中心点C_i的第d个属性的值；p取1时表示1-范数，p取2时表示2-范数，即两个类中心的欧式距离。

Images