基于簡化HMM和時間分段的非侵入式負荷分解算法

劉 凱，符 玲，楊金剛，熊思宇，蒿保龍，劉麗娜

（1.西南交通大學 唐山研究院，河北 唐山 063000；2.西南交通大學 電氣工程學院，四川 成都 611756；3.國網四川省電力公司計量中心，四川 成都 610045）

0 引言

以高級量測體系（advanced metering infrastruc?ture，AMI）等技術作為支撐，家庭用戶與供電公司之間的關系由原來的家庭用戶單方向用電發展為家庭用戶與供電公司雙向靈活互動。負荷監測是AMI的一種重要功能，可為供電公司提供詳細的用戶用電數據信息，同時也可使用戶更清晰地了解自己的用電行為，從而幫助用戶／供電公司選擇更節能高效的用電／供電策略。

負荷監測按照各種設備級負荷數據的獲取方式分為侵入式負荷監測（intrusive load monitoring，ILM）和非侵入式負荷監測（non-intrusive load monitoring，NILM）。相比于直接獲取負荷數據的ILM，NILM 更依賴于算法，通過分解在用戶電力入口端測量的集總數據得到用戶家庭內部各種設備的用電詳情［1］。NILM 由于具有更好的經濟性和隱私保護性，受到了大量學者的關注［2］。NILM 的求解需要以負荷特征作為支撐，根據數據的不同采樣頻率，可以將負荷特征分為高頻（頻率大于1 Hz）特征和低頻（頻率不大于1 Hz）特征［3］。高頻特征主要包括負荷的電流諧波、電壓-電流軌跡、高頻暫態波型等，低頻特征包括負荷的有功功率、無功功率、電流有效值等。雖然高頻特征包含更多的負荷信息，但是現有AMI 中的智能電表一般僅提供低頻數據的采集支持［4］，要獲得高頻特征就需要額外添加采集設備，這會增加經濟和人工成本，因此，利用智能電表提供的低頻特征數據來實現居民家庭的負荷分解，是一種具有實際應用價值的方案［5］。

負荷分解算法可以分為2 類。一類是組合優化，如：文獻［6］提出利用時間概率因子來改進目標函數，再利用遺傳優化算法進行負荷分解，該算法的準確率高但求解時間長；文獻［7］結合改進的交叉熵算法和懲罰函數，通過迭代更新概率求取最優解并將其作為分解結果。另一類是模式識別，如：文獻［8］利用各種電器間的關聯規則，基于信息熵和k近鄰算法實現狀態辨識，但僅考慮兩兩負荷之間的關聯；文獻［9］提出利用模體挖掘的方式，結合事件檢測與調和函數進行事件匹配和設備識別，但算法依賴于運行窗的選取。近年來，隨著神經網絡的發展，利用神經網絡來實現負荷分解也受到了研究人員的關注，如：文獻［10］提出利用深度神經網絡，綜合電器的時間和功率信息進行負荷分解，但是模型參數的訓練需要大量的歷史數據和時間。在負荷分解的求解算法中經常采用隱馬爾可夫模型（hidden Markov model，HMM），如：文獻［11］為每種負荷均建立一條HMM 鏈，再結合二次整數規劃進行求解，通過假設初始分布概率相等簡化了部分運算，但是算法的計算復雜度較高；文獻［12］為一天中的每個小時單獨建立HMM，再結合行為影響因子來實現負荷辨識?？傮w而言，現有的低頻負荷分解算法模型復雜，求解時間長，在部署監測設備方面壓力較大。此外，大部分算法需要將前一個或多個時刻的分解結果作為當前時刻分解的依據，若前面的分解結果出錯，則會造成誤差累積。

基于以上分析，本文提出一種利用有功功率，并在每個時刻獨立求解的負荷分解算法。首先，在訓練階段，利用改進的親和力傳播聚類算法自適應地獲取負荷功率模板，并將數據按照時間分段；然后，對HMM 進行簡化，構建負荷分解算法，并在不同的時間段內分別訓練模型參數，在負荷分解階段，結合時間信息實時求取用戶負荷功率組合；最后，基于2 個國外的低頻公開數據集的測試結果表明，本文算法準確率高且所需的求解時間短。

1 負荷的超狀態編碼

1.1 負荷分解

假設居民家庭內部共包含M個負荷，任意負荷i的有功功率Pi≥0，則在t時刻用戶總的有功功率可以表示為：

