1 Introduction

Electricity markets around the world focus on demand response (DR), energy efficiency and the integration of renewable sources. There are several reasons for this: technology, economics and environment. Regarding the foreseeable success of European Union (EU) policies in 2020, some independent reports [1, 2] show that DR alone could reach 25%–50% of the EU’s 2020 targets on energy saving and CO2 emission reductions. This could also pre-empt the need for up to 200 medium size power plants in the EU. A recent survey data [3] indicates that smart meter penetration rates continue to increase in the developed countries or unions (USA, Canada, Australia or the EU). In the USA, the penetration in residential sectors (37.8%) was slightly greater than in commercial or industrial segments (about 36.1% and 35.2% in 2013). For example, Edison Foundation’s Institute for Electric Innovation gauges that more than half of all USA households will have advanced meters (estimation for 2015). According to European Commission (EC) estimates, its energy policy represents the installation of about 195 million smart meters for electricity by 2020 (ca. 72% of EU-28 consumers) and an accumulated investment of €35 billion [4].

The problem is that the potential of smart metering may not benefit consumers if they are not well informed about the advantages and possibilities of these technologies. Moreover, customers need to be interested in developing the potential benefits that are possible thanks to the deployment of these “enabling” tools in the current decade. Moreover, customers can be engaged in DR in response to an economic signal (e.g. energy price or incentives) or in response to the need for sustainability, but to do so they need information about their possibilities and potential (i.e. they need information about their load, costs and their potential for response).

On the other hand, electricity companies and commercializers in the EU have limited experience in selling complex services to their smaller customers (e.g. DR products) [5–7]. Moreover, household segments account for 30% of the overall energy demand in 2010 in the EU, with an increasing trend of 1.6% between 2008 and 2009 [8]. Therefore, more knowledge and experience need to be acquired and new information made available through the deployment of smart meters and “enabling” technologies in the EU and USA [9]. Unfortunately, it is difficult to evaluate the individual load behavior; especially in the case of loads such as space heating, air conditioning and water heating, which sometimes are the first candidates to be used in DR, and especially for new products in new markets such as capacity and ancillary services, where interesting benefits could arise for the demand-side.

This work explains and proposes a method to achieve the disaggregation of a customer’s overall demand from the information supplied by a meter, to ascertain the main characteristics and behaviors of the elemental loads. One difference of the proposed method is that the integral transform (Hilbert) and the average frequency, filter the individual square/pulse components, according to their decreasing amplitudes (i.e., it is not based on frequency octaves usually applied in wavelets or Fourier transforms). The proposed methodology could be very useful to boost smart meter potential for demand aggregators and small customers [2]. There are several reasons why the aggregator needs to determine the load composition and response: to choose the most suitable DR policy for a customer load profile depending on the season or the time of the day (critical peak pricing or price-DR tariffs); to evaluate and verify (a system operator concern in DR markets) the customer response to an event; to determine billing and expenses that are shared-out between the customer and its aggregator; and finally, customer education and information on energy use, DR products and possibilities.

This paper is structured as follows: Sect. 2 briefly presents the need for the load disaggregation analysis; Sect. 3 deals with the technology and environment being used for data monitoring and validation purposes; Sect. 4 presents the Hilbert transform (HT) as a tool to analyze non-stationary signals, introduces the filtering process and discusses the information that can be extracted from load measurements and end-uses filtering (besides, some limitations of the proposed method are discussed); Sect. 5 shows how concepts explained in Sect. 4 are applied to real demand profiles and customers. Finally, Sect. 6 presents the conclusion.

2 Load disaggregation

Energy planning of DR programs first involves the determination of the percentage of energy by end-uses for each customer, time of use and the demand levels for each of these uses. Second, their response capacity and the time in which the load effectively changes from the “on” to “off” state are required to build a model for the customer (i.e. which policies are the most effective) and also to verify the DR when the system operator sends regulation signals (for example in ancillary services). Utilities and aggregators usually have aggregated power measurements from customers. The easiest way to identify final end-uses is to insert an intermediate monitoring device between the plug and the appliance. This method is known as “intrusive monitoring”. It is still an expensive method and it should be studied in detail (cost vs. effectiveness) for large scale deployment [10]. For small customer segments, it is very complex and costly to deploy such a system (for example, it costs about $50–60 for plug energy meters, $80–200 for the gateway and at least $100 for a basic control and measurement software license). The idea of the methodology proposed in this work, is to explore and recover the investments already done in information and communication technologies (ICTs) that have been deployed on the demand side, i.e. smart meters. Other solutions, proposed by some evaluation protocols for energy efficiency, suggest (for instance, due to the seasonality of energy usage) that some loggers may be placed on energy-efficiency appliances in a random customer sample [11].

The evolution of technology and public authority support of smart metering infrastructure have added value to the measurement in a single centralized point without having to access to individual sockets. This is usually called “non-intrusive load monitoring (NILM)”, resulting in lower implementation costs and less invasive and bothersome solutions for the end-users. In the last decades, new assessment tools have been proposed and tuned to obtain load features [12] or the so-called “load signatures (LS)” [10, 13] from power and waveform measurements at the customer level.

