宋軍英，崔益?zhèn)ィ钚廊?，?偉，劉韜文，李培強

（1. 國網湖南省電力有限公司,湖南省長沙市410077；2. 湖南大學電氣與信息工程學院,湖南省長沙市410082）

0 引言

變電站綜合負荷是由數(shù)量眾多、特性各異的電力用戶構成的,通過在線掌握用戶的負荷特性,可實現(xiàn)變電站綜合負荷構成特性的在線解析。其基本思路是：運用聚類分析的方法,對用戶的日負荷曲線進行分類與綜合,獲得典型用電行業(yè)的分類負荷特性,進而在線解析得到綜合負荷的用電行業(yè)構成比例,最終實現(xiàn)在線負荷建模［1-2］。

隨著智能電網建設的深入,用戶日負荷曲線數(shù)據(jù)量及維度大幅提高,對電網大數(shù)據(jù)平臺負荷特性提取的準確性、魯棒性及運算效率提出了更高的要求［3-4］。負荷曲線采樣頻率的提升雖然能更完整地反映用戶的日用電特性,但高維數(shù)據(jù)集應用于聚類算法時會導致運算效率下降,且由于采樣點增多,某一時間段的小幅噪聲（以某一采樣值水平線為基準上下波動的采樣點）會影響曲線動態(tài)特性的刻畫,從而導致相似度衡量誤差。此時如果依然采用傳統(tǒng)算法（如K 均值（K-means）算法）,直接應用原始數(shù)據(jù)進行聚類,就會存在以下3 個缺點［5］：①聚類數(shù)目需事先劃定,初始聚類中心曲線的選取過程完全隨機；②相似度衡量方法難以準確估量高維曲線的動態(tài)特性；③魯棒性較差,擾動點對算法聚類質量負面影響較大。因此,對于基于電網大數(shù)據(jù)平臺的在線負荷建模,如何從海量的日負荷曲線中提取有價值的信息,高質量、高效率地進行實時準確聚類,完成變電站綜合負荷解析,就成為亟須解決的重要問題［6］。

解決上述問題的有效方法之一就是對用戶日負荷曲線集進行降維處理,提取能準確表征用戶用電特性的特征點或特征指標進行聚類。文獻［7］采用峰谷期負載率、最大最小負荷對應時刻等特征指標對日負荷曲線降維；文獻［8］通過對日負荷曲線進行奇異值分解,以提取負荷曲線的負荷特征；文獻［9］提出一種基于離散小波變換（discrete wavelet transformation,DWT）的模糊聚類方法；文獻［10］基于子空間聚類算法進行負荷曲線集的特征提取,嘗試在相同數(shù)據(jù)集的不同子空間上發(fā)現(xiàn)聚類?？v觀上述文獻,大多采用等分辨率降維方法,強制使數(shù)據(jù)集統(tǒng)一降至某一維度,一定程度上忽略了某些負荷曲線的斜率、極值點等關鍵負荷特性［11］。同時,以歐氏距離作為相似度衡量方法,無法準確衡量負荷曲線的動態(tài)特性,也不適用于自適應降維所構成的不等維時間序列集的相似度衡量［12-13］。

文獻［13］研究表明,以動態(tài)時間彎曲（dynamic time warping,DTW）距離代替歐氏距離作為時間序列相似度衡量指標,能更充分地反映序列的整體動態(tài)特性,但在效率上存在一定的劣勢,因此對于高維度曲線,難以滿足在線、實時的應用要求。

本文在文獻［13］研究的基礎上,提出一種改進分 段 線 性 表 示（improved piecewise linear representation,IPLR）與DTW 距離相結合的基于Canopy 的K-means（CK-means）日負荷曲線聚類方法。首先,該方法以相鄰及間隔采樣點變化量為依據(jù),對原始曲線集進行基于IPLR 方法的自適應重構,得到一組不等維度的降維數(shù)據(jù)組；然后,采用Canopy 算法獲取聚類數(shù)目及初始聚類中心；最后,利用以DTW 距離作為相似度衡量手段的K-means算法對其進行聚類處理。算例結果表明,本文方法所采取的IPLR 自適應降維方法與DTW 距離相似度衡量手段相契合,所得聚類結果與實際相符,且在聚類質量、魯棒性及運算效率上較傳統(tǒng)方法均具有一定的優(yōu)越性,滿足基于電網大數(shù)據(jù)平臺的實時在線負荷建模的要求。

