基于密度峰值聚類的高斯混合模型核電運行工況劃分

崔文浩,鄭勝,楊森權,楊珊珊,曾曙光,羅驍域

(1.三峽大學電氣與新能源學院,宜昌 443002;2.三峽大學理學院,宜昌 443002;3.中核核工業(yè)仿真技術重點實驗室,武漢 443074)

推進能源電力清潔低碳轉型與高質量發(fā)展是貫徹落實中國碳達峰、碳中和戰(zhàn)略目標[1-2]的核心要義。核電作為清潔能源因具有污染排放少、發(fā)電效率高等優(yōu)點,而擁有良好的發(fā)展前景。隨著核電數(shù)字化發(fā)展[3],核電廠中各系統(tǒng)的傳感器數(shù)據(jù)被采集和存儲,這些數(shù)據(jù)蘊含了系統(tǒng)設備的狀態(tài)信息,為工況劃分模型的建立提供了數(shù)據(jù)基礎。核電廠對各工況數(shù)據(jù)進行處理分析,便可實現(xiàn)對運行過程的在線監(jiān)控[4]、實時預警[5]和故障診斷[6],為核電運維人員的安全維護,經(jīng)濟運行提供有效的數(shù)據(jù)基礎和分析手段,因此核電的工況劃分是一項非常重要的工作。但是由于核電廠內部過程的復雜性,無法得知核電運行工況的準確個數(shù),因此需要采用聚類方法來確定工況數(shù)。

目前,有幾種合適的聚類算法可以被應用至工況劃分之上。Hallac等[7]提出了TICC(toeplitz inverse covariance-based clustering)算法,該算法通過考察信號各維度之間的相關性,利用時間窗口來完成對多元時序數(shù)據(jù)的聚類。早前有相關核電研究人員將該方法應用至工況劃分,但是由于時間窗的大小需要靠先驗知識確定,而核電數(shù)據(jù)較為復雜,難以通過人工的方式確定時間窗口的大小,使得該方法的適用性不強。K-means聚類算法[8]收斂速度快,聚類高效且效果好,從理論上看比較適合核電運行數(shù)據(jù)這種維度高體量大的大規(guī)模數(shù)據(jù)集。秦緒華等[9]通過改進K-means算法并應用至火電的外部工況劃分上,使工況劃分結果更合理。然而K-means聚類算法對初始聚類中心的依賴非常嚴重[10],對非凸數(shù)據(jù)集的聚類效果欠佳[11],且會出現(xiàn)局部最小值的情況[12],而核電廠的內部過程復雜,難以獲知準確的工況個數(shù)和聚類中心,這也使得將K-means方法應用于核電工況劃分時所得到的結果并不理想。

為了解決上述算法在應用于核電運行工況劃分時所出現(xiàn)的問題。提出了基于密度峰值聚類的高斯混合模型核電運行工況劃分算法。密度峰值聚類算法[13]的聚類速度快、能夠快速發(fā)現(xiàn)任意形狀的類簇[14]、魯棒性強,且不受聚類對象嵌入的空間維數(shù)的影響,已廣泛用于圖像識別[15]、天文數(shù)據(jù)處理[16]等領域。該算法通過引入決策圖的方式,使得操作人員可以直觀的確定聚類的個數(shù),該算法可以幫助確定核電運行工況的個數(shù)。高斯混合模型(Gaussian mixture model,GMM,)[17-18]能適用于線性與非線性數(shù)據(jù),在異常檢測[19]、工況劃分[20]等方面均得到了廣泛的應用。核電廠在穩(wěn)定運行狀態(tài)下的數(shù)據(jù)可認為是服從高斯分布的時序數(shù)據(jù),因此可利用高斯混合模型來完成工況劃分。但是核電廠數(shù)據(jù)維度高,直接導入模型中會產(chǎn)生計算難度大的問題,且高斯混合模型需要確定聚類初值,不合適的聚類初值會給聚類結果帶來很大的影響。

鑒于此,利用主成分分析(principal component analysis,PCA)算法[21]篩選出高維核電數(shù)據(jù)的主元,完成數(shù)據(jù)降維,降低高斯混合模型的計算難度,然后采用密度峰值聚類算法來確定工況個數(shù),為確定高斯混合模型的初值提供參考。最后,通過模型響應度剔除響應度低的子模型,提高混合高斯模型的準確性,實現(xiàn)準確的工況劃分。

1 算法設計

基于密度峰值聚類的高斯混合模型核電運行工況劃分算法流程圖如圖1所示,具體步驟如下。

圖1 算法流程Fig.1 Algorithmic flow

步驟1 對核電數(shù)據(jù)進行預處理,去除原始數(shù)據(jù)中的噪聲。

