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