1 IPLR 算法

1.1 基本思路

傳 統(tǒng) 分 段 線 性 表 示（piecewise linear representation,PLR）算法基于相鄰采樣點的采樣值變化量對原始時間序列進行重構。首先,計算每個采樣點的相鄰采樣值變化量；然后,根據(jù)變化量是否超過給定閾值,以判斷其是否為特征點；最后,將所有特征點按先后順序依次連接,即完成時間序列的降維重構［14］。但是該算法的特征點提取條件過于寬松,僅關注時間序列的局部特性,受擾動影響較大［15-16］。當時間序列中某一時間段的采樣點維持某一水平采樣值進行近似的等幅小額振動,且此一系列頻繁振動導致此時間段的采樣點相鄰變化量皆超過閾值,則降維重構序列中將包含大量擾動點。為解決這一問題,本文在傳統(tǒng)算法的基礎上加入對采樣點間隔變化量的限制,以達到通過提取少量采樣點即可反映序列的關鍵特性,并增強算法的抗干擾能力的目的。算法流程如圖1 所示。

圖1 IPLR 算法流程圖Fig.1 Flow chart of IPLR algorithm

1.2 算法實現(xiàn)

對于具有m條時間序列的n維數(shù)據(jù)集G=[G1,G2,…,Gm],其中第i條時間序列為Gi=[gi1,gi2,…,gin]。通過IPLR 算法對Gi進行自適應降維得到u維的降維序列Hi=[hi1,hi2,…,hiu]的具體步驟如下。

步驟1：首先,取原始序列Gi的首尾兩端的采樣點作為降維序列的首尾兩端,即hi1=gi1,hiu=gin；然后,輸入相鄰采樣點變化量閾值參數(shù)λ1和λ2［17］及間隔采樣點變化量閾值參數(shù)σ和ε。其中σ通過設定間隔采樣點變化量差閾值,選取間隔變化量絕對值變化較大的特征點；ε通過設定采樣點的間隔變化量乘積閾值,選取左右間隔變化量絕對值相近但值較大的特征點。λ1,λ2,σ和ε的取值,由數(shù)據(jù)集G中隨機選取m條時間序列集的每條時間序列的相鄰采樣點的平均變化量與間隔采樣點的平均變化量決定,如式（1）至式（3）所示。

步驟2：計算第j個采樣點的相鄰變化量k1=xj+1-xj,k2=xj-xj-1,以及間隔變化量k3=xj+2-xj,k4=xj+1-xj-1,其中2 ≤j≤n-2。

步驟3：對于時間序列點gin,若|k1+k2|＞λ1或|k1-k2|＞λ2,則進入步驟4,否則j=j+1,進入步驟5。

步驟4：若||k3|-|k4||＞σ或|k3k4|＞ε,則提取該點,進入步驟5,否則j=j+1,進入步驟5。

步驟5：若j=n-2,輸出降維重構曲線,否則進入步驟2。

1.3 降維效果分析

現(xiàn)對一條維度為96 的日負荷曲線A進行降維分析,該曲線在0～24 和48～65 采樣時段分別維持0.5 和0.85 的負荷水平進行小額波動。若應用基于SEEP 序列的PLR 算法,曲線維度降至49,得到降維曲線A',如圖2 所示。不難看到,該方法將負荷低谷期及雙峰間谷端的處于波動時間段的采樣點采納為特征點,導致重構的降維曲線維數(shù)較高,且含有大量擾動點,無法準確反映原始曲線的動態(tài)特性。若應用本文的IPLR 算法對日負荷曲線進行降維重構,得到降維曲線A″,如圖3 所示,維度降至39,大部分波動時間段的擾動點被剔除。

圖2 PLR 降維示意圖Fig.2 Schematic diagram of dimension reduction of PLR

圖3 IPLR 降維示意圖Fig.3 Schematic diagram of dimension reduction of IPLR

