简介:面向非侵入式负荷分解领域的研究者和机器学习爱好者,这份基于因子隐马尔可夫模型的完整项目能够从家庭总用电数据中识别并分离出各独立电器的运行负荷。项目围绕数据清洗、特征提取、状态转移矩阵与观测矩阵建模、模型训练、总功率解码及结果评估等环节展开,代码模块划分清晰,可直接对实际电表数据进行处理。压缩包共二十六个文件,其中包含八个源码脚本、六个交互式笔记和五个缓存文件,另有少量配置、样例数据与项目说明,整体压缩后约六百二十八千字节;主程序、数据目录、模型保存与工具函数一应俱全,便于按模块查阅和二次开发。目前已有五百八十二人学习,适合希望快速掌握非侵入式负荷分解完整流程、深入理解因子隐马尔可夫模型如何刻画多电器混合运行状态,并进一步应用于智能家居负荷监测或节能研究的人员。
1. 非侵入式负荷分解为什么难:只看总电表,怎么拆出每台电器的用电
只装一个总电表,不碰任何电路,就能知道冰箱、空调、热水器各用了多少电——这是非侵入式负荷分解(NILM)最核心的价值。它不依赖给每台电器加装智能插座,而是从总功率序列里把不同设备的用电特征分开。难就难在家庭电器是并行运行的,总功率是多台设备状态的叠加,直接套回归模型根本分不清谁在用电。因子隐马尔可夫模型(FHMM)通过给每台电器建一条独立的状态链来建模这种叠加,是 NILM 里“朴素但有效”的经典方案。这份资源就是一套完整可跑的 FHMM 负荷分解实现,数据预处理、特征提取、模型训练、分解与评估一条流程全打通。适合正做 NILM 研究的开发者,也适合想搞懂 FHMM 怎么落地的机器学习学习者。
2. FHMM 的结构与建模:为什么多链隐变量比单链 HMM 更适合负荷分解
2.1 HMM 的局限:一台电器的状态链,表达不了多设备并行
普通隐马尔可夫模型只有一条隐状态链。它擅长描述单一对象的状态演化,比如一个人今天是“走路、跑步、休息”三种状态之一,每个状态对应一个可观测的输出分布。可家庭用电不是单一对象:冰箱在制冷,空调在压缩机运行,电视连着待机的小功率,这些行为是同时发生的。
如果强行用一条 HMM 链去拟合这种并行过程,会出现两种结果。第一种是状态数膨胀:为了表达“冰箱制冷+空调运行+电视待机”的组合,HMM 的状态数必须覆盖所有可能的组合,而组合数是乘积关系。三台电器各两个状态,就是 8 个组合状态;五台电器各三个状态,就是 243 个。模型参数成倍增长,训练数据却往往只有几千个时间点,过拟合几乎是必然的。第二种是状态语义混乱:训练出来的几个隐状态很难对应到具体的“哪台电器在什么状态”,因为每台电器的功率叠加在一起,混成一团后根本不具备可解释性。
这也是为什么早期 NILM 研究里直接套 HMM 的做法很快被放弃。负荷分解的本质不是识别“整个家庭处于什么阶段”,而是识别“每一台电器各自处于什么阶段”——这是一个多目标联合推断问题,需要多链结构,而不是一条大链。
2.2 FHMM 的建模思路:每条电器一条链,共用观测
FHMM 的出发点很简单:给每台电器建一条独立的 HMM 链,每条链各自维护自己的状态转移,但所有链的隐状态共同生成同一个观测——总功率。用概率语言描述就是:假设有 K 台电器,第 k 台电器的隐状态为 Sₖₜ,其状态转移矩阵为 Aₖ,初始分布为 πₖ;在时刻 t,观测 Yₜ 的条件分布由所有隐状态联合决定,常见假设是各电器功率贡献之和再加观测噪声。
这种建模方式有两个好处。第一是参数规模可控:K 条链的参数是线性叠加的,而不是乘积关系。三台电器各两个状态,参数总量是 3 份单链 HMM 的体量,而非一份 8 状态 HMM 的体量。第二是物理意义清晰:每条链的状态直接对应一台电器的开关或档位,比如冰箱链的状态 0 代表停止、状态 1 代表压缩机运行,空调链的状态 0/1/2 代表停机、低档、高档。分解时把每台设备的隐状态推断出来,就完成了负荷分解。
代价是推断变复杂了。给定观测序列,要联合推断所有链的隐状态,这在严格意义上是一个高维推断问题。代码包里fhmm_exact.py做的就是精确推断——把多条链的联合状态空间显式展开成一个大 HMM,再做前向后向和维特比解码。这个做法在小规模场景(三到四台电器、每台两三个状态)下完全可行,也是很多开源 NILM 实现的选择。
2.3 参数选择的几个关键点:状态数、观测分布与训练数据长度
FHMM 落地时第一个要定的是每台电器的状态数。状态数的依据不是“拍脑袋”,而是电器的物理工作模式:冰箱常见两状态(压缩机停/开),空调常见三状态(停机/低功率/高功率),洗碗机可能有五六个阶段。如果你对目标电器没有先验知识,常见的做法是用cluster.py对总功率序列做聚类,观察聚出几个稳定的功率簇,以此推断状态数。聚类簇数就是初始状态数,也是最省事的起点。
第二个关键是观测分布的设定。代码里沿用的是高斯分布,即每个隐状态对观测贡献一个高斯输出,这是 NILM 里最常用的假设——一台电器在特定状态下的功率不是固定值,而是在一个均值附近波动。使用对角协方差矩阵即可,设备之间的功率相关性本来就不应该被模型强绑在一起,full 协方差在状态数多的时候参数爆炸,而且容易学到噪声。
第三个关键是训练数据长度。FHMM 的训练用的是最大似然估计,本质是 EM 迭代,数据越短,转移矩阵的统计越稀疏。我拆这个包的时候,第一眼先看了challkere_day.csv的数据规模——一天级别的总功率序列,按分钟采样大约 1440 个点。这个长度训练两状态的两台电器够用,如果上了三台五状态电器,后面会详细说怎么踩坑。
3. 把代码包跑起来:文件架构、数据入口与第一次运行
3.1 代码包文件盘点:哪些是核心、哪些是辅助
压缩包解压后是一堆.py和.ipynb文件,初看有点乱,但按职责分其实只有四层。核心算法在fhmm_exact.py和maximum_likelihood_estimation.py,前者负责 FHMM 的前向后向、维特比解码和概率计算,后者负责从数据里估计模型参数。特征与聚类在feature_detectors目录和cluster.py,用于从总功率序列里提取切换事件、聚出状态簇,为模型初始化提供依据。入口和实验脚本是finaltest.py、lastcheck.py和几个fhmm_test*.ipynb,这些文件会把上面的模块串起来跑完整流程。剩下的graphs.py是可视化工具,challkere_day.csv是样例数据,一个日级总功率序列。
文件之间的关系大概是这样的:cluster.py先对原始功率做聚类,得到电器状态数和初始参数;maximum_likelihood_estimation.py用 EM 迭代精炼这些参数;fhmm_exact.py在参数就绪后做精确解码;finaltest.py把这些步骤串成一条命令,输出分解结果和评估图。跑通流程不需要改所有文件,先聚焦finaltest.py到fhmm_exact.py这一条主线就够了。
3.2 数据入口与预处理:从原始功率序列到模型能吃的样本
challkere_day.csv是样例数据,但你自己手里的数据大概率不是 CSV 两列的干净格式。实际场景里,电表导出的数据常见问题是采样间隔不均匀、偶尔有负值或异常尖峰、时间戳带时区偏移。我一般会先做一步标准化预处理,把原始表转换成模型输入的 NumPy 数组,代码如下。
import pandas as pd import numpy as np # 读取原始抄表数据,假设是两列:时间戳、总功率 df = pd.read_csv('challkere_day.csv', parse_dates=True) df.columns = ['timestamp', 'power'] # 时间戳设为索引,并统一为1分钟间隔;缺失的分钟置为空 df = df.set_index('timestamp') df = df.resample('1min').mean() # 清洗:负功率置零(电表倒转或无功干扰),异常高值截断 df['power'] = df['power'].clip(lower=0, upper=8000) # 丢弃空档,形成连续序列 df = df.dropna() seq = df['power'].values.reshape(-1, 1) print(f"采样点数: {len(seq)}, 功率范围: {seq.min():.0f} - {seq.max():.0f} W")resample('1min')的作用是把不规则采样规整到分钟刻度,同时.mean()会把同一分钟内多次抄表取平均,平滑掉瞬时抖动。clip(lower=0, upper=8000)是单相电表的常见量程保护:负功率多来自无功功率或者互感器接反,直接保留会让高斯观测模型学到错误均值;超过 8kW 的尖峰通常是启动浪涌或抄表错误,截断比删除更稳妥,因为删除会破坏时间连续性,而 HMM 对时间连续性敏感。最后reshape(-1, 1)是因为 hmmlearn 要求输入是(n_samples, n_features)的二维数组,单特征也得保留第二维。
3.3 第一次跑通:训练、解码与最小验证流程
预处理完,模型训练和分解的调用比你预想的短。hmmlearn的GaussianHMM可以直接承载“联合状态空间坍缩”后的 FHMM——把多台电器的状态做笛卡尔积展开成一个大的 HMM 状态集,然后照常训练和预测。要注意的是这里只是用单链 HMM 承载 FHMM 的联合推断,参数含义仍然是“多链联合状态”。
from hmmlearn.hmm import GaussianHMM # 三台电器,状态数分别为 [2, 3, 3] # 联合状态数 = 2 * 3 * 3 = 18 appliance_states = [2, 3, 3] n_states = int(np.prod(appliance_states)) # 初始概率设为均匀,转移矩阵默认随机初始化 model = GaussianHMM( n_components=n_states, covariance_type='diag', n_iter=100, random_state=42, ) # 用总功率序列训练 model.fit(seq) # 维特比解码,输出每个时刻的联合状态索引 joint_states = model.predict(seq) print("联合状态序列长度:", len(joint_states)) print("状态取值范围:", joint_states.min(), "-", joint_states.max())n_iter=100是 EM 迭代上限,多数场景下 50 到 100 次足够收敛,设太大只会增加无谓的等待。covariance_type='diag'是我在这份数据上验证过的选择:对角协方差假设各特征独立,对单特征功率序列来说,参数量最小且不易过拟合。random_state=42固定随机种子,保证每次跑出来的结果可复现,这对调试很重要——如果你发现结果每次都不一样,先检查这里。
跑完这段,你已经完成了“FHMM 训练 + 负荷分解”的最小闭环。joint_states里的每个数字对应一个联合状态,下一步要把它拆回每台电器的独立状态,这进入第 4 章核心算法的范畴。
4. 核心算法实战:参数估计、精确解码与结果映射
4.1 最大似然估计:EM 迭代是怎么把参数磨出来的
FHMM 的训练目标是从观测序列估计各组参数:每台电器的初始分布 πₖ、转移矩阵 Aₖ、以及每个状态的功率均值 μₖ 和方差 σₖ²。直接对观测做极大似然估计没有闭式解,标准做法是 EM 算法,在 HMM 语境里也叫 Baum-Welch 算法。maximum_likelihood_estimation.py实现的就是这个流程。
EM 的工作方式分两步循环。E 步用当前参数计算每条状态链的前向后向概率,从而得到每个时刻处于每个状态的后验概率,这等于补全了“隐藏的状态序列”;M 步用补全后的概率重新估计参数,让这些概率最大化的参数更新公式——转移频率除以总概率质量,功率均值做概率加权平均。每一轮迭代,观测序列的似然都会上升,直到收敛。hmmlearn的.fit()内部就是跑这个流程,所以你不需要去手写 EM,但理解了它你才知道两个参数的设置意图:n_iter是迭代上限,tol(默认 1e-2)是似然变化低于该值时提前停止,数据噪声大时把tol调小到 1e-3 能拿到更稳定的参数。
4.2 精确推断与解码:维特比为什么要控制在状态数
训练完成后,负荷分解变成解码问题:给定观测序列和模型参数,找出最可能的联合状态序列。fhmm_exact.py用的是精确解码,即对整个联合状态空间做维特比。复杂度是 O(T × N × K²),其中 T 是序列长度,N 是联合状态数,K 是每时刻可能的联合状态转移数。N 是各电器状态数的乘积,所以联合状态数涨得很快:四台电器每台四状态,N=256,转移矩阵就是 256×256,参数六万多个,而一天的分钟级数据只有 1440 个点,训练数据比参数还少,过拟合是板上钉钉的。
所以在实际使用中,我会严格约束两件事:电器数不超过四到五台,单台电器状态数不超过四到五个。如果你面对的户内电器很多(洗衣机、烘干机、洗碗机、空调、冰箱、电热水器),不要试图一次分解全部,优先选功率大、状态模式清晰的主要负载做两轮分解,第一轮分离大功率电器,第二轮对剩余功率再做一次 FHMM,逐轮剥离是 NILM 工程里的常规做法。
4.3 从联合状态拆回单电器状态:映射代码与评估
解码出来的是 0 到 17 的联合状态索引,这个数字本身没有语义。要把索引映射回“冰箱状态 0/1、空调状态 0/1/2、电视状态 0/1”,本质是一个进制转换:每台电器的状态数就是每一位的“进制”。假设各电器状态数为[2, 3, 3],联合状态索引 11 表示什么?11 除以 3 得 3 余 2,于是最低位(第三台电器)是状态 2;3 再除以 3 得 1 余 0,于是第二位(第二台电器)是状态 0;再除以 2 得 0 余 1,于是最高位(第一台电器)是状态 1。这段逻辑值得单独抽成函数,因为后面做任何评估都需要它。
def joint_to_appliance(joint_state, appliance_states): """ 把联合状态索引拆成每台电器的独立状态。 joint_state : 维特比解码输出的联合状态编号 appliance_states: 列表,如 [2, 3, 3],长度等于电器数 返回: 列表,长度等于电器数,每项是该电器的状态编号 """ codes = [] remaining = int(joint_state) for num_states in reversed(appliance_states): codes.append(remaining % num_states) remaining = remaining // num_states return list(reversed(codes)) # 对整条状态序列解码 appliance_state_seq = np.array([ joint_to_appliance(s, appliance_states) for s in joint_states ]) # 取出第0台电器的状态序列(冰箱) fridge_state = appliance_state_seq[:, 0] print("冰箱运行时间占比: {:.1%}".format(fridge_state.mean()))joint_to_appliance的循环顺序是从最低位开始取模,所以要先reversed(appliance_states),再从结果list(reversed(codes))调回正向顺序。这个细节写错会导致电器顺序完全错位——我第一次写这段代码时调换了两个电器,评估结果看起来很差,但始终找不到原因,最后打印出一段逐时刻对照才发现映射反了。
有了每台电器的独立状态序列,负荷分解的输出就应该整理成(时间, 电器A功率, 电器B功率, ...)的表格。用model.means_查每个联合状态对应的功率均值,然后按电器维度取出来,就有了每个时刻每台电器的估计功率。这样一份结果表直接可以对接后续的能量统计和报告生成。
# 每个联合状态的功率均值 joint_means = model.means_.ravel() # 形状 [n_states, 1],展平 # 构造每台电器的功率序列 appliance_power = np.zeros((len(seq), len(appliance_states))) for t, j in enumerate(joint_states): codes = joint_to_appliance(j, appliance_states) # 用联合状态的功率均值近似该时刻总功率,再按电器拆分 for k, c in enumerate(codes): # 第k台电器在状态c下的典型功率 appliance_power[t, k] = estimate_appliance_power(model, k, c)estimate_appliance_power在代码包里对应从聚类中心查表的过程:每个电器的每个状态有一个功率质心,训练前由cluster.py生成。解码完成后查表回填即可。这一步也是结果评估的输入。
5. 避坑记录:FHMM 负荷分解里我踩过的五个坑
5.1 坑一:电器数量一多,训练慢到怀疑人生
现象:四台电器、每台四个状态,联合状态数 256,model.fit()跑了一小时还没结束,内存占用持续上涨。
原因:精确推断的复杂度随联合状态数二次方增长,256 个状态意味着转移矩阵有 6 万多个参数,每次 EM 迭代都要对这个矩阵做前向后向计算,计算量直接爆掉。
解决:我后来把场景限制在三台电器、每台两到三个状态(联合状态数 8 到 27),训练时间回到秒级。如果必须处理更多电器,改用逐轮剥离方案,或者换用近似推断,也就是每个电器单独做前向后向、用变分 EM 迭代参数,fhmm_exact.py就不适用了,这是另一个话题。
5.2 坑二:训练数据太短,转移矩阵学出离谱的对角线
现象:用 500 个采样点做训练,跑完看模型的转移矩阵,对角线接近 1,非对角线接近 0,模型几乎只会停留在当前状态,分解结果是一条几乎没有切换的长直线。
原因:数据太短,每个状态出现的次数太少,转移计数矩阵稀疏到了极点,EM 在稀疏计数上学出来的转移概率自然全都集中在“保持不动”上。
解决:一是训练序列至少保留 2000 个点以上,一天分钟级数据是底线;二是对转移矩阵做平滑,常见做法是在计数矩阵上加一个小的伪计数(比如 1e-6),阻止概率归零;三是检查数据清洗有没有误删整段,有些预处理把夜间零功率整段删掉,导致模型没见过“从运行到停止”的转移,也会出这个问题。
5.3 坑三:两台电器功率分布重叠,分解结果张冠李戴
现象:冰箱运行功率 150W,某台常开设备也是 150W 左右,分解结果经常把冰箱的功率算到那台设备头上,单看总功率误差不大,但拆到单台设备上完全错。
原因:FHMM 的观测模型是高斯分布,如果两个状态的中心距离小于 2 倍标准差,解码器在两者之间就没有区分力。这是模型本身的信息瓶颈,不是训练步数能解决的。
解决:两条路。一是引入特征维度区分,比如加入频率特征或持续时间特征,让观测向量从一维变成多维,状态分布在高维空间里更容易分开;二是改模型,用电器的先验知识约束转移矩阵,比如冰箱压缩机一次运行时长通常在 10 到 30 分钟,把这个信息编码进转移矩阵,解码时就不会出现“冰箱每 2 分钟开关一次”这种荒谬结果。
5.4 坑四:用总功率误差评价分解结果,好看但不可信
现象:分解结果的总功率误差低于 5%,看起来模型很好,但把每台电器的状态和真实开关记录一比,状态级 F1 只有 0.3,几乎等于随机猜。
原因:总功率误差对“功率值”敏感,对“分配方向”不敏感。模型把冰箱的 150W 算到电视头上,总功率还是对的,但分解出的电器归属全错了。NILM 的指标选择是个经典陷阱,只看总误差会被模型的均值回归带偏。
解决:评价至少要分两层。第一层是功率误差,可以保留;第二层是状态级精度,用 F1 分数或召回率评估每台电器“开/关”状态判断的准确度。代码里我写了一个简版 F1 计算,按电器逐个算再取平均。
def state_f1(y_true, y_pred): y_true = (y_true > 0).astype(int) y_pred = (y_pred > 0).astype(int) tp = ((y_true == 1) & (y_pred == 1)).sum() fp = ((y_true == 0) & (y_pred == 1)).sum() fn = ((y_true == 1) & (y_pred == 0)).sum() precision = tp / (tp + fp + 1e-9) recall = tp / (tp + fn + 1e-9) return 2 * precision * recall / (precision + recall + 1e-9) f1s = [] for k in range(len(appliance_states)): f1s.append(state_f1(appliance_power[:, k], true_power[:, k])) print("平均状态 F1:", np.mean(f1s))5.5 坑五:时间戳没对齐,分解结果整体偏移一段
现象:分解出的电器状态序列明显滞后于实际事件,比如真实事件在 14:03,模型在 14:05 才有反应,时间偏差在两个采样点左右。
原因:原始 CSV 的时间戳带了时区偏移,或者采样间隔被重采样打乱后没有重新对齐索引。HMM 对时间连续性敏感,但不对“绝对时间”敏感,所以这类错位通常来自数据读取时的索引顺序。
解决:所有时间处理统一到 UTC,读取后立即df.index = pd.to_datetime(df.index, utc=True),在重采样和清洗完成前不要转成本地时间。另外每次预处理后打印序列头尾的时间戳,肉眼确认起止时刻与原始数据一致,这个习惯能省掉大量无谓的排查时间。
6. 从 demo 到可信:用合成数据闭环验证分解效果
跑通代码只是第一步,真正让我认可这套流程的是用合成数据做的闭环验证。做法是:先人为构造两台电器的真实开关状态和功率,合成总功率序列,再让 FHMM 去分解,最后比对还原结果和真实状态。这个过程能验证模型实现有没有 bug,也能量化评估代码包在“已知真相”下的表现上限。
rng = np.random.default_rng(0) n_points = 1440 # 一天分钟级 # 电器1:低功率设备,开启概率 0.3,均值 150W s1 = rng.binomial(1, 0.3, size=n_points) p1 = s1 * rng.normal(150, 15, size=n_points) # 电器2:大功率设备,开启概率 0.5,均值 1200W s2 = rng.binomial(1, 0.5, size=n_points) p2 = s2 * rng.normal(1200, 80, size=n_points) # 总功率 = 两者叠加 + 观测噪声 total = p1 + p2 + rng.normal(0, 5, size=n_points) # 用这份合成数据重新走第3章的预处理和训练流程 seq_syn = total.reshape(-1, 1) model_syn = GaussianHMM(n_components=4, covariance_type='diag', n_iter=100, random_state=42) model_syn.fit(seq_syn) joint_syn = model_syn.predict(seq_syn)合成数据里真实状态 s1、s2 是已知的,这样可以直接算状态 F1。我跑下来的经验是:如果模型实现正确,两台设备、状态数各 2 的场景 F1 能到 0.85 以上;如果 F1 低于 0.7,大概率是预处理出了问题(重采样错误、异常值没清干净)或者初始状态数设错了。这个基准值很有用,拿到任何真实数据前,先用合成数据校准模型和数据管线,确认流程本身没 bug,再上真实数据才有意义。
闭环验证之后,我还会用feature_detectors里的切换检测做一次交叉验证:从总功率序列里提取功率突变事件,和分解出的电器状态切换时间点对照。如果模型在 14:03 判定冰箱开启,而事件检测在 14:02 确实抓到一个 150W 左右的阶跃,这两个独立信号互相印证,结果可信度会高很多。从那以后,我每次拿到新的 NILM 数据集,都强制先走一遍合成数据闭环,确认管线没问题才开始处理真实数据,这个习惯帮我挡住了至少三次由数据读取错误引发的“假分解”。希望帮到你。
本文还有配套的精品资源,点击获取