## An L 0 -Norm Constrained Non-Negative Matrix Factorization Algorithm for the Simultaneous Disaggregation of Fixed and Shiftable Loads

Ahmad Khaled Zarabie

Sanjoy Das

Electrical &amp; Computer Engineering Dept. Kansas State University Manhattan, KS 66506, USA

Abstract -Energy disaggregation refers to the decomposition of energy use time series data into its constituent loads. This paper decomposes daily use data of a household unit into fixed loads and one or more classes of shiftable loads. The latter are characterized by ON/OFF duty cycles. A novel algorithm based on non-negative matrix factorization (NMF) for energy disaggregation is proposed, where fixed  loads  are  represented  in  terms  of  real-valued  basis vectors,  whereas  shiftable  loads  are  divided  into  binary  signals. This  binary  decomposition  approach  directly  applies L 0-norm constraints  on  individual  shiftable  loads.  The  new  approach obviates the need for more computationally intensive methods (e.g. spectral  decomposition  or  mean-field  annealing)  that  have  been used  in  earlier  research  for  these  constraints.  A  probabilistic framework for the proposed approach has been addressed. The proposed  approach's  effectiveness  has  been  demonstrated  with real consumer energy data.

Keywords-binary  matrix;  energy  disaggregation;  fixed  load; non-negative matrix factorization; Frobenius norm, shiftable load; sparsity constraint; unsupervised learning

## I. INTRODUCTION

Energy disaggregation refers to the decomposition of energy usage into multiple components in a physically meaningful way [1],[2].  For instance, daily energy consumption of a household unit can be disaggregated into various loads, such as refrigerator, air  conditioner,  lighting,  pool  pump,  and  other  loads  that  are typically  present  in  homes.  In  disaggregation  tasks,  the  data samples are of the form of energy use over a fixed duration of time, at regularly spaced intervals. Until the last decade, energy disaggregation had met with little success. However, in recent years, NMF [3] has emerged as powerful tool for this purpose and has met with remarkable success.

The  classical  NMF  algorithm  decomposes  an  input  data matrix 𝐗 whose  columns  are 𝐷 × 1 sample vectors,  into  two factors, 𝐖 and 𝐇 , so that their product equals 𝐗 . Usually, 𝐗 has a very large number of columns, which are independent samples. Although there exists an abundance of classical matrix methods to factorize a given matrix in such a manner, in NMF there is the additional  constraint that 𝐖≥𝟎 and 𝐇 ≥ 𝟎 [4],[5].  This  nonnegativity  requirement  placed  on  both 𝐖 and 𝐇 render  NMF suitable for many applications. For instance, NMF can be used to  decompose  image  sequences  or  audio  power  spectra  into factors 𝐖 and 𝐇 that can only have non-negative values, since

Electrical &amp; Computer Engineering Dept. Kansas State University Manhattan, KS 66506, USA

neither pixel values nor power components can be negative. This is also the situation in the present research. The columns of 𝐖 , which are relatively lesser in number, serve as basis vectors so that  each  sample 𝐱(𝑛), 𝑛 ∈ 𝒩 which are columns of 𝐗 can be represented  as  a  weighted combination of  the bases,  with  the non-negative weights being the corresponding column vectors 𝐡(𝑛) of 𝐇 . In  load  disaggregation,  the  basis  vectors  may correspond to individual appliances [1],[2].

Due  to  the  non-negativity  constraints, 𝐗 is not  exactly factorizable into 𝐖 and 𝐇 ; whence the goal of NMF is to seek an approximate solution, so that 𝐗 ≈ 𝐖𝐇 . NMF is an ill-posed problem  since  if 𝐗 ≈ 𝐖𝐇 is  one  solution,  then  so  are  other factorizations of the form 𝐗 ≈ 𝐖𝐏𝐏 T 𝐇 where 𝐏 is any rotation matrix that preserves non-negativity, that can be considered to be other solutions.

Furthermore, NMF involves seeking an optimum of a nonconvex objective function. For example, although the Frobenius norm, ‖⋅‖ 𝐹 , is convex in a single argument, this is not the case with  NMF,  where  the  objective  function  is ‖𝐗 - 𝐖𝐇‖𝐹 2 ,  the square  normed difference between 𝐗 and 𝐖𝐇 (a  very  popular choice in many applications). This is due to the presence of the matrix product 𝐖𝐇 in the argument. Two other very commonly used objectives are the (generalized) Kullback-Leibler divergence and the Itakura-Saito divergence [6]. The latter is a special case of the 𝛽 -divergence that is also used in some NMF applications  [7]. From  the  non-convexity  of  the  objective function, it follows that the NMF problem would encounter local minima. Although most NMF algorithms are first-order descent methods, providing local minima as the outputs, it has been seen that NMF consistently produces good results for a wide variety of applications.

Several  approaches  for  NMF  have  been  proposed,  the earliest of them being the multiplicative method [3],[4]. Further details of this method are postponed until the following section. Another  successful  method  is  the  alternating  least  squares algorithm  [8].  In  this  approach,  the  KKT  conditions  at  the optimum  provide  a  set  of  equations  through  which  iterative updates on the matrices 𝐖 and 𝐇 are carried out. Each iteration typically involves matrix pseudoinversion that is accomplished indirectly  through  LU  decomposition.  The  projected  gradient descent algorithm has also been applied for NMF [9]. This is a straightforward extension of gradient descent, with the matrices 𝐖 and 𝐇 thresholded  at  the  end  of  each  step  to  non-negative

values.  More  recent  NMF  algorithms  include  the  hierarchical NMF  [10]  that  is  designed  to  deal  with  unduly  large  sized matrices, and the randomized NMF algorithm [11].

Earlier  NMF  algorithms  that  restrict  either 𝐖 or 𝐇 to  be sparse matrices tended to replace the non-differentiable 𝐿0 norm constraint with an approximate measure of sparsity based on the 𝐿1 and 𝐿2 norms  [12].  More  recently,  algorithms  that  make explicit use of the 𝐿0 norm have also been suggested [13]. The algorithm in [13] applies k -SVD [14] based rank-1 updates on either  of  the  two  matrices.  Therefore,  its  use  is  confined  to Frobenius norm NMF.

