Energy Consumer Model Estimation Using Machine Learning

1. Introduction

The overarching goal of this research is to develop a data-driven modeling framework of a residential unit’s daily electricity usage. A nonlinear model is proposed here that takes into account time-of-use (TOU) costs and ambient temperatures; it can be useful for a wide variety of tasks [1]. The model can be used by the consumer to detect excessive energy use times, predict monthly budget, in home energy management systems (HEMS), and modifying household activities. It can benefit architects and construction engineers in designing and retrofitting homes, and in the integration of renewable energy. Utility companies can adopt the proposed model to plan in advance installation capacity, analyze TOU pricing policies, improve preparedness for electrical outages, manage infrastructure costs, analyze demand response (DR), incentivize households’ electricity consumption in order to minimize grid stress, and to carry out large scale simulations of neighborhoods and towns for various energy market studies over large time periods.

Unlike a (fully) parametric model which relies on physical parameters of each electrical device (e.g., water heater, space heater, and clothes drier [2]), the model parameters do not reflect appliance level physics. At the same time, it differs from a non-parametric model whose parameters (however many) cannot be interpreted in any meaningful manner. Many such black box models have been proposed in current literature, based on deep neural networks (DNNs) and other machine learning methods (cf. [3]). Ergo, a semi-parametric model is proposed herein, whose parameters are firmly grounded in established econometric and game theoretic principles [4,5,6].

Two broad categories of household electric appliances are parametrized separately. These are, (i) fixed loads, and (ii) shiftable loads. In situations where these loads cannot be monitored separately, a novel algorithm for partial disaggregation of the aggregate load, has been proposed. This step can be used prior to parameter estimation of the utility model. Similar strategies, which do not disaggregate energy usage into individual appliance loads, have been adopted elsewhere [7,8,9,10,11].

Fixed loads are appliances such as lamps, LEDS, heating units, fans, stoves, ovens, refrigerators, freezers, and medical equipment, which can be operated at any time of the day, independently of the TOU cost. For instance, refrigerators and freezers must run throughout the day to prevent food spoilage, inevitably consuming electrical energy even during expensive peak periods. Likewise, lamps, LEDS, and fans are scheduled based on human occupancy and sunlight. Medical devices like CPAP machines or home dialysis units, must run according to patients’ health needs, regardless of TOU cost.

Shiftable loads are appliances such as EV charging, HVAC, and electric water heaters, whose energy usage is intermittent. These include high impact appliances that consume a large amount of energy that can be discretized into a few levels (e.g., OFF/ON or OFF/LO/MED/HI). Smart homes schedule these appliances according to the TOU cost of electricity. For instance, EV charging – often the largest single shiftable load, can be scheduled during overnight “super off-peak” periods, resulting in significant reduction in cost. While HVAC are required for comfort, “pre-cooling” or “pre-heating” a home before peak hours, allows the system to run less when rates are highest.

The schematic in Figure 1 shows how a typical load profile may change over time. Both classes of loads are shown. While fixed load appliances (top) are used throughout, shiftable load appliances (bottom) are operated intermittently. The figure also illustrates temperature dependent load variation. Detailed explanations are in the caption below.

In order to incorporate TOU dependence of shiftable loads, the proposed consumer model incorporates an extension of the Cobb-Douglas preference function [12,13]. The consumer’s preference for scheduling such an appliance in any hourly period of the day is regarded as a set of mutually competing tradeoffs. The log function is treated as the consumer’s utility [14,15,16]. In this manner, the latter is a saturating nonlinearity whose derivative is the unit cost of electricity [15,17].

A novel machine learning-based approach to train the consumer utility model, is proposed in this research. Using historical time-series data pertaining to residential energy consumption, ambient temperature, and local TOU policy, the model parameters can be adapted to reproduce a real consumer’s energy usage patterns during all 24-hourly periods of the day.

The overall scheme consists of one or two steps of unsupervised learning, depending on the availability of load data, (i) partial disaggregation, and (ii) parameter estimation. The latter step, which constitutes the core aspect of this research, is able to provide reliable estimates of the model parameters. It requires separate shiftable and fixed load data, which can directly be obtained from multiple sensors that monitor each load category. A purely non-intrusive alternative is to apply partial disaggregation to load data from single point sensor to monitor the aggregate load over time. In this manner, no prior knowledge of the appliances of a residential unit is assumed [9,18,19,20].

A novel hybrid algorithm has been developed for partial disaggregation. The fixed load component is extracted from the total (i.e., aggregate) load by expressing it as the weighted sum of independent basis vectors by means of Non-negative Matrix Factorization (NMF), an unsupervised learning algorithm that is used in many applications [21,22,23]. NMF uses multiplicative updates of the basis vectors and their weights, that can be revised as new load data arrives. The shiftable load component is represented as a Gaussian Mixture Model (GMM), which is extracted from the aggregated load data using the Expectation Maximization (EM) algorithm. The use of this algorithm in GMMs is a well-established method in machine learning [24,25]. The EM algorithm is flexible enough to allow changes in the GMM parameters’ update rules when new load data is used to replace or appended to the extant dataset. This step will be referred to as the Hybrid Unsupervised Partial Disaggregation Algorithm (HUPDA). To the best of the authors’ knowledge, HUPDA is the first effort to successfully hybridize EM-GMM with NMF as a standalone unsupervised learning algorithm.

The main step of this research involves training the parameters of the consumer utility model using the partially disaggregated load data that is either obtained directly from multi-point sensor monitoring or by means of the previous partial disaggregation step. The proposed consumer model’s parameters are trained so that its estimated hourly energy usage time series closely match the real time series. The difference between the two sequences is the error vector; its squared Euclidean norm serves as the quadratic objective function for minimization. A quadratic regularization term is added to the objective to ensure that there were no dramatic variations in the model parameters across two consecutive hourly periods. Linear inequality constraints ensure that relatively small discrepancies lying within a certain threshold are not penalized. In this manner, the overall problem has been formulated as a constrained quadratic program [26,29], which will be referred to as Semi-Parametric Unsupervised Model Algorithm (SPUMA) for the remainder of this article.

It is shown that there exists a unique global minimum of the underlying optimization problem in SPUMA. The Lagrangian dual is derived from the constrained optimality conditions, whence the training is carried out in an iterative manner [28,29]. Dual decomposition allows the parallel implementation of SPUMA as well as real-time updates as the residential unit’s load usage patterns slowly shift with time.

The overall mathematical framework is identical to that used to train Support Vector Machines (SVMs) for nonlinear regression [30], where estimation errors and model parameters are the primal variables (instead of the SVM’s weights). However, to the best of the authors’ knowledge, the literature on consumer parameter estimation in terms of a SPUMA in this manner, is sparse (if not non-existent).

1.1. Related Research

A comprehensive review of all significant residential consumer models has been published in [31]. In this article, the various modeling schemes have been divided into three basic categories, (i) bottom-up methods, (ii) top-down methods, and, (iii) hybrid methods. Bottom-up schemes use a different model for individual appliances. The appliance model parameters are referred to as micro-variables (e.g., thermal properties, load profiles, etc.). Hence, they are parametric models. These models are then grouped together into the household. The rate in such schemes may be as high as one sample per second. On the other hand, top-down methods represent entire homes as wholistic models [31]. The parameters of these models are referred to as macro-variables (e.g., aggregate loads, occupancy patterns, user preferences, etc.). Therefore, these are semi-parametric models that do not require appliance-level load profiles. These methods only need sampling rates that can be as low as once per hour. Hybrid modeling schemes are intermediate between bottom-up and top-down schemes. This research entails such an approach. It should be noted that non-parametric models are not considered in [31].

Only a limited number of studies have used fully parametric models of residential consumers in electrical grid applications. Such an approach has been proposed in [32] for a HEMS application. In a manner that bears a strong resemblance to this research, appliances have been categorized into curtailable and shiftable loads. A parametric model has been suggested for another HEMS application [33]. The trained model can be applied to schedule HVAC operation based on TOU price and automatically adjust smaller loads according to occupant comfort. Energy management in a community of residential units with solar panels and/or energy storage has been investigated in [34]. Consumer models of residential units are used to simulate their interactive behavior within the wholesale market. Only peak load specifications of home appliances have been used as model parameters. This is because this research is directed at long term planning applications, stretched over multiple years with seasonal changes. A parametric model for smart homes has been proposed in [35], to optimize energy consumption based on TOU pricing.

A variety of semi-parametric models have been investigated recently that render themselves to meaningful interpretation in terms of game theory, mechanism design, and/or machine learning. Fuzzy inferencing is used in a HEMS model for DR in [36]. Q-learning is applied to train the parameters. In the same manner as in this article, two categories of appliances are considered, shiftable loads and non-shiftable loads. In the approach taken in [37], multiple semi-parametric models based on consumer profiles are divided into clusters, in order to plan separate pricing policies for consumer clusters. Another HEMS model for DR applications has been suggested in [38] to obtain optimal thermostat settings of residential heating/cooling devices that minimizes cost. Relevant appliances are modeled separately. Since quadratic curve fits are used in place of physical parameters, they have been treated as semi-parametric models here. The approach in [4], is for aggregator day-ahead scheduling, and market planning applications, where modeling is based on prior consumer preferences to identify consumption patterns of different classes of industrial consumers. Prospect theory has been used to in [40] for semi-parametric consumer modeling. The approach incorporates psychological factors (i.e., risk seeking, risk neutral, and risk averse) into consumer models during variable energy pricing. The semi-parametric model in [41] uses Cobb-Douglas function along with perfect substitution, perfect complement to incorporate residential consumer’s preference tradeoff between energy consumption and battery storage. Different battery discharging rates have been studied in this study.

Machine learning algorithms such as NMF, GMMs, hidden Markov models (HMMs), and decision trees are included under the class of semi-parametric models as they offer some degree of visibility. A modeling method that hybridizes GMMs with Dempster-Shaffer theory has been investigated in [42]. It can be used to monitor unusual activity of patients in homecare facilities. The approach in [43] derives GMM parameters from multiple consumer usage time sequence patterns in buildings. A GMM based disaggregation method has also been used in [44]. An approach that uses a combination of K-means clustering and NMF to estimate consumer profiles has been proposed in [45]. NMF methods are used in [18,46] for load disaggregation. A matrix decomposition method similar to NMF is examined in [19]. HMMs are used to model individual appliances in [5] which preserves consumer privacy by means of a non-cooperative game theoretic framework. HMMs for load disaggregation have been extensively explored (cf. [8,20]). Decision trees and random forests have been used in [47] to model smart homes, which are trained using gradient boosting.