應用文獻［14］中的時間序列壓縮度計算方法,對2 種降維方法關于原始日負荷曲線A的壓縮率進行計算對比分析,可以發(fā)現(xiàn)IPLR 降維方法的壓縮率較PLR 降維方法高25.64%。定義Ddtw(X,Y)為時間序列X與Y之間基于DTW 距離的相似度。將原始曲線A分別與降維曲線A'與A″進行擬合度分析,則 可 得Ddtw(A,A')=0.024 8,Ddtw(A,A″)=0.017 8。不難發(fā)現(xiàn),IPLR 降維方法所得降維曲線與原始曲線擬合程度更高,更能反映負荷曲線的負荷特性。

現(xiàn)從數(shù)據(jù)集中隨機選取負荷曲線B（B'與B″分別為曲線B經PLR 算法與IPLR 算法所得的降維曲線）與負荷曲線A作相似度分析,以DTW 距離為相似度衡量指標［18］,結果如表1 所示。其中,Ⅰ,Ⅱ,Ⅲ分別對應不降維處理、PLR 降維處理及IPLR 降維處理。

表1 相似度衡量方法對比Table 1 Comparison of similarity measuring methods

由表1 可知,通過對負荷曲線降維,可以顯著提升DTW 距離應用于負荷曲線間相似度衡量的運算效率,使運算時間縮減至接近原來的1/4。同時,IPLR 降維方法與原始曲線的相似度指標值及運算時間分別比PLR 降維方法降低近11%和10%。

綜上可知,本文所提IPLR 降維方法通過提取負荷曲線的特征點對原始曲線進行自適應降維,能以盡可能低維度的降維曲線反映原始負荷曲線的關鍵動態(tài)特性。該方法所得降維曲線擬合度高、運算時間短、抗干擾能力強,適用于用戶日負荷曲線的聚類運算。

2 基于IPLR 與DTW 距離的CK-means算法

2.1 CK-means 算法

CK-means 算法為基于Canopy 聚類結果的Kmeans 聚類算法,相比于K-means 算法,其首先應用Canopy 算法選取聚類數(shù)目及初始聚類中心,然后再采用K-means 算法進行聚類。

Canopy 算法是一種快速簡單但精準度一般的聚類算法,其最大特點是無須確定聚類數(shù)目,適用于預處理階段對原始數(shù)據(jù)集進行粗聚類處理,將此聚類結果作為初始聚類中心,雖同樣具有一定的隨機性,但代表性更強。算法基本運算流程如下［19］。

步驟1：設定閾值T1與T2（一般T2取所有曲線平均距離的2 倍,且T1=2T2）,并滿足T1＞T2。

步驟2：從數(shù)據(jù)集中任取一點,作為第1 個Canopy（作為Canopy 的數(shù)據(jù)點應從數(shù)據(jù)集中刪去）。

步驟3：計算數(shù)據(jù)集中其他每個點Zi與所有Canopy（Rj）的距離Ui-j。若Ui-j＜T2,則將其歸入此Rj類別；若其關于所有Canopy 的距離滿足Ui-j＞T1,則將其當作一個新的Canopy；若該點到某個Canopy 距離Ui-j＜T1,并在其與所有Canopy距離計算完成后依然未加入任何Canopy,則將其作為一個新的Canopy。

步驟4：重復步驟3,直至數(shù)據(jù)集為空。

2.2 聚類有效性指標

聚類的有效性一般通過以下2 個方面反映：一是同一種類的對象間相似度較高；二是不同種類的對象間差異性較大。

由文獻［20］可知,DBI（Davies-Bouldin index）指標同時考慮聚類結果的類內及類間的聚類效果,如式（4）所示。IDBI為類內距離之和與類外距離的比值,其計算公式簡單且指標值變化范圍小,能直觀反映聚類質量。因此,相比于其他指標,DBI 指標更適用于評定電力用戶日負荷曲線聚類的有效性。

式中：K為聚類個數(shù)；Si為第i個類中曲線集與該類聚類中心曲線的平均距離,用于反映第i個類中曲線集的聚攏程度；Mij為第i類中心與第j類聚類中心曲線的距離,以反映類間第i類與第j類曲線集的分散程度。

