负荷监测是一种节能技术, 可以给消费者和电网公司提供住宅内每个电器的详细用电信息. 该技术按照实现方式被分为侵入式和非侵入式. 1992年, Hart以稳态有功功率和无功功率为特征, 提出了最初的NILM框架^[1], 相对于侵入式负荷监测, NILM采用家庭的总功率水平作为输入来估计每个设备对总能耗的贡献程度^[2], 大大降低了对用户生活的影响以及装设成本.

为提高能源利用效率和智能电网的发展^[3], 响应科学、合理用电的需求, 越来越多的研究者致力于提高NILM的准确性和适用性, 由此衍生出诸多模型与算法. 通过使用NILM技术, 电力部门可以为用户提供详细的电力使用情况和个性化的节能建议, 保证电网的合理运行. 因此, NILM是一种高效且经济的能源消耗技术框架, 对电力用户、电网和整个社会都具有重要意义^{[4, 5]}. 负荷分解算法可以分为数学优化方法^{[6, 7]}和机器学习方法^[8]. 前者尝试将观察到的功率测量值与 (设备特征库中) 设备功率信号的可能组合进行匹配, 以减少匹配误差作为优化策略^[7], 达到分解目的, 但效果并不理想. 后者则先处理电器的电气信息, 然后通过机器学习算法从数据中学习设备的特征模式, 从而实现对负荷的辨识与分解, 本文的落脚点在于机器学习方法. 根据学习过程是否使用标签信息, 可以进一步划分为监督学习^[8-13]和无监督学习^[14-17]两种. 通常, 监督学习需要已知负载的标签数据, 适用于负载类型较少的情况. 如文献[8, 9]采用支持向量机对负荷功率进行分解, 但其训练时间较长, 且对噪声敏感. 文献[10]通过深度稀疏编码实现负荷分解, 文献[11]通过基于簇的稀疏编码以提取多事件场景下所涉及的潜在单个事件, 但它们无法被用于实时分解. 文献[12]提出了利用域间知识转移的负荷分解方法. 文献[13]将随机森林作为负荷分解的学习算法; 无监督学习并不需要已知负载的标签数据, 只需从数据中挖掘特征的相似性. 隐马尔可夫模型 (HMM) 是一种常用的无监督学习方法, 在时间序列数据的分析方面具有很大优势, 被广泛应用于电器的能耗序列模拟, 其扩展形式, 因子隐马尔可夫模型(FHMM)常被应用于NILM领域. 文献[15]提出了一种基于分段整数二次约束规划的FHMM负荷分解模型. 文献[16]提出加性因子隐马尔可夫模型, 具有较高的分解精度. 文献[17]利用迭代模糊C均值算法以自适应的方式确定HMM的隐状态数. 文献[18]对处于启状态的设备建立FHMM进行负荷分解. 上述文献从不同的角度对基于FHMM的负荷模型进行改进, 但HMM的模型训练通常采用的是Baum-Welch算法^[19]. 这种方法虽然经典且容易实现, 但由于其初值选取的随机性, 模型参数估计时容易受到局部最优的限制, 影响分解结果的精确性.

为克服上述缺点, 一些学者将粒子群优化算法^[20]同HMM相结合, 提升学习的效率. 文献[21]将粒子群算法(PSO)引入HMM的参数训练过程, 并用于识别和预测报警信息, 文献[22]将自适应粒子群算法融入HMM中, 并用于识别地面装备的退化情况. 本文提出一种基于改进粒子群算法优化FHMM的非侵入式负荷分解方法, 通过线性递减策略动态控制惯性权重, 提升PSO算法的全局搜索能力, 并利用改进的PSO算法对因子隐马尔可夫模型的参数学习过程进行优化, 使之跳出局部极值点, 有效地提升了FHMM在负荷分解任务中的表现.

1 负荷分解模型的建立 1.1 单电力负载模型

在不同的运行条件下, 电力负载的电流幅值、有功功率和无功功率等电气参数是不同的, 但在相同的运行条件下总是表现出固有的特性. 通常, 电力负荷的运行过程可以看作是一系列的状态转移, 这与隐马尔可夫模型是一致的^[23].