The approach in [15] applies 𝐿0 norm constraints (referred to  as sum-to-k constraints  in  [15])  indirectly  by  including  a penalty term to a Frobenius normed objective. It has been used for  energy  disaggregation  of  HVAC  load  components  in  an industrial building and in a smart home setting. Another method has been proposed in [2],[16] to impose 𝐿0 constraints, which uses  a  softmax  distribution  for  the  elements  in 𝐇 to  assign weights to them in such a manner that those with higher values are likelier to improve the objective function. Semi-supervised NMF  using  prior  knowledge  of  the  usage  time  profiles  of individual appliances has been proposed in [17]. NMF has been applied  for  data  over  larger  periods  in  time  to  glean  seasonal trends in usage profiles in [18].

In this research, a novel NMF  algorithm for load disaggregation has been proposed. It uses the Frobenius norm as objective,  although  the  approach  is  generalized  enough  to  be extended to others. The novelty of this approach is the manner in which the load is divided into two classes, ( ) i fixed loads , and, ( ii ) shiftable loads . This is a fundamental distinction between the different  appliances  in  a  typical  household  that  has  not  been hitherto considered.

Fixed  loads  are  associated  with  appliances  that  are  in continual use throughout the day. Examples of shiftable loads include lighting, air-conditioners, and refrigerators. On the other hand, shiftable loads pertain to appliances that are used intermittently,  such  as  washers,  dryers,  and  ovens.  The  latter class  of  loads  are  characterized  by  duty  cycles,  with  typical temporal profiles. Our approach tries to exploit this feature in shiftable loads to obtain improved load disaggregation. Moreover,  direct 𝐿0 norm  constraints  are  imposed  on  the number of duty cycles, separately for each shiftable load.

## II. PROBLEM STATEMENT

## A. Load Models

Each input 𝐱(𝑛) is a vector of 𝐷 regularly spaced samples of  energy  usage  of  a  single  household  unit  over  a  24-hour period. Therefore, the index 𝑛 can be regarded as the energy use of the 𝑛 th day in the sample set 𝒩 . The purpose of the proposed NMF  algorithm  is  to  express  each  sample  in  the  following manner.

$$\mathbf x ( n ) \approx \sum _ { k \in \mathcal { F$$

In  (1), ℱ is  the  fixed  load  basis  set  and 𝒮𝑗 is  that  of  the 𝑗 th shiftable load. Each 𝐷 × 1 vector 𝐰𝑘 𝑠 𝑓 is a fixed load basis vector 𝑘 ( 𝑘 ∈ ℱ ).  Likewise,  each 𝐷 × 1 𝐰 𝑗,𝑘 is  a  basis  vector 𝑘 ( 𝑘 ∈ 𝒮𝑗 ) of the 𝑗 th shiftable load. In terms of basis matrices,

$$X \approx \widetilde { X } = W ^ { f } \mathbf H ^$$

The quantities 𝐖 𝑓 and 𝐖𝑗 𝑠 are basis matrices of dimensionalities 𝐷 × |ℱ| and 𝐷 × |𝒮 𝑗 | . The  matrices 𝐇 𝑓 = [𝐡 𝑓 (𝑛)]𝑛∈𝒩 and 𝐇𝑗 𝑠 = [𝐡 𝑗 𝑠 (𝑛)]𝑛∈𝒩 are |ℱ| × 𝑁 and |𝒮 | × 𝑁 𝑗 dimensional arrays of weights ℎ 𝑘 𝑓 (𝑛) and ℎ 𝑗,𝑘 𝑠 (𝑛) ∈ {0,1} . This decomposition is illustrated in Fig. 1.

The  shiftable  loads  are  characterized  by  OFF-ON  duty Fig. 1 . Schematic of NMF factorization with separate basis and weight matrices for the fixed and shiftable loads.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00142v1-figure-ac2249a2-97dd-45c1-99c0-3e02dae70e49.jpg)

cycles. It is assumed that during ON intervals, the 𝑗 th shiftable load draws a constant amount of energy 𝑝𝑗 and that it can stay ON for a maximum of 𝐿𝑗 time intervals. Hence,

$$\left \| \mathbf h _ { j } ^ { s } ( n ) \right \$$



Fig. 2 shows the time profile of a shiftable load for the duration of a single day. It should be noted that as 𝐿𝑗 is the upper limit on the number of duty cycles, the actual number of such cycles may be less than 𝐿𝑗 .

Fig. 2 . Illustrative usage profile of a shiftable load showing a maximum  of 𝐿𝑗 duty  cycles,  with  each  such  cycle  being  a rectangular pulse of uniform peak 𝑝𝑗 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00142v1-figure-1c87a775-5246-4abc-8ea4-6c360d19c138.jpg)

## B. Objective

The squared Frobenius norm of the difference between the real data matrix 𝐗 and its approximation 𝐗 ̃ ,

$$\tau \ u u u u u u u u u u u u u u u u u u$$



The expression for the approximation 𝐱 ̃ (𝑛) in  (4)  has  already been provided in (1). The summation in (4) is carried out over all  samples  in  the  set 𝒩 .  A  probabilistic  justification  for  this choice  is  provided  in  Section  III,  which  also  serves  as  a theoretical basis for obtaining the weights associated with the shiftable loads in the proposed approach.

## A. Probabilistic Framework

Consider  the  following  probabilistic  interpretation.  The joint probability of 𝐗 , assuming that all samples are statistically independent is given by the following expression,