NILM methods can be classified in terms of data sampling and the level of the method supervision. From the point of view of data sampling, these methods are classified into micro and macro-level signatures [13]. Any measurement with a pacer triggering faster than 1 sample per cycle (20 ms at 50 Hz) is considered to be at a micro level, whereby an attempt is made to identify the characteristics of the load waveform (for instance, its harmonics). A measurement that is slower than 1 sample per cycle is considered to be at a macro level. This last level is the one studied by Hart in [14]. With regard to the level of supervision of the method, they can be classified into supervised and unsupervised methods [15]. Taking into account that customers face educational barriers (the understanding of energy and markets), the number of loads needed to participate in markets (the minimum size of energy bids), and the reasonable performance of most of unsupervised methods, this work proposes a methodology which requires limited supervision.

One objective of this paper is to find LS for the end-uses with potential for DR purposes in small customer segments, according to the literature [12]. A DR waveform with own and forced duty cycles embedded at the customer metering point is shown in Fig. 1. It is not trivial to identify this load working at the same time as other loads.

Fig. 1
figure 1

Aggregated load for a residential customer

Full size image

The first approaches at solving the problem were based on the edge detection of power draw levels (real and reactive) and subsequent clustering [14]. However, they were not able to detect appliances consuming similar power or variable-load appliances. Unlike the previous methods, there have been techniques that relied only on real power measurements (easier and cheaper technology). A technique to disaggregate large loads is shown in [15]; the 15 min sampling requiring limits is due to former pacer triggering capacity of monitoring devices (phased out by ICT). In [16, 17] end-uses, very close to the ones treated in this paper, were analyzed; however, appliance-specific decision rules and excessive training need to be developed. There are also other methods that do not require training [18, 19], instead they create a frequency analysis of power changes, and use an optimization algorithm to match the “on/off” events of a large set of appliances. As a drawback, the algorithm nature of the optimization may or may not provide the best solution. A widely-used method for extracting appliances demand is the hidden Markov model (HMM) [20, 21]. HMM defines a number of states in which the model can be moved, representing the operating “on” and “off” states of the load (see Fig. 1), and perhaps possible intermediate states. It requires a previous statistical study of the individual loads in the household, which is not suitable for the analysis of individual loads that change their behavior due to DR policies. About micro level analysis, there are examples of Fourier analysis [22] or noise [23] in literature. They require a faster sampling and do not provide some of the macro level features we require for our analysis, such as pulse width or cycling frequency. A more exhaustive review of the NILM analysis can be checked in [12, 24]. Some example of data set for NILM research can be found in [25]. Finally, a better performance can be achieved by introducing a feature extraction stage: i.e., the application of a transform operator (Fourier, wavelet, S Transform, Hilbert) [24].

Thus, our objective is not only to disaggregate each end-use, but also to analyze the loads under DR policies by establishing the key variables such as: power level, cycling frequency and change of “on” states which can fluctuate due to DR orders.

A well-known mathematical tool, the HT has been used in this study. In power systems, a refinement of HT, Hilbert-Huang Transform (HHT), has helped to solve other problems [25, 26]. In [27], HT is used to extract the individual demand for water heater and electric heater from aggregated demand for energy efficiency and DR uses through the use of the analytic signal s(t). Some problems of this approach are that errors from 10% to 15% in power are reported, and moreover, the fact that only the two loads with highest amplitudes are extracted.

However, HHT methods use the so called empirical mode decomposition (EMD) and masked empirical mode decomposition (M-EMD), which are not effective with pulse-square signals. Besides, EMD does not provide other important information, such as the pulse width. This paper will demonstrate how HT (without conventional HHT tools) can be used for NILM purposes.

3 Datasets and monitoring system prototype

Every problem to be solved through NILM methods requires the availability of data for algorithm parameterization and validation. In most cases, the datasets are recorded in the USA (for example, the REDD and BLUED database [28]), with some examples for Canada and EU countries, for example UK-DALE and Tracebase [20]. These database are of interest for NILM but in some cases, some problems arise: the duration of tests (days or months) and the availability of data while load is under DR policies. The problem is that the appliance patterns or power levels can change due to weather variations during the year or season (solar radiation or temperature), internal loads (inhabitants patterns, other appliances in service) or due to the response to DR policies. In other database, the number of houses (a single customer) or the appliance and aggregated measurement resolution (1–15 min) which makes it difficult to achieve an “ideal” dataset. For these reasons, a specific dataset has been developed for this project and will be licensed in the near future for research purposes.

3.1 Monitoring system prototype

The future of smart metering will be bright, but from a practical point of view it is quite possible that the customer still does not have a smart meter, or does not have access to detailed meter data, or if they exist, there may be an issue of data confidentiality between the commercializer of energy and the aggregator. In this way, the access to demand data is a cornerstone in DR systems.

Taking these problems into account, the proposed monitoring system has the following abilities: enables an optical output reading. The sensor readings have been carried out by a microcontroller board based on Arduino Uno. This board can read other data of interest such as temperature or humidity. The exchange of data from the Arduino to the PC (customer or aggregator) was performed using XBee modules (Digi International). In the event that a customer has his own meter, the procedure is similar but in these cases the system directly records the electric pulse output of the meter (in our case, 500 impulse/(kWh)). If the customer meter is not accessible for measuring the optical output, or an electronic meter is not fitted in his home, the easiest and cheapest option (around €150 for a three-phase unit) is to install a meter with current clamps and signal transmission capability. In our case, the power meter has an emitter, which sends information to a PC through an USB gateway (Z-stick) using Z-wave protocols (with a cost about €50). Three households have been monitored for one year using the three methods.