(1)

(2)

步驟3 采用密度峰值聚類,該算法以決策圖的方式將分類結果可視化。密度峰值聚類算法是由Rodriguze等[13]于2014年提出的聚類算法,其思想是將高密度區(qū)域從低密度區(qū)域的包圍中剝離出來,所以該算法也就滿足以下兩點假設:①聚類中心的密度大于周圍區(qū)域的密度;②聚類中心點與其他高密度點的距離較大。

由此可知,密度峰值聚類算法的核心在于如何計算局部密度ρi與聚類中心點的距離δi。可將ρi定義為

(3)

(4)

式中:dc為截斷距離;dij的含義為任意一點j與i點的距離;χ(x)為判斷點是否處于聚類圓之中的函數(shù)。

式(3)可理解為與數(shù)據(jù)點i的距離dij小于截斷距離dc的點的個數(shù)。對于聚類中心點的距離δi可定義為

(5)

利用該算法對PCA降維后的數(shù)據(jù)進行聚類,在計算得出局部密度ρi和距離δi之后,可通過如圖2所示的決策圖,將同時滿足具有較大密度和較大距離的點認定為類簇中心。對于具有較大距離的但是局部密度較低的點,可認定為噪聲點。對于剩余的點可將其分配到最近鄰的且密度更大點的所在簇中。所找到的類簇中心個數(shù)即為該數(shù)據(jù)內存在的工況個數(shù)。

■為人工選擇出的聚類中心圖2 密度峰值聚類決策圖Fig.2 Density peak clustering decision plot

通過該方法可以得到核電的運行工況類數(shù),為高斯混合模型的聚類初值提供參考,從而提高了高斯混合模型的精度。

步驟4 高斯模型是一種常用的變量分布模型,廣泛應用于數(shù)理領域。一維高斯分布的概率密度函數(shù)可定義為

(6)

式(6)中:μ為均值;σ為標準差;σ2為方差。

高斯混合模型可以理解為多個單高斯模型按一定的權重組合起來,每一個不同的類都分別代表一個高斯分布,因此可把高斯混合模型[24]的概率分布描述為

(7)

關于第k個高斯模型的響應度γk的計算公式為

(8)

為了使高斯混合模型得到最優(yōu)參數(shù),選擇使用EM算法來優(yōu)化參數(shù),EM算法可分為兩步,在進行這兩步工作之前需要設定子分布模型的參數(shù)的初值。E步:求取第k個分模型的對當前觀測數(shù)據(jù)的響應度γ′jk,M步:迭代求新一輪的模型參數(shù)(期望、方差、權重),當?shù)鷿M足|αk-αk-1|≤ε時即可終止迭代,其中ε為閾值。

(9)

(10)

(11)

高維核電數(shù)據(jù)在利用PCA完成降維后,降低了高斯混合模型的計算難度,參考密度峰值聚類的結果,確定高斯混合模型的初值K。將降維數(shù)據(jù)導入模型,迭代至模型均值和方差趨于穩(wěn)定的同時,利用模型響應度γk去排除模型響應度低的類,提高分類合理性。如圖3所示,為高斯混合模型的聚類結果圖,能夠找到聚類中心的聚類結果即為一類工況,未找到聚類中心的結果會根據(jù)模型響應度的大小被剔除掉。降維數(shù)據(jù)與原始數(shù)據(jù)的索引一致,因此在高斯混合模型劃分出工況以后,根據(jù)各工況中所屬數(shù)據(jù)點的索引,依次還原至各個傳感器的原始數(shù)據(jù)當中,找到原始數(shù)據(jù)中各點所屬的工況,實現(xiàn)最終的工況劃分。

圖3 高斯混合模型聚類結果圖Fig.3 Clustering results plot of Gaussian mixture model

2 實驗過程及結果分析

采用國內某核電廠的歷史運行數(shù)據(jù)來完成工況的劃分和算法有效性的驗證。該數(shù)據(jù)的記錄時間為2017年1月1日—2017年3月31日。如表1所示,本次實驗使用了其中12個傳感器的數(shù)據(jù),其中包括了液位、壓力、溫度和流量等傳感器。每個傳感器記錄了2 759 638條運行數(shù)據(jù),可以清晰還原整個系統(tǒng)的運行工況。

表1 參數(shù)分布Table 1 Parameters distribution