$$\begin{array} {$$

It  is  assumed  that  each  sample 𝐱(𝑛) follows  a  Gaussian distribution  around  its  expected  value 𝐱 ̃ (𝑛) .  In  this  case,  the probability  of  each  such  sample  can  be  expressed  in  the following manner,

$$\dots \cdots \circ \dots \circ \dots \circ, & \quad \vdots _ { r } \cdots \\ p [ \mathbf x ( n ) | & \mathbf W ^ { f }, \mathbf h ^ { f } ( n ), \mathbf W ^ { s } _ { j }, \mathbf h ^ { s } _ { j } ( n ) ] \\ & = \frac { 1 } { \sigma ^ { D } ( 2 \pi ) ^ { \frac { D } { 2 } } \prod _ { d } } e ^ { - \frac { 1 } { 2 \sigma ^ { 2 } } \left ( x _ { d } ( n ) - \mathcal { X } _ { d } ( n ) \right ) ^ { 2 } }. \quad \quad ( 6 ) \quad \sin _ { \substack { \text{subjec} \\ \text{deriva} \\ \text{show} } }$$

The negated log probability of 𝐱(𝑛) is,

$$\text{$1$ he negated log probability of $x(n)$ is}, & & \\ - \log p \left [ x ( n ) | W ^ { f }, \mathbf h ^ { f } ( n ), W ^ { s } _ { j }, \mathbf h ^ { s } _ { j } ( n ) \right ] & & \nabla _ { \mathbf v } \\ & = \frac { 1 } { 2 \sigma ^ { 2 } } \sum _ { d } \left ( x _ { d } ( n ) - \tilde { x } _ { d } ( n ) \right ) ^ { 2 } \ + \log \sigma ^ { D } ( 2 \pi ) ^ { \frac { D } { 2 } } & & \\ & = \frac { 1 } { 2 \sigma ^ { 2 } } \| \mathbf x ( n ) - \tilde { x } ( n ) \| _ { 2 } ^ { 2 } \ + \log \sigma ^ { D } ( 2 \pi ) ^ { \frac { D } { 2 } }. & & \nabla _ { \mathbf v }$$

Hence the negative log probability of 𝐗 is,

$$- \log p \left [ X | W ^ { f }, H ^ { f }, W ^ { s } _ { j }, H ^ { s } _ { j } \right ]$$

$$& = \prod _ { n } p [ \mathbf x ( n ) | \mathbf W ^ { f }, \mathbf h ^ { f } ( n ), \mathbf W ^ { s } _ { j }, \mathbf h ^ { s } _ { j } ( n ) ] & \quad \text{ for this} \\ & = \frac { 1 } { 2 \sigma ^ { 2 } } \left \| \mathbf X - \widetilde { \mathbf X } \right \| _ { F } ^ { 2 } \, + N \log \sigma ^ { D } ( 2 \pi ) ^ { \frac { D } { 2 } }, \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \ \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \end{cases}$$

where the Frobenius norm of the difference between 𝐗 and 𝐗 ̃ is given by the earlier expression in (4).

Applying the maximum likelihood criterion to (8),

$$\left [ W ^ { f }, \mathbf H ^ { f }, \mathbf H ^ { s } _ { j } \right ] _ { M L } = \underset { W ^ { f }, \mathbf H ^ { f }, \mathbf H ^ { s } _ { j } } { \arginf } \frac { 1 } { 2 } \left \| \mathbf X - \widetilde { \mathbf X } \right \| _ { F } ^ { 2 }. \quad \text{ (9)} \quad \text{ norr}$$

The expression in (9) above provides a theoretical justification for  the  choice  of  objective  function  in  (4)  in  the  present approach.

The component of the objective function associated with the 𝑛 th sample is,

$$n \text{ samples } n, & \sup _ { \substack { sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma s alpha sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma salpha sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma selta sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigsigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigmasigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma$$

## B. Multiplicative Updates

Multiplicative updates are used for the matrices 𝐖 𝑓 and 𝐇 𝑓 that are associated with the fixed load. Although this method is quite routine in the existing literature on NMF, it is described here for the sake of completeness of this paper. More details can be found in [2],[3],[4],[5].

Consider any parameter 𝐏 (which can be either 𝐖 𝑓 or 𝐇 𝑓 ). Let ∇𝐏Φ be the gradient of tan arbitrary objective function Φ. Its  gradient  can  be  expressed  in  terms  of  its  positive  and negative  components,  so  that ∇𝐏Φ = ∇𝐏 + -∇𝐏 -.    In  gradient descent, the update rule would have been of the form,

$$P \leftarrow P - \eta \nabla _ { P } ^ { + } + \eta \nabla _ { P } ^ { - }.$$

In the multiplicative  method,  the  following  multiplicative update replaces the gradient ascent step,

$$P \leftarrow P \circ \nabla _ { P } ^ { - } \emptyset \, \nabla _ { P } ^ { + }.$$

In  the  proposed  approach,  the  parameters 𝐖 𝑓 and 𝐇 𝑓 are subject  to multiplicative  updates. In order  to  do  so,  the derivatives of the objective must be first computed. It can be shown that,

$$& \text{shown that,} \\ & \nabla _ { W ^ { f } } \frac { 1 } { 2 } \left \| \mathbf X - \widetilde { \mathbf X } \right \| _ { F } ^ { 2 } = \\ & - \mathbf X \circ \left ( W ^ { f } \mathbf H ^ { f } + \sum _ { j } W _ { j } ^ { s } \mathbf H _ { j } ^ { s } \right ) ^ { \circ - 1 } \mathbf H ^ { f } ^ { T } + \mathbf 1 _ { D \times N } \mathbf H ^ { f } ^ { T }. \quad ( 1 1 ) \\ \right ) \quad \nabla _ { \mathbf H ^ { f } } \frac { 1 } { 2 } \left \| \mathbf X - \mathbf X \right \| _ { F } ^ { 2 } = \\ & - / \quad \longrightarrow \quad \bigwedge ^ { \circ - 1 } \quad -$$

$$- W ^ { f ^ { T } } X \circ \left ( W ^ { f } { \mathbf H } ^ { f } + \sum _ { j } W _ { j } ^ { s } { \mathbf H } _ { j } ^ { s } \right ) \ + W ^ { f ^ { T } } { \mathbf 1 } _ { D \times N }. \quad ( 1 2 ).$$

$$& - X \circ \left ( W ^ { f } \mathbf H ^ { f } + \sum _ { j } W ^ { s } _ { j } \mathbf H ^ { s } _ { j } \right ) ^ { \dashv - \cdot \cdot } \mathbf H ^ { f } + 1 _ { D \times N } \mathbf H ^ { f }. \quad ( 1 1 ) \\ & \nabla _ { \mathbf H ^ { f } } \frac { 1 } { 2 } \| \mathbf X - \tilde { X } \| _ { F } ^ { 2 } = \\ & - W ^ { f } \mathbf X \circ \left ( W ^ { f } \mathbf H ^ { f } + \sum _ { j } W ^ { s } _ { j } \mathbf H ^ { s } _ { j } \right ) ^ { \dashv - 1 } + W ^ { f } \mathbf 1 _ { D \times N }. \quad ( 1 2 )$$

Using the above expressions for the gradients, the update rules for the matrices pertaining to the fixed loads are,

$$W ^ { f } \leftarrow W ^ { f } \circ \frac { ( X \circ \widetilde { X } ^ { \circ - 1 } ) \mathbf H ^ { f } } { 1 _ { D \times N } \mathbf H ^ { f } } \,. \quad \quad ( 1 3 )$$

$$\mathbf H ^ { f } \leftarrow \mathbf H ^ { f } \circ \frac { \mathbf W ^ { f } ^ { \mathrm T } ( X \circ \widetilde { X } ^ { \circ - 1 } ) } { \mathbf W ^ { f } \mathbf 1 _ { D \times N } } \,. \quad \quad ( 1 4 )$$

Since  the  columns  of 𝐖 𝑓 must  be  unit  vectors,  they  are normalized  at  the  end  of  each  multiplicative  update  step  as shown below,

$$w _ { k } ^ { f } \leftarrow \frac { w _ { k } ^ { f } } { \| w _ { k } ^ { f } \| }.$$

## C. Sparsity Constrained Binary Updates

In order to train 𝐇𝑗 𝑠 , arginf 𝐇 𝑗 𝑠 Φ is obtained using the binary updating  heuristic  proposed  in  this  research.  This  heuristic directly  imposes  the 𝐿0 norm  constraint.  The  shiftable  basis matrices 𝐖𝑗 𝑠 remain fixed throughout the training process.

It should be noted that since ℎ 𝑗,𝑘 𝑠 (𝑛) ∈ {0,1} and ‖𝐡 𝑗 𝑠 (𝑛)‖ 0 ≤ 𝐿𝑗 , computing the optimal weight matrix 𝐇𝑗 𝑠 is an

NP-hard problem. Therefore a hill-climbing heuristic is proposed. The component of the total error due to any sample 𝑛 is given by the following expression,

$$n \text{ is given by the following expression, } \quad \text{$each$} \quad \text{$cores$} \quad \text{$implen} \\ \quad \text{$= W^{\prime}h^{\prime}(n)+\sum_{j} W_{j}^{s}h_{j}^{s}(n)-x(n)} \quad \text{$Howev subrou} \\ \quad \text{$= W^{\prime}h^{\prime}(n)+W_{j}^{s}h_{j}^{s}(n)} + \sum_{j^{\prime}\neq j} W_{j^{\prime}h_{j^{\prime}}^{s}(n)-x(n)} \quad \text{$r_{j}(i}$} \\ \quad \text{$= W_{j}^{s}h_{j}^{s}(n)-r_{j}(n).} \quad \text{$(16)} \quad \text{$r_{j}(i}$} \\ \text{In the above equality, $r_{j}(n)$ is the residual approximation of the} \quad \text{$h_{j}$} ($$

In the above equality, 𝐫 𝑗 (𝑛) is the residual approximation of the 𝑛 th sample when ignoring the 𝑗 th shiftable load as shown below,

$$n \text{ samples when running use } J \text{ samples low as shown over}, \\ \mathbf r _ { j } ( n ) & = \mathbf x ( n ) - W ^ { f } \mathbf h ^ { f } ( n ) - \sum _ { j ^ { \prime } \neq j } W _ { j ^ { \prime } } ^ { s } \mathbf h _ { j ^ { \prime } } ^ { s } ( n ). & \text{Th vector} \\ & \approx \sum _ { k \in \mathcal { S } _ { j } } w _ { j, k } ^ { s } h _ { j, k } ^ { s } ( n ). & \text{(17)} \text{ \quad \text{the set} } \\ \text{return}$$

The  objective  function  for  any  given  sample 𝑛 provided earlier in (10) can be expressed as given below,

$$\text{earlier in} \, ( \check { 1 0 } ) \, can be expressed as \, \underbrace { \cdot } \, \times \, \cdot \, \cdot \, \cdot \\ \varphi ( n ) & = \frac { 1 } { 2 } \, \| \text{e} ( n ) \| _ { 2 } ^ { 2 } \\ & = \frac { 1 } { 2 } \sum _ { d } e _ { d } ^ { 2 } ( n ) \\ & = \frac { 1 } { 2 } \sum _ { d } \left ( \sum _ { k \in \real s _ { j } } w _ { j, d, k } ^ { s } h _ { j, k } ^ { s } ( n ) - r _ { j, d } ( n ) \right ) ^ { 2 } \,. \, ( 1 8 ) \, \quad \text{(} \, \text{argu} \, \text{consist} \, \text{$w_{j}$}, \, \text{th} \\ \text{In} \, ( 1 8 ), \, \text{consider any term} \, d \, of \, the \, \text{outer summation} \,. \, The \, \text{inner} \, \quad \text{shown}$$

In (18), consider any term 𝑑 of the outer summation. The inner summation  is carried out over all columns 𝑘 ∈ 𝒮 𝑗 of 𝐖𝑗 𝑠 such that 𝑤𝑗,𝑘,𝑑 𝑠 ≠ 0 . Noting that any such non-zero 𝑤𝑗,𝑘,𝑑 𝑠 = 𝑝𝑗 , the error 𝜑(𝑛) can be re-expressed in the following manner,

$$\text{en} \varphi ( n ) = \frac { 1 } { 2 } p _ { j } ^ { 2 } & \sum _ { d \in \mathcal { D } _ { k } } \left ( \sum _ { k \in \mathcal { S } _ { j } } h _ { j, k } ^ { S } ( n ) - \frac { r _ { j, d } ( n ) } { p _ { j } } \right ) ^ { 2 } & \text{$T$} \\ + \frac { 1 } { 2 } \sum _ { d \notin \mathcal { D } _ { k } } r _ { j, d } ^ { 2 } ( n ) \,. & \text{$(19)} \quad \text{there} \quad \text{there} \quad \text{$a Boo} \\ \text{The summation in (19) above is carried out over all elements in} \quad \text{(step}$$

The summation in (19) above is carried out over all elements in the set of indices 𝒟𝑘 defined as, 𝒟𝑘 = {𝑑|𝑤 𝑗,𝑑,𝑘 𝑠 ≠ 0}. The set 𝒟𝑘 can be obtained easily by examining 𝑘 th column of 𝐖𝑗 𝑠 and including  all  indices 𝑑 , that  have  nonzero  entries  in  that column.  After  scaling  appropriately and ignoring  the  term  in (19) not containing ℎ 𝑗,𝑘 𝑠 (𝑛) , the objective is,

$$\varphi ^ { \prime } ( n ) = \sum _ { d \in \mathcal { D } _ { k } } \left ( \sum _ { k \in \mathcal { S } _ { j } } h _ { j, k } ^ { s } ( n ) - \, r _ { j, d } ^ { \prime } ( n ) \right ) ^ { 2 }. \quad ( 2 0 ) \quad \stackrel { \widetilde { \text{iterat} } } { \text{iterat} } } \, \underset { \text{error} } { \text{error} }$$

The proposed algorithm begins with 𝐡𝑗 𝑠 (𝑛) = 𝟎 |𝒮 𝑗 | , updating it in a stepwise manner, one element in each step. In each step 𝑙 a new index 𝑘 ∈ 𝒮 𝑗 is selected and the corresponding ℎ 𝑗,𝑘 𝑠 (𝑛) updated  to  1.  The  algorithm  can  be implemented  using  a  separate  binary  heuristic  subroutine. However,  the  correct  arguments  need  to  be  passed  to  the subroutine. The steps to do so are shown below.

$$& \mathbf r _ { j } ( n ) \leftarrow \mathbf x ( n ) - W ^ { f } \mathbf h ^ { f } ( n ) - \sum _ { j ^ { \prime } \neq j } W ^ { s } _ { j ^ { \prime } } \mathbf h ^ { s } _ { j ^ { \prime } } ( n ) \\ & \mathbf r _ { j } ( n ) \leftarrow \frac { 1 } { p _ { j } } \mathbf r _ { j } ( n ) \\ & \mathbf h ^ { s } _ { j } ( n ) \leftarrow \, \text{hililClimb} \left ( W ^ { s } _ { j }, \mathbf r _ { j } ( n ) \right )$$

The arguments of subroutine hillClimb() are a 𝐷 × 1 vector 𝐫 and  either  a 𝐷 × 𝑆 binary  matrix 𝐖 ,  or  equivalently the sets of indices 𝒟𝑘 for each 𝑘 ∈ {1,2, … , 𝑆} . The subroutine returns an updated 𝐡 that minimizes,

$$\varphi ( \mathbf h | \mathbf w, \mathbf r ) = \sum _ { d } \left ( r _ { d } - \sum _ { k = 1 } ^ { s } w _ { k, d } h _ { k } \, \right ) ^ { 2 }. \quad \ ( 2 1 )$$

Alternately, the expression in (21) can be written as,

$$\varphi ( \mathbf h | \mathcal { D } _ { 1 }, \dots, \mathcal { D } _ { D }, \mathbf r ) = \sum _ { d \in \mathcal { D } _ { k } } \left ( r _ { d } - \sum _ { k = 1 } ^ { S } h _ { k } \right ) ^ { 2 } + \frac { 1 } { 2 } \sum _ { d \notin \mathcal { D } _ { k } } r _ { d } ^ { 2 } \,, \, ( 2 2 )$$

The  proposed  subroutine  is  outlined  in Fig.  4 .  The  two arguments  that  are  passed  to  the  subroutine hillClimb() consist  of  a  basis  matrix 𝐖 ,  and  a  residual  vector 𝐫 .  If 𝐖≡ 𝐖𝑗 𝑠 , the basis matrix of the 𝑗 th shiftable load, then 𝐫 ≡ 𝐫 (𝑛) 𝑗 as shown  in  (17)  but  after  normalization  so  that 𝑝𝑗 = 1 .  It  is assumed that two constants, 𝑆 ,  the  number  of  columns in 𝐖 , and 𝐿 ,  the  maximum  allowable  value  of ‖𝐡‖0 are  implicitly accessible to the subroutine.

The  subroutine hillClimb() maintains  a  set ℒ of  all indices 𝑘 such  that ℎ𝑘 = 0 .  Since  the  subroutine  begins  with 𝐡 = 𝟎𝑆 (step  1),  the  set ℒ is  initialized  to  include  all  basis vectors in 𝐖 (step 2). The quantity 𝑙 stores the value of ‖𝐡‖0 ; therefore it is initialized to 1 (step 3). The subroutine maintains a Boolean variable terminate to indicate if the termination condition of the algorithm is satisfied; it is initialized to FALSE (step 4). Since 𝐖 is  a  sparse  binary  matrix,  the  sets 𝒟𝑘 ( 𝑘 = 1,2, … , 𝑆 are initialized to indicate the elements in column 𝐰𝑘 in 𝐖 that contain 1s (step 5).

During each iteration of the while loop (step 6), an ℎ𝑘 is updated to unity and the corresponding index 𝑘 removed from ℒ .  In  other  words,  the  norm ‖𝐡‖0 is  increased  by  unity  per iteration.  As  a  hill-climbing  procedure,  in  each  iteration  the subroutine  picks  an  index 𝑘 from ℒ that  lowers  the  error ‖𝐫 - 𝐖𝐡‖ 2 2 by  the  maximum  amount.  In step 6a, the  current error is obtained and stored in the variable 𝜑0 . In the for loop that  follows  in  step  5,  for  every 𝑘 in ℒ ,  the  error 𝜑1(𝑘) that would result if the corresponding ℎ𝑘 were to be incremented to unity,  is  computed.  However,  before  updating  the ℎ𝑘 , the termination condition  is  evaluated.  The  indicator  variable terminate is set to TRUE if ‖𝐡‖0 is equal to 𝐿 , in which case no further  updates are possible. The variable  becomes TRUE also if incrementing  any other ℎ𝑘 will only increase the error, that is, if the smallest entry of the vector 𝛗1 exceeds 𝜑0 .  The latter situation arises when all elements of the residual 𝐫 are less than 0.5 . This is shown in step 6e.

If the termination condition is not satisfied, the subroutine proceeds  by  obtaining  the  index 𝑘 that  corresponds  to  the smallest 𝜑1(𝑘) (step 6f) and sets that ℎ𝑘 to 1 (step 6g). Next, the  same  index 𝑘 is  removed  from  the  set ℒ ,  in  step  6f.  The variable   in incremented to indicate the new value of 𝑙 ‖𝐡‖0 .

Fig. 3. Subroutine for shiftable load weight updates.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00142v1-figure-13b37e30-26fe-462f-92ba-c4ce45b4c460.jpg)

It  is  clear  that  the  subroutine  is  a  hill-climbing  heuristic approach  to  update  the  weights  of  the  shiftable  load  basis vectors. Therefore, upon termination, hillClimb() returns the binary vector 𝐡 ≡ 𝐡 𝑗 𝑠 (𝑛) such that,

$$h = \underset { \| h \| _ { 0 } \leq L } { \arginf } \| r - W h \| _ { 2 } ^ { 2 } \,. \quad \quad ( 2 3 )$$

To  the  best  of  the authors' knowledge,  there  are  no established upper bound estimates on the computational complexity of mean field annealing, as has been used in [15]. Although  it  is  reasonable  to  consider 𝑂(|𝒮 𝑗 |) as the per iteration complexity in extracting the largest singular value of a sparse matrix (c.f. [19]), this is only an estimate; SVD algorithms  as used in [13] for NMF, do not have strict upper bounds. In  contrast,  the  hill  climbing  procedure  described in this  research  requires 𝑂(𝐿0 ln|𝒟|) steps  per  iteration  and 𝐿0 can never exceed the total number of basis vectors |𝒮 | 𝑗 .

## D. NMF Algorithm

The proposed NMF algorithm follows a block coordinated descent  method,  with  the  matrices 𝐖 𝑓 , 𝐇 𝑓 ,  and  each 𝐇𝑗 𝑠 . Whereas the  fixed  load  matrices  are  treated  as  in  [3],[4],[5], updating  the  weight  matrices  of  the  shiftable  loads,  a  novel feature  of  this  research,  is  done  by  calling  the  subroutine hillClimb() described above multiple times in an iterative fashion. The overall algorithm is outlined in Fig. 4.

For simplicity, the initialization steps are omitted. At first, the matrices 𝐖 𝑓 , 𝐇 𝑓 are updated in accordance with (13) and (14). The updates are carried out in steps 1 and 3, with step 2 being the implementation of the normalization step in (15). It may be noted that the informal notation used in step 2 should, in reality, be implemented iteratively to normalize each column of 𝐖 𝑓 separately.

In each iteration of step 4, a sample 𝑛 ∈ 𝒩 is picked, either randomly or sequentially, so that all shiftable load weights can be updated using the hill-climbing procedure. The appearance of a second level for loop in step 5 is due to the presence of more than a single shiftable load. As before, the inner loop may proceed in a sequential manner, or in any random order of the shiftable loads.

Given sample 𝑛 and a shiftable load  , the residual vector is 𝑗 extracted in the manner shown in (17). This is implemented in step 5a. The residual is normalized in step 5b, with respect to its peak 𝑝𝑗 . Finally in step 5c, the column 𝐡𝑗 𝑠 (𝑛) of the weight matrix 𝐇𝑗 𝑠 is updated by calling hillClimb() .

$$\text{is peak $p_{j}$ family in step $c$, the column $n_{j}(n) or the weight} \\ \text{matrix $H_{j}$ is updated by calling hillClimb()} \,. \\ \text{while converged  ==  FALSE  do} \\ 1. W ^ { f } \leftarrow W ^ { f } \frac { ( X \circ \tilde { X } ^ { \circ - 1 } ) H ^ { f } } { 1 _ { D \times N } H ^ { f } } \\ 2. W ^ { f } \leftarrow \frac { W ^ { f } } { \| W ^ { f } \| } \\ 3. H ^ { f } \leftarrow H ^ { f } \frac { W ^ { f } ( X \circ \tilde { X } ^ { \circ - 1 } ) } { W ^ { f } 1 _ { D \times N } } \\ 4. for  each n \in N \\ \text{is $5.for  each j (in  any  order)} \\ \text{is $5a.r_{j}(n) \leftarrow x(n) - W ^ { f } h ^ { f } ( n ) - \sum _ { j ^ { \prime } \neq j } W ^ { s } h ^ { s } ( n ) } \\ \text{runs} \\ \text{$5b.r_{j}(n) \leftarrow \frac { 1 } { p _ { j } } h ^ { r } ( n ) \\ \text{$5c.h_{j}(n) \leftarrow \, \text{hillClimb} ( W ^ { s }, r _ { j } ( n ) ) } \\ \text{  end  for} \\ \text{no} \\ \text{end  while} \\ \text{per} \\ \text{$Fig.4.Section of the proposed NMF algorithm containing the} \\ \text{of} \\ \text{main loop.} Proceeding initialization steps have been omitted.}$$

Fig. 4. Section of the proposed NMF algorithm containing the main loop. Preceding initialization steps have been omitted.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00142v1-figure-e368016a-ef2a-4a8c-a8ac-02cecb326952.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00142v1-figure-c5b91687-7e98-4fe2-ac30-b79bcc87ad1a.jpg)

Fig. 5. Disaggregated shiftable loads (top) and aggregate loads (bottom) for Day-1.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00142v1-figure-2207f17e-b31d-4741-93e0-baaf1b106e93.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00142v1-figure-70a4570d-5ff0-4776-9de9-669c5a3c28e3.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00142v1-figure-b218dc6c-6996-44cd-9bac-a3545cac0698.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00142v1-figure-7785fe54-3ef7-4df5-a391-c8849e04f49d.jpg)

Fig. 6. Disaggregated shiftable loads for Day-1.

Fig. 7. Disaggregated shiftable loads (top) and aggregate loads (bottom) for Day-2.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00142v1-figure-90f34820-4938-4a11-b2e9-a7a1d34780c0.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00142v1-figure-cbd9b0ad-a8a9-41a2-9d77-14e926c8586b.jpg)

Fig. 8. Disaggregated shiftable loads for Day-2.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00142v1-figure-fc2af7df-0e41-4ec3-9a8e-b0a47c009f79.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00142v1-figure-1757e2e3-3615-408f-934f-ef855dd89426.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00142v1-figure-be72ee4f-0933-48a8-bf6f-a49addf21c51.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00142v1-figure-d52ea449-4d48-4e21-ab7a-d50631300b9e.jpg)

