Online Learning and Distributed Control for Residential Demand Response

Gaussian process is a non-parametric supervised machine learning method [22] that has been widely used to model nonlinear system dynamics [27]. A formal definition of GP over a function $f(\bm{x})$ is that any finite number of function realizations $(f(\bm{x}_{1}),f(\bm{x}_{2}),f(\bm{x}_{3}),\cdots)$ are random variables and follow a joint Gaussian distribution, which is fully specified by the mean function $m(\bm{x})$ and the (kernel) covariance function $k(\bm{x},\bm{x}^{\prime})$. Consider learning an unknown function $y\!=\!f(\bm{x})$ based on a training dataset $\mathcal{D}$ of $n$ noisy observations, i.e. $\mathcal{D}:=\{(\bm{x}_{i},\hat{y}_{i})|\,i\!=\!1,\cdots,n\}$. It aims to infer the function value $f(\bm{x}_{*})$ for a new point $\bm{x}_{*}$. Denote $\hat{\bm{y}}_{\mathcal{D}}\!:=\!(\hat{y}_{i})_{i=1}^{n}$ and $\bm{f}_{\mathcal{D}}\!:=\!(f(\bm{x}_{i}))_{i=1}^{n}$. By the GP definition, $(\bm{f}_{\mathcal{D}},f(\bm{x}_{*}))$ are assumed to be random variables and follow a joint Gaussian distribution

where vector $\bm{k}_{*,\mathcal{D}}\!:=\![k(\bm{x}_{*},\bm{x}_{1}),\cdots,k(\bm{x}_{*},\bm{x}_{n})]^{\top}$ and $\bm{m}_{\mathcal{D}}\!:=\!(m(\bm{x}_{i}))_{i=1}^{n}$. $K_{\mathcal{D},\mathcal{D}}\!\in\!\mathbb{R}^{n\times n}$ is the covariance matrix, whose $ij$-component is $k(\bm{x}_{i},\bm{x}_{j})$. Conditioning on the given observations $\mathcal{D}$, it is known that the posterior distribution of $f(\bm{x}_{*})|(\bm{f}_{\mathcal{D}}\!=\!\hat{\bm{y}}_{\mathcal{D}})$ is also Gaussian, i.e. $\mathcal{N}(\mu_{*|\mathcal{D}},\sigma^{2}_{*|\mathcal{D}})$, with the closed form

Then the mean value $\mu_{*|\mathcal{D}}$ (1a) can be used as the prediction on $f(\bm{x}_{*})$, and the variance $\sigma^{2}_{*|\mathcal{D}}$ (1b) provides a confidence estimate for this prediction. The merits of GP include 1) GP is a non-parametric method that avoids the bias of model selection; 2) GP works well with small datasets; 3) GP can incorporate prior domain knowledge by defining priors on hyperparameters or using a particular covariance function. The major issue of GP is the computational complexity, which scales cubically in the number of observations, i.e., ${O}(n^{3})$.

Thompson sampling [24] is a prominent Bayesian learning framework that was originally developed to solve the multi-armed bandit (MAB) problem [25] and can be extended to tackle other online learning problems. Consider a classical $T$-periods MAB problem where an agent selects an action (called “arm”) $a_{t}$ from the action set $\mathcal{A}$ at each time $t\in\{1,2,\cdots,T\}$. After taking action $a_{t}$, the agent observes an outcome $z_{t}$ that is randomly generated from a conditional distribution $\mathcal{P}_{\theta}(\cdot|a_{t})$, and then obtains a reward $r_{t}=r(z_{t})$ with known reward function $r(\cdot)$. The agent is initially uncertain about the parameter $\theta$ in $\mathcal{P}_{\theta}$ but aims to maximize the total expected reward using the observation feedback. To achieve good performance, it is generally required to take actions with an effective balance between 1) exploring poorly-understood actions to gather new information that may improve future reward and 2) exploiting what is known for decision to maximize the immediate reward.

TS is a straightforward online learning algorithm that strikes an effective balance between exploration and exploitation. As shown in Algorithm 1, TS treats the unknown $\theta$ as a random variable and represents the initial belief on $\theta$ with a prior distribution $\mathcal{P}$. At each time $t$, TS draws a random sample $\hat{\theta}$ from $\mathcal{P}$ (Step 3), then takes the optimal action that maximizes the expected reward based on the sample $\hat{\theta}$ (Step 4). After outcome $z_{t}$ is observed, the Bayesian rule is applied to update the distribution $\mathcal{P}$ over $\theta$ and obtain the posterior (Step 5).

The main features of the TS algorithm are listed below:

• •

  As outcomes accumulate, the predefined prior distribution will be washed out and the posterior converges to the true distribution or concentrates on the true value of $\theta$.