式中：s為t時刻負荷i是否運行在狀態n的標志，若運行在狀態n，則其值為1，否則為0；pi，n為負荷i在狀態n下的功率模板。通常，在正常工作情況下，負荷功率模板不會發生變化，因此，在確定負荷功率模板后，負荷功率分解問題就可以表示為t時刻所有負荷狀態的求解問題。

1.2 負荷狀態編碼

可以將居民家庭中的負荷按照運行時可能出現的狀態分為只存在開關狀態的Ⅰ型負荷、具有有限運行狀態數的Ⅱ型負荷、狀態可以連續變化的Ⅲ型負荷以及一直處于某一功率運行的Ⅳ型負荷［13］。當負荷為Ⅰ型負荷時，以“0”表示關閉狀態，以“1”表示開啟狀態；當負荷為Ⅱ型負荷時，以有限位數的“0”“1”組成的向量表示負荷狀態；當負荷為Ⅲ型負荷時，負荷功率特征在有限的范圍內連續變化，一般通過聚類的方式對負荷狀態進行離散化，從而將其轉化為有限種狀態［8］；當負荷為Ⅳ型負荷時，可以將其等效為缺少“0”狀態的前3 種負荷。由于通過聚類可將所有負荷均等效為有限狀態機，因此，采用獨熱碼來表示各種負荷狀態［14］。對于擁有Ni種狀態的負荷，其獨熱碼位數為Ni，且對于任意給定的狀態，獨熱碼向量的各元素中僅有1 個元素置1。若將整個家庭所有負荷的獨熱碼按順序拼接，則可以將家庭中所有的負荷狀態表示為1個超狀態向量S［15］，如式（3）所示。

式中：sk(k=1，2，…，K)的取值為0 或1，K為所有負荷的運行狀態總數。例如，假設用戶家庭中只有電燈和洗衣機，其中電燈為Ⅰ型負荷，只有開關狀態，洗衣機為Ⅱ型負荷，有浸泡、洗滌、脫水3 種開啟狀態，則可以在超狀態向量中用第1 個元素表示電燈的狀態，第2 個元素表示洗衣機的浸泡狀態，第3 個元素表示洗衣機的洗滌狀態，第4 個元素表示洗衣機的脫水狀態，若S=[0，1，0，0]，則表示電燈處于關閉狀態，洗衣機處于浸泡狀態。

2 負荷功率模板的獲取

獲取負荷功率模板是進行負荷分解的重要前提，考慮到每種負荷的運行狀態數未知，K-means 和模糊C均值聚類算法需提前確定聚類中心數［16］，而親和力傳播（affinity propagation，AP）算法不需要提前確定最終的聚類中心數，因此，本文利用AP 算法來獲取負荷功率模板。

2.1 聚類數據的預處理

AP 算法需要人為調整的參數少，降低了人為因素對結果的干擾［17］。AP 算法的具體計算流程參考文獻［17］，該算法主要包括相似度矩陣E、歸屬度矩陣A和吸引度矩陣R這3 個參數矩陣。算法所需時間和計算結果受到樣本點數的影響：當聚類的樣本點過少時，聚類結果會不準確；當聚類的樣本點過多時，算法的計算時間會增加。因此，本文提出對大量的樣本進行分層抽樣，在滿足較長時間跨度的同時，減少樣本點。

本文采用的分層抽樣規則如下：

1）由于實際運行時負荷處于關閉狀態的時間遠長于處于開啟狀態的時間，同時考慮到噪聲干擾，先剔除負荷功率小于10 W的部分；

2）將10 W 與最大功率之間的部分等分為5 層，在每層隨機抽取α個樣本點，對于樣本點數不足的層，保留該層的全部樣本點；

3）將樣本點數小于5 的層中所有數據剔除，消去異常樣本點的干擾；

4）對剩余的樣本點進行聚類。

通過分層抽樣可將盡可能多地反映負荷運行狀態的數據作為聚類依據，同時剔除了異常點，減少了用于聚類算法的樣本點。

2.2 AP算法的改進

在傳統AP 算法中，當歸屬度矩陣A和吸引度矩陣R的計算結果在一定的迭代次數內保持不變或者迭代達到設定次數時，停止迭代，并且聚類中心c由式（4）確定。

在確定各點所屬的簇中心時，會出現由式（4）得到節點ci和節點cj均為聚類中心，但節點cj所屬簇中心的計算結果卻為節點ci的情況。傳統AP 算法直接將節點cj所屬簇中心強制確定為節點cj，這會造成2 個相近的點被劃分至2 個簇的情況。為了避免產生相近的簇中心，本文將節點cj和歸屬于節點cj的所有點都重新劃歸為屬于節點ci。