## A. Data Preprocessing

The proposed approach was tested on actual energy usage profiles of a single residential customer that was obtained from the  Pecan  Street  Inc.,  Dataport  database  [20]  sampled  at  one minute intervals, and for the first three weeks in April, 2019.  As energy consumption patterns on weekdays differ significantly from those in weekends, the latter was discarded, yielding a total of 𝑁 = 15 samples.  The  data  was  arranged  as  a 1440 × 15 input  matrix 𝐗 whose  columns  were  the  1440-dimensional sample vectors 𝐱(𝑛), 𝑛 ∈ 𝒩 .

The database in [20] also included individual energy usages of the following four appliances, ( ) a furnace, ( i ii ) a washer/dryer unit, ( iii ) an oven, and ( iv ) kitchen appliances. These measurements provided the basis to evaluate the quality of the disaggregation  obtained  by  the  proposed  algorithm.  All  four appliances were classified as shiftable loads for this purpose and indexed in the above order. The rest of the aggregate load was treated as fixed loads.

The  peak  values 𝑝𝑗 of  the  duty  cycle  of  each  load 𝑗 ∈ {1, 2, 3, 4} ,  as well as the maximum number of its ON cycles, 𝐿𝑗 , of a typical day (denoted as 𝐿 𝑗 max ) was determined through visual inspection. These parameters are shown in Table I.

