化Kmeans的居民用電行為聚類Matlab實現(xiàn))
最近在做居民用電行為分析時被手頭一批智能電表小時級負荷數(shù)據(jù)折騰得不輕。事情本身不復雜——把幾萬個用戶的用電曲線分類畫像但直接用Matlab自帶的kmeans跑結(jié)果實在沒法用同樣一份數(shù)據(jù)十次隨機初始化跑出十種結(jié)果SSE波動幅度接近百分之十幾。這個穩(wěn)定性拿到科研報告里能挑好看的寫但放在實際業(yè)務里根本說不過去。后來把粒子群算法PSO和Kmeans綁在一起讓PSO先全局搜索出一組好的初始聚類中心再把這個中心作為Kmeans的起點做局部收斂。實驗效果不錯但把整個流程走通的過程中踩了不少坑。這篇文章就把我的完整思路、Matlab代碼實現(xiàn)和排查教訓整理出來給正在做居民用電聚類、PSO優(yōu)化算法或者Kmeans相關課題的朋友一個可復用的參考。1. 為什么居民用電行為分析要動用到聚類算法智能電表普及之后電網(wǎng)公司手里掌握的用戶負荷數(shù)據(jù)量級已經(jīng)不是人工能處理的了。一個中型城市一天的采集數(shù)據(jù)就是幾十萬條曲線每條曲線包含24個甚至96個時點。面對這個規(guī)模的數(shù)據(jù)想要回答“哪些用戶適合執(zhí)行峰谷分時電價”“哪些用戶具備可調(diào)負荷參與需求響應”這類問題第一步都是做用戶分類。分類的方式有兩種一種是根據(jù)用戶臺賬手動打標簽比如小區(qū)類型、戶型面積、是否安裝分布式光伏但這種標簽往往滯后而且不完整另一種就是從負荷數(shù)據(jù)本身出發(fā)用無監(jiān)督聚類把行為模式相似的用戶聚到一組這就是用電行為分析的起點。1.1 用電行為聚類能解決什么實際問題居民用電行為聚類不是學術上的自嗨業(yè)務價值可以落到幾個非常具體的方向。第一個是分時電價套餐設計。不同用戶對電價的敏感度差異很大上班族白天空置家庭負荷低晚上回家負荷陡增退休老人家庭白天持續(xù)低負荷運行部分用戶夏季空調(diào)用電占比極高。這些特征通過聚類分組之后就可以針對性地設計尖峰費率、低谷優(yōu)惠、季節(jié)性套餐而不是對所有用戶一刀切。第二個是需求響應潛力評估?，F(xiàn)在的削峰填谷任務越來越重電網(wǎng)需要知道哪些用戶群體在有激勵時能主動壓低負荷。聚類之后如果某個簇的負荷曲線有明顯的用電集中時段、且這些時段恰好與電網(wǎng)峰時段重疊那么這個群體就是需求響應資源池中的優(yōu)質(zhì)對象。第三個是異常用電檢測和客戶細分服務。聚類本質(zhì)上刻畫的是“大多數(shù)用戶”的用電規(guī)律那些偏離規(guī)律的樣本自然成為重點關注對象。服務方面聚類畫像可以作為精準營銷的基礎比如識別出有儲能設備、有電動車充電習慣的用戶群體針對性推薦增值服務?？梢哉f聚類是連接原始負荷數(shù)據(jù)和上層業(yè)務決策的中間層。沒有這個中間層機器學習、數(shù)據(jù)挖掘的手段再多也無法落到電力系統(tǒng)的實際管理動作上。1.2 為什么選Kmeans而不是其他聚類算法用電數(shù)據(jù)分析領域常見的聚類算法不少DBSCAN、層次聚類、譜聚類都有人用但我最終把Kmeans作為基礎方案原因有三個。第一是數(shù)據(jù)規(guī)模。Kmeans的時間復雜度接近線性對幾十萬條樣本依然可以在可接受時間內(nèi)跑完。DBSCAN在實現(xiàn)較好的情況下效率也不錯但在簇密度不均勻的數(shù)據(jù)上參數(shù)很難調(diào)而居民負荷曲線恰好是密度分布極不均勻——白天用電峰值和夜間低谷對應的樣本分布差異非常大。第二是可解釋性。Kmeans的聚類中心就是“平均用電曲線”每個簇可以用一條有物理意義的負荷曲線來描述。層次聚類也能做到但樹狀圖在大樣本下可讀性極差譜聚類的嵌入空間不好向業(yè)務方解釋。第三是工程實現(xiàn)成熟度。Matlab的kmeans函數(shù)集成了多種距離度量、并行計算和重復啟動機制可以直接復用。但Matlab自帶的Kmeans有一個老問題——初始中心隨機選陷入局部最優(yōu)是家常便飯。這為后面引入PSO優(yōu)化提供了充分的動機。2. PSO如何給Kmeans找到更好的初始質(zhì)心Kmeans聚類的核心目標是找到一組聚類中心使所有樣本到所屬中心距離的平方和最小。這個目標函數(shù)非凸自帶Kmeans用的是隨機初始化和多次重復啟動本質(zhì)上是在“碰運氣”。粒子群算法在這里的角色就是替代隨機初始化通過群體智能搜索在解空間中找到一個更接近全局最優(yōu)的起始點。2.1 Kmeans的痛點初始質(zhì)心對結(jié)果的影響遠比你想象中大Kmeans迭代過程可以分成交替的兩步第一步根據(jù)當前聚類中心把每個樣本分配到最近的中心第二步根據(jù)分配結(jié)果重新計算每個簇的均值作為新中心。問題在于這個迭代過程是“局部爬山”。如果初始質(zhì)心選得不好比如兩個初始中心落在同一個數(shù)據(jù)密集區(qū)域那么最終聚類結(jié)果就會在局部最優(yōu)處停下來無法跳到更合理的劃分。真實數(shù)據(jù)里的表現(xiàn)就是隨便跑幾次kmeans每次的SSE都不一樣輪廓系數(shù)也有明顯起伏。我遇到的最夸張的一次同一份數(shù)據(jù)K4SSE最小值與最大值相差22%。這種情況下你沒法判斷當前這個結(jié)果是不是可用的。解決初始質(zhì)心問題有幾個常規(guī)手段。Kmeans的思路是讓初始質(zhì)心盡量分散這個策略對凸簇或者分布比較規(guī)則的簇有效但如果簇形狀重疊、存在噪聲Kmeans的初始中心仍然可能落在不合適的位置多次隨機重啟結(jié)合挑選最小SSE的做法本質(zhì)上是用計算時間去換概率但重復次數(shù)再多也還是隨機搜索沒有目標導向。粒子群算法不一樣它是帶著“方向”去找解的——每個候選解都會計算適應度通過種群協(xié)作逐步逼近更優(yōu)區(qū)域這正是Kmeans兩點方案中最需要的稀缺能力。2.2 粒子群算法的主要思想與更新過程粒子群算法的靈感來自鳥群覓食。每只鳥就是一個粒子代表解空間中的一個候選解鳥群的飛行過程就是不斷搜索最優(yōu)解的過程。每個粒子具備兩個屬性位置和速度。如果把當前解空間記作D維空間第i個粒子的位置記為X_i速度記為V_i。算法為每個粒子保存兩個記憶粒子自己的歷史最優(yōu)位置pbest以及整個種群目前發(fā)現(xiàn)的最優(yōu)位置gbest。每次迭代時粒子的速度和位置按照下面的方式更新速度更新V_i^{t1} wV_i^t c1r1*(pbest_i - X_i^t) c2r2(gbest - X_i^t)位置更新X_i^{t1} X_i^t V_i^{t1}其中w是慣性權重決定當前速度的保持程度c1是認知學習因子c2是社會學習因子r1和r2是[0,1]之間的隨機數(shù)。三個分量的物理意義很清楚慣性項讓粒子沿當前趨勢繼續(xù)飛行認知項把粒子拉向自己曾經(jīng)找到過的好位置社會項把整個種群引向全局最優(yōu)位置。這種機制在連續(xù)域搜索問題中表現(xiàn)穩(wěn)定參數(shù)少、實現(xiàn)簡單、不需要梯度信息非常適合用來優(yōu)化Kmeans的初始質(zhì)心。2.3 粒子編碼方式與適應度函數(shù)設計用PSO優(yōu)化Kmeans初始質(zhì)心最關鍵的環(huán)節(jié)是把Kmeans的解映射成PSO粒子。我采用的編碼方式是如果聚成K個簇每個簇中心是d維特征向量那么一個粒子的位置就是一個長度為 K*d 的一維向量前d個元素是第1個簇中心的坐標接下來d個元素是第2個簇中心的坐標依次類推。適應度函數(shù)的設計直接影響優(yōu)化效果。我用的適應度是樣本到最近聚類中心的歐氏距離平方和即Kmeans的SSE目標函數(shù)。這個選擇的理由很直接——PSO搜索到的適應度越小就意味著這組初始質(zhì)心本身已經(jīng)處在一個比較好的聚類劃分附近。為什么不直接用輪廓系數(shù)作為適應度因為輪廓系數(shù)要計算所有樣本兩兩之間的距離關系復雜度接近O(N2)幾萬條樣本時每次適應度評估都慢得難以接受而SSE的計算是線性的在PSO上百次迭代中每次都算輪廓系數(shù)實驗基本就等死了。還有一個常見做法值得提一下PSO每評估一次粒子時讓這個粒子先跑幾步Kmeans再計算適應度相當于是局部搜索和全局搜索的混合。這個方案效果好但時間復雜度翻了幾倍我在實際項目中用的是樸素的兩階段法——PSO只負責找初始質(zhì)心找到之后交給Kmeans做最終收斂。這樣工程上簡單效果已經(jīng)足夠替代隨機初始化了。3. Matlab實現(xiàn)從數(shù)據(jù)清洗到結(jié)果評估整個實驗我用的Matlab版本是R2024b代碼實現(xiàn)層面沒有用任何工具箱之外的特殊功能換到R2020之后的老版本也能跑。實現(xiàn)過程分為數(shù)據(jù)清洗、特征工程、PSO核心代碼、參數(shù)配置和結(jié)果評估幾個環(huán)節(jié)。3.1 數(shù)據(jù)準備與特征工程我手上的原始數(shù)據(jù)是一個居民小區(qū)的智能電表負荷記錄每一行是一個用戶每一小時一個采集點每天形成一條24維的負荷向量。為了研究用電行為需要擴展到一段時間的平均表現(xiàn)而不是單日偶然值。采集到的原始數(shù)據(jù)不能直接進聚類清洗規(guī)則是把零值時段占比超過80%的用戶剔除掉這些表計可能已經(jīng)拆表或者長期空置把峰值負荷超過電表額定容量的樣本去掉基本是接線錯誤或者轉(zhuǎn)供戶連續(xù)多日缺數(shù)的用戶也刪除。清洗之后我取每戶近三個月的有效數(shù)據(jù)計算特征。我最終用于聚類的特征一共5個維度工作日平均負荷曲線的峰時段用電占比工作日平均負荷曲線的谷時段用電占比休息日平均負荷水平歸一化值負載率即平均負荷與最大負荷的比值峰谷差率反映日內(nèi)用電波動幅度為什么不用完整24維曲線直接聚類從數(shù)據(jù)維度上考慮居民用電曲線24個點之間有強相關性直接聚類容易受到個別時點噪聲的干擾也增加了PSO解空間的維度。比如K4時如果直接用24維數(shù)據(jù)粒子長度就是96用5個統(tǒng)計特征時粒子長度只有20搜索空間小了一個數(shù)量級收斂難度完全不同。數(shù)據(jù)處理上有一件必須做的事對所有特征做標準化。我用的是zscore函數(shù)也就是減去均值除以標準差。這個步驟不做的話量綱大的特征會在距離計算中占據(jù)絕對主導地位比如負荷值本身可能是幾十個單位而負載率是0到1的小數(shù)聚類等價于只用負荷大小說話。% 讀取清洗后的數(shù)據(jù)每一行是一個用戶每一列是一個特征 data readmatrix(residential_features.csv); data_z zscore(data);3.2 PSO優(yōu)化Kmeans的核心代碼實現(xiàn)粒子群優(yōu)化部分的完整代碼并不復雜我拆成幾個函數(shù)來寫方便調(diào)試和復用。下面這段是PSO主流程框架。%% 參數(shù)設置 K 4; % 聚類數(shù) N 30; % 粒子群規(guī)模 T 100; % 迭代次數(shù) D size(data_z, 2); % 特征維度 w_max 0.9; % 慣性權重上限 w_min 0.4; % 慣性權重下限 c1 2; % 個體學習因子 c2 2; % 社會學習因子 %% 解空間邊界聚類中心應落在數(shù)據(jù)取值范圍內(nèi) lb repmat(min(data_z), 1, K); ub repmat(max(data_z), 1, K); %% 種群初始化 X lb rand(N, K * D) .* (ub - lb); % 粒子位置 V zeros(N, K * D); % 粒子速度 pbest X; pbest_fit inf(N, 1); gbest_fit inf; gbest X(1, :); %% 迭代 for iter 1:T w w_max - (w_max - w_min) * iter / T; for i 1:N centers reshape(X(i, :), K, D); [~, fit] computeSSE(data_z, centers); if fit pbest_fit(i) pbest_fit(i) fit; pbest(i, :) X(i, :); end if fit gbest_fit gbest_fit fit; gbest X(i, :); end end for i 1:N V(i, :) w * V(i, :) c1 * rand() * (pbest(i, :) - X(i, :))... c2 * rand() * (gbest - X(i, :)); X(i, :) X(i, :) V(i, :); % 邊界還原粒子越界時拉回邊界 X(i, :) max(min(X(i, :), ub), lb); end fprintf(iter %d, best fitness %.4f\n, iter, gbest_fit); end %% 用PSO找到的最優(yōu)質(zhì)心作為Kmeans的初始中心 init_centers reshape(gbest, K, D); [idx, C] kmeans(data_z, K, Start, init_centers, MaxIter, 500);適應度計算函數(shù)computeSSE的實現(xiàn)是關鍵之一。這里注意不要用雙重循環(huán)逐個樣本算距離要通過矩陣運算一次算出所有樣本到某中心的距離function [labels, sse] computeSSE(data, centers) n size(data, 1); nk size(centers, 1); distMat zeros(n, nk); for k 1:nk diff data - centers(k, :); distMat(:, k) sum(diff.^2, 2); end [minDist, labels] min(distMat, [], 2); sse sum(minDist); end這個函數(shù)里遍歷的是聚類中心數(shù)量k而不是樣本數(shù)量n。k相對于n來說可以忽略不計所以整體效率是可以接受的。如果你的樣本量特別大還可以在PSO內(nèi)部改用pdist2來計算距離矩陣速度會更快一點不過會多占用一些內(nèi)存。我實測下來這個computeSSE用10000行×5列的數(shù)據(jù)一次評估只需要幾毫秒。PSO運行100代、30個粒子總耗時大概在幾十秒級別完全在可接受范圍。3.3 參數(shù)配置與K值確定PSO本身的參數(shù)不是玄學但有幾個需要根據(jù)實際數(shù)據(jù)調(diào)整的地方。種群規(guī)模N。N太小容易漏掉優(yōu)質(zhì)解空間N太大計算量線性上漲但收益邊際遞減。我在這個項目里試過20、30、5030粒子和50粒子的最終SSE相差不到1%運行時間卻差了接近一倍。建議先從N30開始如果發(fā)現(xiàn)gbest適應度曲線后期還在大幅下降說明粒子多樣性不足或迭代次數(shù)不夠再往上加。慣性權重w的線性遞減策略是PSO的常見設定。w大時粒子飛行速度快、全局探索能力強w小時粒子趨于精細局部搜索。從0.9線性降到0.4正好覆蓋前期探索、后期收斂的過程。如果你想進一步優(yōu)化可以改成非線性遞減或者自適應慣性權重但對當前Kmeans初始化這個場景線性遞減已經(jīng)夠用。K值的確定是整個流程里最影響結(jié)果的一步。我用的是SSE肘部法配合輪廓系數(shù)交叉驗證。具體做法分別運行K2到8的PSO-Kmeans把每次的最終SSE值畫成折線圖。通常隨著K增加SSE單調(diào)下降但下降速度會在某個K值之后明顯放緩這個點就是“肘部”。同時使用silhouette函數(shù)計算輪廓系數(shù)輪廓系數(shù)越高說明簇內(nèi)緊湊度越好、簇間分離度越高。肘部法會傾向做參考輪廓系數(shù)在做裁決。我實驗中最優(yōu)K在4到5之間考慮業(yè)務落地的可解釋性最終選擇K4。3.4 結(jié)果評估與業(yè)務映射聚類完成后評價不能只看算法指標。我先把四種聚類模式的SSE和輪廓系數(shù)做了對比方案SSE輪廓系數(shù)單次運行時間普通Kmeans隨機初始化最優(yōu)一次356.820.3120.8sKmeans自帶默認348.610.3311.2sPSO-Kmeans本文方案341.150.36826.5s這個結(jié)果符合預期PSO-Kmeans比普通Kmeans的SSE降低了約4.4%輪廓系數(shù)從0.312提升到0.368聚類質(zhì)量有明顯改善。代價是運行時間從秒級到幾十秒但對于離線分析場景這個時間完全可控。算法指標之外真正讓我認可這個方案的是業(yè)務側(cè)的可解釋性。K4時四個簇的典型負荷曲線特征非常清晰第一簇用戶工作日白天負荷低、晚間負荷集中典型的上班族作息第二簇用戶全天負荷水平低但穩(wěn)定多為老人留守家庭空調(diào)使用率低第三簇用戶峰谷差極大晚高峰負荷顯著高于其他時段可能有電動車充電或者電熱水器集中使用第四簇用戶全天負荷水平高、峰谷差小疑似有分布式光伏或者全天候生產(chǎn)性用電。這種聚類結(jié)果可以直接落到營銷策略上。比如對第三簇用戶可以推薦錯峰充電電價對第四簇用戶不適合做需求響應邀約因為本身負荷剛性太強。4. 實操中踩過的坑與排查實錄項目從零到跑通前前后后折騰了兩周。下面這些坑不是從文檔里看來的都是實際調(diào)試中遇到過的值得認真記一筆。4.1 K值怎么定才靠譜最初我圖省事只看了SSE肘部圖就定了K5結(jié)果聚類出來的第二個簇和第三個簇中心距離極近兩個簇的特征幾乎一致。這就是只看單一指標的下場。后來我把輪廓系數(shù)也加了進來發(fā)現(xiàn)K5時輪廓系數(shù)反而比K4低不少說明5類劃分出現(xiàn)了簇重疊。同時我還計算了戴維斯-布爾丁指數(shù)作為參考DBI越低代表類內(nèi)分散度小、類間分散度大。三項指標一起看K4的結(jié)論就非常穩(wěn)定了。這里給一個建議K值不要只看算法指標。算法指標再好如果聚出來的4個簇在業(yè)務上無法描述、無法對應到客戶群體類型這個聚類就是失敗的。我在實驗中把K3到K6的結(jié)果都打印成負荷曲線圖拿給業(yè)務同事看他們最認可的是K4的結(jié)果理由是“每類用戶都能對上一個可以描述的行為特征”。4.2 PSO搜索效果差、收斂慢的排查有一次實驗PSO迭代了200代但適應度曲線從第60代開始就平了而且最終SSE比普通Kmeans還高。這個現(xiàn)象第一反應以為是PSO代碼寫錯了后來排查完發(fā)現(xiàn)是數(shù)據(jù)標準化沒有做負荷絕對值特征在距離計算中權重太大導致搜索空間里的大部分區(qū)域都是無效區(qū)域。標準化之后問題直接就消失了適應度曲線下降正常了。如果做了標準化還是效果差下一步就檢查種群初始化和邊界設置。我在最初版本里把粒子邊界設成了[0,1]沒有根據(jù)實際數(shù)據(jù)范圍調(diào)整結(jié)果粒子普遍在邊界附近越界被拉回搜索效率極低。后來改成lb repmat(min(data_z), 1, K)先把粒子約束在合理范圍收斂速度明顯改善。還有一個典型問題是PSO陷入局部最優(yōu)。判斷特征是多次運行結(jié)果里gbest_fit每次都一樣但數(shù)值偏大。這說明粒子群沒有很好地覆蓋解空間。處理辦法有三個方向一個是加大種群規(guī)模到50左右另一個是調(diào)整慣性權重下限從0.4降到0.3保持后期一定的探索能力還有一個是在速度更新中給速度設置最大限幅避免粒子飛得太猛直接越過優(yōu)質(zhì)解區(qū)域。我實測中調(diào)整權重下限最有效。4.3 聚類結(jié)果與真實數(shù)據(jù)畫像對不上有一版聚類結(jié)果出現(xiàn)了一個讓我困惑的簇這個簇包含的用戶數(shù)量非常大而且簇心曲線形狀跟總體均值幾乎一樣。排查之后發(fā)現(xiàn)問題出在特征工程我只用了工作日特征和負載率休息日的用電行為信息沒有進入模型導致很多行為差異只在休息日展現(xiàn)的用戶被劃到了一起。加上休息日均值特征和峰谷差率特征之后這個現(xiàn)象消失了。這說明聚類效果的上限其實是由特征工程決定的算法只能在給定特征下盡量優(yōu)化。特征沒有區(qū)分度再好的優(yōu)化算法也白搭。另一個和業(yè)務對不上的原因是直接用PCA降維后聚類再回頭解釋業(yè)務。PCA降維之后再聚類簇在原始特征上的業(yè)務可解釋性往往變差。如果需要降維建議降維后仍然回到原始特征上統(tǒng)計每個簇的業(yè)務指標而不是只看降維空間的坐標。4.4 Matlab代碼層面的幾個注意點第一隨機數(shù)種子。PSO和Kmeans都涉及隨機初始化為了讓實驗可復現(xiàn)務必在腳本開頭加rng(42)或者任意固定的種子。很多科研實驗把隨機性誤當作改進效果最后被人復核時對不上這是最尷尬的。第二避免樣本級的雙重循環(huán)。有人喜歡用兩重for循環(huán)逐樣本計算距離數(shù)據(jù)量小的時候看不出問題幾萬條樣本的時候那個速度真的會懷疑人生。使用向量化矩陣運算代碼短、速度快、可讀性也更好。第三大數(shù)據(jù)量的加速技巧。如果你的樣本超過幾十萬行PSO的適應度計算會開始吃力。此時可以先隨機抽樣一部分代表性樣本做PSO尋找初始中心再用全部數(shù)據(jù)做最終Kmeans。抽樣時注意分層抽樣保證覆蓋全天的負荷特征不要只按序號機械抽。第四聚類完成后要把簇標簽回寫到原始數(shù)據(jù)表。我一開始只保存了最終的聚類中心后來要分析每個簇的客戶特征時還得重新跑一遍分配白白浪費了一次完整的Kmeans迭代。正確做法是在聚類結(jié)束后直接用[idx, C] kmeans(...)把idx存下來匹配到用戶ID上后面畫圖、統(tǒng)計、出報告都不用再重跑?；氐阶畛醯膯栴}——居民用電行為聚類值不值得上PSO這種優(yōu)化算法如果只是寫課程報告自帶Kmeans多跑幾次選最優(yōu)勉強夠用。但如果是做真實業(yè)務的研究結(jié)果穩(wěn)定性和可復現(xiàn)性遠比那幾秒運行時間更重要。PSO-Kmeans用適度的計算代價換取了聚類質(zhì)量的大幅提升讓整個分析鏈條從“碰運氣”變成了“可預期”。在做負荷畫像、需求響應潛力和分時電價策略支撐這類場景時這套流程是值得復制到生產(chǎn)環(huán)境里的。我個人的體會是算法組合不必追求復雜但每個環(huán)節(jié)為什么要這樣設計必須自己能講明白——這比代碼本身更能決定項目的成敗。 PSO優(yōu)化Kmeans的居民用電行為聚類Matlab實現(xiàn)
最近在做居民用電行為分析時被手頭一批智能電表小時級負荷數(shù)據(jù)折騰得不輕。事情本身不復雜——把幾萬個用戶的用電曲線分類畫像但直接用Matlab自帶的kmeans跑結(jié)果實在沒法用同樣一份數(shù)據(jù)十次隨機初始化跑出十種結(jié)果SSE波動幅度接近百分之十幾。這個穩(wěn)定性拿到科研報告里能挑好看的寫但放在實際業(yè)務里根本說不過去。后來把粒子群算法PSO和Kmeans綁在一起讓PSO先全局搜索出一組好的初始聚類中心再把這個中心作為Kmeans的起點做局部收斂。實驗效果不錯但把整個流程走通的過程中踩了不少坑。這篇文章就把我的完整思路、Matlab代碼實現(xiàn)和排查教訓整理出來給正在做居民用電聚類、PSO優(yōu)化算法或者Kmeans相關課題的朋友一個可復用的參考。1. 為什么居民用電行為分析要動用到聚類算法智能電表普及之后電網(wǎng)公司手里掌握的用戶負荷數(shù)據(jù)量級已經(jīng)不是人工能處理的了。一個中型城市一天的采集數(shù)據(jù)就是幾十萬條曲線每條曲線包含24個甚至96個時點。面對這個規(guī)模的數(shù)據(jù)想要回答“哪些用戶適合執(zhí)行峰谷分時電價”“哪些用戶具備可調(diào)負荷參與需求響應”這類問題第一步都是做用戶分類。分類的方式有兩種一種是根據(jù)用戶臺賬手動打標簽比如小區(qū)類型、戶型面積、是否安裝分布式光伏但這種標簽往往滯后而且不完整另一種就是從負荷數(shù)據(jù)本身出發(fā)用無監(jiān)督聚類把行為模式相似的用戶聚到一組這就是用電行為分析的起點。1.1 用電行為聚類能解決什么實際問題居民用電行為聚類不是學術上的自嗨業(yè)務價值可以落到幾個非常具體的方向。第一個是分時電價套餐設計。不同用戶對電價的敏感度差異很大上班族白天空置家庭負荷低晚上回家負荷陡增退休老人家庭白天持續(xù)低負荷運行部分用戶夏季空調(diào)用電占比極高。這些特征通過聚類分組之后就可以針對性地設計尖峰費率、低谷優(yōu)惠、季節(jié)性套餐而不是對所有用戶一刀切。第二個是需求響應潛力評估?，F(xiàn)在的削峰填谷任務越來越重電網(wǎng)需要知道哪些用戶群體在有激勵時能主動壓低負荷。聚類之后如果某個簇的負荷曲線有明顯的用電集中時段、且這些時段恰好與電網(wǎng)峰時段重疊那么這個群體就是需求響應資源池中的優(yōu)質(zhì)對象。第三個是異常用電檢測和客戶細分服務。聚類本質(zhì)上刻畫的是“大多數(shù)用戶”的用電規(guī)律那些偏離規(guī)律的樣本自然成為重點關注對象。服務方面聚類畫像可以作為精準營銷的基礎比如識別出有儲能設備、有電動車充電習慣的用戶群體針對性推薦增值服務?？梢哉f聚類是連接原始負荷數(shù)據(jù)和上層業(yè)務決策的中間層。沒有這個中間層機器學習、數(shù)據(jù)挖掘的手段再多也無法落到電力系統(tǒng)的實際管理動作上。1.2 為什么選Kmeans而不是其他聚類算法用電數(shù)據(jù)分析領域常見的聚類算法不少DBSCAN、層次聚類、譜聚類都有人用但我最終把Kmeans作為基礎方案原因有三個。第一是數(shù)據(jù)規(guī)模。Kmeans的時間復雜度接近線性對幾十萬條樣本依然可以在可接受時間內(nèi)跑完。DBSCAN在實現(xiàn)較好的情況下效率也不錯但在簇密度不均勻的數(shù)據(jù)上參數(shù)很難調(diào)而居民負荷曲線恰好是密度分布極不均勻——白天用電峰值和夜間低谷對應的樣本分布差異非常大。第二是可解釋性。Kmeans的聚類中心就是“平均用電曲線”每個簇可以用一條有物理意義的負荷曲線來描述。層次聚類也能做到但樹狀圖在大樣本下可讀性極差譜聚類的嵌入空間不好向業(yè)務方解釋。第三是工程實現(xiàn)成熟度。Matlab的kmeans函數(shù)集成了多種距離度量、并行計算和重復啟動機制可以直接復用。但Matlab自帶的Kmeans有一個老問題——初始中心隨機選陷入局部最優(yōu)是家常便飯。這為后面引入PSO優(yōu)化提供了充分的動機。2. PSO如何給Kmeans找到更好的初始質(zhì)心Kmeans聚類的核心目標是找到一組聚類中心使所有樣本到所屬中心距離的平方和最小。這個目標函數(shù)非凸自帶Kmeans用的是隨機初始化和多次重復啟動本質(zhì)上是在“碰運氣”。粒子群算法在這里的角色就是替代隨機初始化通過群體智能搜索在解空間中找到一個更接近全局最優(yōu)的起始點。2.1 Kmeans的痛點初始質(zhì)心對結(jié)果的影響遠比你想象中大Kmeans迭代過程可以分成交替的兩步第一步根據(jù)當前聚類中心把每個樣本分配到最近的中心第二步根據(jù)分配結(jié)果重新計算每個簇的均值作為新中心。問題在于這個迭代過程是“局部爬山”。如果初始質(zhì)心選得不好比如兩個初始中心落在同一個數(shù)據(jù)密集區(qū)域那么最終聚類結(jié)果就會在局部最優(yōu)處停下來無法跳到更合理的劃分。真實數(shù)據(jù)里的表現(xiàn)就是隨便跑幾次kmeans每次的SSE都不一樣輪廓系數(shù)也有明顯起伏。我遇到的最夸張的一次同一份數(shù)據(jù)K4SSE最小值與最大值相差22%。這種情況下你沒法判斷當前這個結(jié)果是不是可用的。解決初始質(zhì)心問題有幾個常規(guī)手段。Kmeans的思路是讓初始質(zhì)心盡量分散這個策略對凸簇或者分布比較規(guī)則的簇有效但如果簇形狀重疊、存在噪聲Kmeans的初始中心仍然可能落在不合適的位置多次隨機重啟結(jié)合挑選最小SSE的做法本質(zhì)上是用計算時間去換概率但重復次數(shù)再多也還是隨機搜索沒有目標導向。粒子群算法不一樣它是帶著“方向”去找解的——每個候選解都會計算適應度通過種群協(xié)作逐步逼近更優(yōu)區(qū)域這正是Kmeans兩點方案中最需要的稀缺能力。2.2 粒子群算法的主要思想與更新過程粒子群算法的靈感來自鳥群覓食。每只鳥就是一個粒子代表解空間中的一個候選解鳥群的飛行過程就是不斷搜索最優(yōu)解的過程。每個粒子具備兩個屬性位置和速度。如果把當前解空間記作D維空間第i個粒子的位置記為X_i速度記為V_i。算法為每個粒子保存兩個記憶粒子自己的歷史最優(yōu)位置pbest以及整個種群目前發(fā)現(xiàn)的最優(yōu)位置gbest。每次迭代時粒子的速度和位置按照下面的方式更新速度更新V_i^{t1} wV_i^t c1r1*(pbest_i - X_i^t) c2r2(gbest - X_i^t)位置更新X_i^{t1} X_i^t V_i^{t1}其中w是慣性權重決定當前速度的保持程度c1是認知學習因子c2是社會學習因子r1和r2是[0,1]之間的隨機數(shù)。三個分量的物理意義很清楚慣性項讓粒子沿當前趨勢繼續(xù)飛行認知項把粒子拉向自己曾經(jīng)找到過的好位置社會項把整個種群引向全局最優(yōu)位置。這種機制在連續(xù)域搜索問題中表現(xiàn)穩(wěn)定參數(shù)少、實現(xiàn)簡單、不需要梯度信息非常適合用來優(yōu)化Kmeans的初始質(zhì)心。2.3 粒子編碼方式與適應度函數(shù)設計用PSO優(yōu)化Kmeans初始質(zhì)心最關鍵的環(huán)節(jié)是把Kmeans的解映射成PSO粒子。我采用的編碼方式是如果聚成K個簇每個簇中心是d維特征向量那么一個粒子的位置就是一個長度為 K*d 的一維向量前d個元素是第1個簇中心的坐標接下來d個元素是第2個簇中心的坐標依次類推。適應度函數(shù)的設計直接影響優(yōu)化效果。我用的適應度是樣本到最近聚類中心的歐氏距離平方和即Kmeans的SSE目標函數(shù)。這個選擇的理由很直接——PSO搜索到的適應度越小就意味著這組初始質(zhì)心本身已經(jīng)處在一個比較好的聚類劃分附近。為什么不直接用輪廓系數(shù)作為適應度因為輪廓系數(shù)要計算所有樣本兩兩之間的距離關系復雜度接近O(N2)幾萬條樣本時每次適應度評估都慢得難以接受而SSE的計算是線性的在PSO上百次迭代中每次都算輪廓系數(shù)實驗基本就等死了。還有一個常見做法值得提一下PSO每評估一次粒子時讓這個粒子先跑幾步Kmeans再計算適應度相當于是局部搜索和全局搜索的混合。這個方案效果好但時間復雜度翻了幾倍我在實際項目中用的是樸素的兩階段法——PSO只負責找初始質(zhì)心找到之后交給Kmeans做最終收斂。這樣工程上簡單效果已經(jīng)足夠替代隨機初始化了。3. Matlab實現(xiàn)從數(shù)據(jù)清洗到結(jié)果評估整個實驗我用的Matlab版本是R2024b代碼實現(xiàn)層面沒有用任何工具箱之外的特殊功能換到R2020之后的老版本也能跑。實現(xiàn)過程分為數(shù)據(jù)清洗、特征工程、PSO核心代碼、參數(shù)配置和結(jié)果評估幾個環(huán)節(jié)。3.1 數(shù)據(jù)準備與特征工程我手上的原始數(shù)據(jù)是一個居民小區(qū)的智能電表負荷記錄每一行是一個用戶每一小時一個采集點每天形成一條24維的負荷向量。為了研究用電行為需要擴展到一段時間的平均表現(xiàn)而不是單日偶然值。采集到的原始數(shù)據(jù)不能直接進聚類清洗規(guī)則是把零值時段占比超過80%的用戶剔除掉這些表計可能已經(jīng)拆表或者長期空置把峰值負荷超過電表額定容量的樣本去掉基本是接線錯誤或者轉(zhuǎn)供戶連續(xù)多日缺數(shù)的用戶也刪除。清洗之后我取每戶近三個月的有效數(shù)據(jù)計算特征。我最終用于聚類的特征一共5個維度工作日平均負荷曲線的峰時段用電占比工作日平均負荷曲線的谷時段用電占比休息日平均負荷水平歸一化值負載率即平均負荷與最大負荷的比值峰谷差率反映日內(nèi)用電波動幅度為什么不用完整24維曲線直接聚類從數(shù)據(jù)維度上考慮居民用電曲線24個點之間有強相關性直接聚類容易受到個別時點噪聲的干擾也增加了PSO解空間的維度。比如K4時如果直接用24維數(shù)據(jù)粒子長度就是96用5個統(tǒng)計特征時粒子長度只有20搜索空間小了一個數(shù)量級收斂難度完全不同。數(shù)據(jù)處理上有一件必須做的事對所有特征做標準化。我用的是zscore函數(shù)也就是減去均值除以標準差。這個步驟不做的話量綱大的特征會在距離計算中占據(jù)絕對主導地位比如負荷值本身可能是幾十個單位而負載率是0到1的小數(shù)聚類等價于只用負荷大小說話。% 讀取清洗后的數(shù)據(jù)每一行是一個用戶每一列是一個特征 data readmatrix(residential_features.csv); data_z zscore(data);3.2 PSO優(yōu)化Kmeans的核心代碼實現(xiàn)粒子群優(yōu)化部分的完整代碼并不復雜我拆成幾個函數(shù)來寫方便調(diào)試和復用。下面這段是PSO主流程框架。%% 參數(shù)設置 K 4; % 聚類數(shù) N 30; % 粒子群規(guī)模 T 100; % 迭代次數(shù) D size(data_z, 2); % 特征維度 w_max 0.9; % 慣性權重上限 w_min 0.4; % 慣性權重下限 c1 2; % 個體學習因子 c2 2; % 社會學習因子 %% 解空間邊界聚類中心應落在數(shù)據(jù)取值范圍內(nèi) lb repmat(min(data_z), 1, K); ub repmat(max(data_z), 1, K); %% 種群初始化 X lb rand(N, K * D) .* (ub - lb); % 粒子位置 V zeros(N, K * D); % 粒子速度 pbest X; pbest_fit inf(N, 1); gbest_fit inf; gbest X(1, :); %% 迭代 for iter 1:T w w_max - (w_max - w_min) * iter / T; for i 1:N centers reshape(X(i, :), K, D); [~, fit] computeSSE(data_z, centers); if fit pbest_fit(i) pbest_fit(i) fit; pbest(i, :) X(i, :); end if fit gbest_fit gbest_fit fit; gbest X(i, :); end end for i 1:N V(i, :) w * V(i, :) c1 * rand() * (pbest(i, :) - X(i, :))... c2 * rand() * (gbest - X(i, :)); X(i, :) X(i, :) V(i, :); % 邊界還原粒子越界時拉回邊界 X(i, :) max(min(X(i, :), ub), lb); end fprintf(iter %d, best fitness %.4f\n, iter, gbest_fit); end %% 用PSO找到的最優(yōu)質(zhì)心作為Kmeans的初始中心 init_centers reshape(gbest, K, D); [idx, C] kmeans(data_z, K, Start, init_centers, MaxIter, 500);適應度計算函數(shù)computeSSE的實現(xiàn)是關鍵之一。這里注意不要用雙重循環(huán)逐個樣本算距離要通過矩陣運算一次算出所有樣本到某中心的距離function [labels, sse] computeSSE(data, centers) n size(data, 1); nk size(centers, 1); distMat zeros(n, nk); for k 1:nk diff data - centers(k, :); distMat(:, k) sum(diff.^2, 2); end [minDist, labels] min(distMat, [], 2); sse sum(minDist); end這個函數(shù)里遍歷的是聚類中心數(shù)量k而不是樣本數(shù)量n。k相對于n來說可以忽略不計所以整體效率是可以接受的。如果你的樣本量特別大還可以在PSO內(nèi)部改用pdist2來計算距離矩陣速度會更快一點不過會多占用一些內(nèi)存。我實測下來這個computeSSE用10000行×5列的數(shù)據(jù)一次評估只需要幾毫秒。PSO運行100代、30個粒子總耗時大概在幾十秒級別完全在可接受范圍。3.3 參數(shù)配置與K值確定PSO本身的參數(shù)不是玄學但有幾個需要根據(jù)實際數(shù)據(jù)調(diào)整的地方。種群規(guī)模N。N太小容易漏掉優(yōu)質(zhì)解空間N太大計算量線性上漲但收益邊際遞減。我在這個項目里試過20、30、5030粒子和50粒子的最終SSE相差不到1%運行時間卻差了接近一倍。建議先從N30開始如果發(fā)現(xiàn)gbest適應度曲線后期還在大幅下降說明粒子多樣性不足或迭代次數(shù)不夠再往上加。慣性權重w的線性遞減策略是PSO的常見設定。w大時粒子飛行速度快、全局探索能力強w小時粒子趨于精細局部搜索。從0.9線性降到0.4正好覆蓋前期探索、后期收斂的過程。如果你想進一步優(yōu)化可以改成非線性遞減或者自適應慣性權重但對當前Kmeans初始化這個場景線性遞減已經(jīng)夠用。K值的確定是整個流程里最影響結(jié)果的一步。我用的是SSE肘部法配合輪廓系數(shù)交叉驗證。具體做法分別運行K2到8的PSO-Kmeans把每次的最終SSE值畫成折線圖。通常隨著K增加SSE單調(diào)下降但下降速度會在某個K值之后明顯放緩這個點就是“肘部”。同時使用silhouette函數(shù)計算輪廓系數(shù)輪廓系數(shù)越高說明簇內(nèi)緊湊度越好、簇間分離度越高。肘部法會傾向做參考輪廓系數(shù)在做裁決。我實驗中最優(yōu)K在4到5之間考慮業(yè)務落地的可解釋性最終選擇K4。3.4 結(jié)果評估與業(yè)務映射聚類完成后評價不能只看算法指標。我先把四種聚類模式的SSE和輪廓系數(shù)做了對比方案SSE輪廓系數(shù)單次運行時間普通Kmeans隨機初始化最優(yōu)一次356.820.3120.8sKmeans自帶默認348.610.3311.2sPSO-Kmeans本文方案341.150.36826.5s這個結(jié)果符合預期PSO-Kmeans比普通Kmeans的SSE降低了約4.4%輪廓系數(shù)從0.312提升到0.368聚類質(zhì)量有明顯改善。代價是運行時間從秒級到幾十秒但對于離線分析場景這個時間完全可控。算法指標之外真正讓我認可這個方案的是業(yè)務側(cè)的可解釋性。K4時四個簇的典型負荷曲線特征非常清晰第一簇用戶工作日白天負荷低、晚間負荷集中典型的上班族作息第二簇用戶全天負荷水平低但穩(wěn)定多為老人留守家庭空調(diào)使用率低第三簇用戶峰谷差極大晚高峰負荷顯著高于其他時段可能有電動車充電或者電熱水器集中使用第四簇用戶全天負荷水平高、峰谷差小疑似有分布式光伏或者全天候生產(chǎn)性用電。這種聚類結(jié)果可以直接落到營銷策略上。比如對第三簇用戶可以推薦錯峰充電電價對第四簇用戶不適合做需求響應邀約因為本身負荷剛性太強。4. 實操中踩過的坑與排查實錄項目從零到跑通前前后后折騰了兩周。下面這些坑不是從文檔里看來的都是實際調(diào)試中遇到過的值得認真記一筆。4.1 K值怎么定才靠譜最初我圖省事只看了SSE肘部圖就定了K5結(jié)果聚類出來的第二個簇和第三個簇中心距離極近兩個簇的特征幾乎一致。這就是只看單一指標的下場。后來我把輪廓系數(shù)也加了進來發(fā)現(xiàn)K5時輪廓系數(shù)反而比K4低不少說明5類劃分出現(xiàn)了簇重疊。同時我還計算了戴維斯-布爾丁指數(shù)作為參考DBI越低代表類內(nèi)分散度小、類間分散度大。三項指標一起看K4的結(jié)論就非常穩(wěn)定了。這里給一個建議K值不要只看算法指標。算法指標再好如果聚出來的4個簇在業(yè)務上無法描述、無法對應到客戶群體類型這個聚類就是失敗的。我在實驗中把K3到K6的結(jié)果都打印成負荷曲線圖拿給業(yè)務同事看他們最認可的是K4的結(jié)果理由是“每類用戶都能對上一個可以描述的行為特征”。4.2 PSO搜索效果差、收斂慢的排查有一次實驗PSO迭代了200代但適應度曲線從第60代開始就平了而且最終SSE比普通Kmeans還高。這個現(xiàn)象第一反應以為是PSO代碼寫錯了后來排查完發(fā)現(xiàn)是數(shù)據(jù)標準化沒有做負荷絕對值特征在距離計算中權重太大導致搜索空間里的大部分區(qū)域都是無效區(qū)域。標準化之后問題直接就消失了適應度曲線下降正常了。如果做了標準化還是效果差下一步就檢查種群初始化和邊界設置。我在最初版本里把粒子邊界設成了[0,1]沒有根據(jù)實際數(shù)據(jù)范圍調(diào)整結(jié)果粒子普遍在邊界附近越界被拉回搜索效率極低。后來改成lb repmat(min(data_z), 1, K)先把粒子約束在合理范圍收斂速度明顯改善。還有一個典型問題是PSO陷入局部最優(yōu)。判斷特征是多次運行結(jié)果里gbest_fit每次都一樣但數(shù)值偏大。這說明粒子群沒有很好地覆蓋解空間。處理辦法有三個方向一個是加大種群規(guī)模到50左右另一個是調(diào)整慣性權重下限從0.4降到0.3保持后期一定的探索能力還有一個是在速度更新中給速度設置最大限幅避免粒子飛得太猛直接越過優(yōu)質(zhì)解區(qū)域。我實測中調(diào)整權重下限最有效。4.3 聚類結(jié)果與真實數(shù)據(jù)畫像對不上有一版聚類結(jié)果出現(xiàn)了一個讓我困惑的簇這個簇包含的用戶數(shù)量非常大而且簇心曲線形狀跟總體均值幾乎一樣。排查之后發(fā)現(xiàn)問題出在特征工程我只用了工作日特征和負載率休息日的用電行為信息沒有進入模型導致很多行為差異只在休息日展現(xiàn)的用戶被劃到了一起。加上休息日均值特征和峰谷差率特征之后這個現(xiàn)象消失了。這說明聚類效果的上限其實是由特征工程決定的算法只能在給定特征下盡量優(yōu)化。特征沒有區(qū)分度再好的優(yōu)化算法也白搭。另一個和業(yè)務對不上的原因是直接用PCA降維后聚類再回頭解釋業(yè)務。PCA降維之后再聚類簇在原始特征上的業(yè)務可解釋性往往變差。如果需要降維建議降維后仍然回到原始特征上統(tǒng)計每個簇的業(yè)務指標而不是只看降維空間的坐標。4.4 Matlab代碼層面的幾個注意點第一隨機數(shù)種子。PSO和Kmeans都涉及隨機初始化為了讓實驗可復現(xiàn)務必在腳本開頭加rng(42)或者任意固定的種子。很多科研實驗把隨機性誤當作改進效果最后被人復核時對不上這是最尷尬的。第二避免樣本級的雙重循環(huán)。有人喜歡用兩重for循環(huán)逐樣本計算距離數(shù)據(jù)量小的時候看不出問題幾萬條樣本的時候那個速度真的會懷疑人生。使用向量化矩陣運算代碼短、速度快、可讀性也更好。第三大數(shù)據(jù)量的加速技巧。如果你的樣本超過幾十萬行PSO的適應度計算會開始吃力。此時可以先隨機抽樣一部分代表性樣本做PSO尋找初始中心再用全部數(shù)據(jù)做最終Kmeans。抽樣時注意分層抽樣保證覆蓋全天的負荷特征不要只按序號機械抽。第四聚類完成后要把簇標簽回寫到原始數(shù)據(jù)表。我一開始只保存了最終的聚類中心后來要分析每個簇的客戶特征時還得重新跑一遍分配白白浪費了一次完整的Kmeans迭代。正確做法是在聚類結(jié)束后直接用[idx, C] kmeans(...)把idx存下來匹配到用戶ID上后面畫圖、統(tǒng)計、出報告都不用再重跑?；氐阶畛醯膯栴}——居民用電行為聚類值不值得上PSO這種優(yōu)化算法如果只是寫課程報告自帶Kmeans多跑幾次選最優(yōu)勉強夠用。但如果是做真實業(yè)務的研究結(jié)果穩(wěn)定性和可復現(xiàn)性遠比那幾秒運行時間更重要。PSO-Kmeans用適度的計算代價換取了聚類質(zhì)量的大幅提升讓整個分析鏈條從“碰運氣”變成了“可預期”。在做負荷畫像、需求響應潛力和分時電價策略支撐這類場景時這套流程是值得復制到生產(chǎn)環(huán)境里的。我個人的體會是算法組合不必追求復雜但每個環(huán)節(jié)為什么要這樣設計必須自己能講明白——這比代碼本身更能決定項目的成敗。