此外，最終的聚類中心數受到參考度取值的影響。在得到聚類中心后，可以計算各點到其聚類中心的距離之和D，D的大小能夠反映聚類中心數的合理性。聚類中心過多的極限情況是每個點均為聚類中心，此時D為0；聚類中心過少的極限情況是聚類中心只有1 個，此時D的最大值為樣本點到其余樣本點距離之和的最大值，即相似度矩陣中絕對值最大的列元素之和。為了使不同分布狀況的負荷均能獲得適宜的簇中心，求取D與相似度矩陣中列元素之和最小值的比值γ，再限定該比值的范圍，并以該范圍為依據，采用迭代更新的方式自適應選擇每種負荷的參考度，具體方法如下。

1）由于參考度絕對值越接近0，最終的聚類中心數越多，因此，初始參考度選取為所有樣本點兩兩之間距離中值的1/4，通過聚類得到各點所屬的簇，并計算D。

2）將γ與閾值進行比較。通過實驗，本文選取閾值為0.01 和0.1，即當γ小于0.01 時，減小參考度，當γ大于0.1 時，增大參考度，并重新聚類獲取聚類結果，直到滿足范圍要求。

最終，負荷功率模板的獲取過程包括聚類前的樣本點分層抽樣、聚類時相近簇中心的剔除以及參考度的自適應迭代更新，具體如附錄A 圖A1 所示。在通過聚類得到各負荷功率模板后，可以計算超狀態向量S的元素數，再將S與負荷功率模板Ptem位置一一對應。式（2）可以表示為：

式中：St為t時刻的超狀態向量。

3 HMM的簡化與分時段負荷的分解

3.1 HMM的簡化

HMM 主要包括負荷的初始概率分布矩陣π、負荷的狀態轉移概率矩陣Z以及表示負荷狀態與功率關系的發射矩陣B這3 個參數矩陣?；贖MM 的負荷分解算法滿足馬爾可夫鏈的齊次假設，即馬爾可夫鏈的一步轉移概率與馬爾可夫鏈所處時間無關。但是，考慮到用戶的用電習慣，各種負荷在不同時間段內運行的概率不同，這導致π、Z、B也不同。此外，HMM 需要根據前一時刻的負荷狀態分解結果，結合狀態轉移概率矩陣Z來求解當前時刻的狀態，若前一時刻負荷的狀態識別結果與真實情況不符，則會影響后續的負荷分解結果。

本文將HMM進行簡化：①不考慮狀態轉移概率矩陣，只保留π和B，從而在算法求解時減少運算量，縮減計算時間；②在參數學習階段，通過用戶家庭的歷史用電數據減少超狀態種類數，并通過時間分段進一步限制在不同時間段內可能出現的超狀態種類數，同時，分別在不同時間段通過訓練獲取π和B，從而在這2個矩陣中保留各負荷在不同時間段的運行差異信息。

3.2 時間段的劃分

按照一般用戶家庭生活習慣，可以將一天劃分為睡覺時間段、與飲食相關的時間段、晚飯后的休閑娛樂時間段以及上述各時間段之間的過渡時間段。本文將一天劃分為11 個時間段，其中01:00 — 05:00為第1 個時間段，將其余時間段按照2 h 為1 個時間段進行劃分，如附錄A圖A2所示。

3.3 超狀態參數的獲取

根據負荷功率模板，將訓練集中各負荷功率序列標記為負荷狀態，將訓練集數據按照時間段劃分為各子數據集，在每個時間段，統計超狀態的種類數，計為Jλ。初始概率分布矩陣πλ可以表示為：

式中：V（λλ=1，2，…，11）為第λ個時間段內的樣本總數；VSj(j=1，2，…，Jλ)為超 狀態向 量Sj下的樣本點數。

在每個時間段計算超狀態的發射矩陣Bλ∈RJλ×Oλ，Oλ=P-P+1 為在第λ個時間段總功率的取值區間，P、P分別為第λ個時間段總功率的最大值和最小值。為了減小發射矩陣的維數，將功率值取整。發射矩陣的元素bjo表示當用戶家庭整體狀態構成超狀態向量Sj時功率觀測值為o的概率，計算公式為：

式中：Hjo為超狀態向量Sj下功率觀測值為o的頻數；Hj為超狀態向量Sj下的總頻數。

3.4 負荷功率分解

算法的輸入為總功率的觀測值ot及其對應的時間信息?；诜謺r段初始概率分布矩陣和發射矩陣的功率求解步驟如下。