TABLE I. DUTY CYCLE PARAMETERS

|         |   furnace |   washer/dryer |   oven |   kitchen apps |
|---------|-----------|----------------|--------|----------------|
| 𝒑 𝒋     |     0.465 |            2.5 |      5 |           0.37 |
| 𝑳 𝒋 𝐦𝐚𝐱 |   150     |           20   |     10 |          60    |

## B. Load Disaggregation

The observed 𝐿𝑗 maximums (denoted as 𝐿 𝑗 max ) of the loads in Table I served as upper limits on the number of ON cycles, so that constraints of the form ‖𝐡 𝑗 𝑠 (𝑛)‖ 0 ≤ 𝐿 𝑗 max were applied to each shiftable load. The total number of basis vectors were chosen to be |ℱ| = 1 for the fixed loads, and |𝒮 | = 1440 𝑗 for each shiftable load  . The disaggregated outcomes of the days 𝑗 indexed 𝑛 = 1 and 𝑛 = 5 (indexed  after  dropping  weekends) were picked for illustrative purposes. These will be referred to hereafter as Day-1 and Day-2.

The plots in Figs. 5 and 6 show the results obtained by the NMF  algorithm  for  Day-1.  The  disaggregated  fixed  load 𝐖 𝑓 𝐡 𝑓 (1) (red) and  the  sum  of  all  four  shiftable loads, ∑𝐖𝑗 𝑠 𝐡𝑗 𝑠 (1) (blue) for the entire 24 hour period appears in Fig. 5.  (top).  For  the  sake  of  comparison,  the  expected  aggregate load 𝐱 ̃ (1) from the algorithm (blue) is plotted alongside the real load profile 𝐱(1) (red) in Fig. 5 (bottom). The effectiveness of the proposed approach is evident from the similar patterns of both plots. The peak consumption in both cases occur in the late evening hours (1320 -1440 mins). Additionally the real data shows increased energy usage in the morning hours (540 -600 mins), which is effectively reproduced by the NMF algorithm.