將上述12個參數(shù)的數(shù)據(jù)整理得到一個2 759 638×12的原始數(shù)據(jù)矩陣,利用PCA算法對該數(shù)據(jù)矩陣進行降維,選取總貢獻率達到55%的兩個向量作為投影向量(實際貢獻率為58.80%),得到了一個2 759 638×2的降維數(shù)據(jù)矩陣。如圖3所示,高斯混合模型的結果以兩個維度(第一主成分、第二主成分)的形式呈現(xiàn),即為降維數(shù)據(jù)的可視化呈現(xiàn)。使用密度峰值聚類算法對降維后的數(shù)據(jù)進行聚類。該算法將每個數(shù)據(jù)點的ρ值(密度)和δ值(距離)表示在一個二維的決策圖上。用戶根據(jù)決策圖的分布情況,對聚類中心點進行選擇,這是一個人工操作,無法自動完成。圖2為聚類所得的決策圖,在決策圖中靠近δ(縱軸)的點屬于噪聲點,這一類點密度小且距離其他點也遠,不能選為聚類核心,靠近ρ(橫軸)的屬于正常點,但是周圍有更為合適的點可以選作聚類核心。圖2所示紅圈標記出的3個點,即為人工選擇出的3個聚類中心,因為這3個點具有較大的局部密度以及距離,符合作為聚類中心的標準,因此將這3個點選擇為聚類中心點,每一個聚類中心點就代表著一類工況。由此可知實驗數(shù)據(jù)中所存在的工況個數(shù)為3類。該值為后續(xù)的高斯混合模型在選擇聚類初值時提供參考。

由于密度峰值聚類所得到的結果為3類工況,可以獲知在該數(shù)據(jù)類存在的工況種類應該為3類,其他不屬于這3類的點,會以散點的形式呈現(xiàn)。如圖3所示,將降維數(shù)據(jù)以散點形式呈現(xiàn)時,可以看出,有3個點簇以及其他的散點。為了不把散點強行劃入某一類工況之中,根據(jù)峰值密度聚類所得到的3類工況參考值,可在其參考值上加1,將高斯混合模型的聚類初值定為4。迭代運算100次后,高斯混合模型的方差和均值均趨于穩(wěn)定,此時可輸出最終的聚類結果。從圖3可以看出,Cluster0、Cluster1、Cluster3找到聚類核心。Cluster2包括其他的點,但是沒有找到聚類中心。高斯混合模型的聚類結果以點集的形式呈現(xiàn),處于不同工況的點所屬的集合也不相同,但每個散點所對應的時間點與實驗數(shù)據(jù)的時間點依然保持一致。

獲得了高斯混合模型的聚類結果后,可找到每類工況中所包含點所對應的時間點(索引),按時間點將所對應的工況情況還原到原始的傳感器數(shù)據(jù)當中。因此可以獲得各類傳感器在原始數(shù)據(jù)中的工況分布情況。以RCP012MN液位傳感器為例,該傳感器的數(shù)據(jù)占主成分的40.45%,為占比最大的傳感器,最能代表實驗數(shù)據(jù)的特點,其還原效果也最為優(yōu)秀。在還原過程中使用的是傳感器的原始數(shù)據(jù),利用高斯混合模型所得到的聚類結果,根據(jù)各個類別點所對應的時間點,依次還原至原數(shù)據(jù)點上,即可達到工況分類的可視化。其聚類初值均選擇為4時,利用本文算法(基于密度峰值聚類的高斯混合模型核電運行工況劃分算法)與K-means聚類算法、TICC聚類算法以及DBSCAN(density-based spatial clustering of applications with noise)[25]算法在RCP012MN傳感器數(shù)據(jù)上的工況劃分結果對比如圖4所示。

圖4 RCP012MN液位傳感器的工況劃分結果Fig.4 Working condition division results of RCP012MN liquid level sensor

圖4為各算法對傳感器RCP012MN采集的數(shù)據(jù)的工況劃分結果。該傳感器記錄了穩(wěn)壓器內的液位變化,當液位的處于不同的穩(wěn)定狀態(tài)時,代表著液位所處的不同工況,液位的變化過程不屬于任何工況。從圖4(c)中可以看出,本文算法的劃分結果僅呈現(xiàn)了工況0、工況1、工況3,這3種工況對應圖3中高斯混合模型結果中有聚類中心的3個點簇(Cluster0、Cluster1、Cluster3),未找到聚類中心的點簇中的點就被劃為過渡狀態(tài)(原始點)。如圖4(a)、圖4(b)、圖4(d)所示K-means、TICC和DBSCAN的劃分結果均將過渡態(tài)(變化區(qū)間)做出了錯誤劃分,同時在這3種方法的劃分結果中均出現(xiàn)工況穿插問題。

