化：居民用電負(fù)荷曲線用戶行為分析實(shí)戰(zhàn))
做居民用電行為分析最頭疼的往往不是算法本身而是數(shù)據(jù)背后的規(guī)律看不到。你拿到的是一堆負(fù)荷曲線怎么告訴別人這小區(qū)里哪些用戶是上班族、哪些是全天在家、哪些可能偷偷搞生產(chǎn)經(jīng)營(yíng)聚類就是個(gè)好幫手。但真去用Kmeans的時(shí)候那個(gè)初始質(zhì)心的選擇問題足夠讓人抓狂——同一個(gè)數(shù)據(jù)集換一次初始值跑出來就是另一種分法。我這段時(shí)間正好在Matlab里把粒子群算法和Kmeans拼在一起拿居民用電負(fù)荷數(shù)據(jù)做行為分析整個(gè)過程踩了不少坑也把關(guān)鍵細(xì)節(jié)理清楚了。這篇就把PSO-Kmeans聚類的思路、代碼實(shí)現(xiàn)和實(shí)戰(zhàn)經(jīng)驗(yàn)完整走一遍適合正在做負(fù)荷分析、用戶畫像或者剛接觸聚類優(yōu)化的朋友參考。1. 為什么居民用電分析繞不開聚類優(yōu)化1.1 負(fù)荷曲線背后的用戶行為差異先看數(shù)據(jù)。居民用戶的用電行為最直觀的載體就是日負(fù)荷曲線。把一天24小時(shí)或者96個(gè)采樣點(diǎn)的功率值串聯(lián)起來就構(gòu)成一條代表用戶當(dāng)天用電習(xí)慣的曲線。正常來說上班族用戶的工作日負(fù)荷曲線會(huì)出現(xiàn)明顯的兩峰一谷——早高峰可能是7點(diǎn)到9點(diǎn)晚高峰是18點(diǎn)到22點(diǎn)白天和夜晚則相對(duì)平緩而老人家庭、待業(yè)在家的用戶白天的負(fù)荷往往比上班族高出一截曲線形態(tài)更平穩(wěn)家里有電動(dòng)汽車的用戶可能晚上會(huì)出現(xiàn)一個(gè)持續(xù)的充電功率平臺(tái)。不同行為模式之間的差異在負(fù)荷曲線上是有跡可循的。但問題在于一個(gè)城市或一個(gè)臺(tái)區(qū)往往有幾千幾萬個(gè)用戶靠人工去看曲線分門別類根本不現(xiàn)實(shí)。聚類就是用來做這件事的——它能把相似形態(tài)的負(fù)荷曲線自動(dòng)歸到同一組讓有相同用電習(xí)慣的用戶自然聚到一起。只要聚類方法靠譜分出來的每一類用戶你都能倒推出一套對(duì)應(yīng)的行為描述這對(duì)接下來的需求響應(yīng)、分時(shí)電價(jià)策略、臺(tái)區(qū)負(fù)荷預(yù)測(cè)都特別有價(jià)值。1.2 標(biāo)準(zhǔn)Kmeans的天然短板在負(fù)荷聚類這個(gè)場(chǎng)景里最常用的是Kmeans算法。原理其實(shí)特別直白先在樣本空間里挑K個(gè)點(diǎn)當(dāng)初始質(zhì)心然后把每個(gè)樣本分配給離它最近的質(zhì)心分完以后再重新計(jì)算每個(gè)簇的中心點(diǎn)反復(fù)迭代直到結(jié)果穩(wěn)定。整個(gè)過程就是分配—更新—再分配—再更新很符合直覺代碼也簡(jiǎn)單。但它有個(gè)致命問題——Kmeans對(duì)初始質(zhì)心的選擇極度敏感。初始質(zhì)心選偏了迭代多少次都可能停留在某個(gè)局部最優(yōu)解上。比如有兩類用戶一類是白天用電一類是晚上用電如果你的初始質(zhì)心都落在白天那堆數(shù)據(jù)里晚上那一類很可能被硬生生拆散最終聚類結(jié)果從業(yè)務(wù)角度怎么解釋都不合理。而且Kmeans迭代過程中一旦某個(gè)簇被分配為空算法還會(huì)出現(xiàn)質(zhì)心失效的異常情況。實(shí)際處理負(fù)荷數(shù)據(jù)時(shí)樣本量大、曲線波動(dòng)多Kmeans跑出來的結(jié)果經(jīng)常不穩(wěn)定同一份數(shù)據(jù)跑十次能有七八種分法。這就是我要引入粒子群算法的直接原因。粒子群優(yōu)化算法PSO是一種全局尋優(yōu)方法它的思路是模擬鳥群覓食——每個(gè)粒子代表一個(gè)候選解靠個(gè)體經(jīng)驗(yàn)個(gè)體最優(yōu)和群體經(jīng)驗(yàn)全局最優(yōu)不斷調(diào)整自己的位置逐步逼近全局最優(yōu)解。把PSO和Kmeans結(jié)合簡(jiǎn)單說就是用PSO先把Kmeans的初始質(zhì)心這個(gè)老大難問題解決掉讓聚類從一個(gè)好的起點(diǎn)開始跑。這樣既保留了Kmeans計(jì)算快的優(yōu)點(diǎn)又顯著降低了落入局部最優(yōu)的概率。下面我把這套組合方案的原理和代碼一步步拆開講。2. PSO和Kmeans是怎么配合的2.1 先從Kmeans的目標(biāo)說起要理解PSO-Kmeans的配合邏輯得先把Kmeans在做什么看透。Kmeans本質(zhì)上是求解一個(gè)最小化問題把n個(gè)樣本分到K個(gè)簇中讓所有樣本到所屬簇質(zhì)心的距離平方和也就是簇內(nèi)誤差平方和SSE最小。公式寫出來是[ SSE \sum_{i1}^{K}\sum_{x \in C_i} |x - \mu_i|^2 ]其中(C_i)是第i個(gè)簇(\mu_i)是這個(gè)簇的質(zhì)心。聚類結(jié)果好不好直接看SSE——SSE越小說明簇內(nèi)樣本越緊湊同類用戶之間的相似度越高。但Kmeans采用的是一種貪心式的交替優(yōu)化先固定質(zhì)心分配樣本再固定分配更新質(zhì)心。這種方式求解速度快卻很依賴初始質(zhì)心給得怎么樣。初始質(zhì)心離全局最優(yōu)解太遠(yuǎn)交替優(yōu)化就可能收斂到SSE較大的局部最優(yōu)解。怎么跳出這個(gè)坑我在項(xiàng)目里采用的思路是把Kmeans的初始質(zhì)心當(dāng)作粒子群優(yōu)化算法中的決策變量用PSO去全局搜索一組好的質(zhì)心位置再把搜到的結(jié)果作為Kmeans的起點(diǎn)。換句話說Kmeans負(fù)責(zé)局部精修PSO負(fù)責(zé)全局尋優(yōu)兩者分工合作。這里有一個(gè)技術(shù)路線選擇問題。有些人會(huì)把PSO直接作為聚類工具來用讓每個(gè)樣本以一定概率歸屬某個(gè)簇最后按隸屬度劃分。這個(gè)方案在高維負(fù)荷數(shù)據(jù)上計(jì)算量非常大而且解釋性不如Kmeans清晰。我更推薦的是把PSO定位成質(zhì)心初始化優(yōu)化器后面照常跟Kmeans迭代。實(shí)測(cè)下來這種方式既快又穩(wěn)結(jié)果還好用業(yè)務(wù)語言解釋。2.2 粒子編碼方式與適應(yīng)度函數(shù)設(shè)計(jì)要把PSO用于優(yōu)化Kmeans的初始質(zhì)心第一個(gè)要解決的是粒子怎么編碼的問題。假設(shè)負(fù)荷數(shù)據(jù)經(jīng)過特征處理后每個(gè)樣本是一個(gè)D維向量最簡(jiǎn)單的做法就是24維對(duì)應(yīng)24小時(shí)負(fù)荷值我實(shí)驗(yàn)里也用過96維對(duì)應(yīng)96個(gè)采樣點(diǎn)聚類數(shù)設(shè)為K。那么一組完整的初始質(zhì)心就是K個(gè)D維向量把它們按順序拼接成一個(gè)長(zhǎng)向量這個(gè)長(zhǎng)向量就是一個(gè)粒子的位置。粒子的維度就是(K \times D)。舉個(gè)例子K4D24粒子維度就是96。在Matlab里我習(xí)慣用矩陣來組織蜂群——一個(gè)粒子用一個(gè)(K \times D)的矩陣表示整個(gè)粒子群用一個(gè)三維數(shù)組存儲(chǔ)這樣在計(jì)算距離時(shí)可以避免頻繁的reshape操作。適應(yīng)度函數(shù)的設(shè)計(jì)是整個(gè)算法的靈魂。PSO的迭代方向完全靠適應(yīng)度值牽引。我這里使用的適應(yīng)度函數(shù)就是Kmeans的SSE[ fitness(x) \sum_{j1}^{K}\sum_{x_n \in C_j} |x_n - z_j|^2 ]其中(z_j)是粒子的第j個(gè)質(zhì)心位置。計(jì)算時(shí)先按每個(gè)樣本到各質(zhì)心的歐氏距離做最近鄰分配再累加出SSE。適應(yīng)度值越小代表該粒子對(duì)應(yīng)的質(zhì)心組合越好。有人會(huì)問PSO時(shí)代里的好到底是全局最優(yōu)還是局部最優(yōu)這正是PSO的價(jià)值所在——粒子群中的每個(gè)粒子都在自己的位置附近搜索同時(shí)向全局最優(yōu)粒子靠攏這種信息共享機(jī)制讓種群不容易卡死在單個(gè)局部區(qū)域。配合慣性權(quán)重和學(xué)習(xí)因子PSO能在搜索前期保持較強(qiáng)的全局探索能力后期逐漸收斂到精細(xì)區(qū)域。這比隨機(jī)撒點(diǎn)選初始質(zhì)心要靠譜得多。2.3 算法流程梳理我實(shí)際跑通的PSO-Kmeans完整流程如下讀取并預(yù)處理負(fù)荷數(shù)據(jù)缺失值處理、歸一化。確定聚類數(shù)K用輪廓系數(shù)或肘部法則輔助判斷后面細(xì)說。初始化粒子群體。每個(gè)粒子的位置為K個(gè)隨機(jī)的樣本點(diǎn)速度為全零或小隨機(jī)數(shù)。對(duì)每個(gè)粒子計(jì)算適應(yīng)度SSE更新個(gè)體最優(yōu)pbest和群體最優(yōu)gbest。按標(biāo)準(zhǔn)PSO公式更新粒子速度和位置。速度更新公式為 [ v_{t1} w \cdot v_t c_1 r_1 (pbest - x_t) c_2 r_2 (gbest - x_t) ] 位置更新為簡(jiǎn)單累加。這里w是慣性權(quán)重c1和c2是學(xué)習(xí)因子r1和r2是0到1之間的隨機(jī)數(shù)。檢查是否達(dá)到最大迭代次數(shù)否則返回第4步。將gbest還原為(K \times D)的質(zhì)心矩陣作為Kmeans的初始質(zhì)心。執(zhí)行標(biāo)準(zhǔn)Kmeans迭代分配樣本、更新質(zhì)心直到收斂。輸出聚類標(biāo)簽、質(zhì)心、SSE并做可視化。整個(gè)流程里PSO階段其實(shí)相當(dāng)于在做全局熱身Kmeans階段在做局部沖刺。我用這個(gè)方案對(duì)比過純Kmeans在典型居民負(fù)荷數(shù)據(jù)上SSE能降低15%到25%而且多次運(yùn)行的結(jié)果穩(wěn)定性明顯提升。3. Matlab代碼實(shí)現(xiàn)與參數(shù)配置3.1 數(shù)據(jù)準(zhǔn)備與特征構(gòu)建代碼實(shí)現(xiàn)上我建議先把特征工程做成獨(dú)立腳本不要把數(shù)據(jù)處理和聚類算法混在一起。原始用電數(shù)據(jù)通常是這樣的每15分鐘一個(gè)采樣點(diǎn)一天96個(gè)點(diǎn)連續(xù)若干天。但如果直接用96維做聚類維度較高PSO粒子搜索空間的體積會(huì)指數(shù)增長(zhǎng)不僅慢而且效果不一定好。我的做法是折中做日負(fù)荷曲線特征壓縮。常用的壓縮方式有幾種我列個(gè)表對(duì)比一下特征方案維度優(yōu)點(diǎn)缺點(diǎn)24小時(shí)均值負(fù)荷24直觀、計(jì)算快能保留峰谷形態(tài)丟失了日內(nèi)變化細(xì)節(jié)96點(diǎn)原始負(fù)荷96信息完整維度高PSO粒子維度過大易過擬合峰谷特征峰時(shí)負(fù)荷、谷時(shí)負(fù)荷、峰谷差、日用電量等5-8業(yè)務(wù)解釋性強(qiáng)維度低需要按當(dāng)?shù)胤骞葧r(shí)段定義有一定主觀性統(tǒng)計(jì)特征均值、方差、峰度、偏度、最大負(fù)荷時(shí)間5-8壓縮程度高形態(tài)信息流失多我在項(xiàng)目里最終選的是24小時(shí)均值負(fù)荷幾個(gè)統(tǒng)計(jì)特征的組合總維度約28。折中的原因有兩個(gè)一是24小時(shí)曲線能讓聚類結(jié)果直接畫圖解釋生成工作族居家型這種標(biāo)簽二是維度控制在30以內(nèi)PSO搜索效率高很多。作為補(bǔ)充我也做了96維的對(duì)照實(shí)驗(yàn)后面在問題排查部分會(huì)講這個(gè)方案踩了什么坑。數(shù)據(jù)清洗這一步很關(guān)鍵。居民負(fù)荷數(shù)據(jù)里常見的問題是采集終端偶爾掉線導(dǎo)致整天數(shù)據(jù)是0或者個(gè)別時(shí)段出現(xiàn)異常尖峰。我的處理規(guī)則是連續(xù)3小時(shí)以上全為0的用戶直接剔除非零時(shí)段中超過99.5%分位的數(shù)值視為異常尖峰用前后時(shí)刻的均值替換。這些規(guī)則比單純用是否大于某閾值判斷更魯棒。歸一化也要特別注意。如果不做歸一化用電量大的用戶比如冬夏開空調(diào)日電量幾十度甚至上百度會(huì)在歐氏距離計(jì)算中占據(jù)絕對(duì)主導(dǎo)聚類結(jié)果基本就變成了按用電量分等級(jí)而不是按行為模式分類。我的做法是按特征列做Z-score標(biāo)準(zhǔn)化也就是每列減去均值再除以標(biāo)準(zhǔn)差這樣每個(gè)特征對(duì)距離的貢獻(xiàn)平等。在Matlab里一行代碼就能搞定data_norm zscore(data_raw);處理完以后記得保存一份標(biāo)準(zhǔn)化參數(shù)后面做新用戶分類或者畫原尺度曲線時(shí)要用。3.2 PSO-Kmeans主程序編寫主程序我分了三個(gè)函數(shù)塊粒子初始化、適應(yīng)度計(jì)算、PSO迭代主循環(huán)。這種模塊化寫法方便調(diào)試也便于替換不同的適應(yīng)度函數(shù)或數(shù)據(jù)集。先看粒子初始化% 輸入data為標(biāo)準(zhǔn)化后的樣本矩陣(nxD)K為聚類數(shù)N為種群規(guī)模 % 輸出particle為(N, K, D)的三維數(shù)組 n size(data, 1); D size(data, 2); particle zeros(N, K, D); velocity zeros(N, K, D); for i 1:N idx randperm(n, K); % 隨機(jī)選K個(gè)樣本作為初始質(zhì)心 particle(i, :, :) data(idx, :); velocity(i, :, :) 0.02 * randn(K, D); end初始化方式選擇隨機(jī)取樣本點(diǎn)而不是在整個(gè)搜索空間隨機(jī)撒點(diǎn)。原因是負(fù)荷數(shù)據(jù)做完Z-score標(biāo)準(zhǔn)化后雖然有少數(shù)離群點(diǎn)但絕大多數(shù)樣本都集中在可行區(qū)域內(nèi)。從樣本中選初始質(zhì)心相當(dāng)于一開始就沒有偏離合理區(qū)域能明顯加快收斂。這個(gè)細(xì)節(jié)我建議一定保留。適應(yīng)度函數(shù)我單獨(dú)寫核心邏輯如下function fitness calcFitness(data, particle_i, K) n size(data, 1); distMat zeros(n, K); for j 1:K centroid squeeze(particle_i(j, :)); diff data - centroid; % n x D distMat(:, j) sqrt(sum(diff.^2, 2)); end [~, assign] min(distMat, [], 2); fitness 0; for j 1:K clusterData data(assign j, :); if ~isempty(clusterData) centroid mean(clusterData, 1); fitness fitness sum(sum((clusterData - centroid).^2, 2)); end end end注意這里我在適應(yīng)度計(jì)算中不是用粒子自帶質(zhì)心算SSE而是按分配結(jié)果重新計(jì)算實(shí)際質(zhì)心再算SSE。為什么不直接用粒子里的質(zhì)心因?yàn)榱Ｗ釉赑SO迭代中可能移動(dòng)到遠(yuǎn)離任何樣本的位置用空簇質(zhì)心算距離會(huì)產(chǎn)生虛低的SSE誤導(dǎo)搜索方向。重新計(jì)算簇質(zhì)心相當(dāng)于做了局部投影適應(yīng)度值更真實(shí)。這個(gè)細(xì)節(jié)是我調(diào)試過程中對(duì)比了幾種方案后確定的效果確實(shí)更穩(wěn)。主迭代循環(huán)采用標(biāo)準(zhǔn)的PSO公式慣性權(quán)重w隨迭代次數(shù)線性遞減maxIter 50; N 30; K 4; c1 1.5; c2 1.5; wMax 0.9; wMin 0.4; pbestScore inf(N, 1); pbestParticle particle; gbestScore inf; gbestParticle squeeze(particle(1, :, :)); for t 1:maxIter w wMax - (wMax - wMin) * t / maxIter; for i 1:N fitness calcFitness(data, squeeze(particle(i, :, :)), K); if fitness pbestScore(i) pbestScore(i) fitness; pbestParticle(i, :, :) particle(i, :, :); end if fitness gbestScore gbestScore fitness; gbestParticle squeeze(particle(i, :, :)); end end for i 1:N r1 rand(K, D); r2 rand(K, D); velocity(i, :, :) w * velocity(i, :, :) ... c1 * r1 .* (squeeze(pbestParticle(i, :, :)) - squeeze(particle(i, :, :))) ... c2 * r2 .* (gbestParticle - squeeze(particle(i, :, :))); particle(i, :, :) particle(i, :, :) velocity(i, :, :); end end最后把gbestParticle作為初始質(zhì)心送給Kmeans[clusterIdx, centroid] kmeans(data, K, Start, gbestParticle, MaxIter, 1000);如果Matlab版本較老不支持Start參數(shù)直接傳入矩陣可以先調(diào)用類的靜態(tài)方法設(shè)置選項(xiàng)再執(zhí)行聚類或者自己手寫10-20輪Kmeans迭代。老版本其實(shí)也完全可以用我后面遇到過一次版本兼容問題在常見問題部分會(huì)展開說明。3.3 關(guān)鍵參數(shù)的選擇依據(jù)與調(diào)試建議PSO-Kmeans涉及到的參數(shù)不少我把我實(shí)測(cè)下來比較合適的配置整理一下。種群規(guī)模N我建議取20到40之間。太小了全局搜索能力不足太大了計(jì)算量明顯上升。負(fù)荷曲線的樣本數(shù)通常在幾千到幾萬之間每次適應(yīng)度計(jì)算都要遍歷所有樣本做距離計(jì)算N取30不算大但加上50次迭代在幾千樣本量下Matlab要跑幾十秒可以接受。如果樣本量超過5萬建議先把訓(xùn)練集采樣到1萬規(guī)模做粒子搜索再用跑出來的質(zhì)心初始化全量Kmeans。最大迭代次數(shù)maxIter50次通常夠了。我在調(diào)試時(shí)觀察過適應(yīng)度收斂曲線大約在30次以后下降曲線就趨于平緩50次屬于留有余量。如果追求速度25到30次也能得到差不多的結(jié)果差別在2%以內(nèi)。但首次實(shí)驗(yàn)我建議還是跑到50次先把算法的穩(wěn)定基線摸清楚。慣性權(quán)重w采用0.9到0.4線性遞減。前期w大粒子飛得快、探索范圍廣不容易陷進(jìn)局部最優(yōu)后期w小粒子精細(xì)琢磨加速收斂。這個(gè)區(qū)間是粒子群算法的經(jīng)典經(jīng)驗(yàn)值實(shí)測(cè)在聚類問題上效果穩(wěn)定。學(xué)習(xí)因子c1和c2取1.5是比較均衡的組合。也有文獻(xiàn)推薦c1c22我試過收斂快一些但偶爾會(huì)跳過好的質(zhì)心區(qū)域。1.5加上0.9到0.4的慣性權(quán)重搭配探索和開發(fā)平衡得更舒服。如果你發(fā)現(xiàn)結(jié)果波動(dòng)大可以嘗試把c1降到1.2、c2提到1.8增強(qiáng)向群體最優(yōu)靠攏的趨勢(shì)。聚類數(shù)K用輪廓系數(shù)輔助判斷。輪廓系數(shù)綜合考慮了簇內(nèi)緊密度和簇間分離度取值范圍-1到1越大代表聚類效果越好。我在項(xiàng)目里對(duì)K2到K8分別跑PSO-Kmeans計(jì)算每個(gè)K下的平均輪廓系數(shù)選峰值對(duì)應(yīng)的K。實(shí)際業(yè)務(wù)上K取4或5比較常見這樣每一類用戶都有足夠明確的畫像不會(huì)分得過細(xì)而失去解釋力。4. 實(shí)驗(yàn)效果分析與聚類結(jié)果解讀4.1 與標(biāo)準(zhǔn)Kmeans的對(duì)比實(shí)驗(yàn)我拿來驗(yàn)證的數(shù)據(jù)是某市一個(gè)臺(tái)區(qū)3000戶居民用戶30天的用電記錄按前文方法清洗和特征化后得到3000×28的特征矩陣聚類目標(biāo)K4PSO種群取30迭代50次。為了控制變量標(biāo)準(zhǔn)Kmeans我用Matlab自帶的kmeans函數(shù)跑100次隨機(jī)初始化取SSE最小的一次作為參照這種多次隨機(jī)取最優(yōu)本身就是實(shí)踐中應(yīng)對(duì)Kmeans不穩(wěn)定的常見手段但計(jì)算開銷遠(yuǎn)高于PSO輔助。最終實(shí)驗(yàn)數(shù)據(jù)如下表方案平均SSE最優(yōu)SSE波動(dòng)范圍SSE單次運(yùn)行耗時(shí)標(biāo)準(zhǔn)Kmeans單次1846.71752.3160.40.8秒標(biāo)準(zhǔn)Kmeans100次取最優(yōu)1635.21635.2024秒PSO-Kmeans單次1658.11641.533.218秒PSO-Kmeans3次取最優(yōu)1642.01641.53.254秒幾個(gè)結(jié)論很直觀。PSO-Kmeans單次結(jié)果明顯優(yōu)于Kmeans單次SSE從1846.7降到1658.1下降了大約10.2%即使對(duì)比Kmeans跑100次取最優(yōu)的1635.2PSO-Kmeans的最優(yōu)SSE 1641.5也非常接近差了不到0.4%。更關(guān)鍵的是穩(wěn)定性——PSO-Kmeans三次運(yùn)行的最優(yōu)與最差只差33.2幾乎都在同一水平線上這說明算法已經(jīng)不太受隨機(jī)初始化的影響而標(biāo)準(zhǔn)Kmeans單次運(yùn)行的波動(dòng)范圍高達(dá)160以上這在工程上非常致命。當(dāng)然PSO-Kmeans也不是免費(fèi)的午餐18秒的處理時(shí)間比標(biāo)準(zhǔn)Kmeans單次0.8秒慢得多。但對(duì)離線用戶畫像分析這種場(chǎng)景18秒完全可接受。4.2 聚類結(jié)果如何映射到用電行為聚類跑完只是第一步更重要的工作是把每一類用戶的行為模式描述出來。我是這樣做的拿到聚類標(biāo)簽后把原始負(fù)荷數(shù)據(jù)未標(biāo)準(zhǔn)化按類分組計(jì)算每類用戶的平均24小時(shí)負(fù)荷曲線然后結(jié)合日用電量、峰谷比等業(yè)務(wù)指標(biāo)做解讀。在我的實(shí)驗(yàn)里K4時(shí)的四類用戶畫像如下第一類工作日早、晚雙峰特別突出白天負(fù)荷很低午間有小幅回落周末曲線相對(duì)平緩。結(jié)合日用電量處于中低水平可以判定為典型的上班族家庭工作日只有早晚在家用電。第二類白天負(fù)荷較高曲線全天相對(duì)平穩(wěn)夜晚略降但不會(huì)降到很低日用電量處于中上水平。這是全天居家型用戶可能是老人、家庭主婦或自由職業(yè)者。第三類夜間和凌晨負(fù)荷異常偏高白天反而較低日用電量也比較大。結(jié)合當(dāng)?shù)仉妰r(jià)政策這類用戶很可能是有意將洗衣機(jī)、熱水器等大功率設(shè)備挪到夜間使用甚至可能有電動(dòng)汽車充電行為。第四類整體負(fù)荷水平低曲線平緩無峰長(zhǎng)時(shí)間維持很小的用電功率。這種通常是空心戶或者出租率較高的房屋用電行為不活躍。每類用戶對(duì)應(yīng)的策略建議也不一樣第一類適合宣傳分時(shí)電價(jià)引導(dǎo)削峰填谷第二類可以推薦節(jié)能設(shè)備第三類可以作為需求響應(yīng)的重點(diǎn)對(duì)象第四類則需要在臺(tái)區(qū)管理上排查是否有空置房或者表計(jì)異常。這些業(yè)務(wù)層面的延伸才是分析工作真正產(chǎn)生價(jià)值的地方。4.3 可視化技巧如何把聚類結(jié)果畫得讓業(yè)務(wù)方看懂聚類結(jié)果可視化我踩過不少坑。最開始我直接用plot畫所有用戶的原始曲線3000條線疊在一起密密麻麻根本看不出差異。后來改成每個(gè)類畫一條平均曲線標(biāo)準(zhǔn)差帶效果立刻不一樣。Matlab里用fill可以畫帶meanCurve mean(clusterData, 1); stdCurve std(clusterData, 1); t 1:24; fill([t fliplr(t)], [meanCurvestdCurve fliplr(meanCurve-stdCurve)], ... [0.9 0.9 0.9], FaceAlpha, 0.4, EdgeColor, none); hold on; plot(t, meanCurve, LineWidth, 2);標(biāo)準(zhǔn)差帶能夠直觀表達(dá)這一類用戶內(nèi)部的波動(dòng)程度。如果某類的帶很窄說明這類用戶的負(fù)荷形態(tài)高度一致聚類可信度高帶很寬則說明這一類內(nèi)部還存在細(xì)分可以考慮是否增加K值。另外一個(gè)可視化技巧是降維散點(diǎn)圖。高維特征矩陣不好直接展示可以用t-SNE或者PCA降到2維再按聚類標(biāo)簽著色。不過我要提醒一句降維后再看聚類是否分得開只能作為輔助參考因?yàn)榻稻S過程會(huì)扭曲真實(shí)距離關(guān)系。業(yè)務(wù)匯報(bào)時(shí)這東西很好看內(nèi)部驗(yàn)證時(shí)別太當(dāng)真。5. 常見問題與排查技巧實(shí)錄5.1 粒子維度爆炸和計(jì)算速度慢怎么辦我在96維特征上嘗試過直接跑PSO-Kmeans粒子維度是(K \times 96)K取4就是384維。粒子群優(yōu)化在這么高的維度上進(jìn)行搜索效果非常差——適應(yīng)度收斂慢、粒子群容易散開、結(jié)果還不穩(wěn)定。因?yàn)楦呔S空間里距離度量變得稀疏隨機(jī)初始化的粒子互相之間差異很小PSO很難通過對(duì)比分辨哪個(gè)方向更好。解決思路有兩個(gè)。第一是在特征層面降維比如用24小時(shí)均值替代96點(diǎn)數(shù)據(jù)或者先用PCA把特征壓到15到20維再做聚類。第二是改變PSO的搜索策略比如將速度初始化設(shè)置為0限制粒子的搜索半徑但這樣又會(huì)犧牲全局搜索能力。我的建議是優(yōu)先做特征降維因?yàn)榫用褙?fù)荷數(shù)據(jù)本身的冗余度很高96個(gè)采樣點(diǎn)之間存在很強(qiáng)的時(shí)序相關(guān)性強(qiáng)行保留全部維度得不償失。計(jì)算速度問題還有另一層來源適應(yīng)度函數(shù)里頻繁的矩陣運(yùn)算。如果循環(huán)寫的效率低幾千樣本都?jí)蜃孧atlab卡上幾分鐘。我把計(jì)算距離的代碼從for循環(huán)改成矩陣廣播后原來45秒一次迭代縮到3秒左右。Matlab效率的關(guān)鍵就是不要讓循環(huán)套循環(huán)多用維度廣播和矩陣運(yùn)算如果還想更快可以把calcFitness寫成mex函數(shù)或者用parfor并行計(jì)算粒子群中不同粒子的適應(yīng)度。5.2 陷入局部最優(yōu)的判斷與處理有一種情況PSO迭代結(jié)束后gbest對(duì)應(yīng)的質(zhì)心組其實(shí)還不是理想解Kmeans再迭代也跳不出來。怎么判斷我會(huì)把PSO-Kmeans的SSE和多次隨機(jī)初始化的Kmeans最優(yōu)SSE做對(duì)比如果前者顯著大于后者基本可以斷定PSO階段早收斂了。處理辦法有這么幾種。一是檢查粒子群初始化如果初始粒子全都擠在樣本集中的區(qū)域多樣性不夠PSO很容易早熟。初始化時(shí)除了隨機(jī)采樣樣本點(diǎn)我還會(huì)刻意加幾個(gè)遠(yuǎn)離中心的點(diǎn)。二是增大慣性權(quán)重或者調(diào)節(jié)學(xué)習(xí)因子如果w從0.9降到0.4太快個(gè)體經(jīng)驗(yàn)權(quán)重過大可以在實(shí)驗(yàn)中把wMax提到1.0wMin提到0.5讓粒子飛得更激進(jìn)一點(diǎn)。三是重啟策略如果一個(gè)粒子連續(xù)N代都沒有改進(jìn)自己的pbest給它重新初始化到隨機(jī)位置這是個(gè)簡(jiǎn)單但很有效的辦法。我再強(qiáng)調(diào)一次PSO-Kmeans不是銀彈它只能顯著降低落入局部最優(yōu)的概率不能完全消除。所以在項(xiàng)目落地時(shí)我通常跑3次PSO-Kmeans取SSE最小的那次。由于單次已經(jīng)很穩(wěn)定3次取最優(yōu)帶來的額外收益也有限更多是買個(gè)心理保險(xiǎn)。5.3 K值怎么選最合理選擇K值最常見的是肘部法則畫SSE隨K變化的折線圖找那個(gè)拐點(diǎn)。但實(shí)際數(shù)據(jù)里肘部往往不明顯SSE下降曲線保持平滑你很難說出3和4哪個(gè)是肘。我用輪廓系數(shù)配合業(yè)務(wù)可解釋性一起判斷。輪廓系數(shù)對(duì)第i個(gè)樣本的定義是[ s_i \frac{b_i - a_i}{\max(a_i, b_i)} ]其中(a_i)是樣本i與同簇其他樣本的平均距離(b_i)是樣本i與最近其他簇的平均距離。把全部樣本的輪廓系數(shù)平均就是總體輪廓系數(shù)。我一般要求總體輪廓系數(shù)大于等于0.5如果某個(gè)K下只有0.3說明簇內(nèi)不夠緊湊或者簇間分得不清楚這個(gè)K值基本不可用。但我也要說業(yè)務(wù)可解釋性有時(shí)候比數(shù)值指標(biāo)更關(guān)鍵。比如K5時(shí)輪廓系數(shù)最高但其中有一類用戶曲線形態(tài)和另一類非常接近業(yè)務(wù)上完全無法區(qū)分和應(yīng)對(duì)那K5就沒有實(shí)際意義。我的習(xí)慣是先選2到3個(gè)候選K輪廓系數(shù)比較高的然后把這幾個(gè)K下的聚類結(jié)果拿給業(yè)務(wù)同事看問哪一版最容易講故事通常答案很明確。5.4 版本兼容和Matlab環(huán)境的坑我在實(shí)驗(yàn)過程中遇到過一次運(yùn)行環(huán)境導(dǎo)致的怪問題在Matlab R2021b上能正常運(yùn)行的腳本換到老版本后kmeans的Start參數(shù)傳矩陣就報(bào)錯(cuò)。Matlab每個(gè)版本對(duì)聚類函數(shù)輸入?yún)?shù)的校驗(yàn)機(jī)制不一樣如果公司或?qū)嶒?yàn)室的Matlab版本不統(tǒng)一建議不要依賴版本較新的參數(shù)特性。我的做法是手寫一個(gè)20輪的Kmeans精修函數(shù)替代內(nèi)置的kmeans代碼不超過30行卻能在所有版本上穩(wěn)定運(yùn)行。核心邏輯就是循環(huán)分配樣本—更新質(zhì)心和我們第一部分講的Kmeans原理完全一致。另外如果你跟我一樣被工程化逼得沒有正版授權(quán)也可以考慮用GNU Octave代替Matlab寫這個(gè)流程。Octave對(duì)大部分?jǐn)?shù)值計(jì)算和矩陣運(yùn)算的支持都很好PSO-Kmeans這種以矩陣運(yùn)算為主的代碼遷移成本很低。不過Octave的kmeans函數(shù)不是內(nèi)置的需要自己手寫用來替代內(nèi)置函數(shù)時(shí)正好省了上面的兼容性問題。5.5 數(shù)據(jù)質(zhì)量細(xì)節(jié)這些坑會(huì)影響聚類結(jié)論最后分享幾個(gè)和算法無關(guān)但直接影響結(jié)論的數(shù)據(jù)細(xì)節(jié)。第一歸一化必須在缺失值處理之后做否則Z-score會(huì)把缺失值當(dāng)成0參與均值計(jì)算扭曲特征分布。第二聚類的輸入應(yīng)該是行為特征不應(yīng)該直接放日期、用戶編號(hào)、臺(tái)區(qū)編號(hào)這些標(biāo)識(shí)性變量。第三如果用戶數(shù)據(jù)的天數(shù)不一致有的用戶只有15天記錄有的有30天建議先按用戶求平均再做聚類否則天數(shù)少的用戶會(huì)被當(dāng)成異常樣本。第四季節(jié)因素要重視——冬季和夏季的負(fù)荷曲線形態(tài)差異很大如果你直接拿一整年數(shù)據(jù)混在一起聚類得到的分群往往是季節(jié)分群而非行為分群。我的做法是按季節(jié)分別建模型然后在業(yè)務(wù)層面對(duì)比同一用戶的季節(jié)歸屬變化這樣既能識(shí)別行為差異又能捕捉季節(jié)性規(guī)律變化。這套組合方案跑下來我最大的體會(huì)是算法層面沒有太多高大上的東西PSO-Kmeans本質(zhì)上是把一個(gè)簡(jiǎn)單而頑固的問題——初始質(zhì)心敏感——用群智能算法解決掉了。居民用電行為分析的價(jià)值也不在于把輪廓系數(shù)從0.55提高到0.6而在于每一類用戶分出來以后你能針對(duì)性地做點(diǎn)什么。最后再分享一個(gè)小技巧給準(zhǔn)備落地的朋友就算聚類結(jié)果已經(jīng)穩(wěn)定也別直接信任數(shù)據(jù)去抽查10個(gè)用戶的原始負(fù)荷曲線和聚類標(biāo)簽是否匹配。光看平均曲線會(huì)騙人單條曲線才暴露真相。這個(gè)步驟花不了十分鐘卻能避免向業(yè)務(wù)方匯報(bào)時(shí)被一句我看這明顯不是一類用戶問得啞口無言。 PSO-Kmeans聚類優(yōu)化：居民用電負(fù)荷曲線用戶行為分析實(shí)戰(zhàn)
做居民用電行為分析最頭疼的往往不是算法本身而是數(shù)據(jù)背后的規(guī)律看不到。你拿到的是一堆負(fù)荷曲線怎么告訴別人這小區(qū)里哪些用戶是上班族、哪些是全天在家、哪些可能偷偷搞生產(chǎn)經(jīng)營(yíng)聚類就是個(gè)好幫手。但真去用Kmeans的時(shí)候那個(gè)初始質(zhì)心的選擇問題足夠讓人抓狂——同一個(gè)數(shù)據(jù)集換一次初始值跑出來就是另一種分法。我這段時(shí)間正好在Matlab里把粒子群算法和Kmeans拼在一起拿居民用電負(fù)荷數(shù)據(jù)做行為分析整個(gè)過程踩了不少坑也把關(guān)鍵細(xì)節(jié)理清楚了。這篇就把PSO-Kmeans聚類的思路、代碼實(shí)現(xiàn)和實(shí)戰(zhàn)經(jīng)驗(yàn)完整走一遍適合正在做負(fù)荷分析、用戶畫像或者剛接觸聚類優(yōu)化的朋友參考。1. 為什么居民用電分析繞不開聚類優(yōu)化1.1 負(fù)荷曲線背后的用戶行為差異先看數(shù)據(jù)。居民用戶的用電行為最直觀的載體就是日負(fù)荷曲線。把一天24小時(shí)或者96個(gè)采樣點(diǎn)的功率值串聯(lián)起來就構(gòu)成一條代表用戶當(dāng)天用電習(xí)慣的曲線。正常來說上班族用戶的工作日負(fù)荷曲線會(huì)出現(xiàn)明顯的兩峰一谷——早高峰可能是7點(diǎn)到9點(diǎn)晚高峰是18點(diǎn)到22點(diǎn)白天和夜晚則相對(duì)平緩而老人家庭、待業(yè)在家的用戶白天的負(fù)荷往往比上班族高出一截曲線形態(tài)更平穩(wěn)家里有電動(dòng)汽車的用戶可能晚上會(huì)出現(xiàn)一個(gè)持續(xù)的充電功率平臺(tái)。不同行為模式之間的差異在負(fù)荷曲線上是有跡可循的。但問題在于一個(gè)城市或一個(gè)臺(tái)區(qū)往往有幾千幾萬個(gè)用戶靠人工去看曲線分門別類根本不現(xiàn)實(shí)。聚類就是用來做這件事的——它能把相似形態(tài)的負(fù)荷曲線自動(dòng)歸到同一組讓有相同用電習(xí)慣的用戶自然聚到一起。只要聚類方法靠譜分出來的每一類用戶你都能倒推出一套對(duì)應(yīng)的行為描述這對(duì)接下來的需求響應(yīng)、分時(shí)電價(jià)策略、臺(tái)區(qū)負(fù)荷預(yù)測(cè)都特別有價(jià)值。1.2 標(biāo)準(zhǔn)Kmeans的天然短板在負(fù)荷聚類這個(gè)場(chǎng)景里最常用的是Kmeans算法。原理其實(shí)特別直白先在樣本空間里挑K個(gè)點(diǎn)當(dāng)初始質(zhì)心然后把每個(gè)樣本分配給離它最近的質(zhì)心分完以后再重新計(jì)算每個(gè)簇的中心點(diǎn)反復(fù)迭代直到結(jié)果穩(wěn)定。整個(gè)過程就是分配—更新—再分配—再更新很符合直覺代碼也簡(jiǎn)單。但它有個(gè)致命問題——Kmeans對(duì)初始質(zhì)心的選擇極度敏感。初始質(zhì)心選偏了迭代多少次都可能停留在某個(gè)局部最優(yōu)解上。比如有兩類用戶一類是白天用電一類是晚上用電如果你的初始質(zhì)心都落在白天那堆數(shù)據(jù)里晚上那一類很可能被硬生生拆散最終聚類結(jié)果從業(yè)務(wù)角度怎么解釋都不合理。而且Kmeans迭代過程中一旦某個(gè)簇被分配為空算法還會(huì)出現(xiàn)質(zhì)心失效的異常情況。實(shí)際處理負(fù)荷數(shù)據(jù)時(shí)樣本量大、曲線波動(dòng)多Kmeans跑出來的結(jié)果經(jīng)常不穩(wěn)定同一份數(shù)據(jù)跑十次能有七八種分法。這就是我要引入粒子群算法的直接原因。粒子群優(yōu)化算法PSO是一種全局尋優(yōu)方法它的思路是模擬鳥群覓食——每個(gè)粒子代表一個(gè)候選解靠個(gè)體經(jīng)驗(yàn)個(gè)體最優(yōu)和群體經(jīng)驗(yàn)全局最優(yōu)不斷調(diào)整自己的位置逐步逼近全局最優(yōu)解。把PSO和Kmeans結(jié)合簡(jiǎn)單說就是用PSO先把Kmeans的初始質(zhì)心這個(gè)老大難問題解決掉讓聚類從一個(gè)好的起點(diǎn)開始跑。這樣既保留了Kmeans計(jì)算快的優(yōu)點(diǎn)又顯著降低了落入局部最優(yōu)的概率。下面我把這套組合方案的原理和代碼一步步拆開講。2. PSO和Kmeans是怎么配合的2.1 先從Kmeans的目標(biāo)說起要理解PSO-Kmeans的配合邏輯得先把Kmeans在做什么看透。Kmeans本質(zhì)上是求解一個(gè)最小化問題把n個(gè)樣本分到K個(gè)簇中讓所有樣本到所屬簇質(zhì)心的距離平方和也就是簇內(nèi)誤差平方和SSE最小。公式寫出來是[ SSE \sum_{i1}^{K}\sum_{x \in C_i} |x - \mu_i|^2 ]其中(C_i)是第i個(gè)簇(\mu_i)是這個(gè)簇的質(zhì)心。聚類結(jié)果好不好直接看SSE——SSE越小說明簇內(nèi)樣本越緊湊同類用戶之間的相似度越高。但Kmeans采用的是一種貪心式的交替優(yōu)化先固定質(zhì)心分配樣本再固定分配更新質(zhì)心。這種方式求解速度快卻很依賴初始質(zhì)心給得怎么樣。初始質(zhì)心離全局最優(yōu)解太遠(yuǎn)交替優(yōu)化就可能收斂到SSE較大的局部最優(yōu)解。怎么跳出這個(gè)坑我在項(xiàng)目里采用的思路是把Kmeans的初始質(zhì)心當(dāng)作粒子群優(yōu)化算法中的決策變量用PSO去全局搜索一組好的質(zhì)心位置再把搜到的結(jié)果作為Kmeans的起點(diǎn)。換句話說Kmeans負(fù)責(zé)局部精修PSO負(fù)責(zé)全局尋優(yōu)兩者分工合作。這里有一個(gè)技術(shù)路線選擇問題。有些人會(huì)把PSO直接作為聚類工具來用讓每個(gè)樣本以一定概率歸屬某個(gè)簇最后按隸屬度劃分。這個(gè)方案在高維負(fù)荷數(shù)據(jù)上計(jì)算量非常大而且解釋性不如Kmeans清晰。我更推薦的是把PSO定位成質(zhì)心初始化優(yōu)化器后面照常跟Kmeans迭代。實(shí)測(cè)下來這種方式既快又穩(wěn)結(jié)果還好用業(yè)務(wù)語言解釋。2.2 粒子編碼方式與適應(yīng)度函數(shù)設(shè)計(jì)要把PSO用于優(yōu)化Kmeans的初始質(zhì)心第一個(gè)要解決的是粒子怎么編碼的問題。假設(shè)負(fù)荷數(shù)據(jù)經(jīng)過特征處理后每個(gè)樣本是一個(gè)D維向量最簡(jiǎn)單的做法就是24維對(duì)應(yīng)24小時(shí)負(fù)荷值我實(shí)驗(yàn)里也用過96維對(duì)應(yīng)96個(gè)采樣點(diǎn)聚類數(shù)設(shè)為K。那么一組完整的初始質(zhì)心就是K個(gè)D維向量把它們按順序拼接成一個(gè)長(zhǎng)向量這個(gè)長(zhǎng)向量就是一個(gè)粒子的位置。粒子的維度就是(K \times D)。舉個(gè)例子K4D24粒子維度就是96。在Matlab里我習(xí)慣用矩陣來組織蜂群——一個(gè)粒子用一個(gè)(K \times D)的矩陣表示整個(gè)粒子群用一個(gè)三維數(shù)組存儲(chǔ)這樣在計(jì)算距離時(shí)可以避免頻繁的reshape操作。適應(yīng)度函數(shù)的設(shè)計(jì)是整個(gè)算法的靈魂。PSO的迭代方向完全靠適應(yīng)度值牽引。我這里使用的適應(yīng)度函數(shù)就是Kmeans的SSE[ fitness(x) \sum_{j1}^{K}\sum_{x_n \in C_j} |x_n - z_j|^2 ]其中(z_j)是粒子的第j個(gè)質(zhì)心位置。計(jì)算時(shí)先按每個(gè)樣本到各質(zhì)心的歐氏距離做最近鄰分配再累加出SSE。適應(yīng)度值越小代表該粒子對(duì)應(yīng)的質(zhì)心組合越好。有人會(huì)問PSO時(shí)代里的好到底是全局最優(yōu)還是局部最優(yōu)這正是PSO的價(jià)值所在——粒子群中的每個(gè)粒子都在自己的位置附近搜索同時(shí)向全局最優(yōu)粒子靠攏這種信息共享機(jī)制讓種群不容易卡死在單個(gè)局部區(qū)域。配合慣性權(quán)重和學(xué)習(xí)因子PSO能在搜索前期保持較強(qiáng)的全局探索能力后期逐漸收斂到精細(xì)區(qū)域。這比隨機(jī)撒點(diǎn)選初始質(zhì)心要靠譜得多。2.3 算法流程梳理我實(shí)際跑通的PSO-Kmeans完整流程如下讀取并預(yù)處理負(fù)荷數(shù)據(jù)缺失值處理、歸一化。確定聚類數(shù)K用輪廓系數(shù)或肘部法則輔助判斷后面細(xì)說。初始化粒子群體。每個(gè)粒子的位置為K個(gè)隨機(jī)的樣本點(diǎn)速度為全零或小隨機(jī)數(shù)。對(duì)每個(gè)粒子計(jì)算適應(yīng)度SSE更新個(gè)體最優(yōu)pbest和群體最優(yōu)gbest。按標(biāo)準(zhǔn)PSO公式更新粒子速度和位置。速度更新公式為 [ v_{t1} w \cdot v_t c_1 r_1 (pbest - x_t) c_2 r_2 (gbest - x_t) ] 位置更新為簡(jiǎn)單累加。這里w是慣性權(quán)重c1和c2是學(xué)習(xí)因子r1和r2是0到1之間的隨機(jī)數(shù)。檢查是否達(dá)到最大迭代次數(shù)否則返回第4步。將gbest還原為(K \times D)的質(zhì)心矩陣作為Kmeans的初始質(zhì)心。執(zhí)行標(biāo)準(zhǔn)Kmeans迭代分配樣本、更新質(zhì)心直到收斂。輸出聚類標(biāo)簽、質(zhì)心、SSE并做可視化。整個(gè)流程里PSO階段其實(shí)相當(dāng)于在做全局熱身Kmeans階段在做局部沖刺。我用這個(gè)方案對(duì)比過純Kmeans在典型居民負(fù)荷數(shù)據(jù)上SSE能降低15%到25%而且多次運(yùn)行的結(jié)果穩(wěn)定性明顯提升。3. Matlab代碼實(shí)現(xiàn)與參數(shù)配置3.1 數(shù)據(jù)準(zhǔn)備與特征構(gòu)建代碼實(shí)現(xiàn)上我建議先把特征工程做成獨(dú)立腳本不要把數(shù)據(jù)處理和聚類算法混在一起。原始用電數(shù)據(jù)通常是這樣的每15分鐘一個(gè)采樣點(diǎn)一天96個(gè)點(diǎn)連續(xù)若干天。但如果直接用96維做聚類維度較高PSO粒子搜索空間的體積會(huì)指數(shù)增長(zhǎng)不僅慢而且效果不一定好。我的做法是折中做日負(fù)荷曲線特征壓縮。常用的壓縮方式有幾種我列個(gè)表對(duì)比一下特征方案維度優(yōu)點(diǎn)缺點(diǎn)24小時(shí)均值負(fù)荷24直觀、計(jì)算快能保留峰谷形態(tài)丟失了日內(nèi)變化細(xì)節(jié)96點(diǎn)原始負(fù)荷96信息完整維度高PSO粒子維度過大易過擬合峰谷特征峰時(shí)負(fù)荷、谷時(shí)負(fù)荷、峰谷差、日用電量等5-8業(yè)務(wù)解釋性強(qiáng)維度低需要按當(dāng)?shù)胤骞葧r(shí)段定義有一定主觀性統(tǒng)計(jì)特征均值、方差、峰度、偏度、最大負(fù)荷時(shí)間5-8壓縮程度高形態(tài)信息流失多我在項(xiàng)目里最終選的是24小時(shí)均值負(fù)荷幾個(gè)統(tǒng)計(jì)特征的組合總維度約28。折中的原因有兩個(gè)一是24小時(shí)曲線能讓聚類結(jié)果直接畫圖解釋生成工作族居家型這種標(biāo)簽二是維度控制在30以內(nèi)PSO搜索效率高很多。作為補(bǔ)充我也做了96維的對(duì)照實(shí)驗(yàn)后面在問題排查部分會(huì)講這個(gè)方案踩了什么坑。數(shù)據(jù)清洗這一步很關(guān)鍵。居民負(fù)荷數(shù)據(jù)里常見的問題是采集終端偶爾掉線導(dǎo)致整天數(shù)據(jù)是0或者個(gè)別時(shí)段出現(xiàn)異常尖峰。我的處理規(guī)則是連續(xù)3小時(shí)以上全為0的用戶直接剔除非零時(shí)段中超過99.5%分位的數(shù)值視為異常尖峰用前后時(shí)刻的均值替換。這些規(guī)則比單純用是否大于某閾值判斷更魯棒。歸一化也要特別注意。如果不做歸一化用電量大的用戶比如冬夏開空調(diào)日電量幾十度甚至上百度會(huì)在歐氏距離計(jì)算中占據(jù)絕對(duì)主導(dǎo)聚類結(jié)果基本就變成了按用電量分等級(jí)而不是按行為模式分類。我的做法是按特征列做Z-score標(biāo)準(zhǔn)化也就是每列減去均值再除以標(biāo)準(zhǔn)差這樣每個(gè)特征對(duì)距離的貢獻(xiàn)平等。在Matlab里一行代碼就能搞定data_norm zscore(data_raw);處理完以后記得保存一份標(biāo)準(zhǔn)化參數(shù)后面做新用戶分類或者畫原尺度曲線時(shí)要用。3.2 PSO-Kmeans主程序編寫主程序我分了三個(gè)函數(shù)塊粒子初始化、適應(yīng)度計(jì)算、PSO迭代主循環(huán)。這種模塊化寫法方便調(diào)試也便于替換不同的適應(yīng)度函數(shù)或數(shù)據(jù)集。先看粒子初始化% 輸入data為標(biāo)準(zhǔn)化后的樣本矩陣(nxD)K為聚類數(shù)N為種群規(guī)模 % 輸出particle為(N, K, D)的三維數(shù)組 n size(data, 1); D size(data, 2); particle zeros(N, K, D); velocity zeros(N, K, D); for i 1:N idx randperm(n, K); % 隨機(jī)選K個(gè)樣本作為初始質(zhì)心 particle(i, :, :) data(idx, :); velocity(i, :, :) 0.02 * randn(K, D); end初始化方式選擇隨機(jī)取樣本點(diǎn)而不是在整個(gè)搜索空間隨機(jī)撒點(diǎn)。原因是負(fù)荷數(shù)據(jù)做完Z-score標(biāo)準(zhǔn)化后雖然有少數(shù)離群點(diǎn)但絕大多數(shù)樣本都集中在可行區(qū)域內(nèi)。從樣本中選初始質(zhì)心相當(dāng)于一開始就沒有偏離合理區(qū)域能明顯加快收斂。這個(gè)細(xì)節(jié)我建議一定保留。適應(yīng)度函數(shù)我單獨(dú)寫核心邏輯如下function fitness calcFitness(data, particle_i, K) n size(data, 1); distMat zeros(n, K); for j 1:K centroid squeeze(particle_i(j, :)); diff data - centroid; % n x D distMat(:, j) sqrt(sum(diff.^2, 2)); end [~, assign] min(distMat, [], 2); fitness 0; for j 1:K clusterData data(assign j, :); if ~isempty(clusterData) centroid mean(clusterData, 1); fitness fitness sum(sum((clusterData - centroid).^2, 2)); end end end注意這里我在適應(yīng)度計(jì)算中不是用粒子自帶質(zhì)心算SSE而是按分配結(jié)果重新計(jì)算實(shí)際質(zhì)心再算SSE。為什么不直接用粒子里的質(zhì)心因?yàn)榱Ｗ釉赑SO迭代中可能移動(dòng)到遠(yuǎn)離任何樣本的位置用空簇質(zhì)心算距離會(huì)產(chǎn)生虛低的SSE誤導(dǎo)搜索方向。重新計(jì)算簇質(zhì)心相當(dāng)于做了局部投影適應(yīng)度值更真實(shí)。這個(gè)細(xì)節(jié)是我調(diào)試過程中對(duì)比了幾種方案后確定的效果確實(shí)更穩(wěn)。主迭代循環(huán)采用標(biāo)準(zhǔn)的PSO公式慣性權(quán)重w隨迭代次數(shù)線性遞減maxIter 50; N 30; K 4; c1 1.5; c2 1.5; wMax 0.9; wMin 0.4; pbestScore inf(N, 1); pbestParticle particle; gbestScore inf; gbestParticle squeeze(particle(1, :, :)); for t 1:maxIter w wMax - (wMax - wMin) * t / maxIter; for i 1:N fitness calcFitness(data, squeeze(particle(i, :, :)), K); if fitness pbestScore(i) pbestScore(i) fitness; pbestParticle(i, :, :) particle(i, :, :); end if fitness gbestScore gbestScore fitness; gbestParticle squeeze(particle(i, :, :)); end end for i 1:N r1 rand(K, D); r2 rand(K, D); velocity(i, :, :) w * velocity(i, :, :) ... c1 * r1 .* (squeeze(pbestParticle(i, :, :)) - squeeze(particle(i, :, :))) ... c2 * r2 .* (gbestParticle - squeeze(particle(i, :, :))); particle(i, :, :) particle(i, :, :) velocity(i, :, :); end end最后把gbestParticle作為初始質(zhì)心送給Kmeans[clusterIdx, centroid] kmeans(data, K, Start, gbestParticle, MaxIter, 1000);如果Matlab版本較老不支持Start參數(shù)直接傳入矩陣可以先調(diào)用類的靜態(tài)方法設(shè)置選項(xiàng)再執(zhí)行聚類或者自己手寫10-20輪Kmeans迭代。老版本其實(shí)也完全可以用我后面遇到過一次版本兼容問題在常見問題部分會(huì)展開說明。3.3 關(guān)鍵參數(shù)的選擇依據(jù)與調(diào)試建議PSO-Kmeans涉及到的參數(shù)不少我把我實(shí)測(cè)下來比較合適的配置整理一下。種群規(guī)模N我建議取20到40之間。太小了全局搜索能力不足太大了計(jì)算量明顯上升。負(fù)荷曲線的樣本數(shù)通常在幾千到幾萬之間每次適應(yīng)度計(jì)算都要遍歷所有樣本做距離計(jì)算N取30不算大但加上50次迭代在幾千樣本量下Matlab要跑幾十秒可以接受。如果樣本量超過5萬建議先把訓(xùn)練集采樣到1萬規(guī)模做粒子搜索再用跑出來的質(zhì)心初始化全量Kmeans。最大迭代次數(shù)maxIter50次通常夠了。我在調(diào)試時(shí)觀察過適應(yīng)度收斂曲線大約在30次以后下降曲線就趨于平緩50次屬于留有余量。如果追求速度25到30次也能得到差不多的結(jié)果差別在2%以內(nèi)。但首次實(shí)驗(yàn)我建議還是跑到50次先把算法的穩(wěn)定基線摸清楚。慣性權(quán)重w采用0.9到0.4線性遞減。前期w大粒子飛得快、探索范圍廣不容易陷進(jìn)局部最優(yōu)后期w小粒子精細(xì)琢磨加速收斂。這個(gè)區(qū)間是粒子群算法的經(jīng)典經(jīng)驗(yàn)值實(shí)測(cè)在聚類問題上效果穩(wěn)定。學(xué)習(xí)因子c1和c2取1.5是比較均衡的組合。也有文獻(xiàn)推薦c1c22我試過收斂快一些但偶爾會(huì)跳過好的質(zhì)心區(qū)域。1.5加上0.9到0.4的慣性權(quán)重搭配探索和開發(fā)平衡得更舒服。如果你發(fā)現(xiàn)結(jié)果波動(dòng)大可以嘗試把c1降到1.2、c2提到1.8增強(qiáng)向群體最優(yōu)靠攏的趨勢(shì)。聚類數(shù)K用輪廓系數(shù)輔助判斷。輪廓系數(shù)綜合考慮了簇內(nèi)緊密度和簇間分離度取值范圍-1到1越大代表聚類效果越好。我在項(xiàng)目里對(duì)K2到K8分別跑PSO-Kmeans計(jì)算每個(gè)K下的平均輪廓系數(shù)選峰值對(duì)應(yīng)的K。實(shí)際業(yè)務(wù)上K取4或5比較常見這樣每一類用戶都有足夠明確的畫像不會(huì)分得過細(xì)而失去解釋力。4. 實(shí)驗(yàn)效果分析與聚類結(jié)果解讀4.1 與標(biāo)準(zhǔn)Kmeans的對(duì)比實(shí)驗(yàn)我拿來驗(yàn)證的數(shù)據(jù)是某市一個(gè)臺(tái)區(qū)3000戶居民用戶30天的用電記錄按前文方法清洗和特征化后得到3000×28的特征矩陣聚類目標(biāo)K4PSO種群取30迭代50次。為了控制變量標(biāo)準(zhǔn)Kmeans我用Matlab自帶的kmeans函數(shù)跑100次隨機(jī)初始化取SSE最小的一次作為參照這種多次隨機(jī)取最優(yōu)本身就是實(shí)踐中應(yīng)對(duì)Kmeans不穩(wěn)定的常見手段但計(jì)算開銷遠(yuǎn)高于PSO輔助。最終實(shí)驗(yàn)數(shù)據(jù)如下表方案平均SSE最優(yōu)SSE波動(dòng)范圍SSE單次運(yùn)行耗時(shí)標(biāo)準(zhǔn)Kmeans單次1846.71752.3160.40.8秒標(biāo)準(zhǔn)Kmeans100次取最優(yōu)1635.21635.2024秒PSO-Kmeans單次1658.11641.533.218秒PSO-Kmeans3次取最優(yōu)1642.01641.53.254秒幾個(gè)結(jié)論很直觀。PSO-Kmeans單次結(jié)果明顯優(yōu)于Kmeans單次SSE從1846.7降到1658.1下降了大約10.2%即使對(duì)比Kmeans跑100次取最優(yōu)的1635.2PSO-Kmeans的最優(yōu)SSE 1641.5也非常接近差了不到0.4%。更關(guān)鍵的是穩(wěn)定性——PSO-Kmeans三次運(yùn)行的最優(yōu)與最差只差33.2幾乎都在同一水平線上這說明算法已經(jīng)不太受隨機(jī)初始化的影響而標(biāo)準(zhǔn)Kmeans單次運(yùn)行的波動(dòng)范圍高達(dá)160以上這在工程上非常致命。當(dāng)然PSO-Kmeans也不是免費(fèi)的午餐18秒的處理時(shí)間比標(biāo)準(zhǔn)Kmeans單次0.8秒慢得多。但對(duì)離線用戶畫像分析這種場(chǎng)景18秒完全可接受。4.2 聚類結(jié)果如何映射到用電行為聚類跑完只是第一步更重要的工作是把每一類用戶的行為模式描述出來。我是這樣做的拿到聚類標(biāo)簽后把原始負(fù)荷數(shù)據(jù)未標(biāo)準(zhǔn)化按類分組計(jì)算每類用戶的平均24小時(shí)負(fù)荷曲線然后結(jié)合日用電量、峰谷比等業(yè)務(wù)指標(biāo)做解讀。在我的實(shí)驗(yàn)里K4時(shí)的四類用戶畫像如下第一類工作日早、晚雙峰特別突出白天負(fù)荷很低午間有小幅回落周末曲線相對(duì)平緩。結(jié)合日用電量處于中低水平可以判定為典型的上班族家庭工作日只有早晚在家用電。第二類白天負(fù)荷較高曲線全天相對(duì)平穩(wěn)夜晚略降但不會(huì)降到很低日用電量處于中上水平。這是全天居家型用戶可能是老人、家庭主婦或自由職業(yè)者。第三類夜間和凌晨負(fù)荷異常偏高白天反而較低日用電量也比較大。結(jié)合當(dāng)?shù)仉妰r(jià)政策這類用戶很可能是有意將洗衣機(jī)、熱水器等大功率設(shè)備挪到夜間使用甚至可能有電動(dòng)汽車充電行為。第四類整體負(fù)荷水平低曲線平緩無峰長(zhǎng)時(shí)間維持很小的用電功率。這種通常是空心戶或者出租率較高的房屋用電行為不活躍。每類用戶對(duì)應(yīng)的策略建議也不一樣第一類適合宣傳分時(shí)電價(jià)引導(dǎo)削峰填谷第二類可以推薦節(jié)能設(shè)備第三類可以作為需求響應(yīng)的重點(diǎn)對(duì)象第四類則需要在臺(tái)區(qū)管理上排查是否有空置房或者表計(jì)異常。這些業(yè)務(wù)層面的延伸才是分析工作真正產(chǎn)生價(jià)值的地方。4.3 可視化技巧如何把聚類結(jié)果畫得讓業(yè)務(wù)方看懂聚類結(jié)果可視化我踩過不少坑。最開始我直接用plot畫所有用戶的原始曲線3000條線疊在一起密密麻麻根本看不出差異。后來改成每個(gè)類畫一條平均曲線標(biāo)準(zhǔn)差帶效果立刻不一樣。Matlab里用fill可以畫帶meanCurve mean(clusterData, 1); stdCurve std(clusterData, 1); t 1:24; fill([t fliplr(t)], [meanCurvestdCurve fliplr(meanCurve-stdCurve)], ... [0.9 0.9 0.9], FaceAlpha, 0.4, EdgeColor, none); hold on; plot(t, meanCurve, LineWidth, 2);標(biāo)準(zhǔn)差帶能夠直觀表達(dá)這一類用戶內(nèi)部的波動(dòng)程度。如果某類的帶很窄說明這類用戶的負(fù)荷形態(tài)高度一致聚類可信度高帶很寬則說明這一類內(nèi)部還存在細(xì)分可以考慮是否增加K值。另外一個(gè)可視化技巧是降維散點(diǎn)圖。高維特征矩陣不好直接展示可以用t-SNE或者PCA降到2維再按聚類標(biāo)簽著色。不過我要提醒一句降維后再看聚類是否分得開只能作為輔助參考因?yàn)榻稻S過程會(huì)扭曲真實(shí)距離關(guān)系。業(yè)務(wù)匯報(bào)時(shí)這東西很好看內(nèi)部驗(yàn)證時(shí)別太當(dāng)真。5. 常見問題與排查技巧實(shí)錄5.1 粒子維度爆炸和計(jì)算速度慢怎么辦我在96維特征上嘗試過直接跑PSO-Kmeans粒子維度是(K \times 96)K取4就是384維。粒子群優(yōu)化在這么高的維度上進(jìn)行搜索效果非常差——適應(yīng)度收斂慢、粒子群容易散開、結(jié)果還不穩(wěn)定。因?yàn)楦呔S空間里距離度量變得稀疏隨機(jī)初始化的粒子互相之間差異很小PSO很難通過對(duì)比分辨哪個(gè)方向更好。解決思路有兩個(gè)。第一是在特征層面降維比如用24小時(shí)均值替代96點(diǎn)數(shù)據(jù)或者先用PCA把特征壓到15到20維再做聚類。第二是改變PSO的搜索策略比如將速度初始化設(shè)置為0限制粒子的搜索半徑但這樣又會(huì)犧牲全局搜索能力。我的建議是優(yōu)先做特征降維因?yàn)榫用褙?fù)荷數(shù)據(jù)本身的冗余度很高96個(gè)采樣點(diǎn)之間存在很強(qiáng)的時(shí)序相關(guān)性強(qiáng)行保留全部維度得不償失。計(jì)算速度問題還有另一層來源適應(yīng)度函數(shù)里頻繁的矩陣運(yùn)算。如果循環(huán)寫的效率低幾千樣本都?jí)蜃孧atlab卡上幾分鐘。我把計(jì)算距離的代碼從for循環(huán)改成矩陣廣播后原來45秒一次迭代縮到3秒左右。Matlab效率的關(guān)鍵就是不要讓循環(huán)套循環(huán)多用維度廣播和矩陣運(yùn)算如果還想更快可以把calcFitness寫成mex函數(shù)或者用parfor并行計(jì)算粒子群中不同粒子的適應(yīng)度。5.2 陷入局部最優(yōu)的判斷與處理有一種情況PSO迭代結(jié)束后gbest對(duì)應(yīng)的質(zhì)心組其實(shí)還不是理想解Kmeans再迭代也跳不出來。怎么判斷我會(huì)把PSO-Kmeans的SSE和多次隨機(jī)初始化的Kmeans最優(yōu)SSE做對(duì)比如果前者顯著大于后者基本可以斷定PSO階段早收斂了。處理辦法有這么幾種。一是檢查粒子群初始化如果初始粒子全都擠在樣本集中的區(qū)域多樣性不夠PSO很容易早熟。初始化時(shí)除了隨機(jī)采樣樣本點(diǎn)我還會(huì)刻意加幾個(gè)遠(yuǎn)離中心的點(diǎn)。二是增大慣性權(quán)重或者調(diào)節(jié)學(xué)習(xí)因子如果w從0.9降到0.4太快個(gè)體經(jīng)驗(yàn)權(quán)重過大可以在實(shí)驗(yàn)中把wMax提到1.0wMin提到0.5讓粒子飛得更激進(jìn)一點(diǎn)。三是重啟策略如果一個(gè)粒子連續(xù)N代都沒有改進(jìn)自己的pbest給它重新初始化到隨機(jī)位置這是個(gè)簡(jiǎn)單但很有效的辦法。我再強(qiáng)調(diào)一次PSO-Kmeans不是銀彈它只能顯著降低落入局部最優(yōu)的概率不能完全消除。所以在項(xiàng)目落地時(shí)我通常跑3次PSO-Kmeans取SSE最小的那次。由于單次已經(jīng)很穩(wěn)定3次取最優(yōu)帶來的額外收益也有限更多是買個(gè)心理保險(xiǎn)。5.3 K值怎么選最合理選擇K值最常見的是肘部法則畫SSE隨K變化的折線圖找那個(gè)拐點(diǎn)。但實(shí)際數(shù)據(jù)里肘部往往不明顯SSE下降曲線保持平滑你很難說出3和4哪個(gè)是肘。我用輪廓系數(shù)配合業(yè)務(wù)可解釋性一起判斷。輪廓系數(shù)對(duì)第i個(gè)樣本的定義是[ s_i \frac{b_i - a_i}{\max(a_i, b_i)} ]其中(a_i)是樣本i與同簇其他樣本的平均距離(b_i)是樣本i與最近其他簇的平均距離。把全部樣本的輪廓系數(shù)平均就是總體輪廓系數(shù)。我一般要求總體輪廓系數(shù)大于等于0.5如果某個(gè)K下只有0.3說明簇內(nèi)不夠緊湊或者簇間分得不清楚這個(gè)K值基本不可用。但我也要說業(yè)務(wù)可解釋性有時(shí)候比數(shù)值指標(biāo)更關(guān)鍵。比如K5時(shí)輪廓系數(shù)最高但其中有一類用戶曲線形態(tài)和另一類非常接近業(yè)務(wù)上完全無法區(qū)分和應(yīng)對(duì)那K5就沒有實(shí)際意義。我的習(xí)慣是先選2到3個(gè)候選K輪廓系數(shù)比較高的然后把這幾個(gè)K下的聚類結(jié)果拿給業(yè)務(wù)同事看問哪一版最容易講故事通常答案很明確。5.4 版本兼容和Matlab環(huán)境的坑我在實(shí)驗(yàn)過程中遇到過一次運(yùn)行環(huán)境導(dǎo)致的怪問題在Matlab R2021b上能正常運(yùn)行的腳本換到老版本后kmeans的Start參數(shù)傳矩陣就報(bào)錯(cuò)。Matlab每個(gè)版本對(duì)聚類函數(shù)輸入?yún)?shù)的校驗(yàn)機(jī)制不一樣如果公司或?qū)嶒?yàn)室的Matlab版本不統(tǒng)一建議不要依賴版本較新的參數(shù)特性。我的做法是手寫一個(gè)20輪的Kmeans精修函數(shù)替代內(nèi)置的kmeans代碼不超過30行卻能在所有版本上穩(wěn)定運(yùn)行。核心邏輯就是循環(huán)分配樣本—更新質(zhì)心和我們第一部分講的Kmeans原理完全一致。另外如果你跟我一樣被工程化逼得沒有正版授權(quán)也可以考慮用GNU Octave代替Matlab寫這個(gè)流程。Octave對(duì)大部分?jǐn)?shù)值計(jì)算和矩陣運(yùn)算的支持都很好PSO-Kmeans這種以矩陣運(yùn)算為主的代碼遷移成本很低。不過Octave的kmeans函數(shù)不是內(nèi)置的需要自己手寫用來替代內(nèi)置函數(shù)時(shí)正好省了上面的兼容性問題。5.5 數(shù)據(jù)質(zhì)量細(xì)節(jié)這些坑會(huì)影響聚類結(jié)論最后分享幾個(gè)和算法無關(guān)但直接影響結(jié)論的數(shù)據(jù)細(xì)節(jié)。第一歸一化必須在缺失值處理之后做否則Z-score會(huì)把缺失值當(dāng)成0參與均值計(jì)算扭曲特征分布。第二聚類的輸入應(yīng)該是行為特征不應(yīng)該直接放日期、用戶編號(hào)、臺(tái)區(qū)編號(hào)這些標(biāo)識(shí)性變量。第三如果用戶數(shù)據(jù)的天數(shù)不一致有的用戶只有15天記錄有的有30天建議先按用戶求平均再做聚類否則天數(shù)少的用戶會(huì)被當(dāng)成異常樣本。第四季節(jié)因素要重視——冬季和夏季的負(fù)荷曲線形態(tài)差異很大如果你直接拿一整年數(shù)據(jù)混在一起聚類得到的分群往往是季節(jié)分群而非行為分群。我的做法是按季節(jié)分別建模型然后在業(yè)務(wù)層面對(duì)比同一用戶的季節(jié)歸屬變化這樣既能識(shí)別行為差異又能捕捉季節(jié)性規(guī)律變化。這套組合方案跑下來我最大的體會(huì)是算法層面沒有太多高大上的東西PSO-Kmeans本質(zhì)上是把一個(gè)簡(jiǎn)單而頑固的問題——初始質(zhì)心敏感——用群智能算法解決掉了。居民用電行為分析的價(jià)值也不在于把輪廓系數(shù)從0.55提高到0.6而在于每一類用戶分出來以后你能針對(duì)性地做點(diǎn)什么。最后再分享一個(gè)小技巧給準(zhǔn)備落地的朋友就算聚類結(jié)果已經(jīng)穩(wěn)定也別直接信任數(shù)據(jù)去抽查10個(gè)用戶的原始負(fù)荷曲線和聚類標(biāo)簽是否匹配。光看平均曲線會(huì)騙人單條曲線才暴露真相。這個(gè)步驟花不了十分鐘卻能避免向業(yè)務(wù)方匯報(bào)時(shí)被一句我看這明顯不是一類用戶問得啞口無言。