The  disaggregated  loads  of  the  individual  appliances  for Day-1  are  provided  in  Fig.  6.  The  plots  for  the  furnace, washer/dryer unit, oven, and kitchen appliances appear in order (top  through  bottom).  The  actual  load  data  of  this  24  hour period appears in red.

Upon  close  observation,  it  is  apparent  that  the  NMF algorithm accurately reproduces the ON and OFF periods for each load (blue). During Day-1, there were no ON periods for the oven in either case. The ON periods for the washer/dryer unit occur in the interval 1340 -1400 mins. The effectiveness of the NMF  algorithm  is  obvious  from  the  very  strong resemblance of the disaggregated washer/dryer usage to the real data.  Likewise,  the  NMF  algorithm  faithfully  reproduces  the real furnace usage profile, albeit to a somewhat lesser extent than  before.  In  comparison  to  the  others,  there  are  some discrepancies  in  the  usage  profile  that  the  NMF  algorithm yields and the real data for the kitchen appliances.

Figs. 7 and 8 pertain to the energy usage occurring in Day2. As the results are remarkably similar to those in Day-1, we focus only on a few key observations. In contrast to the previous results, from the actual data in Fig. 7 (bottom) it is seen that the increased usage of shiftable loads during the evening is more spread out. This feature is reflected in the disaggregated signals. Unlike before, there are a few ON cycles for the oven, which is again captured by the algorithm. As in the previous case, it can be seen that the algorithm is not able to pick the ON cycles of the  kitchen  appliance  with  the  same  precision  as  with  other shiftable loads. We attribute the difference to the smaller peak value 𝑝4 = 0.37 relative to those of the other appliances.