根据负载稳态特征的不同, 负荷类型可分为 ON/OFF型负荷、FSM型负荷、CVD型负荷. 其中, ON/OFF型负荷与FSM型负荷并无本质区别, 都具有相对稳定的有限工作状态, 以上负荷类型可以涵盖大多数家用电器的情况^[1]. 单个电器的HMM模型示例如图1所示.

电力负载的工作状态通常采用有限状态集: $S = \{ {s_1}, {s_2}, \cdots, {s_K}\} $ 来描述, 其中, K指设备工作状态的个数. 离散时间序列 $Q = \{ {q_1}, {q_2}, \cdots, {q_t}, \cdots, {q_T}\} $ 表示负载的运行过程, 其中, ${q_t} \in S$ 表示设备在时间t的状态, T是运行时间的长度. 负载状态的输出有功功率可以用 $O$ $= \{ {o_1}, {o_2}, \cdots, {o_t}, \cdots, {o_T}\}$ 表示, 其中, ${q_t}$ 表示负载在t时刻的输出有功功率, 或者是一个包含多维负载特性的向量. 在该模型中, 状态序列S是马尔可夫链, 状态转移具有Markov性质, 不可观测. 因此, 单负载模型参数可以表示为 $\lambda = (\pi , A, B)$ .

$\pi $ 代表负载处于初始状态的概率, 即: ${\pi _i} = p({q_1}|{s_i})$ , $1 \leqslant i \leqslant K$ .

A代表状态转移概率矩阵, 即 $A = \{ {a_{ij}}, 1 \leqslant i, j \leqslant K\} $ ; ${a_{ij}} = p({q_t} = {s_j}|{q_{t - 1}} = {s_i})$ . 表示该负载从状态i到状态j的转移概率.

B指的是观测概率矩阵, 即 $B = \{ {b_k}({o_l}), 1 \leqslant k \leqslant M, $ $1 \leqslant l \leqslant T\} $ ; ${b_k}({o_l}) = p({o_t} = {n_k}|{q_t} = {s_i})$ 表示t时刻 ${s_i}$ 状态的输出观测值为 ${n_k}$ 的概率.

1.2 基于FHMM的总负荷模型

FHMM^[24]可用于建立具有多电力负载的负荷分解模型并进行分析^[25]. 在这里, 每个设备都可基于HMM建立负荷模型: 对于负荷 $i \in N$ , 采用两个序列来对模型进行描述, 其中 ${Q^i} = \{ q_1^i, q_2^i, \cdots, q_t^i, \cdots, q_T^i\} $ 表示设备运行状态, ${O^i} = \{ o_1^i, o_2^i, \cdots, o_t^i, \cdots, o_T^i\}$ 表示观测序列. 较之单个HMM, FHMM中的 ${O^i}$ 无法被观测. 电力入口处的总功率序列O则可以观测. 总负荷模型如图2所示. 总负荷模型的参数可定义为 $\Psi = (\Pi , A, B)$ , 其中, $\Pi $ 代表负载处于初始状态的概率, 由所有电器的隐含状态在初始时刻的概率决定, 即: ${\Pi _i} = p(q_1^1, q_1^2, \cdots, q_1^N)$ . A代表状态转移概率矩阵, 其中 $p({q_t}|{q_{t - 1}}) = p($ $ q_t^1, q_t^2, \cdots, q_t^N|q_{t - 1}^1, q_{t - 1}^2, \cdots, q_{t - 1}^N) $ . B指的是观测概率矩阵, 表示在某状态下对应总输出的概率, 其中 $p({O_t}|{q_t}) = p({O_t}|q_t^1, q_t^2, \cdots, q_t^N)$ , ${O_t}$ 指的是t时刻的观测序列.

2 模型求解 2.1 基于GMM的负荷状态聚类

