非侵入式负荷分解实战:用FHMM从总电表数据拆分电器功率
简介面向从事非侵入式负荷分解NILM研究的学生、工程师与数据科学人员该压缩包提供基于因子隐马尔可夫模型FHMM的完整Python实现通过分析总用电量数据识别并分离各电器负荷适用于智能家居、能源管理与节能策略等场景。包内共有26个文件包括8个Python脚本、6个Jupyter Notebook、1个CSV数据文件以及配置文件等。Python脚本覆盖数据预处理、特征提取、模型构建、最大似然估计训练、负荷分解与结果评估等关键环节并附有主程序和可调用的辅助模块Notebook则便于分步运行、直观观察分解效果。目录中data存放原始电力数据models保存训练后的模型utils提供数据读取与特征工程函数整体层次清晰便于二次开发。资源包大小约628KB已有582人学习浏览。借助源码与演示数据可深入理解FHMM如何建模多电器的开关状态与功率消耗快速复现分解实验为后续研究或工程应用提供坚实基础。1. 非侵入式负荷分解总电表数据里怎么变成每台电器的功率一套两居室电表箱里只有一块总表。你想知道空调到底耗了多少电、冰箱每天有没有偷偷频繁启动却不能给每个回路装独立计量。非侵入式负荷分解NILM解决的就是这个从总功率信号里用算法把各个电器的功率和运行状态分离出来。因子隐马尔可夫模型FHMM是这类问题里历史最久、也最容易被低估的方案。它假设每台电器的运行状态是一条隐马尔可夫链所有链并行运行叠加起来就是你读到的总功率。FHMM参数少、解释性强、适合1 Hz甚至更低频率的数据在智能电表普及的背景下特别实用。适合谁做能耗分析的、做智能家居的、接智慧园区项目的开发者。这篇笔记给你一条可落地的复现路径包括数据怎么清洗、状态怎么定、训练怎么迭代以及那些不调必翻车的坑。2. FHMM为什么适合拆负荷模型拆解与三个选型理由2.1 因子隐马尔可夫模型在NILM里到底做了什么非侵入式负荷分解的核心是一个逆问题。你只观测到总功率Y_t而Y_t是屋里所有电器功率之和。如果用单个隐马尔可夫模型你得在状态空间里表达所有电器的组合这个组合状态数是各设备状态数的乘积。冰箱两种状态、空调四种、微波炉两种组合起来16个状态转移矩阵16×16还能接受但家里经常有8个目标电器组合状态动辄上万参数根本估不出来。因子隐马尔可夫模型从结构上化解了组合爆炸每台电器是独立的隐链状态s_t^(m)转移概率A_m独立观测方程只有一个Y_t等于各链发射的高斯随机变量之和。也就是说我们不再为一个超级状态估转移矩阵而是分别估M个小转移矩阵参数从指数级降到线性级。这个“独立隐链线性叠加观测”的结构正好对应家用电路并联的物理含义。从数学表达看FHMM的似然函数是P(Y_1:T | θ) ∑_{s_1} … ∑_{s_T} ∏_t N(Y_t | ∑_m μ_{m, s_t^m}, ∑_m σ²_{m, s_t^m}) × ∏_{m,t} A_m(s_t^m → s_{t1}^m)这里每个设备链的隐状态不用合并成大状态向量训练时也可以分别维护各自的转移矩阵。不过真正做前向-后向推理时为了计算方便很多人会把所有设备状态组合成一个超状态。假设设备m有K_m个状态联合状态空间最大是∏K_mFHMM的训练参数仍然是∑(K_m² K_m)而不是超状态的平方。这个差别在三五台设备时还不明显到8台以上时联合状态矩阵根本存不下。这也是FHMM在NILM里替代普通HMM的核心原因。2.2 三个选型理由低频采样、数据量小、可解释第一低频采样下依然有效。深度序列模型在30秒甚至分钟级数据上性能掉得很快因为时间上下文被压缩但FHMM本质上把设备建模成有限状态机低频只影响状态转移概率的分辨率不影响模型的表达。第二训练数据需求量小。FHMM的参数量大约是∑(K_m² K_m)一个3设备、每个2状态的小模型只有20个参数几十小时的子表计数据就能稳定收敛。相比之下LSTM或transformer类方法动辄上万参数很容易在小数据集上过拟合。第三可解释性拉满。训练完成后每个设备链的状态集合可以直接对应“关 / 低档 / 高档”均值就是该档下的平均功率。做能效报告的时候特别方便不用像深度模型那样再做一轮SHAP分析去解释黑匣子。我一开始做NILM也试过LSTM模型调得再好客户问“你这0.8 kW到底从哪来的”解释不了换FHMM后至少能说清。2.3 FHMM的边界不是所有电器都能拆干净FHMM最大的问题是发射模型假设各设备功率叠加后是一个高斯总和。实际上很多电器功率并不稳定变频空调功率随温差连续变化洗衣机滚筒启动瞬间有冲击功率而且同一设备在不同档位的方差也可能很大。高斯发射模型只能抓住一个中心值所以FHMM对“稳定多状态设备”拆得好对“连续变化设备”拆不好。另一个边界是同时工作的设备数量。联合状态空间随设备数指数增长8个状态数3以上的设备有6561个联合状态前向-后向每次迭代都要遍历耗时长而且状态间的可分性变差。因此我在做项目时通常把目标设备限制在5个以内优先选择冰箱、微波炉、热水壶这类启停边界清晰的负载。空调、取暖器这种大功率连续负载需要单独做一条回归通路不直接塞进FHMM。2.4 我用的落地路线独立HMM初始化联合状态推理严格训练FHMM需要对联合状态空间做EM下面会给代码。但工程上有个更稳的落地路线先用带子表计的数据分别训练每台设备的HMM用hmmlearn就能做再把训练好的参数填进FHMM的联合状态模型里最后在总功率观测上重新做一轮轻量微调。这样做的好处是子表计训练收敛快又能利用总功率修正设备之间的串扰。有些开源代码直接跳过了联合微调只用单表模型预测效果会差不少。我一般会至少跑50轮联合EM再收手。3. 数据准备从原始功率记录到FHMM的训练样本3.1 采样频率的选择1 Hz够用为什么不能迷信高频FHMM是状态转移模型对采样频率并不敏感。但采样频率决定了你能捕捉到的最小状态时长。常见做法是把高频数据先计算有功功率再降采样到1 Hz。1 Hz能分辨持续1秒以上的状态冰箱压缩机一次启动约5分钟微波炉30秒以上完全够用。降到0.1 Hz10秒一个点也可以但状态切换发生在采样间隔内时转移概率估计会被稀释分解准确率会下降。所以我一般默认1 Hz除非现场数据只有分钟级聚合。如果原始数据是更高频用resample(1s).mean()直接降采样即可。这里有一个取舍升采样到10 Hz没有意义FHMM的前向-后向计算时间会变长但信息量没有增加因为设备动态并没有那么快。3.2 数据清洗缺失值、跳点、负功率真实的总功率数据远没有论文里漂亮。最常见的问题是开始和结束时段有NaN、智能电表偶尔上报错误值、以及负功率光伏反送或电表接线问题。直接删除会造成时间轴断裂影响前后向算法的连续性。我一般会先把时间戳重采样到等间隔再对缺失段做线性插值但连续缺失超过10分钟就直接切断为独立序列不要强行拼接。对于单点跳变用中值滤波剔除超过3σ的点。下面是一段可用的清洗代码import pandas as pd import numpy as np def clean_power_series(s, freq1s, max_gap_minutes10, zscore_thresh3): # s 是带时间戳的 pandas Series s s.resample(freq).mean() # 记录连续缺失的分组 is_na s.isna() group (is_na ! is_na.shift()).cumsum() na_size s.groupby(group).transform(size) # 超过 max_gap_minutes 的缺失段不做插值保留 NaN max_gap_points max_gap_minutes * 60 # 在 1Hz 下 s s.where(~(is_na (na_size max_gap_points)), np.nan) # 短缺失用线性插值补上 s s.interpolate(methodlinear, limitmax_gap_points) # 3 sigma 去尖峰超过阈值的点设为 NaN 再插值 mu, sigma s.mean(), s.std() s s.where((s - mu).abs() zscore_thresh * sigma, np.nan).interpolate() return s这段代码里的max_gap_points在1 Hz下等于秒数如果你把freq改成10s要同步把换算改成max_gap_minutes * 6。我见过有人在10秒采样时还写* 60结果把10分钟缺失当成6小时缺失处理短缺失全部被切断。清洗后还有一个容易忽略的动作连续两个有效段之间如果隔了较长NaN模型会把它们当成连续时间序列训练前向算法在缺失处强行跳变。所以要在长缺失处手动拆分成多个子序列后续训练和推理都按子序列独立算。做法很简单判断is_na连续长度超阈值的地方用np.split或者直接按索引切片。3.3 子表计数据对齐训练FHMM必须有的“标准答案”FHMM训练时如果完全只有总功率是个一鸡生蛋问题你不知道状态序列也没法估发射参数。常见做法是使用带子表计的公开数据集如REDD、UK-DALE训练或者自己用智能插座采集几天的设备功率。子表计数据也要与总功率在时间上对齐。我踩过的一个大坑是时区没统一总表和子表一个UTC一个本地时间导致状态与功率错位。建议统一转成UTC时间戳并严格对齐到相同秒网格。下面是对齐代码def align_submeters(total, sub_meter: dict, freq1s): # total 是总功率 Seriessub_meter 是 {设备名: Series} df pd.DataFrame({total: total}) for name, ser in sub_meter.items(): df df.join(ser.rename(name), howleft) # 统一重采样到相同频率 df df.resample(freq).mean() # 总功率缺失的时段直接丢弃因为主要观测必须可靠 df df.dropna(subset[total]) # 子表计缺失的极短段用 0 填充 df df.fillna(0.0) return df参数说明sub_meter是一个字典键是设备名值是带时间戳的功率Series。如果用join而不是merge能保持总功率的时间索引不动避免索引错位。填充0是为了让训练样本完整但前提是子表计缺失时间很短。如果某个设备缺失了半小时这段数据实际上可能是开着但没记录填0会告诉模型“设备关机”污染发射均值。所以我会额外打印每个设备缺失占比超过5%就缩短训练时段或者干脆剔除该设备。3.4 设备筛选与状态粗定义拿到对齐数据后不要急着一股脑把8台设备全塞进FHMM。先用可视化或简单统计筛选。第一个标准是状态数不超过4看子表计功率直方图有多少个峰峰数就是初始K_m。第二个标准是开关状态功率差要明显待机功率5W和运行功率1000W的冰箱就很好拆但两个档位功率差只有20W的净化器就很难拆因为总功率里还有别的设备噪声。第三个标准是设备启停次数不能太少否则转移矩阵只有对角线有值前向-后向算法退化成静态分类。筛选后为每台设备记录一个初始状态列表状态0关状态1低档状态2高档等。这部分工作直接决定后续EM收敛速度和分解精度我一般会在这里花掉整个项目40%的时间。没有干净的原始数据再好的模型也是白搭。4. 用Python把FHMM展开成联合HMM训练核心代码与参数设置4.1 初始化从功率直方图或KMeans生成设备状态进入联合EM前先用子表计数据独立初始化每台设备的发射参数。最简单的办法是对每台设备的功率序列做KMeans聚类簇中心就是状态均值簇内标准差就是发射方差。下面的函数返回均值和一组自环偏大的转移矩阵from sklearn.cluster import KMeans import numpy as np def init_from_power_series(power, n_states, random_state0): x power.values.reshape(-1, 1) km KMeans(n_clustersn_states, random_staterandom_state).fit(x) means np.sort(km.cluster_centers_.ravel()) stds [] for k in range(n_states): vals power.values[km.labels_ k] stds.append(vals.std() 1e-6) stds np.array(stds) # 转移矩阵初始化自环 0.6其余均分 trans np.full((n_states, n_states), 0.2) np.fill_diagonal(trans, 0.6) trans / trans.sum(axis1, keepdimsTrue) return means, stds, trans为什么转移矩阵自环要偏大因为FHMM里设备一旦开启通常会持续一段时间比如冰箱压缩机一次运行1020分钟如果初始化为均匀分布模型会把设备切成每秒钟随机闪烁。这里n_states是你自己定的我通常先画直方图数峰峰不明显的就取2宁可漏档位也不要引入一个过拟合的中间状态。4.2 组装联合超状态空间与发射参数FHMM要落地最常见的工程做法是把所有设备状态组合成超状态再用标准高斯HMM做训练和推理。虽然联合状态数是指数级的但在设备数不超过5、每台状态数不超过4的场景下组合数最多1024完全能放进内存。下面的FHMM类负责构造联合超状态from itertools import product class FHMM: def __init__(self, means_list, stds_list, trans_list, starts_list): # means_list[i] 是第 i 台设备各状态的均值数组 self.M len(means_list) self.K [len(m) for m in means_list] # 所有设备状态组合 self.combos list(product(*[range(k) for k in self.K])) self.n_combo len(self.combos) # 每个超状态的发射均值与方差高斯叠加 self.mu np.array([ sum(means_list[m][s[m]] for m in range(self.M)) for s in self.combos ]) self.var np.array([ sum(stds_list[m][s[m]]**2 for m in range(self.M)) for s in self.combos ]) # 超状态间的转移矩阵各设备转移概率相乘 self.A np.ones((self.n_combo, self.n_combo)) for i, si in enumerate(self.combos): for j, sj in enumerate(self.combos): for m in range(self.M): self.A[i, j] * trans_list[m][si[m], sj[m]] # 初始分布同理 self.pi np.array([ np.prod([starts_list[m][s[m]] for m in range(self.M)]) for s in self.combos ])这段代码里有几个值得说的点。combos保存的是每个超状态对应的设备状态元组比如(0, 1, 0)表示设备1关、设备2低档、设备3关。mu和var是总功率下该组合的期望均值和方差因为独立高斯变量的和仍然是高斯。A是超状态转移矩阵形状是(n_combo, n_combo)它没有直接学因子化参数但初始值来自各设备转移矩阵的乘积已经带上了“设备独立性”的先验。如果设备数继续增大这里会发生组合爆炸所以我在类里加了self.n_combo方便你提前判断规模。4.3 用hmmlearn训练联合HMM搭建好联合空间后直接用hmmlearn的GaussianHMM来处理前向后向和EM更新。这里的关键是手动把初始参数塞进模型再调用fit。观测必须是二维数组形状为(T, 1)代表每个时刻的总有功功率。from hmmlearn.hmm import GaussianHMM def train_fhmm(fhmm: FHMM, total_power: np.ndarray, n_iter50, tol1e-3, random_state42): # total_power 必须是 (T, 1) 的二维数组 X total_power.reshape(-1, 1) model GaussianHMM( n_componentsfhmm.n_combo, covariance_typediag, n_itern_iter, toltol, random_staterandom_state ) # 手动塞入初始化参数 model.startprob_ fhmm.pi model.transmat_ fhmm.A model.means_ fhmm.mu.reshape(-1, 1) model.covars_ fhmm.var.reshape(-1, 1) model.fit(X) return modeln_iter和tol是训练收敛的两个关键参数。我一般先设n_iter50如果model.monitor_.converged为False再加到100。这里有一个常见误解n_iter不是越大越好FHMM在高斯方差比较小的时候EM很容易在几十轮内就把转移矩阵推成“全0或全1”的极端形态继续迭代反而过拟合。tol控制对数似然增量设1e-3是经验值太小会多跑很多轮但精度不提升。random_state一定要固定否则不同运行结果可能差异巨大。如果你不想用hmmlearn也可以自己写前向-后向但要注意使用scipy.special.logsumexp防止数值下溢。总功率的联合状态数量上千时直接算概率密度会在长序列上乘出一个趋近0的数日志空间是必须的。自写版本的优势是可以直接对因子参数做投影更新但代码复杂度高对付实际项目用hmmlearn足够了。4.4 分解解码超状态并映射到设备功率训练完成后用Viterbi解码出每个时刻的超状态索引再通过combos映射回每台设备的状态。得到状态后用各设备在该状态下的均值作为该时刻的分解功率。def decompose(model: GaussianHMM, fhmm: FHMM, total_power: np.ndarray): X total_power.reshape(-1, 1) hidden model.decode(X, algorithmviterbi)[1] # hidden 里是超状态编号映射回 (T, M) 的设备状态矩阵 device_states np.array(fhmm.combos)[hidden] # shape (T, M) M fhmm.M power_est np.zeros((len(hidden), M)) for m in range(M): power_est[:, m] fhmm.mu_by_device[m][device_states[:, m]] return device_states, power_est注意这里还需要在FHMM类里保存一份mu_by_device否则解不出来。我在实际代码里初始化时顺便保存了原始的means_list和stds_list便于最后按设备取均值。Viterbi给出的是全局最优状态路径适合离线分析如果你想做实时分解应该改用predict_proba取每个时刻后验概率最高的状态再加权得出功率估计。5. 避坑与排查FHMM负荷分解最常见的5个翻车现场5.1 现象模型把所有负荷都分到同一台设备上我刚开始做联邦表项目时训练完FHMM分解结果里冰箱那一列几乎等于总功率其他设备全是0。查了三天原因不是模型结构而是初始化。我用均匀分布初始化所有设备的发射均值结果EM陷入了一个局部最优那台设备的状态数多后验概率把大部分观测都解释掉了。解决方法是改用第4.1节的KMeans初始化并且把每台设备的状态均值严格限制在子表计数据的真实范围里不要让模型自由漂移太远。另外如果有一台设备的初始方差设得特别大它也容易吞掉所有负荷所以发射方差的下限要设为正常设备功率波动的1/10。5.2 现象训练Loss不降反升hmmlearn打印的monitor_.history里对数似然在训练中途掉头向上通常不是模型问题而是数据没有分段。联合状态HMM的前向-后向在长序列上会有概率值下溢到0虽然用了日志空间但序列里如果有长时间缺失段还是会算成NaN。你会在训练日志里看到一条NaN。解决方法是训练前按3.2节切分长缺失段并且把每个子序列长度控制在12小时以内。另外重新检查一下total_power里是不是有绝对值极大的尖峰一个错误值会让发射概率的方差被拉到100倍之后的log-likelihood当然不稳定。5.3 现象分解出的设备功率序列有锯齿状抖动正常冰箱运行一次至少5分钟但分解结果里它每隔几秒就开关一次锯齿非常明显。原因是转移矩阵里“保持原状态”的概率太低。FHMM的转移矩阵如果不加约束会更倾向于用频繁切换来解释总功率变化因为这样能在每个时刻都拟合得更准。解决思路有两个。一是初始化时把自环概率调到0.7以上让模型在训练前期就看到“设备不会秒开关”的先验。二是后处理加中值滤波窗口长度取3060个采样点形态学上消除单点毛刺。但中值滤波只能救场景不能救参数真正有效的还是调转移矩阵的初始化。5.4 现象预测的启停时刻和实际偏差很大有的设备待机功率接近0比如手机充电器待机0.5W而总功率波动有几十瓦模型很容易忽略这个状态。反过来如果某台设备实际没开但总功率里有一台微波炉在加热微波炉的发射均值可能被分配给了另一台同功率等级的电器导致错位。我的排查顺序是先看每个设备子表计实测的功率直方图和模型学到的发射均值是否一致偏差超过20%就重新检查KMeans聚类时的状态数。很多时候你漏掉了设备的一个待机档位模型就会用另一台设备去补。5.5 现象相同参数不同初始化结果差异巨大FHMM的似然面不是凸的不同的随机初始化会收敛到完全不同的分解方案。这不是代码bug而是模型的固有特性。我的习惯是用5个不同的随机种子各训练一遍取对数似然最高的那个模型同时检查每台设备的发射均值是否落在合理范围内。如果5个结果差异很大说明数据里设备重叠时段太多应该减少同时建模的设备数量。不要试图用单个模型跑一次就下结论FHMM不是那种“一劳永逸”的模型多跑几个种子才能看出稳定性。6. 验证与进阶把分解精度抠到能落地6.1 验证指标与后处理组合模型训练完不能只看分解曲线对不对要用三个指标验证。第一个是状态级F1把每台设备的分档预测和真实子表计分档结果对比F1要超过0.85才算能用于能效报告。第二个是功率分解误差用估计的设备功率均值与子表计实际功率逐样本求RMSE除以设备的平均功率得到归一化误差控制在20%以内。第三个是能量占比误差统计一段时间的累积电量FHMM如果长期偏移后面两个指标再好都白搭。后处理我固定做两步第一步对设备状态序列施加最小开关持续时间约束比如冰箱最短持续开启5分钟把短于这个时间的孤立脉冲直接反转到前一状态。第二步用设备功率均值乘以状态序列得到功率估计之后做一个30秒滑动平均消除Viterbi路径在状态切换瞬间的毛刺。这两步能把F1从0.7拉到0.9代价只是让分解结果延迟30秒离线分析完全能接受。6.2 一个实用的进阶技巧用少量插座校准FHMM如果你只训练一次就交付第二个星期可能就跑偏了。电器是活的冰箱制冷周期会随季节变化空调的功率档位也会因为电压波动漂移。我现在的习惯是保留一个低成本校准流程在关键设备上临时装一个智能插座采3天空调或冰箱的功率然后用这3天数据重新初始化该设备的KMeans参数冻结其他设备的参数对FHMM单独做20轮微调。这样做的计算量很小但能让分解精度在季节交替时不掉下去。FHMM的优势恰恰在这里——参数有物理含义你可以只动一台设备的状态均值不必把整个网络重新训练一遍。做这件事时要注意校准数据里不要包含全屋总功率的异常大事件否则那台设备的方差会被一次微波炉加热拉大。我在微调空调参数时会先把总功率里超过空调满负荷1.5倍的点全部剔除只保留空调单独运行或待机的时段。这个习惯帮我避过无数次“空调开一整夜但分解出来只有2小时”的现场事故。6.3 收尾FHMM不是最花哨的NILM方法但它的可复现性和可维护性在工程现场很难替代。你只需要一个稳定的总功率序列、几路子表计数据、几十行初始化代码就能跑出一套能解释的负荷分解结果。踩过的那些坑最后都变成了项目交付前的检查清单先看数据缺失再看状态初始化然后固定随机种子最后用F1和误差率收口。希望这些经验能帮你少走一段弯路把更多时间花在看设备行为、调业务策略上。希望帮到你。本文还有配套的精品资源点击获取