1）利用輸入的時間信息選擇不同時間段的初始概率分布矩陣πλ和發射矩陣Bλ。

2）由ot確定發射矩陣中的列，t時刻估計的超狀態矩陣為：

式中：argmax 為求參函數。由于Bλ中每行對應一種超狀態，因此，式（8）實際是求取使得概率值最大的行數，從而確定對應的超狀態。

3）由超狀態和負荷功率模板矩陣解碼求得各負荷在t時刻的功率值。

4 算例分析

4.1 算法評估指標

負荷識別結果準確率指標主要包括負荷狀態識別準確率指標和負荷功率分解準確率指標。

為了便于與現有方法進行對比，本文選擇文獻［18］定義的有限狀態識別準確率指標F1，計算公式為：

式中：p為識別為開啟狀態樣本的分解結果中正確結果的比例；r為真實開啟狀態樣本的分解結果中正確結果的比例；Tp為正確識別負荷處于開啟狀態的次數；Fp為負荷實際處于關閉狀態但識別結果為開啟狀態的次數；Fn為負荷實際處于開啟狀態但識別結果為關閉狀態的次數；Tinacc為附加懲罰項，引入該變量是由于對于有多種開啟狀態的負荷，存在雖然識別結果和實際都是開啟狀態但識別的具體開啟狀態與實際不符的情況，其計算公式如式（12）所示。

式中：T為測試集的總時長包含的時刻數；、st，i分別為t時刻負荷i的估計狀態和真實狀態，當估計狀態和真實狀態不同時，兩者差值為1。

負荷功率分解準確率apower也常用于評估負荷識別的效果［19］，計算公式為：

4.2 分解結果

為了測試算法效果，本文選取2 個由真實家庭采集的低頻數據集REDD（reference energy disag?gregation dataset）［20］和AMPds（almanac of minutely power dataset）［21］，其中REDD 來自美國波士頓地區的多個家庭，約每3 s 進行1 次負荷功率數據的采集和儲存。AMPds 記錄了溫哥華地區1 戶居民2 a 的功率數據，每1 min進行1次數據記錄。同時，這2個數據集均以時間戳為時標對每個數據進行標注，結合地理時區劃分能夠得到每個數據在當地的具體時間，方便進行時間段的劃分。

1）在REDD上的測試。

將文獻［11］基于二次規劃的因子HMM 算法和文獻［12］利用行為影響因子的HMM 算法與本文算法進行對比，選取相同的家庭1、2、3、6 進行測試。為了獲取負荷在各時間段內的參數以及功率模板，將前2 周的數據作為訓練集，將其余數據作為測試集，負荷功率分解準確率對比如表1所示。

表1 負荷功率分解準確率對比Table 1 Comparison of load power decomposition accuracy rate

由表1 可知，本文算法的負荷功率分解準確率均高于文獻［11］算法，本文算法對家庭1和家庭6的負荷功率分解準確率高于文獻［12］算法，對家庭2和家庭3 的負荷功率分解準確率與文獻［12］算法相當。

為了便于驗證本文算法對單個負荷的識別準確率，選擇與文獻［7］中場景5 相同的實驗設置，并與該文獻中的改進交叉熵（modified cross-entropy，MCE）算法和分段整數約束規劃（segmented integer quadratic constraint programming，SIQCP）算法進行對比，分解結果如附錄A 表A1 所示。由表可知：本文算法對大多負荷的識別效果均更具優勢，且整體識別效果優于文獻［7］算法；雖然本文算法對個別負荷（如烤箱）的識別效果與文獻［7］算法差距很大，但是整體識別效果和F1是將整個測試階段的所有負荷分解結果作為整體進行計算的，烤箱等負荷使用較少，對整體結果的影響很小。

2）在AMPds上的測試。

相較于REDD，AMPds 的時間跨度更大，采樣頻率更低，與實際的智能電表數據更接近。選擇文獻［7］中的場景3對MCE算法和本文算法進行對比，同樣選擇600 h 的測試數據，選擇相同的7 種運行頻次較高的負荷，即地下室、干衣機、洗碗機、暖通空調、冰箱、熱泵、壁爐。用于訓練的數據時長為40 d，識別結果準確率如表2所示。

表2 各負荷識別結果準確率Table 2 Accuracy rate of identification results for each load