DNNs have been used to model consumers in [48], which takes into account the impact of weather in occupant comfort while minimizing HVAC load. Although DNNs are black boxes, the research makes use of the Shapley value, a cooperative game theoretic concept to garner meaningful explanations of DNNs; for this reason, the research has been labeled as an example of semi-parametric modeling in this article. A non-parametric residential model using Long Short-Term Memory (LSTM), which is a DNN has been proposed in [49]. DNNs are also applied in [50] for disaggregation, which are trained using federated learning. The article also provides a review of federated learning. For other detailed reviews of residential consumer modeling methods, the reader is referred to [1,2].

Concepts from algorithmic game theory have seen widespread use in the energy grid. Nash equilibrium is applied in [5] for a scheduling application. Prospect theory has been taken up in [6,40]. Mechanism design is used in bilevel energy trading in [14,15]. The double auction approach in [16] applies fairness to determine the unit energy costs that is charged to residential consumers. Disaggregated household consumption is applied for lower-level online bidding with the retailer in [51]. The bilevel approach is used to establish unit energy costs in a manner similar to [16]. Double auction energy trading algorithms have also been proposed in [14,15,16,52].

The Cobb Douglas preference function is routinely used in various application, such as wireless networks [53], governmental decision making [54,55,56], agriculture [57], economic risk assessment in energy systems [58], mining industry [59], community heating systems [60], renewable energy [61] and resource allocation [62,63]. It has also been adopted in energy grid research [41,64]. The Cobb Douglas function is very effective to represent tradeoffs, such as that between energy and spectral efficiency [53], renewable and non-renewable energy [55], capital and labor costs [59], repair and replacement [60], solar and wind energy [61], energy consumption and residential batteries [41] and energy consumption community-level storage [64]. A genetic algorithm is used to obtain suitable coefficients for the Cobb Douglas function [65], whereas the research reported here applies LICQP – a mathematically rigorous convex optimization method.

The utility of a consumer is quantified in econometrics as a saturating function of the quantity of resource consumed, with an asymptotic upper bound. For mathematical convenience, the logarithm function is widely used for this purpose, in various applications such as sensor networks [66], energy harvesting [67,68], surgical robotics [69] and energy systems [14,15,16,52,70,71]. Residential units equipped with solar panels and community level energy storage are modeled in [70]. The logarithm of the aggregate energy load serves as the consumer utility function. The derivative of the utility function is the unit energy cost of energy in trading applications [14,15,16,52,71].

Non-intrusive load disaggregation is a widely studied research topic nowadays. Several reviews and survey articles on the state-of-the-art methods have also been published. Many such reviews focus on HEMS for various applications, such as occupancy estimation, load forecasting, DR [72], appliance detection using simultaneous voltage and current sensor data [73], and their OFF/ON statuses in senior residences [74]. Reviews on disaggregation methods in HEMS and microgrids for various applications such as DR, assisted living, and anomaly detection [75], for DNN methods with sampling rates lower than 50-60 Hz [76], for e-monitoring techniques to encourage consumer participation in demand side management [77], and for fault diagnosis of HVAC appliances in buildings [78] have been published. Future research directions have been outlined in [79].

Since load disaggregation involves separating a single temporal sequence into its multiple constituents. Sans ample supplementary information, this is a very difficult task. Research proposals attempt to circumvent this difficulty in a variety of ways. Approaches in [45,46,80,81,82] apply supervised learning. Prior appliance profiles of all $2^N-1$ combinations of $N$ appliances are used with GMMs in [80] and HMMs [83]. The EM algorithm in [84] uses appliance level time series data to train HMMs. Pretrained HMM are used in [7]. The HMM approaches in [85,86] require prior appliance load profiles. Probabilistic methods are applied in [87] to analyze variation margins in trained machine learning models. Double Fourier spectra of voltage and current time series input data are used in [88]. Spectral graph theory is used in [89] to disaggregate energy usage in only a set of prior known appliances, while the remainder is treated as unknown appliance loads or noise. DNNs are perhaps the most widely used models for disaggregation. Convolutional neural networks (CNNs) are used in [51,88,91,92]. Temporal processing is handled by incorporating LSTM and/or biLSTM blocks within the DNNs in [51,92], as well as by attention mechanism and transformer blocks in [91,93,94,95]. A relational network is used in [96], and a generative adversarial network (GAN), in [97]. A novel pyramidal structured DNN is proposed in [98]. Since DNNs depend on supervised learning, to maintain consumer privacy, current research methods use transfer learning [8,96,99], and federated learning [11,95,98,100,101]. Clustering algorithms are sometimes applied to group together similar appliances in [83,102], or multiple home units [11]. The disaggregation method in [103] exploits strongly correlated PV generation of buildings in a neighborhood. Many approaches detect and/or use switching events of appliances for disaggregation [94,96,104]. Synthetic training data have been used in [10,93,97]. A very limited amount of research do not require prior disaggregated data, pre-trained models, appliance characteristics, simultaneous voltage and current measurements, occupancy information, or correlated usage (cf. [5,9,18,42]). Since single-point monitoring is the only requirement, these methods can be adopted universally.

Problems formulated as linear and quadratic programs have been used in energy grid applications. A mixed integer linear programming problem in [10] is handled using an algorithm that hybridizes branch & bound with genetic search. Disaggregation methods involving quadratic programming have been investigated in [7,84,89]. An algorithm that hybridizes sequential quadratic programming (SQL) with a genetic algorithm has been explored in [105] for a microgrid protection application. SQL has also been considered in [106] for a market-oriented scheduling problem, and in [107] for use in an electric vehicle problem and in [108] for economic dispatch. Karush-Kuhn-Tucker (KKT) conditions for optimality have been explicitly applied in [109] for another economic dispatch problem. The performance of SQL has been compared with swarm intelligence and other metaheuristics in [106,107,108].

1.2. Research Contributions

The significant contributions of this research are along the following main directions, (i) a semi-parametric model of the daily electricity consumption within a residential unit, (ii) a data-driven parameter estimation algorithm that finds rigorous theoretical support, (iii) an unsupervised, hybrid learning algorithm for categorical disaggregation of loads; moreover, the above discussion and extensive references show that this article provides (iv) a concise, state-of-the-art review of recently published research in all relevant topics. These directions have been further elaborated, each dealt individually below.

• Semi-Parametric Model

The semi-parametric model aims to fill a crucial gap in current large-scale simulation-oriented research, which tend to represent energy consuming units with ad hoc models [4,14,15,16,36,37,38,40,52]. This is the case with the macro-variables used within semi-parametric energy consumer models that apply fuzzy logic or Dempster-Shaffer theory; moreover, their interpretabilities are domain-specific in nature. General purpose, game-theoretic models, which can be used in energy market research, are sparse in the existing literature. The methodology proposed herein attempts to address these limitations.

There are other desirable features in this article’s approach. It can be extended to include other influencing factors, such as the occupancy rates at different periods in a day, and the home consumer’s own environmental awareness. It is also upwardly scalable to buildings, offices, apartment complexes, neighborhoods, and even cities; if necessary to avoid computational overloads, the model can be simplified by considering larger time periods, reducing the load sampling rate, scalarizing the model’s vector parameters, or excluding cost and/or temperature data.

• Parameter Estimation Algorithm

SPUMA – the proposed estimation algorithm is of equal significance. Metaheuristics, that are often used in similar research work [65,105,106,107,108] are associated with unavoidable shortcomings. Due to their randomness, multiple runs are needed in order to determine statistically significant performance metrics. Unlike gradient ascent, they use populations of candidate model parameters that not only place unnecessary computation overheads, but also make them slower to re-adapt to shifting electricity usage patterns. In contrast, the proposed algorithm is wholly deterministic, with guaranteed performance metrics.

It will be shown in the next section that the algorithm, being a dual decomposition method, can be readily parallelizable. As a gradient ascent approach, its computational requirements are minimal in nature. Lastly, all aspects of SPUMA are well grounded in convex optimization theory. Despite its close resemblance to iterative training in SVMs, significant differences do exist in how regularization function as well as the inequality constraints have been formulated. These have been provided in the following section, with derivations postponed until the appendix.

• Partial Disaggregation

The article introduces HUPDA as a new algorithm to disaggregate the overall energy load into two broad categories. It represents one of the few algorithms that have successfully hybridized NMF and GMM. Unlike many other non-intrusive load monitoring (NILM) approaches [10,80,83,93,97], this algorithm does not need any prior knowledge of installed appliances, multiple sensors, switching event detection, or model pre-training. It requires uniformly sampled aggregate load data, which can be obtained directly from non-intrusive, single-point monitoring alone.

• Literature Survey

The previous subsection of this article provides a summary of current research in the all topics closely related to this research with appropriate references. It could be of benefit to the energy grid community as a brief review of recent research advances.

2. Materials and Methods

The symbol $\mathcal{T}$ represents the ordered sets of time instants covering an entire day. Likewise, $\mathcal{N}$ is the ordered set of samples, where each sample pertains to a full day. At time instant $t\in \mathcal{T}$ in sample day $n\in \mathcal{N}$, the aggregate, shiftable, and fixed loads are $x_t\left(n\right)$, $x_t^s\left(n\right)$, and $x_t^f\left(n\right)$. As they must be non-negative quantities, $x_t\left(n\right),x_t^s\left(n\right),x_t^f\left(n\right)\in \left[0,\infty \right)$, and,

$$\begin{array}{c}{x}_t\left(n\right)=x_t^s\left(n\right)+x_t^f\left(n\right)+{\nu }_t\left(n\right).\end{array}$$

(1)