由于不同相似度衡量手段側重于衡量曲線的不同特性,故本文以歐氏距離作為相似度的指標定義為I1,以DTW 距離作為相似度的指標定義為I2,分別以負荷曲線的整體分布特性與整體動態(tài)特性評估聚類質量［13］。

2.3 算法實現(xiàn)過程

首先,本文算法基于IPLR 對原始數(shù)據(jù)集進行降維處理；然后,應用Canopy 算法確定算法的聚類數(shù)目及初始聚類中心；最后,應用基于DTW 距離的K-means 算法對降維數(shù)據(jù)集進行聚類運算。具體步驟如下。

步驟1：預處理。首先,對于部分采樣點缺失的日負荷曲線進行插值補全,并基于曲線相鄰采樣點波動量篩除異常波動曲線［21］；然后,對日負荷曲線集采取極值歸一化處理［22］,得到含有m條曲線的歸一化數(shù)據(jù)集Y=[y1,y2,…,ym]T。

步驟2：初始化。對數(shù)據(jù)集Y進行IPLR 降維處理,得到不等維數(shù)據(jù)集P=[p1,p2,…,pm]T,并對原始數(shù)據(jù)組進行基于Canopy 算法的粗聚類運算,得到聚類數(shù)目L與初始聚類中心曲線集C=[C1,C2,…,CL]T。

步驟3：相似度的衡量。首先,對不等維數(shù)據(jù)組P=[p1,p2,…,pm]T中的曲線關于每類聚類中心曲線進行基于DTW 距離的相似度計算,之后,將每條負荷曲線分至與其最相似（即DTW 距離值最小）的類別中［23］。

步驟4：對聚類中心矩陣C=[C1,C2,…,CL]T進行更新,如式（5）所示。當算法進行至滿足以下任一條件時,結束運算,本次循環(huán)的聚類中心曲線即為最終聚類結果：①2 次迭代所得成本損失函數(shù)Serror滿 足|Serror(I+1)-Serror(I)| ＜e（I表 示 算 法 迭 代次數(shù),e為收斂閾值）,成本損失函數(shù)如式（6）所示［24］；②2 次迭代所得聚類中心曲線的誤差在閾值范圍內［25］。否則,繼續(xù)進行步驟3 與4,直至完成規(guī)定最大迭代次數(shù)或滿足以上條件之一。

式中：CL(I)為算法進行第I次迭代后的第L類聚類中心曲線；Ddtw(Yn,CL(I))為Yn與CL(I)之間基于DTW 距離的相似度。

3 算例檢驗

3.1 實例檢驗與比較

本文隨機選取某省區(qū)電網110 kV 變電站下屬1 200 個典型用戶（包含工商業(yè)及居民用戶）的不同采樣頻率（采樣頻率分別為15 min/點和30 min/點）的一天的日負荷曲線數(shù)據(jù)集作為實驗對象。為了驗證本文所提降維方法與所采用相似度衡量手段的合理性及優(yōu)越性,本文對基于歐氏距離的CK-means算法（方法1）、基于DTW 距離的CK-means 算法［13］（方法2）、基于PLR 與DTW 距離的CK-means 算法（方法3）和基于IPLR 與DTW 距離的CK-means 算法（本文方法）展開聚類分析,并對這4 種方法的聚類結果（聚類中心曲線）、聚類質量（DBI 指標）及聚類效率（運算時間及迭代次數(shù)）進行綜合比較。實驗所用機器：單臺計算機,配置為i5-4570s CPU@2.90 GHz,1050 Ti 4 GB,操作系統(tǒng)為Windows 7,內存為16 GB。

對48 點日負荷曲線數(shù)據(jù)集展開聚類分析。在預處理步驟中,本算例認為若曲線相鄰負荷點變化超過20%,即為異常波動曲線,需排除。本算例降維算法閾值參數(shù)取值為：λ1=0.055,λ2=0.035,σ=0.025,ε=0.000 45。