• •

  The TS algorithm encourages exploration by the random sampling (Step 3). As the posterior distribution gradually concentrates, less exploration and more exploitation will be performed, which leads to an effective balance.

• •

  The key advantage of the TS algorithm is that a complex online problem is decomposed into a Bayesian learning task (Step 5) and an offline optimization task (Step 4) [26], while the optimization remains the original formulation without being complicated by the learning task.

The practical implementation of AC control for residential DR is generally achieved through smart plugs or smart AC units with built-in controllers and sensors. These smart devices are able to measure, store, and communicate the temperature and AC power data in real time. Hence, the thermal dynamics model (4) can be estimated based on fine-grained historical measurement data. To this end, we provide the following two thermal model estimation schemes.

Since the customer opt-out status $z_{i,t}$ is binary, logistic regression [23] is used to model the transition probability function $p_{i}(\bm{w}_{i,t})$ in (3b). Because the output of logistic regression is naturally a probability value within $[0,1]$, and it is easy to implement and interpret. Moreover, logistic regression is compatible with the online Bayesian learning framework with efficient posterior update approaches (see Section IV-C for details). Accordingly, we formulate the transition probability function $p_{i}(\bm{w}_{i,t})$ as the logistic model (11):

where $\bm{\beta}_{i}$ is the weight vector describing how customer $i$ reacts to the environmental factors $\bm{w}_{i,t}$, and $\alpha_{i}$ depicts the individual preference. Define $\hat{\bm{w}}_{i,t}:=(1,\bm{w}_{i,t})$ and $\bm{\theta}_{i}:=(\alpha_{i},\bm{\beta}_{i})$. Then the linear term in (11) becomes $-\hat{\bm{w}}_{i,t}^{\top}\bm{\theta}_{i}$. Without causing any confusion, we use $p_{\bm{\theta}_{i}}(\bm{w}_{i,t})$ and $p_{i}(\bm{w}_{i,t})$ interchangeably.

As a consequence, the unknown information of customer $i$’s behaviors is summarized as vector $\bm{\theta}_{i}$, which can be estimated from the observations of $(\hat{\bm{w}}_{i,t},z_{i,t})$ in DR events. In contrast to the thermal dynamics model learning, the observation data of $(\hat{\bm{w}}_{i,t},z_{i,t})$ are not historically available but can only be obtained along with the real implementation of DR events. This leads to an online customer behavior learning and AC power control problem. Thus we employ the TS framework to develop the online AC control algorithm in Section IV to effectively balance exploration and exploitation.

In a typical DR setting, once a customer opts out, the AC unit will automatically be switched back to the defaulted operation mode with the original customer-set temperature $s_{i}^{\mathrm{set}}$. Taking the opt-out status into account, the actual AC power consumption $\hat{u}_{i,t}$ can be formulated as

which equals $u_{i,t}$ if customer $i$ stays in ($z_{i,t-\!1}\!=\!1$) or $u_{i,t}^{\mathrm{set}}$ if opts out ($z_{i,t-\!1}\!=\!0$). Denote $\bm{u}_{i}\!:=\!(u_{i,t})_{t\in[T]}$, $\bm{s}_{i}\!:=\!(s_{i,t})_{t\in[T]}$, $\bm{d}_{i}\!:=\!(d_{i,t})_{t\in[T]}$, $\bm{r}_{i}\!:=\!(r_{i,t})_{t\in[T]}$, and $\bm{z}:=(\bm{z}_{i})_{i\in[N]}$.

To simplify the expression, we reformulate the AC control constraints (2), the opt-out status transition (3), (8), (11), and the dynamics of environmental factors (5), (6), (9) or (10), for customer $i\in[N]$, as the following compact form (13):

where $\mathcal{X}_{i}$ denotes the corresponding feasible set. Then two system-level optimal control (SOC) models, i.e. (14) and (15), with different DR goals are established for the SA to solve optimal AC power control schemes over $N$ customers.

1) SOC-1 model (14) aims to reduce as much AC load as possible in a DR event, which can be used to flatten the load peaks or mitigate reserve deficiency emergency.

where objective (14a) minimizes the expected total AC energy consumption over the DR event, plus an opt-out penalty term in the last time step. $\rho_{i}$ is the penalty coefficient that can be tuned to balance load reduction and opt-out outcomes. $\mathbb{E}_{\bm{z}}[\cdot]$ denotes the expectation that is taken over the randomness of the customer opt-out status $\bm{z}$. Constraint (14b) collects the counterparts of (13) for all customers.

2) SOC-2 model (15) aims to closely track a regulated power trajectory $(L_{t})_{t\in[T]}$, which is determined by the upper-level power dispatch or the DR market.

Objective (15a) minimizes the expected total squared power tracking deviation from the global target $L_{t}$, plus the same opt-out penalty term defined in (14a). Constraint (15b) is the same as (14b).