The quantity ${\nu }_t\left(n\right)$ in the above expression is a zero mean, iid random variable representing noise, so that $\mathbb{E}\left[{\nu }_t\left(n\right)\right]=0$ and $\mathbb{E}\left[{\nu }_t\left(n\right){\nu }_{t^{\text{'}}}\left(n\text{'}\right)\right]=0$. Since HUPDA and SPUMA have in-built error components that subsume noise, the third term ${\nu }_t\left(n\right)$ in the above expression can be ignored.

For the sake of conciseness, the following convention is adopted hereafter. Suppose $\mathcal{I}$ is some ordered set of indices, such that there exists a scalar $y_i$ for each $i\in \mathcal{I}$. The vector whose $i$^th element is $y_i$ is, $\mathbf{y}={\left[y_1{y}_2\dots y_{\left|\mathcal{I}\right|}\right]}^{\mathrm{T}}$ which is of size $\left|\mathcal{I}\right|\times 1$. The same vector may be represented in a more convenient manner as shown,

$$\mathbf{y}={\left[y_i\right]}_{i\in \mathcal{I}}.$$

(2)

Next, suppose that there exist vectors $\mathbf{y}\left(j\right)$ for all $j\in \mathcal{J}$, of equal size, $L\times 1$. The vectors can be arranged as the columns of a block matrix $\mathbf{Y}=\left[\mathbf{y}\left(1\right)\mathbf{y}\left(2\right)\dots \mathbf{y}\left(\left|\mathcal{I}\right|\right)\right]$, which is of size $L\times \left|\mathcal{I}\right|$. Our convention allows $\mathbf{Y}$ to be expressed more concisely as,

$$\mathbf{Y}={\left[\mathbf{y}\left(j\right)\right]}_{j\in \mathcal{J}}.$$

(3)

Another way to express the same matrix $\mathbf{Y}$ whose rows are in $\mathcal{I}$ and columns in $\mathcal{J}$ is,

$$\mathbf{Y}={\left[y_i\left(j\right)\right]}_{j\in \mathcal{J}}.$$

(4)

This convention, allows us to define three vectors $\mathbf{x}\left(n\right),{\mathbf{x}}^s\left(n\right),{\mathbf{x}}^f\left(n\right)\in {\mathbb{R}}^{\left|\mathcal{T}\right|\times 1}$, that represent sequences of aggregate, shiftable, and fixed energy usage during any day (or sample) $n\in \mathcal{N}$. In the manner of Equation (2), the vectors are,

$$\left\{\begin{array}{c}x\left(n\right)\triangleq {\left[x_t\left(n\right)\right]}_{t\in \mathcal{T}}\\ {\mathbf{x}}^s\left(n\right)\triangleq {\left[x_t^s\left(n\right)\right]}_{t\in \mathcal{T}}\\ {\mathbf{x}}^f\left(n\right)\triangleq {\left[x_t^f\left(n\right)\right]}_{t\in \mathcal{T}}\end{array}\right..$$

(5)

The unit cost of electricity remains fixed throughout a period; based on the TOU, any change can only take place only at the first instant of any period. Furthermore, the temperature data is recorded only once per period. Therefore, for compatibility, the instants in $\mathcal{T}$ are grouped into contiguous, non-overlapping periods of equal duration. If $\mathcal{K}$ is the ordered set of periods, then $\left|\mathcal{K}\right|$ must divide $\left|\mathcal{T}\right|$. For example, if the load is sampled once every min, and each period is an hour, then for a day, $\left|\mathcal{T}\right|=\mathrm{1,440}$ mins, and $\left|\mathcal{K}\right|=24$ hours, so that $L=60$ mins. Therefore, for some $L\in \mathbb{N}$,

$$L={\displaystyle \frac{\left|\mathcal{T}\right|}{\left|\mathcal{K}\right|}}.$$

(6)

It can be readily seen that the $k$^th period ($k\in \mathcal{K}$) of the day must begin at time instant, $t_{L(k-1)+1}\in \mathcal{T}$, and end at another instant, $t_{Lk}\in \mathcal{T}$. Associated with the same period $k$ is a binary vector, ${\mathbf{e}}_k\in {\left\{\mathrm{0,1}\right\}}^{\left|\mathcal{K}\right|\times 1}$ and its $j$^th element ($j\in \mathcal{K}$) is $e_{k,j}$. This indicator vector is given as,

$${\mathbf{e}}_k={\left[e_{k,j}\right]}_{j\in \mathcal{K}},e_{k,j}=\left\{\begin{array}{c}1,j=k\\ 0,j\ne k\end{array}\right..$$

(7)

The quantity $1_L$ represents the $L\times 1$ vector of all ones (i.e., $1_L={\left[1 1\dots 1\right]}^{\mathrm{T}})$. Clearly, the tensor product ${\mathbf{e}}_k⨂1_L$ is a binary vector of size $\left|\mathcal{T}\right|\times 1$, which is the also the size of any load vector defined earlier in Equation (5). All but $L$ elements in this tensor product are zeros. The $L$ elements between the positions $L(k-1)+1$ and $Lk$ are ones. It follows that the inner product of ${\mathbf{e}}_k⨂1_L$ with a daily load vector is the corresponding net load used during period $k$. The net aggregate, shiftable, and fixed loads used during period $k\in \mathcal{K}$ of any given day $n$ can now be defined,

$$\left\{\begin{array}{c}{}^p{x}_k\left(n\right)\triangleq {\left({\mathbf{e}}_k⨂1_L\right)}^{\mathrm{T}}x\left(n\right)\\ {}^p{x}_k^s\left(n\right)\triangleq {\left({\mathbf{e}}_k⨂1_L\right)}^{\mathrm{T}}{\mathbf{x}}^s\left(n\right)\\ {}^p{x}_k^f\left(n\right)\triangleq {\left({\mathbf{e}}_k⨂1_L\right)}^{\mathrm{T}}{\mathbf{x}}^s\left(n\right)\end{array}\right..$$

(8)

The corresponding load vectors, ${}^p\mathbf{x}\left(n\right),{}^p{\mathbf{x}}^s\left(n\right),{}^p{\mathbf{x}}^f\left(n\right)\in {\mathbb{R}}^{\left|\mathcal{K}\right|\times 1}$ covering all periods in $\mathcal{K}$ for the same sample $n\in \mathcal{N}$, are given as,

$$\left\{\begin{array}{c}{}^p\mathbf{x}\left(n\right)\triangleq {\left[{}^p{x}_k\left(n\right)\right]}_{t\in \mathcal{T}}\\ {}^p{\mathbf{x}}^s\left(n\right)\triangleq {\left[{}^p{x}_k^s\left(n\right)\right]}_{t\in \mathcal{T}}\\ {}^p{\mathbf{x}}^f\left(n\right)\triangleq {\left[{}^p{x}_k^f\left(n\right)\right]}_{t\in \mathcal{T}}\end{array}\right..$$

(9)

Temperatures have been recorded periodically each day. The recorded temperature for sample $n$ during period $k$ is ${\theta }_k\left(n\right)$, and the corresponding temperature vector $\mathsf{\theta }\left(n\right)\in {\mathbb{R}}^{\left|\mathcal{K}\right|\times 1}$ is,

$$\mathsf{\theta }\left(n\right)\triangleq {\left[{\theta }_k\left(n\right)\right]}_{k\in \mathcal{K}}.$$

(10)

Based on TOU pricing, the unit cost of electricity for period $k$ is $c_k$, which does not change with day. The vector of daily unit costs, denoted $\mathbf{c}\in {\mathbb{R}}^{L\times 1}$ is given as,

$$\mathbf{c}\triangleq {\left[c_k\right]}_{k\in \mathcal{K}}.$$

(11)

The load vectors defined in Equation (5) can be arranged as the columns of the load profile matrices, $\mathbf{X},{\mathbf{X}}^s,{\mathbf{X}}^f\in {\mathbb{R}}^{\left|\mathcal{T}\right|\times \left|\mathcal{N}\right|}$ as shown below,

$$\left\{\begin{array}{c}X\triangleq {\left[\mathbf{x}\left(n\right)\right]}_{n\in \mathcal{N}}\\ {\mathbf{X}}^s\triangleq {\left[{\mathbf{x}}^s\left(n\right)\right]}_{n\in \mathcal{N}}\\ {\mathbf{X}}^f\triangleq {\left[{\mathbf{x}}^f\left(n\right)\right]}_{n\in \mathcal{N}}\end{array}\right..$$

(12)

In the same manner, three periodic load profile matrices, ${}^p\mathbf{X},{}^p{\mathbf{X}}^s,{}^p{\mathbf{X}}^f\in {\mathbb{R}}^{\left|\mathcal{K}\right|\times \left|\mathcal{N}\right|}$ can be constructed in this manner,

$$\left\{\begin{array}{c}{}^p\mathbf{X}\triangleq {\left[{}^p\mathbf{x}\left(n\right)\right]}_{n\in \mathcal{N}}\\ {}^p{\mathbf{X}}^s\triangleq {\left[{}^p{\mathbf{x}}^s\left(n\right)\right]}_{n\in \mathcal{N}}\\ {}^p{\mathbf{X}}^f\triangleq {\left[{}^p{x}_k^f\left(n\right)\right]}_{n\in \mathcal{N}}\end{array}\right..$$

(13)

2.1. Utility

The resources complement one another when a consumer’s requirements can be met by using both. For instance, a left and a right shoe are complementary resources so long as they have the same color. Suppose a person owns three left shoes (red, green, blue) but two right shoes (red, green). Effectively, it (To remove potential bias, the neuter form of the third person pronoun is used throughout this article.) has only two pairs of shoes to choose from. Conversely, two resources can supplement each other. For example, a person gains equally by consuming either 12 oz. Coke, 12 oz. Pepsi, or 7 oz. Coke and 5 oz. Pepsi. Unlike shoes, drinks are divisible resources. A person who prefers Pepsi over Coke, may gain equally by consuming either 12 oz. Coke, or 10 oz. Pepsi, or 6 oz. Coke and 5 oz. Pepsi. Resource divisibility allows the use of coefficients to quantify relative benefits derived from using resources; in the last example they may be in the 5:6 ratio.

Next, let us examine a more generic situation with two resources. Let $x_1,x_2>0$ be the amounts consumed by a certain consumer. The overall gain that it derives from their combined use is some non-decreasing function, say $y\left(x_1,x_2\right)$. When the two resources exhibit perfect complementarity, the consumer’s gain is $y\left(x_1,x_2\right)=\min \left(a_1{x}_1,a_2{x}_2\right)$, where $a_1,a_2>0$ are two coefficients intrinsic to the consumer. On the other hand, when the resources display perfect substitutability, then its net gain $y\left(x_1,x_2\right)=a_1{x}_1+a_2{x}_2$.

The Cobb-Douglas function is to represent the intermediate situation where the resources are neither complements nor substitutes of one another. The consumer’s net gain can be expressed as $y\left(x_1,x_2\right)=x_1^{a_1}{x}_2^{a_2}$. The schematic in Figure 2 shows the shape of a typical isoquant for the Cobb-Douglas function. For comparison, isoquants for perfect substitution and perfect complement, are also shown in Figure 2.

In general, with $\mathbf{x}$ representing the vector of all resources, the Cobb-Douglas preference, with multiple resources indexed $k=\mathrm{1,2},\dots$, consumed in amounts $x_k$ is,

$$y\left(\mathbf{x}\right)=\prod _k{\left(x_k\right)}^{a_k}.$$

(14)

The parameters $a_k>0$ are referred to as elasticity constants that sum up to one.

The consumer’s utility from using $\mathbf{x}$, is its logarithm, $\log y\left(\mathbf{x}\right)$. Using the logarithm offers a threefold advantage, for it makes the utility: (i) non-decreasing, (ii) strictly concave, and, (iii) Lipschitz continuous. Together, these advantageous features ensure that the shape of the utility function is in conformity with the law of diminishing returns. Based on the above rational, if ${}^p{x}_k^s\left(n\right)$ is the shiftable load used by a residential unit, its occupants’ utility $u\left({}^p{x}_k^s\left(n\right)\right)$ must be $a_k{\log }_{+}\left({}^p{x}_k^s\left(n\right)-b_k-d_k{\theta }_k\left(n\right)\right)$, where $a_k,d_k$ and $d_k$ are model parameters for the period $k$ that do not vary with day $n$. The function ${\mathrm{l}\mathrm{o}\mathrm{g}}_{+}\left(·\right)$ is equal to the natural logarithm of its enclosed argument, but lower bounded at zero. The parameter $b_k$ is the base component of the shiftable load for period $k$, which is governed by extraneous factors (e.g., it is not practical to operate a washer/dryer unit at 2:45 AM). The parameters $a_k$ and $d_k$ are two constants of proportionality such that $d_k{\theta }_k\left(n\right)$ is a temperature-dependent baseline. In general, the energy consumption increases when the ambient temperature deviates in either direction from an occupant-specific optimum ${\theta }^{\mathrm{*}}$ (that maximizes its comfort). Under these circumstances, the absolute difference between the real temperature ${\theta }_k\left(n\right)$ and ${\theta }^{\mathrm{*}}$ (i.e., $\left|{\theta }_k\left(n\right)-{\theta }^{\mathrm{*}}\right|$) should have been incorporated in the argument of the function ${\mathrm{l}\mathrm{o}\mathrm{g}}_{+}\left(·\right)$. However, as the specific data sets used in this study involve only warmer, summer days, it may be conveniently assumed that ${\theta }_k>{\theta }^{\mathrm{*}}$. Thus, the expression for the utility function is,

$$u\left({}^p{x}_k^s\left(n\right)\right)=a_k\log \left(x_k^s\left(n\right)-b_k-d_k{\theta }_k\left(n\right)\right).$$

(15)

Figure 3 illustrates how the utility function (top) and its derivative (bottom) vary with shiftable load. Note that the function’s derivative is also the cost per unit of energy. The figure also shown the fixed load.

2.2. Estimation

The derivative of the utility function in Equation (15) with respect to ${}^p{x}_k^s\left(n\right)$ is the marginal utility, which is also the unit cost $c_k$, so that $c_k={\displaystyle \frac{\partial }{\partial {}^p{x}_k^s\left(n\right)}}u\left({}^p{x}_k^s\left(n\right)\right)$. A few algebraic steps lead to the following model-based approximation of the shiftable load,

$${}^p{\widehat{x}}_k^s\left(n\right)=c_k^{-1}{a}_k+b_k+d_k{\theta }_k\left(n\right).$$

(16)

The above equation splits the shiftable load into cost dependent, constant, and temperature dependent terms. The ‘hat’ ($\widehat{·}$) over any load variable $x$ indicates that the resultant $\widehat{x}$ is an estimate of $x$. Thus, ${}^p{\widehat{x}}_k^s\left(n\right)$ is an estimate of the true load ${}^p{x}_k^s\left(n\right)$.

Fixed loads pertain to appliances such as lighting and refrigeration, which are typically operated continuously and independent of the unit energy cost. It is assumed that the net fixed load at any period $k$ consists of a constant term and a temperature-dependent term. Using $p_k$ and $q_k$ as two other model parameters, the estimated fixed load is given by,

$${}^p{\widehat{x}}_k^f\left(n\right)=p_k+q_k{\theta }_k\left(n\right).$$

(17)

2.3. Hybrid Unsupervised Partial Disaggregation Algorithm

The total shiftable load $x_t^s\left(n\right)$ at instant $t$ in day $n$ follows a GMM distribution. The probability that $x_t^s\left(n\right)=x$ is the sum of Gaussian terms as shown below,

$$\mathrm{p}\mathrm{r}\left[x_t^s\left(n\right)=x\right]=\sum _{j\in \mathcal{G}\cup \left\{0\right\}}{\pi }_{j}N\left[x|{\mu }_j\right].$$

(18)

where $N\left[x|{\mu }_j\right]$ is defined as per the Gaussian distribution,

$$N\left[x|{\mu }_j\right]={\displaystyle \frac{1}{\sigma \sqrt{2\pi }}}{e}^{{{\displaystyle \frac{1}{{\sigma }^2}}(x-{\mu }_j)}^2}.$$

(19)

The mean of the $j$^th Gaussian is ${\mu }_j$ ($j\in \mathcal{G}$), where $\mathcal{G}$ is the set of non-zero mean Gaussians in the GMM and ${\mu }_0=0$ when no shiftable appliance is in use. Figure 4 shows an example of such a GMM with six Gaussians (i.e., $\left|\mathcal{G}\right|=5$). All Gaussians share the same width parameter $\sigma$, a small, positive constant ($0<\sigma \ll 1$) that accounts for random noise. Only a single Gaussian, $j\in \mathcal{G}\cup \left\{0\right\}$ may be active at any instant $t$. The parameter ${\pi }_j$ is the prior probability of the $j$^th Gaussian.

The mean ${\mu }_j$ of a Gaussian $j$ can be equated to the load usage (in kW) of any shiftable appliance. However, as the total number of shiftable appliances present in the household is not known a priori. To address this lack of information, ${\mu }_j$ must also represent the overall load when combinations of two or more appliances are operated simultaneously; the cardinality $\left|\mathcal{G}\right|$ must be high enough to cover all cases. The caption placed immediately below Figure 4 provides further insights that are specific to the GMM example shown in the figure.

Based on NMF, the fixed load vector, ${\mathbf{x}}^f\left(n\right)$ of each sample day $n\in \mathcal{N}$, is represented through of a set $\mathcal{B}$ of basis profiles, ${\mathbf{w}}_b^f\in {\mathbb{R}}^{\left|\mathcal{T}\right|\times 1}$ ($b\in \mathcal{B}$), that are non-negative and orthonormal vectors (i.e., ${\mathbf{w}}_b^f\ge 0$, ${{\mathbf{w}}_{b^{\text{'}}}^f}^{\mathrm{T}}{\mathbf{w}}_b^f={\delta }_{b^{\text{'}}=b}$ where ${\delta }_{.}$ is the Kronecker delta– ${\delta }_{TRUE}=1,$ ${\delta }_{FALSE}=0$). Each basis vector ${\mathbf{w}}_b^f$ has its weight $h_b^f\left(n\right)\ge 0$, and the fixed load ${\mathbf{x}}^f\left(n\right)$, is their weighted sum,

$${\mathbf{x}}^f\left(n\right)=\sum _{b\in \mathcal{B}}{h}_b^f\left(n\right){\mathbf{w}}_b^f.$$

(20)

It may be observed that all sample days are represented through the same set of basis vectors, but weighted differently with day. This is because each basis reflects a specific pattern, whereas actual load profiles are not identical every day $n\in \mathcal{N}$. To further explain this rationale, let us consider a simple scenario where a residential unit’s fixed load consists of lighting equipment and a freezer. Since lighting is used only during evenings, HUPDA may train itself to incorporate this pattern as a basis vector in $\mathcal{B}$. On the other hand, while the freezer’s load directly varies with daily temperature. HUPDA may select one basis vector (“warm”) to model a sunny day in the dataset, and another basis (“cool”) for a day that was overcast at all times. Whereupon, a day that is sunny during the morning periods, and overcast thereafter, may be apply higher weights to the “warm” basis until noon, slowly transitioning to instants where the “cool” basis is weighed more.

The basis vectors make up the columns of a basis matrix, ${\mathbf{W}}^f\in {\mathbb{R}}^{\left|\mathcal{T}\right|\times \left|\mathcal{B}\right|}$,

$${\mathbf{W}}^f\triangleq {\left[{\mathbf{w}}_b^f\right]}_{b\in \mathcal{B}}.$$

(21)

A weight matrix ${\mathbf{H}}^f\in {\mathbb{R}}^{\left|\mathcal{B}\right|\times \left|\mathcal{N}\right|}$ can be constructed whose elements are the weights $h_b^f\left(n\right)$ ($b\in \mathcal{B}$) as follows,

$${\mathbf{H}}^f\triangleq {\left[h_b^f\left(n\right)\right]}_{b\in \mathcal{B},n\in \mathcal{N}}.$$

(22)

Figure 4 illustrates how the columns of the fixed load matrix ${\mathbf{X}}^f$ are captured by NMF with the matrices ${\mathbf{W}}^f$ and ${\mathbf{H}}^f$; as per Equation (20). The caption appearing beneath the figure elucidates this representation. An estimate of the fixed load profile matrix is computed in this manner, whence,

$${\widehat{\mathbf{X}}}^f={\mathbf{W}}^f{\mathbf{H}}^f.$$

(23)

The trainable parameters are the Gaussian priors ${\pi }_{j,t}$ ($j\in \mathcal{G}\cup \left\{0\right\},t\in \mathcal{T}$), and means ${\mu }_j$ ($j\in \mathcal{G}$), as well as the NMF’s basis and weight matrices, ${\mathbf{W}}^f$ and ${\mathbf{H}}^f$. These parameters are incremented using ${\eta }_{\mathrm{G}}\ll 1$ and ${\eta }_{\mathrm{N}}\ll 1$, which are the GMM and NMF learning rates respectively. It must be noted that all these quantities are meant for internal use within HUPDA, whose purpose is to deliver estimates of the load matrices, ${\widehat{\mathbf{X}}}^f$ and ${\widehat{\mathbf{X}}}^s$ from their aggregate $\mathbf{X}$.

Although the sequential order of time instants in $\mathcal{T}$ is preserved throughout the learning process, the samples in $\mathcal{N}$ are randomly rearranged at the onset of an iteration. An iteration hitherto will refer to a single loop through all sample days $n\in \mathcal{N}$.

The EM algorithm involves latent variables $z_{j,t}\left(n\right)$ which are posterior probabilities; more specifically, $z_{j,t}\left(n\right)=\mathrm{p}\mathrm{r}\left[\mathrm{G}\mathrm{a}\mathrm{u}\mathrm{s}\mathrm{s}\mathrm{i}\mathrm{a}\mathrm{n}j\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}|x_{j,t}\left(n\right)\right].$ During the E-step (estimation), the posteriors are approximated in two steps in the manner shown next,

$$Z_t\left(n\right)\leftarrow \sum _{j\in \mathcal{G}\cup \left\{0\right\}}{\pi }_{j}N\left[x_{j,t}\left(n\right)|{\mu }_j\right];$$

(24)

$$z_{j,t}\left(n\right)\leftarrow {\displaystyle \frac{{\pi }_{j}N\left[x_{j,t}\left(n\right)|{\mu }_j\right]}{Z_t\left(n\right)}}.$$

(25)

Equation (24) merely obtains the variables $Z_t\left(n\right)$, which is the denominator in the RHS of Equation (25) to normalize ${\mathrm{z}}_{j,t}\left(n\right)$ so that $\sum _j{z}_{j,t}\left(n\right)=1$.

The prior probabilities are incremented as shown below,

$${\pi }_j\leftarrow \left(1-{\eta }_{\mathrm{G}}\right){\pi }_j+{\displaystyle \frac{{\eta }_{\mathrm{G}}}{\left|\mathcal{N}\right|}}\sum _{n\in \mathcal{N}}{z}_{j,t}\left(n\right).$$

(26)

As the Gaussians are directly linked with the overall load when some combination of (zero, one, or more) shiftable appliances are in use at the same time, their prior probabilities at different time instants should be different. However, it can be seen that the prior ${\pi }_j$ does not have $t$ as subscript. This is because the increment in Equation (26) is applied in their proper order in $\mathcal{T}$, so that at any instant $t$, ${\pi }_j$ in the RHS of Equation (26) is the prior from the previous instant, $t-1$. As long as ${\pi }_j$ is strictly meant for internal use inside HUPDA, incrementing it in this manner ensures that there are no abrupt changes between two consecutive instants.

The Gaussian means are incremented in a similar manner,

$${\mu }_j\leftarrow \left(1-{\eta }_{\mathrm{G}}\right){\mu }_j+{\displaystyle \frac{{\eta }_{\mathrm{G}}}{\left|\mathcal{N}\right|}}\sum _{n\in \mathcal{N}}{\mathrm{z}}_{j,t}{x}_t\left(n\right).$$

(27)

The shiftable load can be estimated in the following manner,

$${\widehat{x}}_t^s\left(n\right)\leftarrow \underset{\mathcal{G}\cup \left\{0\right\}}{\max }{\mu }_j\le x_t\left(n\right).$$

(28)

The estimated shiftable load matrix is removed from the aggregate load matrix to provide the estimated fixed load matrix that is used to update the NMF matrices,

$${\widehat{\mathbf{X}}}^f\leftarrow \mathbf{X}-{\widehat{\mathbf{X}}}^s.$$

(29)

NMF attempts to decompose ${\tilde{\mathbf{X}}}^f$ into matrices ${\mathbf{W}}^f$ and ${\mathbf{H}}^f$, and there are several NMF algorithms that can accomplish this task. The NMF used in this research iteratively converges towards optimal values of ${\mathbf{W}}^f$ and ${\mathbf{H}}^f$ at the minimum of the Frobenius norm of the difference, ${\widehat{\mathbf{X}}}^f-$ ${\mathbf{W}}^f{\mathbf{H}}^f$, while preserving the orthonormality of ${\mathbf{W}}^f$, as well as the non-negativities of ${\mathbf{W}}^f$ and ${\mathbf{H}}^f$,

$${\mathbf{W}}^f,{\mathbf{H}}^f\equiv \underset{\begin{array}{c}W\ge 0,H\ge 0\\ {\mathbf{W}}^{\mathrm{T}}W={\mathbf{I}}_{\left|\mathcal{B}\right|}\end{array}}{\mathrm{arginf}}{‖{\widehat{\mathbf{X}}}^f-\mathbf{W}\mathbf{H}‖}_F^2.$$

(30)

Everywhere in this article, ${\mathbf{I}}_L$ is the $L\times L$ identity matrix. The basis and weight matrices from the previous iteration of HUPDA are incremented next. The weight matrix is incremented in the following manner,

$${\mathbf{H}}^f\leftarrow \left(1-{\eta }_N\right){\mathbf{H}}^f+{\eta }_N{\mathbf{H}}^f\circ {\left({{\mathbf{W}}^f}^{\mathrm{T}}{\mathbf{W}}^f{\mathbf{H}}^f\right)}^{\circ -1}\circ {{\mathbf{W}}^f}^{\mathrm{T}}{\widehat{\mathbf{X}}}^f.$$

(31)

In the above, $\mathbf{X}\circ \mathbf{Y}$ is the Hadamard product of $\mathbf{X}$ and $\mathbf{Y}$, and ${\mathbf{Z}}^{\circ -1}$ is their element-wise inversion of $\mathbf{Z}$; so that if the element in row $i$, column $j$ of $\mathbf{Z}$ is $z_{i,j}$, then the element in the same row and column of ${\mathbf{Z}}^{\circ -1}$ is $z_{i,j}^{-1}$. In this manner if $\mathbf{X}$, $\mathbf{Y}$, $\mathbf{Z},$ are matrices of the same size, then so would be ${\mathbf{Z}}^{\circ -1}$, $\mathbf{X}\circ \mathbf{Y}$ and any other Hadamard product formed from them.

To preserve orthonormality, the basis matrix is incremented in a somewhat different manner as follows,

$${\mathbf{W}}^f\leftarrow \left(1-{\eta }_N\right){\mathbf{W}}^f+{\eta }_N{\mathbf{W}}^f\circ {\left({\widehat{\mathbf{X}}}^f{\mathbf{H}}^f{{\widehat{\mathbf{X}}}^f}^{\mathrm{T}}{\mathbf{W}}^f\right)}^{\circ -1}\circ {\widehat{\mathbf{X}}}^f{{\mathbf{H}}^f}^{\mathrm{T}}.$$

(32)

2.4. Semi-Parametric Unsupervised Model Algorithm

Since unit costs and temperatures are recorded once per period, SPUMA uses a somewhat coarsely grained load data. Its inputs are the shiftable and fixed load matrices ${\widehat{\mathbf{X}}}^s,{\widehat{\mathbf{X}}}^f\in {\mathbb{R}}^{\left|\mathcal{T}\right|\times \left|\mathcal{N}\right|}$, which are transformed into the corresponding periodic matrices ${}^p{\widehat{\mathbf{X}}}^s,{}^p{\widehat{\mathbf{X}}}^f\in {\mathbb{R}}^{\left|\mathcal{K}\right|\times \left|\mathcal{N}\right|}$ as in Equation (8). Although all noteworthy aspects of the algorithm have been explained here, for the full derivations of some of the equations, the reader is referred to the Appendix.

As a convenient reminder, $1_L,0_L$ are two $L\times 1$ vectors of all ones and zeros, ${\mathbf{I}}_L$ is an $L\times L$ identity matrix, and ${\mathbf{O}}_L$ is an $L\times L$ matrix of all zero elements. All of the above quantities are used throughout this article, with subscripts indicating size.

2.4.1. Constraints

SPUMA uses two algorithmic constants ${ϵ}^s$ and ${ϵ}^f$, that establish tolerance bounds. Smaller errors that are within these bounds are not penalized; this is how errors are typically handled in SVM regression.

The quantities $a_k,$ $b_k$ and $d_k$ in Equation (16) are shiftable load model parameters at period $k\in \mathcal{K}$, whose optimal values are to be determined by SPUMA. In the same manner as in Equation (9) and Equation (10), let the vectors $\mathbf{a},\mathbf{b},\mathbf{d}\in {\mathbb{R}}^{\left|\mathcal{K}\right|\times 1}$ be defined as,

$$\mathbf{a}\triangleq {\left[a_k\right]}_{k\in \mathcal{K}};$$

(33)

$$\mathbf{b}\triangleq {\left[b_k\right]}_{k\in \mathcal{K}};$$

(34)

$$\mathbf{d}\triangleq {\left[d_k\right]}_{k\in \mathcal{K}}.$$

(35)

Given any sample $n\in \mathcal{N}$, let us consider the shiftable load consumed at $k\in \mathcal{K}$. The LHS of Equation (16) is the estimated load, ${}^p{\widehat{x}}_k^s\left(n\right)$. The RHS of Equation (16) is the model-based approximation ${\displaystyle \frac{a_k}{c_k}}+b_k+d_k{\theta }_k\left(n\right)$ of the same load. The difference between them should be within established bounds, i.e., $-{ϵ}^s\le {\displaystyle \frac{a_k}{c_k}}+b_k+d_k{\theta }_k\left(n\right)-{}^p{\widehat{x}}_k^s\left(n\right)\le +{ϵ}^s$. Any excess difference is an error that incurs a penalty. With the vector of such errors being represented as ${\mathsf{\xi }}^s\in {\mathbb{R}}^{\left|\mathcal{K}\right|\times 1}$, the inequality constraints for all periods in $\mathcal{K}$ of the sample day $n\in \mathcal{N}$, can be expressed as,

$$-1_{\left|\mathcal{K}\right|}{ϵ}^s-{\mathsf{\xi }}^s\le {\mathbf{c}}^{\circ -1}\circ \mathbf{a}+\mathbf{b}+\mathsf{\theta }\left(n\right)\circ \mathbf{d}-{}^p{\mathbf{x}}^s\left(n\right)\le 1_{\left|\mathcal{K}\right|}{ϵ}^s+{\mathsf{\xi }}^s.$$

(36)

Referring to Equation (16), it is evident that the fixed load must be bounded on both sides as, $-{ϵ}^f\le p_k+q_k{\theta }_k\left(n\right)-{}^p{\widehat{x}}_k^f\left(n\right)\le +{ϵ}^f$. Before introducing the analogous, vectorized constraints for the fixed load, let us define the parameter vectors $\mathbf{p},\mathbf{q}\in {\mathbb{R}}^{\left|\mathcal{K}\right|\times 1}$,

$$\mathbf{p}\triangleq {\left[p_k\right]}_{k\in \mathcal{K}};$$

(37)

$$\mathbf{q}\triangleq {\left[q_k\right]}_{k\in \mathcal{K}};$$

(38)

The inequality constraints for the fixed load vector can be expressed in terms of the error vector ${\mathsf{\xi }}^f\in {\mathbb{R}}^{\left|\mathcal{K}\right|\times 1}$. in the same manner as in Equation (36), as follows,

$$-1_{\left|\mathcal{K}\right|}{ϵ}^f-{\mathsf{\xi }}^f\le \mathbf{p}+\mathsf{\theta }\left(n\right)\circ \mathbf{q}-{}^p{\mathbf{x}}^f\left(n\right)\le 1_{\left|\mathcal{K}\right|}{ϵ}^f+{\mathsf{\xi }}^f.$$

(39)

2.4.2. Error & Regularization

An error within any period in $\mathcal{K}$ can only occur if there is a bound violation; the corresponding scalar element in ${\mathsf{\xi }}^s$ or ${\mathsf{\xi }}^f$ may be either positive or negative. To prevent cancellations, error functions in SVMs and elsewhere use normed error terms. Specifically, the weighted sum of the squared Euclidean error norms ${‖{\mathsf{\xi }}^s‖}_2^2$ and ${‖{\mathsf{\xi }}^f‖}_2^2$ are adopted in SPUMA. With ${\gamma }^s,{\gamma }^f>0$ as weights, the error function $E:{\mathbb{R}}^{2\left|\mathcal{K}\right|\times 1}\to \mathbb{R}$, is given as,

$$E\triangleq {\displaystyle \frac{1}{2}}\left({\gamma }^s{{\mathsf{\xi }}^s}^{\mathrm{T}}{\mathsf{\xi }}^s+{\gamma }^f{{\mathsf{\xi }}^f}^{\mathrm{T}}{\mathsf{\xi }}^f\right).$$

(40)

The parameters of the semi-parametric model are $\mathbf{a},\mathbf{b},\mathbf{d},\mathbf{p},\mathbf{q}\in {\mathbb{R}}^{\left|\mathcal{K}\right|\times 1}$. Since they are vectors of the same size, let $\mathbf{r}\in {\mathbb{R}}^{\left|\mathcal{K}\right|\times 1}$ be a representative parameter vector, which may be any of $\mathbf{a},\mathbf{b},\mathbf{d},\mathbf{p},\mathbf{q}$. At two consecutive periods $k,k+1\in \mathcal{K}$, the components $r_k$ and $r_{k+1}$, must not be significantly different. Putting it in another way, it is desired to keep their squared difference ${\left(r_k-r_{k+1}\right)}^2$ low. This can be accomplished by adding regularization terms for the model parameters. These terms will be formulated using the shift matrix $\mathbf{S}\in {\mathbb{R}}^{\left|\mathcal{K}\right|\times \left|\mathcal{K}\right|}$, whose columns have been defined earlier in Equation (7),

$$\mathbf{S}={\left[{\mathbf{e}}_{{〈k+1〉}_{\left|\mathcal{K}\right|}}^{\mathrm{T}}\right]}_{k\in \mathcal{K}}.$$

(41)

The subscript ${〈k+1〉}_{\left|\mathcal{K}\right|}$ in Equation (41) is a unit increment of index $k$, modulo $\left|\mathcal{K}\right|$. For $k<\left|\mathcal{K}\right|$, ${〈k+1〉}_{\left|\mathcal{K}\right|}=k+1$, but ${〈\left|\mathcal{K}\right|+1〉}_{\left|\mathcal{K}\right|}=1$. The product $\mathbf{S}\mathbf{r}={\left[r_2{r}_3\dots r_{\left|\mathcal{K}\right|}{r}_1\right]}^{\mathrm{T}}$, so that ${‖\mathbf{r}-\mathbf{S}\mathbf{r}‖}_2^2$ is the sum of terms of the form ${\left(r_k-r_{k+1}\right)}^2$. Regularization terms of this form are used in SPUMA. Let $\mathbf{V}\in {\mathbb{R}}^{\left|\mathcal{K}\right|\times \left|\mathcal{K}\right|}$ be the following matrix,

$$\mathbf{V}={\mathbf{S}}^{\mathrm{T}}\mathbf{S}-\left({\mathbf{S}}^{\mathrm{T}}+\mathbf{S}\right)+{\mathbf{I}}_{\left|\mathcal{K}\right|}.$$

(42)

A few algebraic steps will show that ${\mathbf{r}}^{\mathrm{T}}\mathbf{V}\mathbf{r}={‖\mathbf{r}-\mathbf{S}\mathbf{r}‖}_2^2$. The regularization function $R:{\mathbb{R}}^{5\left|\mathcal{K}\right|\times 1}\to \mathbb{R}$ is formulated as shown in the following expression,

$$R\triangleq {\displaystyle \frac{1}{2}}\left({\gamma }^a{\mathbf{a}}^{\mathrm{T}}\mathbf{V}\mathbf{a}+{\gamma }^b{\mathbf{b}}^{\mathrm{T}}\mathbf{V}\mathbf{b}+{\gamma }^d{\mathbf{d}}^{\mathrm{T}}\mathbf{V}\mathbf{d}+{\gamma }^p{\mathbf{p}}^{\mathrm{T}}\mathbf{V}\mathbf{p}+{\gamma }^q{\mathbf{q}}^{\mathrm{T}}\mathbf{V}\mathbf{q}\right).$$

(43)

2.4.3. Problem Formulation

The primal variables in the optimization framework are the error vectors ${\mathsf{\xi }}^s,{\mathsf{\xi }}^f$ and the model parameter vectors $\mathbf{a},\mathbf{b},\mathbf{d},\mathbf{p},\mathbf{q}$. For brevity, the error vectors are consolidated into the following single vector $\mathsf{\xi }\in {\mathbb{R}}^{2\left|\mathcal{K}\right|\times 1}$,

$$\mathsf{\xi }\triangleq \left[\begin{array}{c}{\mathsf{\xi }}^s\\ {\mathsf{\xi }}^f\end{array}\right].$$

(44)

From Equations (16) and (17), it can be seen that $\mathbf{a},\mathbf{b},\mathbf{p}$ are the coefficients of temperature independent terms, while $\mathbf{d},\mathbf{q}$ are those of temperature dependent terms, they are incorporated into two vectors, $\mathsf{\alpha }\in {\mathbb{R}}^{3\left|\mathcal{K}\right|\times 1}$ and $\mathsf{\beta }\in {\mathbb{R}}^{2\left|\mathcal{K}\right|\times 1}$, which are,

$$\mathsf{\alpha }\triangleq \left[\begin{array}{c}a\\ b\\ p\end{array}\right];$$

(45)

$$\mathsf{\beta }\triangleq \left[\begin{array}{c}d\\ q\end{array}\right].$$

(46)

Before presenting a concise framework of the same optimization problem, the matrices ${\mathbf{Q}}_1\in {\mathbb{R}}^{2\left|\mathcal{K}\right|\times 2\left|\mathcal{K}\right|},{\mathbf{Q}}_2\in {\mathbb{R}}^{3\left|\mathcal{K}\right|\times 3\left|\mathcal{K}\right|},{\mathbf{Q}}_3\in {\mathbb{R}}^{2\left|\mathcal{K}\right|\times 2\left|\mathcal{K}\right|}$, which will be incorporated into the objective function are as below,

$${\mathbf{Q}}_1=\left[\begin{array}{cc}{\gamma }^s& 0\\ 0& {\gamma }^f\end{array}\right]⨂{\mathbf{I}}_{\left|\mathcal{K}\right|};$$

(47)

$${\mathbf{Q}}_2=\left[\begin{array}{ccc}{\gamma }^a& 0& 0\\ 0& {\gamma }^b& 0\\ 0& 0& {\gamma }^p\end{array}\right]⨂\mathbf{V};$$

(48)

$${\mathbf{Q}}_3=\left[\begin{array}{cc}{\gamma }^d& 0\\ 0& {\gamma }^q\end{array}\right]⨂\mathbf{V}.$$

(49)

Using Equations (41) and (43) and the above definitions, the loss function can be expressed simply as ${\displaystyle \frac{1}{2}}\left({\mathsf{\xi }}^{\mathrm{T}}{\mathbf{Q}}_1\mathsf{\xi }+{\mathsf{\alpha }}^{\mathrm{T}}{\mathbf{Q}}_2\mathsf{\alpha }+{\mathsf{\beta }}^{\mathrm{T}}{\mathbf{Q}}_3\mathsf{\beta }\right)$.

A similar approach is used in SPUMA, its constraints expressed in terms of three matrices ${\mathbf{A}}_1\in {\mathbb{R}}^{4\left|\mathcal{K}\right|\times 2\left|\mathcal{K}\right|},{\mathbf{A}}_2\in {\mathbb{R}}^{4\left|\mathcal{K}\right|\times 3\left|\mathcal{K}\right|},{\mathbf{A}}_3\left(n\right)\in {\mathbb{R}}^{4\left|\mathcal{K}\right|\times 2\left|\mathcal{K}\right|}$, and a vector $\mathbf{t}\left(n\right)\in {\mathbb{R}}^{4\left|\mathcal{K}\right|\times 1}$ that are introduced next,

$${\mathbf{A}}_1=-1_2⨂{\mathbf{I}}_{2\left|\mathcal{K}\right|};$$

(50)

$${\mathbf{A}}_2=\left[\begin{array}{c}+1\\ -1\end{array}\right]⨂\left[\begin{array}{c}\mathrm{diag}{\mathbf{c}}^{\circ -1}\\ {\mathbf{O}}_{\left|\mathcal{K}\right|}\end{array}\begin{array}{c}{\mathbf{I}}_{\left|\mathcal{K}\right|}\\ {\mathbf{O}}_{\left|\mathcal{K}\right|}\end{array}\begin{array}{c}{\mathbf{O}}_{\left|\mathcal{K}\right|}\\ {\mathbf{I}}_{\left|\mathcal{K}\right|}\end{array}\right];$$

(51)

$${\mathbf{A}}_3\left(n\right)=\left[\begin{array}{c}+1\\ -1\end{array}\right]⨂{\mathbf{I}}_2⨂\mathrm{diag}\mathsf{\theta }\left(n\right);$$

(52)

$$\mathbf{t}\left(n\right)=\left[\begin{array}{c}-1\\ +1\end{array}\right]⨂\left[\begin{array}{c}{}^p{\mathbf{x}}^s\left(n\right)\\ {}^p{\mathbf{x}}^f\left(n\right)\end{array}\right]-1_2⨂\left[\begin{array}{c}{ϵ}^s\\ {ϵ}^f\end{array}\right]⨂1_{\left|\mathcal{K}\right|}.$$

(53)

Using the above quantities, the constraints in Equations (36), (39) can be expressed in a compact manner as ${\mathbf{A}}_1\mathsf{\xi }+{\mathbf{A}}_2\mathsf{\alpha }+{\mathbf{A}}_3\left(n\right)\mathsf{\beta }-\mathbf{t}\left(n\right)\le 0$.

The constrained optimization problem can be formulated with the loss function $f\left(\mathsf{\xi },\mathsf{\alpha },\mathsf{\beta }\right)$ as the mapping $f:{\mathbb{R}}^{7\left|\mathcal{K}\right|\times 1}\to \mathbb{R}$, and the linear constraint function $\mathbf{g}\left(\mathsf{\xi },\mathsf{\alpha },\mathsf{\beta }\right)$ as another mapping $\mathbf{g}:{\mathbb{R}}^{7\left|\mathcal{K}\right|\times 1}\to {\mathbb{R}}^{4\left|\mathcal{K}\right|\times 1}$. The constrained optimization problem in SPUMA is presented in a formal manner below.

Minimize (w.r.t. $\mathsf{\xi },\mathsf{\alpha },\mathsf{\beta }$):

$$f\left(\mathsf{\xi },\mathsf{\alpha },\mathsf{\beta }\right)={\displaystyle \frac{1}{2}}\left({\mathsf{\xi }}^{\mathrm{T}}{\mathbf{Q}}_1\mathsf{\xi }+{\mathsf{\alpha }}^{\mathrm{T}}{\mathbf{Q}}_2\mathsf{\alpha }+{\mathsf{\beta }}^{\mathrm{T}}{\mathbf{Q}}_3\mathsf{\beta }\right).$$

(54)

Subject to:

$$\mathbf{g}\left(\mathsf{\xi },\mathsf{\alpha },\mathsf{\beta }\right)={\mathbf{A}}_1\mathsf{\xi }+{\mathbf{A}}_2\mathsf{\alpha }+{\mathbf{A}}_3\left(n\right)\mathsf{\beta }-\mathbf{t}\left(n\right)\le 0_{4\left|\mathcal{K}\right|}.$$

(55)

Appendix A.1 and A.2 show detailed derivations of the above expressions for the loss and constraint functions, which are also Equation (A1) and Equation (A2).

Theorem 1.

In SPUMA: (i) The loss $f\left(\mathit{\xi },\mathit{\alpha },\mathit{\beta }\right)$, is strictly convex with respect to$\mathit{\xi },\mathit{\alpha },\mathit{\beta }$; (ii) The feasible region given as$\left\{\mathit{\xi },\mathit{\alpha },\mathit{\beta }\in {\mathbb{R}}^{7\left|\mathcal{K}\right|\times 1}|\mathit{g}\left(\mathit{\xi },\mathit{\alpha },\mathit{\beta }\right)\le 0_{4\left|\mathcal{K}\right|}\right\}$ is a convex polytope; (iii) There is a unique global minimum$\left({\mathit{\xi }}^{*},{\mathit{\alpha }}^{*},{\mathit{\beta }}^{*}\right)=\underset{\mathit{g}\left(\mathit{\xi },\mathit{\alpha },\mathit{\beta }\right)\le 0_{4\left|\mathcal{K}\right|}}{\mathit{arginf}}f\left(\mathit{\xi },\mathit{\alpha },\mathit{\beta }\right)$.

Proof. The matrix ${\mathbf{Q}}_1$ in Equation (47) is obtained as the tensor product of $\mathrm{diag}{\left[{\gamma }^s{\gamma }^f\right]}^{\mathrm{T}}$ with an identity matrix. Both factors in the product are symmetric matrices. Moreover, since the weights ${\gamma }^s,{\gamma }^f$, are strictly positive quantities, and ${\mathbf{I}}_{\left|\mathcal{K}\right|}$ is positive definite, their tensor product ${\mathbf{Q}}_1$ must be a symmetric, positive definite matrix,

$${\mathbf{Q}}_1^{\mathrm{T}}={\mathbf{Q}}_1\succ {\mathbf{O}}_{2\left|\mathcal{K}\right|}$$

The shift matrix $\mathbf{S}$ in Equation (42) is a permutation matrix that has one and only one element in each row and in each column that is one; its remaining elements are zeros.

Let $\mathbf{y}\in {\mathbb{R}}^{\left|\mathcal{K}\right|\times 1},\mathbf{y}\ne 0_{\left|\mathcal{K}\right|}$ be an arbitrary vector. The squared L₂ norm of the difference between the original and shifted vectors is,

$${‖\mathbf{S}\mathbf{y}-\mathbf{y}‖}_2^2={\left(\mathbf{S}\mathbf{y}-\mathbf{y}\right)}^{\mathrm{T}}\left(\mathbf{S}\mathbf{y}-\mathbf{y}\right).$$

Replacing ${\left(\mathbf{S}\mathbf{y}-\mathbf{y}\right)}^{\mathrm{T}}$ with ${\mathbf{y}}^{\mathrm{T}}{\mathbf{S}}^{\mathrm{T}}-{\mathbf{y}}^{\mathrm{T}}$ in the RHS and expanding it yields,

${‖\mathbf{S}\mathbf{y}-\mathbf{y}‖}_2^2={\mathbf{y}}^{\mathrm{T}}{\mathbf{S}}^{\mathrm{T}}\mathbf{S}\mathbf{y}-{\mathbf{y}}^{\mathrm{T}}{\mathbf{S}}^{\mathrm{T}}\mathbf{y}-{\mathbf{y}}^{\mathrm{T}}\mathbf{S}\mathbf{y}+{\mathbf{y}}^{\mathrm{T}}\mathbf{y}.$

Separating the factors ${\mathbf{y}}^{\mathrm{T}}$ and $\mathbf{y}$ in each term in the RHS,

${‖\mathbf{S}\mathbf{y}-\mathbf{y}‖}_2^2={\mathbf{y}}^{\mathrm{T}}\left({\mathbf{S}}^{\mathrm{T}}\mathbf{S}-\left({\mathbf{S}}^{\mathrm{T}}+\mathbf{S}\right)+{\mathbf{I}}_{\left|\mathcal{K}\right|}\right)\mathbf{y}.$

From Equation (41), it can be seen that the second factor (${\mathbf{S}}^{\mathrm{T}}\mathbf{S}-\left({\mathbf{S}}^{\mathrm{T}}+\mathbf{S}\right)+{\mathbf{I}}_{\left|\mathcal{K}\right|}$) in the RHS of the above equality is the matrix $\mathbf{V}$. Replacing it with the latter,

${‖\mathbf{S}\mathbf{y}-\mathbf{y}‖}_2^2={\mathbf{y}}^{\mathrm{T}}\mathbf{V}\mathbf{y}>0.$

Since the squared Euclidean norm in the LHS is necessarily positive, it is evident from the above equality that the RHS ${\mathbf{y}}^{\mathrm{T}}\mathbf{V}\mathbf{y}>0$ for any vector $\mathbf{y}\ne 0_{\left|\mathcal{K}\right|}$. In addition, it can be seen that since ${\mathbf{S}}^{\mathrm{T}}\mathbf{S},\left({\mathbf{S}}^{\mathrm{T}}+\mathbf{S}\right)$ are symmetric matrices, so is $\mathbf{V}$. This shows that $\mathbf{V}$ is a symmetric, positive definite matrix,

${\mathbf{V}}^{\mathrm{T}}=\mathbf{V}\succ {\mathbf{O}}_{\left|\mathcal{K}\right|}.$

The matrix ${\mathbf{Q}}_2$ defined in Equation (48) is the tensor product of the diagonal matrix $\mathrm{diag}{\left[{\gamma }^a{\gamma }^b{\gamma }^p\right]}^{\mathrm{T}}$, and $\mathbf{V}$, both of which are symmetric, positive definite. As the tensor product of two symmetric, positive definite matrices must also be symmetric, positive definite, it follows that,

${\mathbf{Q}}_2^{\mathrm{T}}={\mathbf{Q}}_2\succ {\mathbf{O}}_{3\left|\mathcal{K}\right|}.$

The same rationale can be used with Equation (49) to establish that,

${\mathbf{Q}}_3^{\mathrm{T}}={\mathbf{Q}}_3\succ {\mathbf{O}}_{2\left|\mathcal{K}\right|}.$

The Hessians of $f\left(\mathsf{\xi },\mathsf{\alpha },\mathsf{\beta }\right)$ shown in Equation (54) with respect to its arguments are,

${\nabla }_{\mathsf{\xi }}^{2}f={\mathbf{Q}}_1\succ {\mathbf{O}}_{2\left|\mathcal{K}\right|};$

${\nabla }_{\mathsf{\alpha }}^{2}f={\mathbf{Q}}_2\succ {\mathbf{O}}_{3\left|\mathcal{K}\right|};$

${\nabla }_{\mathsf{\beta }}^{2}f={\mathbf{Q}}_3\succ {\mathbf{O}}_{2\left|\mathcal{K}\right|}.$

As the loss function’s Hessian with respect to $\mathsf{\xi },\mathsf{\alpha }$ and $\mathsf{\beta }$ are symmetric, positive definite matrices, the function $f\left(\mathsf{\xi },\mathsf{\alpha },\mathsf{\beta }\right)$ is strictly convex. This establishes statement (i).

From Equation (55), the inequality constraint is,

$\mathbf{g}\left(\mathsf{\xi },\mathsf{\alpha },\mathsf{\beta }\right)={\mathbf{A}}_1\mathsf{\xi }+{\mathbf{A}}_2\mathsf{\alpha }+{\mathbf{A}}_3\left(n\right)\mathsf{\beta }-\mathbf{t}\left(n\right)\le 0_{4\left|\mathcal{K}\right|}.$Each row of the above expression is a hyperplane in ${\mathbb{R}}^{7\left|\mathcal{K}\right|}$. There are $4\left|\mathcal{K}\right|$ such hyperplanes. It follows that their intersection is a convex polytope (not necessarily bounded). This establishes statement (ii).

A strictly convex function with a convex feasible domain must have a unique global minimum, it follows that the constrained optimization problem in SPUMA also has a unique global minimum. This establishes statement (iii).

All three statements, (i), (ii) and (iii) of this theorem have been established above. This concludes the proof of Theorem 1. □

Theorem 1 shows that starting from an arbitrary initial point, iteratively updating the variables $\mathsf{\xi },\mathsf{\alpha }$ and $\mathsf{\beta }$ along the gradient direction would lead to steady convergence towards the optimum $\left({\mathsf{\xi }}^{*},{\mathsf{\alpha }}^{*},{\mathsf{\beta }}^{*}\right)$. Unlike the other methods discussed earlier, stochastic methods or metaheuristics are not needed. This is the underlying justification for using dual gradient ascent to obtain the optimal model parameters in SPUMA.

2.4.4. Dual Decomposition

For each sample $n\in \mathcal{N}$, there exists the corresponding inequality constraint of the form, ${\mathbf{A}}_1\mathsf{\xi }+{\mathbf{A}}_2\mathsf{\alpha }+{\mathbf{A}}_3\left(n\right)\mathsf{\beta }-\mathbf{t}\left(n\right)\le 0_{4\left|\mathcal{K}\right|}$. Let its dual variable be $\mathsf{\lambda }\left(n\right)\in {\mathbb{R}}^{4\left|\mathcal{K}\right|\times 1}$, so that $\mathsf{\lambda }\left(n\right)\ge 0_{4\left|\mathcal{K}\right|}$. Applying the KKT first order stationary conditions to the problem yields $\left({\mathsf{\xi }}^{\mathrm{*}},{\mathsf{\alpha }}^{\mathrm{*}},{\mathsf{\beta }}^{\mathrm{*}}\right)$, expressed in terms of the aggregate variables, $\mathsf{\nu }\in {\mathbb{R}}^{4\left|\mathcal{K}\right|\times 1},\mathsf{\mu }\in {\mathbb{R}}^{3\left|\mathcal{K}\right|\times 1},\mathsf{\zeta }\in \mathbb{R}$, which are given as,

$$\mathsf{\nu }=\sum _{n\in \mathcal{N}}\mathsf{\lambda }\left(n\right);$$

(56)

$$\mathsf{\mu }=\sum _{n\in \mathcal{N}}{\mathbf{A}}_3^{\mathrm{T}}\left(n\right)\mathsf{\lambda }\left(n\right);$$

(57)

$$\mathsf{\zeta }=\sum _{n\in \mathcal{N}}{\mathbf{t}\left(n\right)}^{\mathrm{T}}\mathsf{\lambda }\left(n\right).$$

(58)

These expressions appear in Appendix A.3 as Equations (A3-A5).

The stationary conditions are applied to yield the following optimal primal vectors,

$${\mathsf{\xi }}^{\mathrm{*}}=-{\mathbf{Q}}_1^{-1}{\mathbf{A}}_1^{\mathrm{T}}\mathsf{\nu }.$$

(59)

$${\mathsf{\alpha }}^{\mathrm{*}}=-{\mathbf{Q}}_2^{-1}{\mathbf{A}}_2^{\mathrm{T}}\mathsf{\nu };$$

(60)

$${\mathsf{\beta }}^{\mathrm{*}}=-{\mathbf{Q}}_3^{-1}\mathsf{\mu }.$$

(61)

These expressions are identical to Equations (A6-A8) in Appendix A.3.

The gradient update rule for the dual variables $\mathsf{\lambda }\left(n\right)$, whose derivation is again relegated to Appendix A.4, is,

$$\mathsf{\lambda }\left(n\right)\leftarrow \mathsf{\lambda }\left(n\right)+{\eta }_D\left({\mathbf{A}}_1\mathsf{\xi }+{\mathbf{A}}_2\mathsf{\alpha }+{\mathbf{A}}_3\left(n\right)\mathsf{\beta }-\mathbf{t}\left(n\right)\right).$$

(62)

This is Equation (A9) in Appendix A4, where $0<{\eta }_D<1$ is the gradient step size.

2.4.5. Dual Gradient Ascent

Equations (56)–(58) are used in each iteration of SPUMA to update values of $\mathit{\nu },\mathit{\mu },\zeta$ using the previous values of $\mathit{\lambda }\left(n\right)$, whereupon Equations (59)–(61) are applied to determine the new values of $\mathit{\xi },\mathit{\alpha },\mathit{\beta }$. The dual $\mathit{\lambda }\left(n\right)$ is initialized to zeroes and subject to iterative increments as per Equation (62). The algorithm terminates when the inequality constraint in Equation (55) is satisfied. In order to keep track of the termination condition, a Boolean variable ($DONE$) is maintained internally in the algorithm.

The steps of the algorithm are outlined in Figure 5. In the flowchart it is assumed that ${\mathit{Q}}_1,{\mathit{Q}}_2,{\mathit{Q}}_3$ in the loss function, have been precomputed in accordance with Equations (47)–(49), that ${\mathit{A}}_1,{\mathit{A}}_2,{\mathit{A}}_3\left(n\right)$ and $\mathit{t}\left(n\right)$ in the constraint bound have also been determined a priori through Equations (50)–(53); furthermore that they are directly accessible to the iterative parameter estimation algorithm.

The schematics in Figure 6 illustrate the steps of the proposed approach,

3. Results

The simulation results reported here used actual energy usage data pertaining to two residential customers, User-1 and User-2, which was obtained from the Pecan Street Inc., Dataport database [110] sampled at the rate of once per minute, and for a total of 61 days during the months of March and April of 2018. The database included individual appliance usage measurements at the circuit level. The following are the main appliances, (i) air condition, (ii) electric car, (iii) washer/dryer, (iv) dishwasher, (v) microwave, (vi) refrigerator, (vii) furnace, (viii) bedroom apps, (ix) light plugs, and (x) kitchen apps. The first five appliances are classified as shiftable loads, and the remainder as fixed loads. It may be noted that some of these appliances were not common to User-1 and User-2. User-1 did not use bedroom apps and furnace, while User-2 did not own an electric car. The shiftable and fixed loads were stored as $1440\times 61$ arrays, separately for User-1 and User-2. In this manner, two raw datasets were obtained, each consisting of $\left|\mathcal{N}\right|=61$ sample vectors, ${\mathbf{x}}^s\left(n\right),{\mathbf{x}}^f\left(n\right)\in {\mathbb{R}}^{\left|\mathcal{T}\right|\times 1}$, where $\left|\mathcal{T}\right|=1440$.

The hourly temperature data for same 61 days was also obtained from [110]. The hourly time of use rates were obtained from the Austin Energy website [111]. The website provides electricity rates in the Austin area for a time of use pilot project. The user data used in this study were part of DER that received pricing information from the utility and were expected to defer their loads during peak or high energy price hours.

Prior to applying SPUMA, the raw load datasets of User-1 and User-2 as well as the temperature data were further divided into two sets – one for Month-1 (March), and the other for Month-2 (April). There were $\left|\mathcal{N}\right|=22$ samples in Month-1, and $\left|\mathcal{N}\right|=21$ in Month-2. Since the unit cost time sequence is the same in all samples, it was not changed.

A total of four datasets were obtained in this manner, uniquely identifiable from the user label (i.e., User-1 or User-2) and the month label (i.e., Month-1 or Month-2). Each dataset contained the load matrices ${\mathbf{X}}^s,{\mathbf{X}}^f\in {\mathbb{R}}^{\left|\mathcal{T}\right|\times \left|\mathcal{N}\right|}$, hourly temperatures of each daily sample ${\left\{\mathsf{\theta }\left(n\right)\right\}}_{n\in \mathcal{N}},\mathsf{\theta }\left(n\right)\in {\mathbb{R}}^{\left|\mathcal{K}\right|}$ and hourly unit cost $\mathbf{c}\in {\mathbb{R}}^{\left|\mathcal{K}\right|}$, where $\left|\mathcal{K}\right|=24$. The periodic load matrices ${}^p{\mathbf{X}}^s,{}^p{\mathbf{X}}^s\in {\mathbb{R}}^{\left|\mathcal{K}\right|\times \left|\mathcal{N}\right|}$ were computed as provided in Equation (13).

All simulations reported in this section used the MATLAB^® software package and the constrained optimization toolbox (fmincon() built-in function). Table 1 shows the values of various algorithmic constants in Equations (36,40,43). The weights of the parameters $\mathbf{a},\mathbf{d},\mathbf{q}$ were assigned relatively small values in order to increase their sensitivities to variations in the daily unit costs and temperatures. The constraint bounds were the sum of a constant baseline and a suitably scaled load norm; ${‖\mathbf{X}‖}_F$ denotes the Frobenius norm of matrix $\mathbf{X}$ in Table 1.

The two plots in Figure 7 depicting the learning algorithm’s steady convergence towards the unique, global optima for each dataset. The log Lagrangian dual function and the log maximum constraint violation are plotted against the log iteration. In all four cases the algorithm terminated when the constraints were tightly bounded. The y-coordinates in both plots are raw values that have not been normalized w.r.t. dataset sample size or the algorithmic constants in Table 1.

Table 2 contains the outcome after termination of the algorithm. The error vectors ${\mathsf{\xi }}^s,{\mathsf{\xi }}^f$ are incorporated in the bound constraints in Equations (36) and (39). The two terms in the quadratic loss are the error function $E$ and the regularizing function $R$ in Equations (40) and (43). The total computer time taken to converge is also shown.

Figure 8 shows the values of the parameters $a_k,b_k,d_k$ in the shiftable load and $p_k,q_k$ of the fixed load in the semi-parametric consumer model for each period in $\mathcal{K}$. Each parameter is shown separately, as a plot of its numerical value vs. period $k$.

Figure 9 and Figure 10 show scatter plots. The x-coordinates of each point are the model parameters obtained using SPUMA with partially disaggregated real load matrices as inputs. The y-coordinates were obtained after HUPDA was used for partial disaggregation which were used by SPUMA to estimate the same parameters. The scatter plots in Figure 9 are model parameters $a_k,b_k$ and $d_k$ of the shiftable loads, whereas those in Figure 10 are $p_k$ and $q_k$ of the fixed loads.

Figure 11 and Figure 12 illustrate how the proposed semi-parametric model can be leveraged to generate the consumers’ electricity usage profiles on typical days. They shaded regions show the ranges of distributions of the raw loads in the datasets. The shiftable load consumption patterns in Figure 11 were obtained as ${}^p{\widehat{x}}_k^s=c_k^{-1}{a}_k+b_k+d_k{\overline{\theta }}_k$. Those of the fixed loads shown in Figure 12 were obtained as ${}^p{\widehat{x}}_k^f=p_k+q_k{\overline{\theta }}_k$. The four datasets are plotted separately in both figures.

4. Discussion

The steady convergence of the iterative dual decomposition algorithm can be readily observed from the plots in Figure 7. Conservative dual increments per iteration successfully avoided overshooting. Only 100 iterations were sufficient to ensure that the optimization algorithm had been attained while using each of the four datasets.

Table 1 provides a summary of the results of the optimization process in SPUMA. The raw loss components $E$ and $R$ in Table 1 quite high as is the sum of quadratic terms and pertain to 21 or 22 days. Additionally, the raw values are lower for User-1 than for User-2. The same pattern is seen after scaling with respect to the weights. While refraining from providing concrete reasons, it is possible that the energy consumption profile of User-1 is more routine than that of User-2. An alternative artifact may lie in the nature of the appliances installed in the home units of User-1 and User-2. It should be noted that the computing time taken by the built-in function (fmincon) in MATLAB^® in this table are when the dual variables were sequentially updated. Running SPUMA on a hardware that supports parallel processing would improve the speed.

Figure 8 shows how the model parameters, $a_k,$ $b_k$, $d_k$ for the shiftable load, and $p_k$, $q_k$ for the fixed load, estimated separately for each dataset and for each hour of the day. The effect of regularizing the loss function in SPUMA was to ‘flatten’ them, so that they did not rise or fall too sharply across hourly periods. It can be observed that the shiftable load parameters $a_k,$ $b_k$ for User-2 follow very similar hourly trends in spite of different scales (specifically the green lines). Since $a_k$ is sensitive to the unit cost ($c_k$) while $b_k$ is a baseline component for hour $k$, it is hypothesized that User-2 is less cost conscious compared to User-1. More experimentation may provide further insights. The weights assigned to the regularization terms may always be adjusted to produce the most desirable outcomes depending on the application.

The scatter plots in Figure 9 are those of the consumer model’s shiftable load parameters $a_k,$ $b_k$ and $d_k$. In all four datasets, it can be seen that the absence of separate sequences for each load category did not produce dissimilar parameter estimates, indicating the effectiveness of HUPDA to preprocess aggregate load prior to SPUMA.

It can be seen that the scatter plots in Figure 10, which are those of $p_k$, $q_k$, show very similar parameter correlations.

Figure 11 shows the high degree of fidelity of the proposed model in capturing a consumer’s intrinsic usage patterns for shiftable loads. Each dataset is plotted separately in the figure. Figure 12 further corroborates this claim.