4 種算法所提取的不同行業(yè)類別的聚類結果如附錄A 圖A1 至圖A4 所示,可以發(fā)現(xiàn),每一行業(yè)的聚類中心曲線相似度較高。其中降維算法（方法3和本文方法）與方法1 和方法2 的分歧點主要集中在第3 類負荷的劃分上。附錄A 圖A1 與圖A2 中第3 類負荷特征體現(xiàn)為白天負荷水平較高,負荷高峰期主要集中于11：30—12：30 與17：00—20：00,且晚高峰負荷遠大于午間負荷；圖A3 與圖A4 的第3 類負荷特征體現(xiàn)為三峰負荷曲線,3 段峰值集中在07：30—09：00,11：30—12：30 與17：00—20：00,同樣,晚高峰負荷遠大于前2 段峰值負荷,這一類負荷的負荷特征屬于市政生活用電類型。定義Ic為類內距離指標,以衡量各算法第3 類聚類中心曲線的聚類質量,公式如式（7）所示,結果如附錄A 圖A5所示,可知本文方法的類內聚攏效果最佳。

式 中：m0為 屬 于 第3 類 的 負 荷 曲 線 數(shù) 目；xi為 第j條曲線第i個采樣點的采樣值；ci為第3 類聚類中心曲線第i個采樣時刻的對應值。

分析另外3 類聚類中心曲線：第1 類負荷曲線白天除午間出現(xiàn)相對低谷期,整體變化較平緩,晚間負荷水平較低,反映的是采用單班制的輕工業(yè)企業(yè)的用電行為；第2 類負荷曲線全天基本保持在一個高負荷水平,反映的是以采礦、煉鋼行業(yè)為代表,一般采用三班制作業(yè)的重工業(yè)企業(yè)的用電行為；第4 類負荷曲線從08：00—22：00 變化較平緩,反映的是以服務業(yè)為代表的第三產業(yè)用電行為。因此,這4 種算法聚類結果均與實際相符,具有一定的實際工程參考價值。

附錄A 圖A6 所示為4 種算法的聚類指標隨迭代次數(shù)變化的曲線；表2 所示為4 種算法對日負荷曲線數(shù)據(jù)進行聚類分析的耗時、迭代次數(shù)及最終聚類指標的性能對比。由附錄A 圖A6 與表2 可知,通過對原始負荷曲線集進行降維重構,本文方法的算法性能得到顯著提升,相比于方法2,聚類指標優(yōu)化明顯,運算時間降低近50%；相比于方法1,雖然犧牲了一定的運算效率,但也因此獲得了更高的聚類質量（I1指標下降16.51%、I2指標下降53.02%）；而相比于同樣對原始負荷曲線集進行了降維重構的方法3,本文方法由于采用了更為嚴格的特征點篩選條件,在運算效率及聚類質量上都得到了進一步提升。因此,本文方法相比于其他3 種方法,具有最優(yōu)的綜合性能。

表2 不同算法的48 點曲線集聚類指標Table 2 Clustering indices of 48-point curves of different algorithms

對96 點日負荷曲線數(shù)據(jù)集展開聚類分析。在預處理步驟中,本算例認為若曲線相鄰負荷點變化超過15%,即為異常波動曲線,需排除。本算例中,降維算法閾值參數(shù)取值為：λ1=0.05,λ2=0.03,σ=0.025,ε=0.000 45。

在該工況下,4 種算法的聚類結果與48 點負荷曲線數(shù)據(jù)集的聚類結果基本一致,限于篇幅,此處不再展示。聚類指標如表3 所示。由表可知,本文方法的綜合性能相比于其他3 種方法依然為最優(yōu),但其相比于對48 點負荷曲線數(shù)據(jù)集進行聚類,算法運算效率出現(xiàn)顯著下滑,耗時增長170.89%。這是因為對于一個維度為n的負荷曲線集,DTW 距離算法復雜度為o（n2）,當數(shù)據(jù)集的維度翻倍時,運算效率將大幅下降。

表3 不同算法的96 點曲線集聚類指標Table 3 Clustering indices of 96-point curves of different algorithms

綜上可知,本文算法得到的聚類結果符合工程實際,且相比于其他方法綜合性能更優(yōu),但受限于算法運算效率,所以更適用于48 點負荷曲線數(shù)據(jù)集的聚類運算。

3.2 算法魯棒性分析