## V. CONCLUSION

The major contributions of this study are as follows.

- ( )  The load disaggregation of real input data at a significantly i high resolution of 1 minute intervals clearly demonstrates the efficacy of the NMF algorithm suggested in this research.
- ( ii ) Categorizing individual appliances into fixed and shiftable loads allowed the proposed NMF based algorithm to leverage the observe d characteristics of each appliance's ON -cycle.
- ( iii ) It has been shown that the NMF algorithm, including the treatment  of  the  shiftable  loads,  can  be  viewed  in  terms  of maximum likelihood. This provides a theoretical justification of the new approach introduced in this research.
- ( iv )   The hill climbing heuristic, requiring up to only 𝐿0 steps per iteration  of  the  outer  loop,  offers  a  significant  computational advantage  over  other  NMF  approaches  if  used  for  a  similar application. (As an example, this algorithm requires up to 𝐿0 = 𝐿1 max = 150 steps for the furnace, whereas an SVD based NMF [13], would need |𝒮 | = 1440 1 steps).

In spite of the high fidelity of the disaggregated loads obtained using the approach introduced in this research, there is ample scope  for  future  research.  Incorporating  a  semi-supervised learning algorithm to obtain best fits of the duty cycles would obviate  the  need  of  prior  knowledge  of  the  duty  cycles  of individual  appliances.  Another  potential  enhancement  is  to allow the algorithm to learn the ON cycles' probabilities of the appliances.  The  algorithm  may  be  extended  to  include  these probabilities within the theoretical framework.