由表2 可知，本文算法總的負荷功率分解準確率達到了89.5 %，且有限狀態識別準確率達到了94.2 %，比MCE 算法的計算結果有了較大幅度的提高，這表明本文算法可得到較好的負荷分解結果，其中，暖通空調主要是在夜晚使用，以保持室內溫度舒適，具有較強的時間段分布特性，因此本文算法對暖通空調的識別準確率提升幅度較大。

為了驗證本文算法在較長時間下的分解準確率以及算法訓練集時長對準確率的影響，設計包含8 種負荷（在前文7 種負荷的基礎上增加電視機）的4 種情景：場景1 — 4 分別利用20、40、60、80 d 的歷史數據訓練參數，再對360 d 的數據進行分解，負荷識別結果準確率如表3所示。

表3 不同時長的負荷識別結果準確率Table 3 Accuracy rate of load identification results for different time lengths

由表3 可知，當用于訓練的數據時長達到40 d后，識別結果準確率趨于穩定，延長訓練時長不會顯著提高識別結果準確率，且當用于訓練的數據時長為40 d時，負荷功率分解準確率能夠達到90.0 %，有限狀態識別準確率能夠達到90.7 %。

將場景2 下本文算法對各負荷的分解結果與真實值進行比較，結果如表4 所示。由表可知：本文算法對電視機功率的分解值與真實值間有較大差別，達到3 %，這主要是由于在AMPds 中電視機的子數據集并非僅記錄了電視機的功率數據，還包括個人錄像機頂盒、音箱等設備的數據，這些數據構成一個合成的待機狀態功率，在整個數據集中的第132 天，待機狀態功率由39 W 降到19 W，并在之后保持不變，而場景2 的負荷功率模板是通過對數據集前40 d的數據進行聚類而獲取的，因此，無法描述該差別，從而給分解帶來了誤差；本文算法對其他負荷功率的分解值與真實值均較接近。

表4 場景2下各負荷消耗的功率Table 4 Power consumption of each load under Scenario 2

4.3 算法復雜度分析

本文算法不考慮HMM中的轉移概率矩陣，且利用獨熱碼對負荷狀態進行編碼縮短了算法所需時間。由于訓練階段可離線進行，且參數訓練完成后可以離線調用，因此，不考慮訓練階段的計算時間。假設在負荷分解時各時間段內最大的超狀態數為G，在求解時輸入觀測值ot后，得到對應的發射矩陣中的列索引，再將發射矩陣的該列向量與初始概率矩陣進行按位乘法運算，即最多進行G次乘法運算，因此，算法的最大復雜度為O（G）。實際上，由于發射矩陣具有稀疏性，某一確定的觀測值對應的超狀態數遠小于G，因此，復雜度也小于O（G）。本文的測試軟件為MATLAB R2021b，硬件為搭載了Intel（R） Core（TM） i5-7300HQ、內存為8 GB 的筆記本電腦。

在長時間數據集AMPds 上進行測試，分別選取2、4、6、8種負荷，待分解的時長為1～12個月，每天的數據點數為1 440。本文算法分解所需時間如附錄A圖A3所示。由圖可知，本文算法可以在短時間內完成負荷分解。為了分析本文算法在每個測試數據點的平均分解時間，將其與文獻［12］添加行為影響因子的HMM算法和文獻［15］的擴展稀疏HMM算法進行對比，結果如表5所示。

表5 平均分解時間Table 5 Average decomposition time

由表5可知，本文算法相比文獻［12］和文獻［15］算法在平均分解時間上均具有明顯優勢。此外，為了分析本文所提簡化HMM 與傳統HMM 在準確率和求解時間上的差異，選擇與表3中場景2相同的負荷種類及訓練時長，并采用相同的負荷功率模板，采用MATLAB 2021b 中集成的hmmestimate 函數和hmmviterbi 函數對模型進行訓練和求解，測試結果如附錄A 圖A4 所示。由圖可知，相較于傳統HMM，簡化HMM在長時間下的分解準確率更高，且在計算時間上具有顯著優勢。

5 結論

由于現有AMI 的低頻采集條件，利用低頻有功功率進行居民家庭NILM 更具可實施性。針對傳統算法對于分解結果的依賴性，本文提出一種在各時刻獨立求解的實時非侵入式負荷分解算法。本文算法在2 種公開數據集上的負荷功率分解準確率和有限狀態識別準確率均不低于90 %，并且算法的計算時間短，對含有8 種負荷場景的平均分解時間僅為0.07 ms。

筆者后續將進一步考慮不同用戶負荷類型以及用電行為的相似性，并進一步簡化模型參數的獲取過程，從而使得算法更具泛化性。

附錄見本刊網絡版（http：//www.epae.cn）。