從劃分結果來看,高斯混合模型劃分的工況個數(shù)為3,工況2這一類劃分結果因為其模型響應度低的原因而遭到剔除,因此在還原圖中并未出現(xiàn)Cluster2的分類。從高斯混合模型的還原結果來看,不同工況的劃分邊界比較清晰,在每一類工況區(qū)域中幾乎未出現(xiàn)明顯的工況穿插情況,在變化態(tài)階段以原始點的形式出現(xiàn),這是因為過渡段的點是不屬于任何工況,這也說明劃分結果的合理性。同時,由于在工況劃分過程中,對于工況變化的判斷是通過變化趨勢來判別的,當變化趨勢持續(xù)出現(xiàn)時才會將這一段判斷為變化態(tài),因此會出現(xiàn)變化態(tài)邊界上的點被劃分到某個工況的情況。

通過人工判斷,屬于工況0中的數(shù)據(jù)點應該為1 630 080個(2017-01-01T00:59:27—2017-01-21T11:12:53和2017-02-14T11:23:02—2017-03-10T22:19:23),工況1中的數(shù)據(jù)點應為158 201個(2017-01-21T23:29:24—2017-02-05T12:13:00),工況3中的數(shù)據(jù)點應為907 447個(2017-02-06T21:13:32—2017-02-14T09:19:21和2017-03-11T08:45:25—2017-03-31T23:59:54)。在通過人工方式去除各算法得出劃分結果中的錯誤分類點后,可知利用本文算法獲得的工況0中的數(shù)據(jù)點為1 618 540個,工況1中的數(shù)據(jù)點為158 201個,工況3中的數(shù)據(jù)點為885 380個,三類工況中的錯分總數(shù)為33 607個。TICC算法的劃分結果為工況0中的數(shù)據(jù)點為1 603 661個,工況1中的數(shù)據(jù)點為153 850個,工況3中的數(shù)據(jù)點為782 568個,三類工況中的錯分總數(shù)為155 649個。K-means算法的劃分結果為工況0中的數(shù)據(jù)點為1 417 213個,工況1中的數(shù)據(jù)點為158 201個,工況3中的數(shù)據(jù)點為519 800個,三類工況中的錯分總數(shù)為700 514個。DBSCAN算法的劃分結果為工況0中的數(shù)據(jù)點為1 581 140個,工況1中的數(shù)據(jù)點為149 000個,工況3中的數(shù)據(jù)點為904 147個,三類工況中的錯分總數(shù)為61 441個。在此工況劃分準確率為劃分出的數(shù)據(jù)點數(shù)與該類工況準確點數(shù)的商。錯誤率為錯誤劃分個數(shù)與三類工況的準確數(shù)據(jù)點數(shù)總和的商。表2為各算法的劃分準確率對比。本文算法對工況0、工況1、工況3的劃分準確率分別達到99.29%,100%,97.57%,且錯誤率僅為1.25%。其他3種算法中DBSCAN的效果最好,但是其在工況0和工況1的劃分準確率上低于本文算法,且錯誤率也更高。

表2 算法準確率對比Table 2 Algorithm accuracy comparison

通過實驗表明,K-means算法、TICC算法以及DBSCAN算法與本文算法相比,本文算法對三類工況的劃分準確率均高于其他三類算法,且錯誤率更低。K-means算法錯誤率高的原因在于K-means在聚類過程中只考慮了數(shù)據(jù)點的空間分布卻忽略了數(shù)據(jù)的時間順序問題,因此造成了很多數(shù)據(jù)點被錯誤歸類。TICC算法錯誤率較高的原因是該算法是基于滑動窗口來做聚類,每一次滑動均會在該窗口內做一次聚類,從而出現(xiàn)了穿插結果的出現(xiàn)。這種穿插現(xiàn)象的出現(xiàn)對劃分精度產(chǎn)生了極大的影響。DBSCAN算法雖然在準確率和錯誤率上的表現(xiàn)均比較好,但在對過渡態(tài)數(shù)據(jù)的處理上,DBSCAN算法與K-means算法和TICC算法一樣錯誤的將過渡段劃分為一類工況,根據(jù)專家經(jīng)驗,過渡狀態(tài)是不屬于任何一類工況的,因此這種劃分方式是不合理的。

3 結論

提出了一種基于密度峰值聚類的高斯混合模型核電工況劃分方法。該方法有效地解決了高維核電數(shù)據(jù)工況類數(shù)難以確定的問題,且劃分結果符合實際情況。將本文方法與K-means方法、TICC方法以及DBSCAN算法應用于真實核電數(shù)據(jù)中,得出如下結論。

(1)本文方法的劃分準確率更高,其三類工況的劃分準確率分別達到99.29%、100%、97.57%,且錯誤率僅為1.25%,減少了類與類之間穿插現(xiàn)象的出現(xiàn)。

(2)在對過渡段的處理上,本文方法較K-means、TICC以及DBSCAN的劃分結果更具有合理性。