The two SOC models (14) and (15) are indeed discrete-time finite-horizon control planning problems, which are in the form of nonconvex stochastic optimization, and the stochasticity results from the probabilistic opt-out status transition (3). We develop the distributed solution methods for the SOC models (14) and (15) in the next section.

Since the opt-out status transition of one customer is assumed to be independent of other customers in (3), objective (14a) in the SOC-1 model has no substantial coupling among different customers. Hence, the SOC-1 model (14) can be equivalently decomposed into $N$ local problems, i.e., model (16) for each customer $i\!\in\![N]$.

The sum of objectives (16a) over all $N$ customers is essentially objective (14a) in the SOC-1 model, and constraint (16b) is the individual version of (14b) for customer $i$.

The local model (16) is a stochastic optimization with the expectation over $\bm{z}_{i}$ in the objective. Since $\bm{z}_{i}$ follows the transition (3) with the probability function $p_{\bm{\theta}_{i}}(\bm{w}_{i,t})$ (11), we can derive the analytic form of the expectation in (16a), which leads to expression (17):

Expression (17) only differs from the expectation in (16a) by a constant term $\Delta t\sum_{t=1}^{T}u_{i,t}^{\mathrm{set}}+\rho_{i}$, thus they are equivalent in optimization. See Appendix A for the detailed derivation.

When $\bm{\theta}_{i}$ is given, we can obtain the optimal AC control schemes $\bm{u}_{i}^{*}$ for each customer $i$ via solving the local model (16), which is implemented in a fully distributed manner. Meanwhile, the aggregation of $\bm{u}_{i}^{*}$ over all $N$ customers is essentially an optimal solution to the SOC-1 model (14).

The objective (15a) in the SOC-2 model has coupling among different customers due to tracking a global power trajectory $(L_{t})_{t\in[T]}$. To solve this problem distributedly, we introduce a local tracking trajectory $\bm{l}_{i}\!:=\!(l_{i,t})_{t\in[T]}$ for each customer $i\in[N]$ with $\sum_{i=1}^{N}l_{i,t}=L_{t}$ for all $t\in[T]$. Then we substitute $L_{t}$ by $\sum_{i=1}^{N}l_{i,t}$ in objective (15a), so that (15a) can be approximated³³3The approximation is made by dropping the term $2\sum_{t=1}^{T}\sum_{i\neq j}\mathbb{E}(\hat{u}_{i,t}-l_{i,t})\mathbb{E}(\hat{u}_{j,t}-l_{j,t})$ from the expansion of (15a). This term is expected to be relatively small and thus neglectable. by the decomposable form (18a), which takes the form of a sum over $N$ customers. As a result, the SOC-2 model (15) is modified as (18):

where $l_{i,t}$ is also a decision variable.

Consequently, the only substantial coupling among customers in the modified SOC-2 model (18) is the equality constraint (18c). Therefore, we can introduce the dual variable $\bm{\lambda}\!:=\!(\lambda_{t})_{t\in[T]}$ for the equality constraint (18c) and employ the dual gradient algorithm [31] to solve the modified SOC-2 model (18) in a distributed manner. The specific distributed solution method is presented as Algorithm 2. The implementation of Algorithm 2 needs the two-way communication of $\bm{\lambda}$ and $\bm{l}_{i}$ between the SA and every customer in each iteration. Due to the simple structure with only one equality coupling constraint (18c), Algorithm 2 can converge quickly with appropriate step sizes, which is verified by our simulations.

Similar to (17), we can derive the equivalent analytic form (21) for the expectation term in (19a):

By substituting the expectations terms in (16a) and (19a) with their analytic forms (17) and (21) respectively, the local AC control models (16) and (19) become deterministic nonconvex optimization problems. Given parameter $\bm{\theta}_{i}$, they can be solved efficiently via available nonlinear optimizer tools, such as the IPOPT solver [32]. For concise expression, we denote the above distributed solution methods together with the optimizer tools as an oracle

which generates optimal $\bm{u}_{i}^{*}$ with the input of parameter $\bm{\theta}_{i}$.

Based on the TS framework, we develop the distributed online DR control algorithm as Algorithm 3 to learn customer opt-out behaviors and optimally control the AC power in a DR event. Since this online algorithm is implemented distributedly, we present Algorithm 3 from the perspective of an individual customer $i\in[N]$. The practical implementation of Algorithm 3 is illustrated as Figure 2.

In Algorithm 3, the unknown customer behavior parameter $\bm{\theta}_{i}$ is treated as a random variable, and we construct a Gaussian prior distribution $\mathcal{N}(\bm{\mu}_{i},\bm{\Sigma}_{i})$ for it based on historical information. At each time $t\in[T]$ of a DR event, $\hat{\bm{\theta}}_{i}$ is randomly sampled from the distribution for decision-making. Two key techniques used in Algorithm 3 are explained as follows.