电器每个工作状态的功率可认为服从均值为 $\;\mu $ 、方差为 ${\sigma ^{\text{2}}}$ 的高斯分布^[26]. 利用K-means进行状态聚类受功率波动的影响较大^[27]. 故本文采用高斯混合聚类模型来确定负荷状态, 并建立其状态集合 $Q = \{ {q_1}, {q_2}, $ $\cdots, {q_t}, \cdots, {q_T}\} $ , 该状态集合服从均值为 $\{ {\mu _1}, {\mu _2}, \cdots, {\mu _T}\} $ , 方差为 $\{ \sigma _1^2, \sigma _2^2, \cdots, \sigma _T^2\} $ 的混合高斯分布. GMM聚类假设数据是由潜在概率分布的混合产生的, 其中每个分量代表不同的类别, 因此, 可基于GMM聚类的方法确定各负荷设备的状态集合. GMM包括M个高斯分布模型的加权和:

$ \begin{gathered} p(x) = \sum\limits_{m = 1}^M {{\omega _m}p(x|m)} = \sum\limits_{m = 1}^M {{\omega _m}{N}({x_m}|{\mu _m}, {\Sigma _m})} \end{gathered} $

(1)

其中, M为模型个数, ${\omega _m}$ 为第m个高斯分布的权重, $p(x|m)$ 为第m个高斯分布的概率密度函数, 其均值为 $\;{\mu _m}$ , 观测向量的协方差矩阵为 ${\Sigma _m}$ .

通常采用EM算法估计GMM的参数. 先假定数据的类别, 然后根据迭代估计出GMM的权值 ${\omega _m}$ 、均值 $\;{\mu _m}$ 、协方差 ${\Sigma _m}$ , 其步骤如下所示.

(1) E步: 计算后验概率:

$ \gamma ({{\textit{z}}_{nm}}) = \frac{{{\pi _m}{ N}({x_n}|{\mu _n}, {\Sigma _n})}}{{\displaystyle\sum\limits_{i = 1}^M {{\omega _i}{ N}({x_n}|{\mu _i}, {\Sigma _i})} }} $

(2)

(2) M 步: 根据 E 步计算的后验概率 $\gamma ({{\textit{z}}_{nm}})$ 再计算 ${\omega _m}$ 、 $\;{\mu _m}$ 、 ${\Sigma _m}$ .

(3) 计算对数似然函数:

$ L(p(x)) = \sum\limits_{m = 1}^M {\log \left\{ {\sum\limits_{m = 1}^M {{\omega _m}{ N}({x_n}|{\mu _m}, {\Sigma _m})} } \right\}} $

(3)

(4) 重复迭代直到似然函数的值收敛.

2.2 负荷模型的参数估计

基于FHMM的负荷分解可以描述为: 在时段T内, 对于包含N个负载的负荷对象, 在已知模型参数 $\Psi $ 和总负荷观测序列 $O = \{ {o_1}, {o_2}, \cdots, {o_t}, \cdots, {o_T}\} $ 的条件下, 求解每个负载的最优状态序列 ${Q^i} = \{ q_1^i, q_2^i, \cdots, q_t^i, \cdots, q_T^i\} $ 和对应的负荷序列 ${O^i} = \{ o_1^i, o_2^i, \cdots, o_t^i, \cdots, o_T^i\} $ , 求解最优状态序列可视为解码问题, 采用Viterbi算法^[28]求解.

模型参数的估计方法可分为无监督和有监督两种. 无监督方法直接利用总负荷功率数据训练模型, 而非利用单个设备的功率数据. 本文中单个设备的工作状态为隐变量, 故采用无监督方法对模型参数进行估计. Baum-Welch算法作为一种典型方法常被用于解决这类问题, 根据初始参数 ${\Sigma ^0}$ 迭代估计出最优的参数 ${\Sigma ^*}$ .

根据HMM的前/后向算法可知:

$ {\gamma _t}(i) = \frac{{p({q_t} = {s_i}, O|\Sigma )}}{{p(O|\Sigma )}} = \frac{{{\alpha _t}(i){\beta _t}(i)}}{{\displaystyle\sum\limits_{j = 1}^T {{\alpha _t}(i){\beta _t}(j)} }} $

(4)