## REFERENCES

- [1] M. Matsumoto, Y. Fujimoto, and  Y.  Hayashi, ' Energy  disaggregation based on semi-binary NMF,' Proc. International Conference on Machine Learning and Data Mining in Pattern Recognition , pp. 401 -414, 2016.
- [2] A.  Miyasawa,  Y.  Fujimoto,  and  Y.  Hayashi ,  ' Energy  disaggregation based  on  smart  metering  data  via  semi-binary  nonnegative  matrix factorization ,' Energy and Buildings , vol. 183, pp. 547 -558, 2019.
- [3] D.  D.  Lee,  and  H.  S.  Seung,  'Learning  the  Parts  of  Objects  by  Non -negative matrix f actorization,' Nature , vol. 401, no. 6755, pp. 788, 1999.
- [4] D.  D.  Lee,  and  H.  S. Seung,  'Algorithms  for  non -negative  matrix f actorization,' Advances in Neural Information Processing Systems , pp. 556 -562. 2001.
- [5] N. Gillis, 'The How and Why of Non -Negative Matrix Factorization ,' ArXiv , arXiv:1401.5226v2 [stat.ML], 2014.
- [6] C. Févotte, N. Bertin, and JL Durrieu, 'Nonnegative matrix factorization with the ItakuraSaito divergence: With application to music analysis,' Neural Computation , vol. 21, no. 3, pp. 793 -830, 2009.
- [7] C. Févotte, and J. Idier , ' Algorithms for nonnegative matrix factorization with the β -divergence ,' Neural Computation ,  vol.  23,  no. 9,  pp. 2421 -2456, 2011.
- [8] H.  Kim,  and  H.  Park , ' Nonnegative  matrix  factorization  based  on alternating nonnegativity constrained least squares and active set method ,' SIAM Journal on Matrix Analysis and Applications , vol. 30, no. 2, pp. 713 -730, 2008.
- [9] C-J Lin, ' Projected gradient methods for nonnegative matrix factorization ,' Neural Computation , vol. 19, no. 10 pp. 2756 -2779, 2007.
- [10] A. Cichocki, R. Zdunek, and SI. Amari, 'Hierarchical ALS algorithms for nonnegative matrix and 3D tensor factorization,' Proc. International Conference on Independent Component Analysis and Signal Separation , pp. 169 -176, Springer, Berlin, Heidelberg, 2007.
- [11] N.  Benjamin  Erichson,  Ariana  Mendible,  Sophie  Wihlborn,  J.  Nathan Kutz, 'Randomized nonnegative matrix f actorization,' Pattern Recognition Letters , vol. 104, pp. 1 -7, 2018.
- [12] P. O. Hoyer, 'Non -negative matrix factorization with sparseness constraints,' Journal of Machine Learning Research ,  vol.  5, pp. 1457 -1469, 2004.
- [13] R. Peharz, and F. Pernkopf, 'Sparse nonnegative matrix factorization with ℓ 0 -Constraints,' Neurocomputing , vol. 80, no. 1, pp. 38 -46, 2012.
- [14] M. Aharon, M. Elad, and A. Bruckstein, 'The K -SVD: An algorithm for designing of overcomplete dictionaries for sparse r epresentation,' IEEE Transactions on Signal Processing , vol. 54, no. 11, pp. 4311 -4322, 2006.
- [15] A. Rahimpour, H. Qi, D. Fugate, and T.Kuruganti , ' Non-intrusive energy disaggregation  using  non-negative  matrix  factorization  with  sum-tok constrain t,' IEEE  Transactions  on  Power  Systems ,  vol.  32,  no.  6,  pp. 4430 -4441, 2017.
- [16] A.  Miyasawa,  M.  Matsumoto,  Y.  Fujimoto,  and  Y.  Hayashi ,  ' Energy disaggregation  based  on  semi-supervised  matrix  factorization  using feedback information from consumers ,' IEEE PES Innovative Smart Grid Technologies Conference-Europe , pp. 1 -6, 2017.
- [17] B. Humala, A. S. N. Nambi, and V. R. Prasad , ' UniversalNILM: A semisupervised  energy  disaggregation  framework  using  general  appliance models ,' Proc. 9 th International Conference on Future Energy Systems , pp. 223 -229, 2018.
- [18] M.  Figueiredo,  Marisa,  B.  Ribeiro,  and  A.  de  Almeida ,  ' Analysis  of trends in seasonal electrical energy consumption via non-negative tensor factorization ,' Neurocomputing , vol. 170, pp. 318 -327, 2015.
- [19] J.  Baglama,  and  L.  Reichel,  ' Augmented  implicitly  restarted  Lanczos bidiagonalization methods, ' SIAM Journal on Scientific Computing , vol 21, no. 1, pp. 9 -42, 2005.
- [20] ' Pecan Street Inc. Dataport, ' [Online]. Available: https://www.pecanstreet.org/dataport/. [Accessed 27 07 2019].
Ahmad Khaled Zarabie Sanjoy Das
An L0-Norm Constrained Non-Negative Matrix Factorization Algorithm for the Simultaneous Disaggregation of Fixed and Shiftable Loads
Energy disaggregation refers to the decomposition of energy use time series
data into its constituent loads. This paper decomposes daily use data of a
household unit into fixed loads and one or more classes of shiftable loads. The
latter is characterized by ON OFF duty cycles. A novel algorithm based on
nonnegative matrix factorization NMF for energy disaggregation is proposed,
where fixed loads are represented in terms of real-valued basis vectors,
whereas shiftable loads are divided into binary signals. This binary
decomposition approach directly applies L0 norm constraints on individual
shiftable loads. The new approach obviates the need for more computationally
intensive methods e.g. spectral decomposition or mean field annealing that have
been used in earlier research for these constraints. A probabilistic framework
for the proposed approach has been addressed. The proposed approach s
effectiveness has been demonstrated with real consumer energy data.