Another problem is the need for synchronized measures from aggregated demand and individual appliances to validate the proposed NILM methodology. For this purpose, four customers have been monitored with an electronic meter and five Z-wave wall plug meter switches (manufactured by Fibaro [29]) which send data to a PC using a USB gateway.

The hardware is managed with the help of software developed on the IP-Symcon platform [30] and allows sampling from seconds to one minute. In our case, samples of 5 s, 10 s and 1 min pacer trigger have been used.

3.2 Customers’ characteristics

Four typical residential customers in Spain have been selected for test purposes (from 3.3 to 10 kW of rated power, single and three phases, with an average energy demand from 400 to 800 kWh/month).

Main end-uses are presented in Table 1 for customer number 4 (customer with individual appliance monitoring) and some representative end-uses for the other customers. It is necessary to state that some end-uses exhibit a few discrete states (on-off states) in demand while others have multiple states or continuous (for instance HPAC) with a wider range of values (see Table 1). Typical household end-uses in the European Union EU-28 member states and the USA [31] are considered when selecting test customers, in order to demonstrate that the validation scenario and the residential customer are standard.

Table 1 Main end-uses and their rated power
Full size table

Figure 2 covers some typical end-uses. Accordingly, it is clear that WH, EH, RF and WM have a cycling demand. The frequency, amplitudes and “on/off” times are quite different, but it is really difficult to extract end-uses, except for certain time periods (late at night, early morning, etc). One feature common to all is that their demand evolves according to square/pulse waves.

Fig. 2
figure 2

Some examples of end-uses demand

Full size image

4 Demodulation and filtering of pulse waveforms

4.1 HT principle

HT is a mathematical transform of interest in signal processing which is valuable due to its ability to analyze non-linear and non-stationary signals. For instance, it has been used to study power and oscillatory transients [25, 26]. For any real time series g(t), its transform H(g(t)) is defined as:

$$H(g(t)) = \frac{1}{\uppi}PV\int\limits_{ - \infty }^{\infty } {\frac{g(u)}{t - u}{\text{d}}u}$$
(1)

where PV is the Cauchy principal value of the integral. After the HT is defined, its analytic signal based on real time series with its transform can be defined:

$$s(t) = g(t) + {\text{j}}H(g(t)) = A(t){\text{e}}^{{{\text{j}}\theta (t)}}$$
(2)

where A(t) and θ(t) are the instantaneous amplitude and phase. The instantaneous frequency f(t) can be defined as:

$$f(t) = \frac{w(t)}{2\pi } = \frac{1}{{2\uppi}}\frac{{{\text{d}}\theta (t)}}{{{\text{d}}t}}$$
(3)

A complete definition can be found in [32]. However, as stated in [33], in order to obtain a meaningful frequency f(t) by HT, g(t) must be “mono-component”. This means it can only have one frequency value at any given time, otherwise the resulting instantaneous frequency would have no significance. As noted by Huang [33], a mono-component signal must be locally symmetric with respect to the local zero mean to achieve a meaningful f(t). This fact highlights the need to extract several components that achieve this. The EMD algorithm extracts these mono-component functions called intrinsic mode functions (IMF). However, the EMD generally does not always provide mono-component signals [26]. To solve this problem, M-EMD has been proposed to force the component extraction [34]. The main problem of the masking process is that it is of heuristic nature and, in this way its practical application is very limited.

The extraction of mono-component (frequency) functions with EMD is not really relevant for our application. The problem is that conventional EMD cannot decompose pulse-square signals (with a physical sense) due to the nature of its algorithm: based on envelopes of extreme points. Thus, the decomposition spoils the square shape of the signal, corrupting the result as it always tries to extract sine-shape components. In a certain way, HHT tries to work in a similar way to that used by Fourier: mathematically expressing a signal as an addition of pseudo-sinusoidal components. This is a problem for pulse-square waves, because it merges signal components of diverse end-uses and makes disaggregation harder. Therefore, to extract the information contained in an aggregated signal, a different and direct approach had to be developed. The proposed method is based on the approach presented in [27] but performs a modified and improved decomposition and filtering method. The advantages and characteristics of the proposed method are as follows.

  1. 1)

    It is based on the properties of HT and not only in the values of s(t): this allows a better performance of signal extraction and filtering (the behavior of HT is better than the one of s(t) for pulse signals).

  2. 2)

    The amplitude of each signal is evaluated through reconstruction and filtering through inverse transform. This allows to bound the errors in amplitude (from 10%–15% in [27] to 5% in the proposed method).

  3. 3)

    The filtering of HT permits the extraction of several components (not only the components with the highest amplitudes) including filtering of loads with a continuous demand, especially when they are engaged in DR response.

For applying HT, the signal must comply with the IMF conditions. Although this is often true, reference [27] demonstrates that the signal does not have to be an IMF.

4.2 Single square signal

First, a square signal g(t) is considered (note that this upholds the IMF conditions [33]):