$ \begin{split} {\zeta _t}(i, j)& = \frac{{p({q_t} = {s_i}, {q_{t + 1}} = {s_j}, O|\Sigma )}}{{p(O|\Sigma )}} \\ & = \frac{{{\alpha _t}(i){a_{ij}}{b_j}{o_{t + 1}}{\beta _t}(j)}}{{\displaystyle\sum\limits_{j = 1}^T {{\alpha _t}(i){\beta _t}(j)} }} \end{split} $

(5)

其中, ${\gamma _t}(i)$ 代表t时刻隐状态 ${s_i}$ 的概率, ${\zeta _t}(i, j)$ 代表t时刻隐状态为 ${s_i}$ , t + 1时刻隐状态为 ${s_j}$ 的概率.

参数估计的过程如下.

(1) 设定初始值: ${\Sigma ^{(0)}} = ({\Pi ^{(0)}}, {A^{(0)}}, {B^{(0)}})$ .

(2) 按以下公式进行迭代:

$ \Pi _i^{n + 1} = {\gamma _1}(i) $

(6)

$ a_{ij}^{n + 1} = \frac{{\displaystyle\sum\limits_{t = 1}^{T - 1} {{\zeta _t}(i, j)} }}{{\displaystyle\sum\limits_{t = 1}^{T - 1} {{\zeta _t}(i)} }} $

(7)

$ b_i^{n + 1}(k) = \frac{{\displaystyle\sum\limits_{t = 1}^{T - 1} {{I_{{o_t} = {n_k}}}{\gamma _t}(i)} }}{{\displaystyle\sum\limits_{t = 1}^{T - 1} {{\gamma _t}(i)} }} $

(8)

设置迭代条件: $\left| {\log p(O|{\Sigma ^{n + 1}}) - \log p(O|{\Sigma ^n})} \right| < \varepsilon $ , 当 $p(O|{\Sigma ^n})$ 收敛时, 得到的 ${\Sigma ^n}$ 就是模型的参数.

由于Baum-Welch算法容易陷入局部最优, 本文考虑将粒子群算法与Baum-Welch算法相结合, 对上文中参数的训练过程进行优化.

PSO算法是Kennedy等人^[20]提出的一种进化优化算法. 它通过对鸟类的社会行为进行模拟来寻找优化问题的全局最优值, 并已应用于诸多科学领域^[29-31]. 在PSO中, 每个粒子就是一个解, 代表着该粒子目前的位置, 种群的解由一个群体组成. 为寻到最优解, 也即粒子的最佳位置, PSO以迭代的方式不断更新群体中所有具有速度向量的粒子的位置, 直到满足终止条件. 假设粒子i的位置由一个K维向量 ${x_i}(t)$ 表示, 其对应的速度向量为 $ {v_i}(t) $ . 每次迭代时记录粒子的最佳位置 ${p_{{\rm{best}}}}$ , 并进行信息共享, 筛选整个群体当前的最好位置 ${g_{{\rm{best}}}}$ , 然后根据以下规则来对速度和位置进行更新:

$ \begin{split} {v_i}(t) = &\omega {v_i}(t - 1) + {c_1}{r_1}({p_{{\rm{best}}}} - {x_i}(t - 1)) \\ & + {c_2}{r_2}({g_{{\rm{best}}}} - {x_i}(t - 1)) \end{split} $

(9)

$ {x_i}(t) = {x_i}(t - 1) + {v_i}(t) $

(10)

其中, $\omega $ 是惯性权重, 其作用是放缩上一时刻的速度, 决定了对当前速度继承程度, ${c_1}$ , ${c_2}$ 是学习因子, ${r_1}$ , ${r_2}$ 是[0, 1]上满足均匀分布的随机数.

用标准的粒子群算法求解优化问题时, 如果开始时 $\omega $ 的取值较大, 则粒子的全局搜索能力会被加强, 代价是局部搜索能力减弱, 不易得到精确解; 如果开始时 $\omega $ 的取值较小, 则可以增强粒子的局部搜索能力, 代价是全局搜索能力削弱, 这样做有利于算法的收敛, 但是也有收敛慢、容易陷入局部最优解的弊端. 选取合适的惯性权重值对提高算法性能而言至关重要. 针对这种情况, 引入一种随着迭代次数增加而线性递减的惯性权重, 在搜索过程中根据迭代次数的增加而对惯性权重不断进行调整, 即:

$ \omega = {\omega _{\max }}{{ - }}\frac{{T({\omega _{\max }} - {\omega _{\min }})}}{{{T_m}}} $

(11)

其中, ${\omega _{\max }}$ 是最大迭代次数时对应的惯性权重, ${\omega _{\min }}$ 指的是初始时刻的惯性权重, T是当前的迭代次数, ${T_m}$ 是最大迭代次数. 依据如此规则更新惯性权重, 使得粒子在初始阶段具有较强的搜索能力, 在后期又能够得到相对精确的结果, 进而使得算法性能提升.

方便起见, 将粒子i的位置坐标 ${x_i}$ 表示为FHMM 的模型参数. 适应度函数由产生观测序列O的概率函数 $P(O|\Sigma )$ 来确定, 即 $f(\Sigma ) = P(O|\Sigma )$ , 由于 $P(O|\Sigma )$ 通常较小, 而且只代表模型输出概率值的相对大小, 故本文将上述函数两边取对数后作为适应度函数. 本文所提出的IPSO-FHMM算法的流程如算法1.

算法1. IPSO-FHMM

1) 初始化: 依据训练规模, 随机生成N组模型参数 $\scriptstyle {\Sigma ^0}$ ;

2) 依据Baum-Welch算法, 基于初始参数 $\scriptstyle {\mathrm{\Sigma }}^{0} $ 训练N组参数 $\scriptstyle {\Sigma ^1}$ ;

3) 将这N组参数作为粒子, 根据式 (9)与式(10), 通过比较每次更新后的适应度值的大小来调整粒子的位置和速度.

4) 根据式 (11), 调整惯性权重 $\scriptstyle \omega $ , 以调整粒子的搜索能力.

5) 确定最大训练次数 $\scriptstyle {T_m}$ , 设定最小允许误差值. 当训练次数达到上限时, 状态估计值对应的参数就是全局最优参数 $\scriptstyle {\Sigma ^*}$ .

2.3 负荷状态的估计

在确定了FHMM的最优参数 ${\Sigma ^*}$ 之后, 可以通过Viterbi算法求解概率最大路径. 观测序列 $O = \{ {o_1}, {o_2}, \cdots,$ ${o_t}, \cdots, {o_T}\}$ 已知, 则目标函数可等价为状态序列和观测序列的联合概率最大化:

$ {\rho _t}(q) = \max [p({q_1}, {q_2}, \cdots, {q_T}, {o_1}, {o_2}, \cdots, {o_T}|\Sigma )] $

(12)

递推可知:

$ {\rho _{t + 1}}(q) = \max [{\rho _t}(p){a_{pq}}{b_q}({o_{t + 1}})] $

(13)

其中, $ {b_q}({o_{t + 1}}) $ 是状态为 ${s_q}$ 条件下观测为 ${o_{t + 1}}$ 的概率. 根据迭代结束时得到的最优状态组合 $q_T^*$ 进行路径回溯, 最终可得最佳状态序列 ${Q^*} = \{ q_1^*, q_2^*, \cdots, q_T^*\} $ .

3 实验结果与分析

本文选用公开数据集AMPds2^[32]进行非侵入式负荷分解, 这是一个低频数据集, 采样周期为1 min, 监测持续时间为2年. 从2个负荷中选择5种负荷作为研究对象, 分别为洗衣机、洗碗机、热力泵、壁挂炉、烘干机, 包含了FSM型负荷、CVD型负荷, 如此选取的研究对象使得研究更加贴合现实生活中的负荷使用状况. AMPds2中包含了11种电力测量参数, 本文选择有功功率作为负荷特征进行研究, 选用F1-score和归一化分解误差(NDE)指标评估负载分解算法的性能. F1-score是准确率和召回率的调和平均值, 其计算方法如下:

$ P = \frac{{TP}}{{TP + FP}} $

(14)

$ R = \frac{{TP}}{{TP + FN}} $

(15)

$ F1 {\textit{-}} score = \frac{{2PR}}{{P + R}} $

(16)