1) To utilize the latest information and take future time-slots into account, we employ the model predictive control (MPC) method in the optimization and action steps. Specifically, at each time $t$, it solves the SOC model (14) or (15) for the rest of the DR event to obtain the optimal AC control scheme $\bm{u}_{i}^{*}$, but only implements the first control action $\bm{u}_{i}^{*}$(1). In addition, the latest predictions or estimations of the environmental factors, the updated thermal dynamics model, and recalculated baseline AC power $u_{i,t}^{\mathrm{set}}$ can be adopted in the optimization step.

2) After observing the outcome pair $(z_{i,t},\hat{\bm{w}}_{i,t})$, the variational Bayesian inference approach introduced in [33] is applied to obtain the posterior distribution on $\bm{\theta}_{i}$ with the update rules (23)-(25). It is well known that Bayesian inference for the logistic regression model (11) is an intrinsically hard problem [34], and the exact posterior $\mathcal{P}(\bm{\theta}_{i}|z_{i,t},\hat{\bm{w}}_{i,t})$ is intractable to compute. Thus we use the variational approach [33] for efficient Bayesian inference, which provides an accurate Gaussian approximation to the exact posterior with a closed form (24). The scalar $\xi_{i}$ in (23)-(25) is an intermediate parameter that affects the approximation accuracy, thus we alternate three times between the posterior update (24) and the $\xi_{i}$ update (25), which leads to an optimal approximated posterior. See reference [33] for details.

For online learning problems, the notion of “regret” and its variants are standard metrics that are defined to measure the performance of online learning and decision algorithms [35]. Accordingly, we denote $\bm{\theta}^{\star}:=(\bm{\theta}^{\star}_{i})_{i\in[N]}$ as the underlying true customer behavior parameter that the LSEs do not know but aim to learn. Then the regret of the proposed online DR control algorithm at $m$-th DR event is defined as

where $C^{\bm{\theta}^{*}}_{m}(\cdot)$ denotes the objective function (i.e., (14a) in the SOC-1 model or (18a) in the modified SOC-2 model) under the true value $\bm{\theta}^{\star}$. $\bm{u}^{\mathrm{online}}_{m}$ denotes the AC control scheme generated by the proposed online algorithm, while $\bm{u}^{\star}_{m}$ is the optimal AC control scheme that minimizes the objective $C^{\bm{\theta}^{*}}_{m}(\cdot)$. Thus $\mathrm{regret}(m)$ in (26) is always non-negative and measures the performance distance between the proposed online algorithm and the underlying best control scheme.

To evaluate the overall learning performance, we further define the cumulative regret until $M$-th DR event as

which is simply the sum of $\mathrm{regret}(m)$ over the first $M$ DR events. Generally, a sublinear cumulative regret over $M$ is desired, i.e., $\mathrm{cumu\_regret}(M)/M\to 0$ as $M\to+\infty$, since it indicates that $\mathrm{regret}(m)\to 0$ as $m\to+\infty$. In other word, it means that the proposed online algorithm can eventually learn the customer behaviors well and make optimal AC control schemes as more and more DR events are experienced. We demonstrate the regret results of the proposed algorithm via numerical simulations in the next section.

In this part, we compare the thermal dynamics prediction performance of the linear model (9) and the GP model (10). Real customer metering data, including indoor temperature and AC power, from ThinkEco⁴⁴4ThinkEco Inc. is a New York City-based energy efficiency and demand response solution company (http://www.thinkecoinc.com/). are used for model training and testing. The outdoor temperature data are procured from an open-access meteorological database⁵⁵5Iowa Environmental Mesonet [online]: https://mesonet.agron.iastate.edu/request/download.phtml?network=MA_ASOS.. Specifically, we use the time series of data in consecutive 5 days with the resolution of 15 minutes to fit the two thermal dynamics models (9) and (10). The GPML toolbox [36] is applied to implement the GP model and optimize the hyperparameters. Then, the fitted models are tested on the time series of real data in the next 3 days for indoor temperature prediction.

The prediction results of one time step ahead (15 minutes) and three time steps ahead (45 minutes) are presented in Figure 3. The average prediction errors of the indoor temperature are $0.632^{\circ}F$ and $0.737^{\circ}F$ for the GP model in these two cases, and $0.846^{\circ}F$ and $1.16^{\circ}F$ for the linear model. It is observed that both the GP model (10) and the linear model (9) work well in the thermal dynamics modelling, and the GP model achieves better prediction accuracy.

We then perform the proposed online DR control algorithm on the SOC-2 model, where the simulation configurations are the same as Section V-B. Since the SOC-2 problem involves the coordination among all customers, we set the customer number as $N=500$, and the global tracking target $L_{t}$ is randomly generated from $\text{Unif}(450,550)$(kW).