$$g(t) = \left\{ \begin{aligned} & A;\; 0 \le t < a + kT_{1} \\ & - A;\; a + kT_{1} < t \le (k + 1)T_{1} \\ \end{aligned} \right.\quad k = 0,1, \ldots ,n$$
(4)

where T 1 is the period; kT 1 and a + kT 1 are the trailing and leading edges of the pulse; A is the pulse amplitude.

Mathematically, g(t) can be developed as the sum of a set of n characteristic functions X [a, b](t) (a pulse with trailing and leading-edges a & b):

$$X_{[a,b]} (t) = \left\{ {\begin{array}{*{20}c} {1;\;a \le t \le b} \\ {0;\;t \notin [a,b]} \\ \end{array} } \right.$$
(5)

With pulse amplitude A:

$$g(t) = \sum\limits_{k = 0}^{n} {(AX_{{[kT_{1} ,\;kT_{1} + a]}} } (t) - AX_{{[a + kT_{1} ,\;T_{1} + kT_{1} ]}} (t))$$
(6)

The transform of X [a, b](t) (see [32]) can be substituted in (6), and then applying linearity properties of HT:

$$H(g(t)) = \frac{1}{\uppi}A\sum\limits_{k = 0}^{n} {\left[ {\ln \left| {\frac{{t - kT_{1} }}{{t - (a + kT_{\;1} )}}} \right| - \ln \left| {\frac{{t - (a + kT_{1} )}}{{t - (T_{1} + kT_{1} )}}} \right|} \right]}$$
(7)

The corresponding analytical signal s(t) in [0, T 1] is:

$$s(t) = \left\{ {\begin{array}{*{20}c} {A + {\text{j}}\frac{1}{\uppi}A\sum\limits_{k = 0}^{n} {\left[ {\ln \left| {\frac{{t - kT_{1} }}{{t - (a + kT_{1} )}}} \right| - \ln \left| {\frac{{t - (a + kT_{1} )}}{{t - (T_{1} + kT_{1} )}}} \right|} \right]} ;\quad 0 < t < a_{{}} } \\ { - A + {\text{j}}\frac{1}{\uppi}A\sum\limits_{k = 0}^{n} {\left[ {\ln \left| {\frac{{t - kT_{1} }}{{t - (a + kT_{1} )}}} \right| - \ln \left| {\frac{{t - (a + kT_{1} )}}{{t - (T_{1} + kT_{1} )}}} \right|} \right]} ;\quad a < t < T_{1} } \\ \end{array} } \right.$$
(8)

for k = 1, and its instantaneous phase θ(t) is given by:

$$\theta (t) = \left\{ {\begin{array}{l} {{\text{arc}}\tan \left( {\sum\limits_{k = 0}^{n} {\left[ {\ln \left| {\frac{{t - kT_{1} }}{{t - (a + kT_{1} )}}} \right| - \ln \left| {\frac{{t - (a + kT_{1} )}}{{t - (T_{1} + kT_{1} )}}} \right|} \right]} } \right);\quad 0 < t < a} \\ {{\text{arc}}\tan \left( {\sum\limits_{k = 0}^{n} {\left[ { - \ln \left| {\frac{{t - kT_{1} }}{{t - (a + kT_{1} )}}} \right| + \ln \left| {\frac{{t - (a + kT_{1} )}}{{t - (T_{1} + kT_{1} )}}} \right|} \right]} } \right);\quad a < t < T_{1} } \\ \end{array} } \right.$$
(9)

The use of the analytic signal changes HT values when g(t) < 0. This property is important to determine the duty cycle of the signal. Figure 3 shows the tangent of θ(t) for a pulse waveform.

Fig. 3
figure 3

Square waveform and its transform

Full size image

From Figure 3, it can be seen that changes in the pulse waveform show discontinuities in its HT. The duty cycle can be obtained through positive to negative transitions of phase (corresponding to trailing and leading edges of g(t)).

4.3 Instantaneous frequency f(t)

The instantaneous frequency f(t) is defined in related literature as the first derivative of θ(t) (see equation (3)). The f(t) has physical sense if and only if g(t) is mono-component and sinusoidal. If g(t) is multi-component or a pulse signal f(t) does not give the frequency of the waveform (see Fig. 4).

Fig. 4
figure 4

Instantaneous frequency of a sinusoidal (mono-component) wave and the frequency of a bi-component wave

Full size image

Nevertheless, it is possible to show that the average of instantaneous frequency (f me) provides information about any waveform [27], in this case, for g(t):

$$f_{\text{me}} = \frac{1}{{T_{1} }}\int\limits_{0}^{{T_{1} }} {f(t){\text{d}}t} = \frac{1}{{2\uppi\,T_{1} }}\int\limits_{0}^{{T_{1} }} {{\text{d}}\theta (t)} = \frac{1}{{T_{1} }}$$
(10)

Therefore, f me(t) equals the frequency of g(t) in certain periods of time. To obtain its f me, the restrictions of the signal s(t) can be reduced: its maximum values are above zero and its minimum values are below zero.

4.4 Analysis of multi-component signals

In real cases, the power waveform at the meter point (or aggregated measurement) is not a single square-pulse signal; it is more complex (see Fig. 2) and involves continuous, multi-state and “on/off” loads (see Table 1). If the customer owns electrical heating, water heating, heat pump (air conditioning), and fridge loads, which are working in the same period, it is foreseeable that two or more square-pulse waveforms could be present in the aggregated load curve.

The average value of f(t), i.e. f me, is seen to provide the cycling frequency of a square signal. How this still applies to multi-component signals is to be demonstrated. Although the proposed method is still appropriate for more than two components, it is sensitive to the amplitude value of the pulse components; so, it extracts the frequency value of the highest power level. Let us consider a signal g(t) with two components: a base (load) level B (DC level), and N-square signals g m (t) defined by decreasing amplitudes and with frequencies 1/T m .

$$g(t) = g_{0} (t) + \sum\limits_{m = 1}^{N} {g_{m} (t)}$$
(11)
$$\left\{ \begin{aligned} & g_{0} (t) = B \\ & g_{m} (t) = \left\{ {\begin{array}{*{20}l} {A_{m} ;\;a_{m} \le t < b_{m} } \hfill \\ { - A_{m} ;\;b_{m} < t \le a_{m} + T_{m} } \hfill \\ {0;\;t = b_{m} } \hfill \\ \end{array} } \right. \\ \end{aligned} \right.$$
(12)

where A m is the amplitude; T m is the period; a m and b m are the trailing and leading edges of the pulse.

Mathematically, g(t) can be developed as the sum of a set of k-characteristic functions \(X_{{[a_{m}}, {b_{m}]}}(t)\) for each square/pulse wave just like in (6):

$$\begin{aligned} g(t) & = g_{0} (t) + \sum\limits_{m = 1}^{N} {g_{m} (t)} = B \\ & + \sum\limits_{m = 1}^{N} {\left( {\sum\limits_{k = 0}^{n} {} A_{k} \left[ {X_{{[a_{k} + kT_{k} ,\;kT_{k} + b_{k} ]}} (t) - X_{{[b_{k} + kT_{k} ,\;T_{k} + kT_{k} ]}} (t)} \right]} \right)} \\ \end{aligned}$$
(13)

The procedure presented at the beginning of the Sect. 4.3 will not be repeated. First, the HT of g(t) should be assessed, then the analytic signal s(t) should be composed, and finally, its phase should be extracted. Afterwards, the tangent of the phase θ(t) for the multi-component signal (i.e. H(g(t)) over g(t) is obtained). The HT of the baseload (DC) signal is zero. This is shown in Fig. 5.

Fig. 5
figure 5

Signal with two pulse components and a base (load) component

Full size image

Note that the values of the HT of g(t) have changed in each side of the points (a 1, b 1 and T 1) of singularity of H(g 1(t)) (for instance, they remain positive in a 1-ε and change to negative in a 1+ε, ε → 0) when s(t) is defined (the waveform with the highest amplitude). Specifically,

$$H(g_{1} (t)) = \left\{ {\begin{array}{*{20}c} { - \infty ;\;t \to a_{1} + kT_{1} } \\ { + \infty ;\;t \to b_{1} + kT_{1} } \\ \end{array} } \right.$$
(14)
$$\tan (\theta (t)) = \left\{ \begin{aligned} \frac{{ -\uppi}}{2};\;t \to b_{1} + kT_{1} + \varepsilon \hfill \\ \frac{{{ + }\uppi}}{2};\;t \to b_{1} + kT_{1} - \varepsilon \hfill \\ \frac{{{ + }\uppi}}{2};\;t \to a_{1} + kT_{1} + \varepsilon \hfill \\ \frac{{ -\uppi}}{2};\;t \to a_{1} + kT_{1} - \varepsilon \hfill \\ \end{aligned} \right.$$
(15)

Whereas the values of singularities of H(g m (t)) (a m , b m and T m ) remain the same on both sides (positive or negative), i.e.,

$$H(g_{m} (t)) = \left\{ {\begin{array}{l} { - \infty;} \\ { + \infty;} \\ \end{array} } \right.\begin{array}{*{20}c} {t \to a_{m} + kT_{m} \pm \varepsilon } \\ {t \to b_{m} + kT_{m} \pm \varepsilon } \\ \end{array} \quad k = 0,1,2, \ldots$$
(16)
$$\tan (\theta (t)) = \left\{ {\begin{array}{*{20}c} {\frac{{ -\uppi}}{2};\;t \to a_{m} + kT_{m} \pm \varepsilon } \\ {\frac{{ +\uppi}}{2};\;t \to b_{m} + kT_{m} \pm \varepsilon } \\ \end{array} } \right.$$
(17)

Note that this result can be applied with several pulse components, if and only if \(\sum\nolimits_{m = 2}^{n} {A_{m} < A_{1} }\), the evaluation of f me is easy because a lot of terms disappear (i.e. the discontinuities of g m (t)). For simplicity, to improve the physical meaning and the understanding of (13), the number of square signals in (12) and (13) is reduced to N = 2 (g 1(t) and g 2(t), see Fig. 5):

$$\begin{aligned} f_{\text{me}} = \frac{1}{{2\uppi\,T_{1} }}\mathop {\lim }\limits_{\varepsilon \to 0} \left( {\int\limits_{0 + \varepsilon }^{{a_{2} - \varepsilon }} {\frac{{{\text{d}}\theta }}{{{\text{d}}t}}{\text{d}}t} + \int\limits_{{a_{2} + \varepsilon }}^{{a_{1} - \varepsilon }} {\frac{{{\text{d}}\theta }}{{{\text{d}}t}}{\text{d}}t} + \int\limits_{{a_{1} + \varepsilon }}^{{b_{2} - \varepsilon }} {\frac{{{\text{d}}\theta }}{{{\text{d}}t}}{\text{d}}t} } \right. \hfill \\ \left. {\quad + \int\limits_{{b_{2} + \varepsilon }}^{{a_{2} + T_{2} - \varepsilon }} {\frac{{{\text{d}}\theta }}{{{\text{d}}t}}{\text{d}}t} \int\limits_{{a_{2} + T_{2} + \varepsilon }}^{{T_{1} - \varepsilon }} {\frac{{{\text{d}}\theta }}{{{\text{d}}t}}{\text{d}}t} } \right) = \frac{1}{{2\uppi\,T_{1} }}\mathop {\lim }\limits_{\varepsilon \to 0} \left[ {\left[ { - \theta (0 + \varepsilon )} \right]} \right. \hfill \\ \quad + \left. {\left[ {\theta (a_{1} - \varepsilon )} \right] + \left[ { - \theta (a_{1} + \varepsilon )} \right] + [\theta (T_{1} - \varepsilon )]} \right] = \frac{1}{{2\uppi\,T_{1} }}2\uppi = \frac{1}{{T_{1} }} \hfill \\ \end{aligned}$$
(18)

This means that the average value of f(t) (i.e. f me) still provides important information about the signal and the frequency of the component of the highest amplitude while masking other components with lower amplitude, i.e. g m (t) (m ≠ 1). Figure 6 shows H(g(t)) and tan(θ(t)) extracted from a non-symmetrical square wave with a continuous component of 5 and two pulse signals: g 1(t) (15.9 Hz), and g 2(t) (39.8 Hz) of amplitudes 30 and 20. Only times with singularities of H(g 1(t)) give a change from positive to negative values in tan(θ(t)), but H(g 1(t)) shows a higher value in these singularities than in times where g 2(t) changes (and H(g 2(t)) exhibits a singularities).

Fig. 6
figure 6

Pulse waveform components

Full size image

4.5 Filtering of pulse components

The disaggregation and identification of elemental load components are based on results presented in paragraphs 4.2 to 4.4. In [28], the instantaneous frequency f(t) and its singularities (Fig. 6) were applied to obtain main characteristics of the appliance with the highest amplitude. In the present case, the information within HT will be considered to filter not only one, but several pulse components, in the aggregated demand curve. To clarify this procedure, an example with a synthetic waveform g(t) with four components has been used (sampling rate 40 kHz). The complex waveform is shown in Table 2.

Table 2 Characteristics of waveform g(t)
Full size table

The software, specifically developed for filtering, considers first the calculation of HT from g(t). This transform is shown in Fig. 7. Theoretically, the local maxima of H(g(t)) are due to H(g 1(t)) (the highest amplitude component). The program computes these maxima, splitting the time-period into subintervals defined considering f(t) behavior and the autocorrelation function (ACF) of H(g(t)), see Fig. 8. With the help of θ(t) (Fig. 6), the condition of maximum amplitude is verified.

Fig. 7
figure 7

Pulse waveform and maxima and minima of its transform

Full size image
Fig. 8
figure 8

Autocorrelation function for the transform of pulse waveform (maximum at lag = 1000 defines period of the first component)

Full size image

With the maximum and minimum values in each subinterval, the time values a and b in equation (7) have been calculated to obtain the characteristic function X [a, b](t) (amplitude is equal to 1). The linear properties of the transform can help the user to get the amplitude of each function g m (t) in (12):

$$A = \frac{H(g(a)) - H(g(b))}{{H(X_{[a,\;b]} (a)) - H(X_{{[a,{\kern 1pt} \;b]}} (b))}}$$
(19)

According to (18), the inverse transform of H(g 1(t)) component gives the filtered component (g F1 (t)) with average frequency f me, irrespective of the values of the other components. Afterwards, this transform is subtracted from HT:

$$\left\{ \begin{aligned} & g(t) = g_{0} (t) + g_{1} (t) + \sum\limits_{m = 2}^{N} {g_{m} (t)} \\ & g_{1}^{\text{F}} (t) = H^{ - 1} (H(g_{1}^{\text{F}} (t))) \\ & H\left( {g_{0} (t) + \sum\limits_{m = 2}^{N} {g_{m} (t)} } \right) = H(g(t)) - H(g_{1}^{\text{F}} (t)) \\ \end{aligned} \right.$$
(20)

The component g 2(t) is now the maximum amplitude component of the residual function and inherits the properties of g 1(t) in the new waveform (average frequency, maxima, minima of HT), see Fig. 9.

Fig. 9
figure 9

Maxima and minima of transform of pulse waveform without the first component

Full size image

The maxima of H(g(t)) due to g 1(t) are shown in Fig. 9 (note that these maxima disappear after the filtering), and the maxima of this transform are due to g 2(t) which drives mean frequency. This process is repeated until the last square-pulse component is extracted (in this case, g 3(t)) (see Fig. 10). Some peaks appear in the extraction process due to errors in the evaluation of time positions (t = a, t = b, see equation (17)), where HT presents singularities. Keep in mind that the sinusoidal component (for example, this case of constant loads in aggregated demand) does not affect the results because HT of continuous and derivable functions do not usually have any singularities (see Fig. 6).

Fig. 10
figure 10

Maxima and minima of transform of pulse waveform without the first and the second component, and the filtered component and measured value (Maxima of HT due to g 2(t) disappear after filtering)

Full size image

4.6 Limitations of the proposed methodology

Obviously, the proposed methodology has some limitations to overcome in the future: firstly, the sampling rate provided by the meter device (from 1 sample/min to 12 sample/min in the proposed dataset) is necessary to capture the behavior of some loads (for example WH-2 in Fig. 2). For this load, 1 sample/min will produce a triangular waveform which transform is more difficult to filter than pulse transforms (this is due to the transform of a pulse signal presents a singularity at leading and trailing edges of square pulses, see Fig. 3, whereas there is not any singularity for a triangular waveform). The second factor to be considered is that the method is not 100% unsupervised at this moment, and a basic but necessary knowledge and expertise of the loads and their dynamics inside the household is needed (the aggregator should develop this expertise). Perhaps the main limitation is the potential presence of several appliances with the same end-uses (for example, some electric heaters to heat several rooms in the house). To explain this problem, three elemental signals with the same frequency (20 Hz) and amplitudes ranging from 10.5 to 10.3 have been considered in Fig. 11 (i.e the square pulse component being analyzed do not fulfil \(\sum\nolimits_{m=2}^n {A_{m} < A_{1} }\)). Here, the problem is that the sum of signals is a new signal (g(t)) with three times the amplitude of the individual components, the same frequency, but with a reduced duty cycle. In this case HT detects the singularities of the “apparent” higher amplitude (g max(t) “signal” in Fig. 7) in g(t) (see the arrows shown in Fig. 11). Thus, the proposed method does not filter the signal with the highest amplitude (g 1(t) in this example) and an error may arise in the disaggregation. An artificial neural network could be used in the future to detect this problem and, in this way, reach a more unsupervised procedure to avoid the aggregator needing a minimum level of experience with modelling and some knowledge of the integral transform tools to evaluate the inconsistency of the results.

Fig. 11
figure 11

Analysis of a multi-component square wave with a three rectangular signals with different duty cycles and amplitudes

Full size image

5 Demand response

5.1 Evaluation of residential end-uses

For the evaluation of end-uses, some measurements (aggregated and individual) from customer 4 were considered. Figure 2 covers some typical end-uses during the winter season for this customer. Accordingly, it is clear that WH, EH, FR, OV and WM have a cycling demand (with two or multiple states). The frequency, amplitudes and “on/off” times are quite different, but it is really difficult to extract end-uses, except in certain time-periods (late at night, early morning, see Figs. 1 and 2). One feature common to these appliances, is that their demand evolves according to square/pulse waves. These patterns of main loads allow advantage to be taken of the results presented in Sect. 4. Moreover, other load patterns (quasi-continuous patterns) arise from filtering (for example HPAC with inverters).

Small customers are becoming increasingly more important for aggregators and system operators, although some feedback is required from the load when DR is deployed or expected to be used, especially in DR for energy and ancillary services to test and probe their performance face to power system events. To explore the possibilities of DR evaluation through information on smart meters and Hilbert tools, some DR policies have been applied to the customers under study. It seems that DR would benefit WH, EH and HPAC the most, because they have high rated power and are often used when a power system is not reliable or the energy price is quite high to persuade customers to become engaged in DR.

For this research, two different days in the winter season have been analyzed: day 0 (a day without DR control at 12 sample/min); and day 1 (where a simple DR policy was applied to HPAC).

The first step is to center the waveform with respect to its maximum and minimum values, or with respect to its “DC value”. This last option can be done numerically or through the use of HT because the transform of a constant (our DC value) is zero. The results (the daily load was previously presented in Fig. 1) can be seen in Fig. 12.

Fig. 12
figure 12

Aggregated load (“day 0”, centered)

Full size image

Afterwards, H(g(t)) and ACF are computed. Results are presented in Fig. 13a, b. The ACF reports to the user (aggregator), that a periodic signal is detected with lag = 60 (note that at 12 sample/min this should be a waveform with a period of 5 min). Then, the instantaneous frequency f(t) and average frequency f me(t) (computed each 5 min or 60 samples) have been computed (see Fig. 13c, d). The frequency f(t) detects that the first component appears first at t = 11 hour and vanishes at t = 15 hour; it appears again at t = 20.5 hour to finish at t = 21.5 hour. Average frequency confirms this statement and shows a different frequency pattern in each period.

Fig. 13
figure 13

Extraction process of the highest amplitude component into the aggregated load

Full size image

With this information, H(g 1(t)) component is filtered and its inverse transform g 1(t) is determined. Some details of results are shown in Fig. 14a, b. Especially in Fig. 14b, the extraction of H(g 1(t)) from the aggregated waveform is quite certain.

Fig. 14
figure 14

Disaggregated load processing: the inverse transform of the component (which corresponds to WH in the morning and OV in the afternoon) and HT residue after the extraction of that component

Full size image

Once the first component has been filtered, the process flows, according to the procedure presented in Sect. 4.5. The autocorrelation function “ACF” of the residual function (H(g(t) − g 1(t)) is processed. In this case, the lag is around 250 (i.e. about 20 min. Note that this pattern also appears in Fig. 15b). After that, the average frequency is calculated. The results are presented in Fig. 15.

Fig. 15
figure 15

Disaggregation process for the second component

Full size image

It is important to pay attention to the patterns present in HT and in its instantaneous frequency f(t), see Fig. 16. From 7 to 12 hour, f(t) positive values are mainly displayed, whereas from 17 to 20 hour, f(t) high negative and positive values appear, see Fig. 16b (negative values in f(t) are responsible for negative values in f me, and this has not any physical sense, this means that the signal has changed or is not centered, see Fig. 12, and this problem is attributable to a change in the amplitude). This gives the user a flag about a change in load pattern, i.e. another load is extracted in the process or the process should be revisited to split the daily demand in two periods to perform the best evaluation of HT and ACF.

Fig. 16
figure 16

Disaggregation of the second component with the help of the instantaneous frequency

Full size image

Day 0 involves a 8/12 natural duty-cycle for EH load from 7 to 14 hour (approximately 8 min ON, 12 min OFF, changing according to outdoor temperature). Figure 17 shows the individual load (recorded with an energy plug) and the load disaggregated with the filtering of components through HT.

Fig. 17
figure 17

Comparison between the second component (extracted through HT) and the real demand of EH (recorded through plug meter)

Full size image

Finally, previous works in NILM have used different metrics to evaluate performance. For example, the proposal in [13] has been widely used; it considers the “detection accuracy (DA)” as the ratio of the correctly classified events to the total detected events (excluding missing detection). In [12, 28] the “percentage of energy correctly classified (ACC)”, considered for evaluation of HMM methodology, defined formally as:

$$ACC = 1 - \frac{{\sum\limits_{t = 1}^{T} {\sum\limits_{k = 1}^{n} {\left| {\bar{e}_{k} (t) - e_{k} (t)} \right|} } }}{{2\sum\limits_{k = 1}^{n} {\left| {e_{k} (t)} \right|} }}$$
(21)

where ē k (t) and e k (t) are the algorithm’s prediction and measurement for the k th load at the t th time step respectively. An important parameter for practical use is the energy consumed by each load k (for example for DR verification). In [21], the performance of individual loads is measured through the “normalized disaggregation error (NDE)”:

$$NDE_{k} = \sqrt {\frac{{\sum\limits_{t = 1}^{T} {\left| {\bar{e}_{k} (t) - e_{k} (t)} \right|^{2} } }}{{\sum\limits_{k = 1}^{n} {\left| {e_{k} (t)} \right|} }}}$$
(22)

Table 3 shows a comparison of performance ranges of the proposed methodology with respect to HMM.

Table 3 Performance ranges of HT and HMM [20, 21]
Full size table

The comparison is difficult because WH, EH or HPAC are not considered in [12, 20, 28]. For example [35] states that a lot of methods decrease performance with air conditioning loads (up to 25%–30%). This is not the case of the proposed approach as demonstrates next paragraph.

5.2 Evaluation of DR: Ancillary services

The proposed methodology can be used to evaluate the compliance of a possible response of loads to regulation or reserve services. For example, in this study (load curve on day 1), the HPAC load is tuned automatically by the aggregator (through the remote change of the thermostat set point at 12 hour) to reduce load and help a non-spinning reserve call commissioned by independent system operator (ISO) in the short term (minutes). To verify and test if this load response was on, EH and WH have to be filtered from the overall load of the customer (the same procedure used in Sect. 5.1). Note that ISO impose high monitoring requirements [36] for ancillary services, and consequently, a statistical solution could be necessary to avoid the use of high performance meters in the entire customer group (i.e. to achieve cost-effectiveness for these customers) [37].

This is a new problem for aggregators in the near future, but it can be solved by a rough iterative calculation of HT for EH and WH loads based on average frequency, instantaneous frequency and filtering. The overall HT is filtered by the subtraction of elemental pulse HT transforms (EH and WH load transforms). The results are shown in Fig. 18.

Fig. 18
figure 18

Aggregated load on day 1 and HT for individual pulse loads (EH, WH) and residual component on day 1

Full size image

The inverse transform (H −1) is necessary to reconstruct any waveform. In the Hilbert domain H −1 is:

$$g(t) = H^{ - 1} (t) = H( - H(g(t))$$
(23)

This equation is applied to the “HT filtered” signal in Fig. 18. This means that the “residual” load profile shown in Fig. 18b can appear by a new iteration of the methodology proposed in paragraph 4.2 to paragraph 4.5. For comparison purposes, the HPAC and residual end-uses are also drawn in Fig. 19a (note that, without a more detailed calculation, the response to the new thermostat set point appears from 12 to 14 hour and from 14 to 18 hour). It can be seen that HPAC load appears (Fig. 19a) together with other end-uses, for instance FR and WM with multiple states of power, but the change of load demand due to DR can be verified. Moreover, those small (residual) demands could be also filtered with a similar procedure. It can be seen that the changes of HPAC load are due to the response to thermostat control.

Fig. 19
figure 19

Inverse transform of residual component (after filtering) and AC inverter load (HPAC) behavior

Full size image

6 Conclusion

This paper presents an alternative method to disaggregate end-uses with cycling “on/off” demand from the data supplied by smart meters. Moreover, the load behavior can be obtained after, during and before a DR policy is implemented to verify if the control related to the suitable cycling load is being applied correctly. Here, HT is shown to provide relevant information through frequency analysis, which can even be used in multi-component signals, to determine pulse-square components. This is done without the traditional use of EMD and IMFs which does not work well with square-pulse signals. It also combines HT properties with auto-correlation analysis to determine other important properties and perform an iterative filtering of pulse components according to decreasing amplitudes. The use of different metrics for the analytical evaluation of results shows that the method has a remarkable performance in comparison with other methodologies, and achieves the disaggregation of important loads for DR, while the accuracy of results remain in similar levels in spite of the presence of quasi-continuous state loads (inverter HPAC). Hence, the information obtained by the proposed methodology from smart meters (aggregated demand) allows the aggregator to verify and evaluate the end-use composition on the customer and the response of selected (representative) customers in each demand segment and in this way, to obtain the maximum available potential from DR.