其中, TP表示预测为正类实际也为正类的样本数, FP表示预测为正类实际为负类的样本数, FN表示预测为负类实际为正类的样本数, 对于ON/OFF设备, 开状态是正例, 关状态是反例; 对于多状态设备, 指定的工作状态为正例, 其余为负例. NDE则描述了实际功率和分解功率之间的重合程度^[9], 其定义如下:

$ NDE = \sqrt {\frac{{{{(y_t^i - \hat y_t^i)}^2}}}{{{{(y_t^i)}^2}}}} $

(17)

其中, $ y_t^i $ 代表电器的实际功率, $ \hat y_t^i $ 代表被分解出来的功率值. NDE值越接近0, 说明分解精度越高.

采用GMM对负荷数据进行聚类, 所选5台设备的聚类结果如表1所示. 洗衣机属于CVM负载, 基于分解精度和计算速度的考量, 将该设备的聚类数设置为 4.

选取8周的采样数据进行负荷分解. 其中前3/4的数据设置为训练样本集, 用以获取设备的工作状态, 并估计出模型参数. 剩余1/4的数据作为测试集, 对算法的性能进行测试. 负荷分解的结果如图3 所示.

F1-score越趋于1, NDE越趋于0则代表算法的性能越理想, 为了验证本文所述方法的准确性, 进行了计算机仿真, 并将结果与Hart方法、Hart w/ MAP方法以及AFAMAP算法^[16]、IFCM 算法^[17]的结果进行对比, 仿真研究过程分为4组, 即有2–5个设备同时工作的情况.

从图3可以看出, 本文提出的模型在两种指标上明显优于Hart及其改进后的模型. 也优于AFAMAP模型的表现. 在与IFCM模型进行比较时, 可以发现当2–4台电器同时工作的状态下, 本文模型表现最佳. 当分解电器数量较多时, 比如5台电器同时工作时, 本文的模型精度略低于IFCM模型. 由此可同样得出, 负荷数增多, F1-score逐渐减小, NDE逐渐增大, 即分解精度逐渐下降, 显然负荷数对分解结果有较大影响.

所提出方法对于各个电器的分解效果如图4所示. 随着同时工作的设备数量的增加, 对于每个电器而言, F1-score降低, NDE值上升. 可见对于单个电器, 分解精度也随电器数目的增加而降低. 烘干机、热力泵、壁挂炉这3种电器都属于典型的FSM型负荷, 工作状态为有限个, 且能按照一定顺序切换状态, 各个工作状态的功率水平相差较大, 分解结果在各个电器数目的情况下都较为理想. 洗碗机的分解结果随电器数目增加而呈现显著下降的趋势. 洗衣机的分解结果始终低于其他设备, 可见连续状态型负荷通过FHMM建模的准确性并不够理想, 进而影响了分解的效果.

4 结论与展望

本文针对目前关于负荷分解研究中存在的问题, 提出了基于改进粒子群优化因子隐马尔可夫模型的非侵入式负荷分解方法. 通过对负荷基于FHMM建立总负荷模型, 将负荷分解问题转变为求解最优状态组合的问题. 根据电器工作状态的功率服从高斯分布的特性, 采用GMM进行电器状态聚类, 避免了K-means聚类中心不稳定问题. 同时, 本文对FHMM的训练算法, 即Baum-Welch算法进行了改进, 引入了粒子群优化算法对模型参数估计过程进行优化, 避免了陷入局部最优问题. 同时考虑到粒子群算法的惯性权重对搜索能力的影响, 采用线性递减的方式来选取惯性权重. 最后通过Viterbi算法求得了最优的状态序列, 实现了对总负荷功率的分解. 最后将本文提出的方法分别与Hart方法、Hart w/ MAP方法、AFAMAP方法IFCM方法在不同设备同时工作的情况下进行性能比较. 结果表明本文的方法在性能上优于上述其他方法. 值得注意的是, 本文所提出的负荷分解方法, 计算成本相对较高, 同时随着负荷数增多, 分解精度会呈现一定程度的下降. 在下一步的研究中需要针对这些问题, 考虑结合深度学习方法^[33-35], 对负荷分解问题进行更深入的讨论.