為檢驗本文方法的魯棒性,對3.1 節(jié)所選曲線集加入大小為r（r=5%,10%,15%,20%,25%）的隨機擾動,以模擬實際用戶負荷曲線采樣過程中因天氣等隨機因素造成的負荷波動。然后,本節(jié)分別采用3.1 節(jié)中的4 種方法對擾動曲線集進行聚類對比分析,并以聚類質量指標DBI 作為魯棒性考量指標。

由表4 可知,隨著對負荷曲線所加擾動的增加,各類方法的聚類質量指標基本呈下降趨勢。方法1與方法2 在小擾動干擾下,聚類質量指標尚可,但當擾動r≥10%時,聚類質量指標出現(xiàn)大幅下降,且方法2 的I2指標下降速度要低于方法1。這是因為,方法1 與方法2 直接應用原始負荷曲線數(shù)據(jù)集進行聚類,導致每個采樣點的采樣值都會對聚類結果產生直接影響,從而使其無法準確提取曲線的動態(tài)特性,產生較大的聚類偏差,但方法2 由于采取DTW 距離作為相似度衡量手段,所以在動態(tài)特性指標I2上表現(xiàn)出更強的魯棒性。方法3 和本文方法由于對原始曲線進行了特征點提取,并以DTW 距離作為相似度衡量方法,綜合考慮負荷曲線的動態(tài)特性,因此受擾動影響較前兩者小,故在中小擾動工況下（r≤20%）都能表現(xiàn)出較好的魯棒性。但在大擾動工況下,方法3 和本文方法魯棒性出現(xiàn)了較大幅度下滑,這主要是因為閾值參數(shù)是按照一般工況進行的設定,取的是日負荷曲線相關特性參數(shù)的平均值,當擾動過大時,部分擾動點將被誤提取為特征點。但是,本文方法在中小擾動及大擾動工況下的魯棒性,依然略優(yōu)于方法3,這主要得益于其更嚴格的特征點提取方法,進一步降低了擾動的負面影響。

表4 不同算法的聚類指標Table 4 Clustering indices of different algorithms

綜上可知,本文方法魯棒性較強,在一般工況（中小擾動工況）下依然能以較高質量提取用戶原始負荷曲線的負荷特性；但在大擾動的極端工況下,魯棒性會出現(xiàn)一定程度的下降?？傮w而言,本文方法能夠滿足基于大數(shù)據(jù)的在線負荷建模要求。

4 結語

為準確提取負荷曲線的特征點,并提升現(xiàn)有聚類算法的聚類質量,本文提出一種基于IPLR 降維與DTW 距離相似度衡量的CK-means 算法。首先,該算法基于IPLR,根據(jù)原始數(shù)據(jù)集中負荷曲線自身負荷特性進行自適應分辨率降維；然后,應用基于DTW 距離相似度衡量方法的CK-means 算法對此不等維降維數(shù)據(jù)組進行聚類運算分析,以準確提取不同用電特性用戶的用電特征。本文算例得到如下結論。

1）IPLR 降維方法采用合理。其可在一定程度上過濾負荷曲線的擾動采樣點,準確提取出負荷曲線的關鍵特征點。

2）DTW 距離相似度衡量手段采用合理。此相似度衡量方法能較準確地衡量負荷曲線間的動態(tài)特性相似度,且其可衡量不等維時間序列間相似度的特點與IPLR 降維方法可對數(shù)據(jù)組進行自適應降維的優(yōu)點相契合。

3）本文所提聚類方法較傳統(tǒng)方法具有更高的綜合性能,滿足基于電網大數(shù)據(jù)平臺的實時在線負荷建模的要求,對變電站綜合負荷構成比例解析具有重要參考價值。

但本文方法對數(shù)據(jù)預處理要求較高,且運算效率有進一步提升的空間。下一步工作將在現(xiàn)有工作基礎上,針對大擾動工況下魯棒性較低及DTW 距離算法計算耗時較長的問題,對本文算法做進一步研究改進。

本文得到湖南省電力公司重點計劃項目(5216A5180018)的資助，特此感謝！

附錄見本刊網絡版（http：//www.aeps-info.com/aeps/ch/index.aspx），掃英文摘要后二維碼可以閱讀網絡全文。