End-to-End Demand Response Model Identification and Baseline Estimation with Deep Learning

Estimation of demand baselines is a fundamental requirement to systematically incorporate demand response into power system operation and market settlements: the system operator must disaggregate the effective change in the demand in response to grid needs from the demand baseline to monitor the demand resource, settle compensations based on market prices or tariffs [10].

Existing baseline estimation approaches can be grouped into top-down or bottom-up approaches. Top-down approaches assume demand response “buyers” (system operators or utilities) estimate baselines, such as use conventional demand forecast or use specialized price-responsive demand models [11], or ex-post baseline re-estimations based on realized system and weather conditions through methods such as clustering [6] and control groups [7]. Similar baseline estimation methods have been implemented in the industry. The Scheduled Load Reduction Program (SLRP) operated by Pacific Gas & Electric Company (PG&E) calculates the baseline by averaging the load demand of the selected time periods in the ten previous normal operating days [12]. Bottom-up approaches, on the other hand, require demand response participants, the “provider” of demand response, to report baselines a-prior to the dispatch period, accompanied with supporting tariff mechanisms that encourage participants to truthfully report baselines base on the best of their knowledge [9, 13].

The decision-making process of an incentive-based DR participant can be generalized as an optimization problem with corresponding objective functions and constraints. The goal of the optimization is to minimize the total cost of the consumer responding to DR incentives which include the DR reward and the disutility. Quadratic objective functions are often used for modeling disutilities based on the intuition that the marginal “sacrifice” increases as a consumer reduces more demand, such as turning off lighting loads or reducing the HVAC comfort level [14, 15, 16]. The approach of representing DR decision-making with a disutility function has been widely adopted in DR coordination studies to design incentives that properly navigate consumers’ behavior [17, 18].

Constraints in the demand response model reflect the limits that a consumer can change its electricity demand regardless of the disutility cost. This could include demand shifting needs to shift the electricity consumption from peak price periods to off-peak periods, including rescheduling laundries, temporarily lower water heater settings, or use behind-the-meter energy storage to purposefully arbitrage electricity price differences. Operation of behind-the-meter energy storage and electric vehicles, subject to more sophisticated state-of-charge evolution and efficiency constraints, can also be represented as linear constraints [19, 20].

We assume a demand response agent conducts a private optimization over time horizon $T$ to determine its responses $\{y_{t}\}_{t=1}^{T}$ to a time-of-use price $\{\lambda_{t}\}_{t=1}^{T}$ on top of its baseline demand $\{D_{t}\}_{t=1}^{T}$. The objective of the demand response agent is to minimize its total electricity bill and the personal discomfort as follows,

where $y_{t}$ is the optimization decision variable, and $\theta=\{\alpha,\underline{P},\overline{P},\underline{E},\overline{E}\}$ includes parameters in the agent model. We model the personal discomfort as a quadratic function with respect to the preferred consumption schedule (no demand response) [21] and $\{\alpha(D_{t})\}_{t=1}^{T}$ is a time-varying comfort deviation coefficient. We model $\alpha$ a function of $D_{t}$ indicating the user comfort deviation coefficient is dependent on the baseline demand. For example, reducing the same amount of demand may incur a higher discomfort to a user during low demand periods. (1b) models the response limit at each time period in which $\overline{P}$ is the maximum demand increase and $\underline{P}$ is the maximum decrease. (1c) is the demand response balance constraint capping the total amount of demand response over the considered time periods between $\underline{E}$ and $\overline{E}$. This constraint loosely models resources such as thermal balancing constraint or energy storage state-of-charge constraints [22]. We also list the dual variables associated with these constraints $\{\underline{\mu}_{t},\overline{\mu}_{t},\underline{\nu}_{t},\overline{\nu}_{t}\}$, these dual variables will later be utilized for designing our gradient-based training algorithm.

In our DR identification problem, we assume parameters in $\theta$ are unknown to the system operator. For example, the system operator does not know what are the discomfort coefficient for their DR participants (i.e., $\alpha$), nor the DR response capacity that a participant can provide (i.e., $\{\underline{P},\overline{P},\underline{E},\overline{E}\}$).

Baseline demand forecast is an important topic in power system planning and operation, and there have been many works proposed for both individual consumer’s and aggregated users’ load forecasting, such as linear models [23], deep learning methods [24, 25] and many more. Our framework is flexible to incorporate different baseline demand forecast methods as a plug-in module, as long as it supports the automatic differentiation [26].

In particular, we use a fully connected multi-layer perceptron (MLP) model for the baseline demand forecast in all our experiments. Inputs to the module include useful features such as weather data (eg. temperature, humidity, wind speed), calendar information (e.g., month, day, weekday/weekend), and effects of human activities. The baseline demand forecast model is represented with the following formulation,

where $\hat{D}=\{\hat{D_{t}}\}_{t=1}^{T}\in\mathbb{R}^{T}$ is the predicted baseline demand, $x\in\mathbb{R}^{m}$ is the demand forecast features and $f(\cdot)$ is the baseline demand prediction network with model parameters $\beta$.

The goal of the system operator is to jointly identify the DR agent model $\theta=\{\alpha,\underline{P},\overline{P},\underline{E},\overline{E}\}$ and baseline demand $\{D_{t}\}_{t=1}^{T}$. However, the challenge is that the system operator cannot directly observe agent’s response $\{y_{t}\}_{t=1}^{T}$ with respect to the incentive price signal $\{\lambda_{t}\}_{t=1}^{T}$. Instead, the system operator only observes agent’s net demand $\{D_{t}+y_{t}\}_{t=1}^{T}$, which is the summation of agent’s baseline demand and demand response. Existing works usually take a two-step approach. In particular, it first estimates the baseline demand using historical demand data, then subtracts baseline estimation from observed net demand for agent DR model identification. Yet, as we have discussed in the literature review, in this approach the system operator must have access to a sufficient amount of “vanilla” demand data to make sure the baseline estimation is fair and accurate. However, such data might not be available if the participant has been constantly responding to DR signals.

Different from the two-step approach, we propose an end-to-end differential learning approach to jointly minimize the baseline prediction and agent model identification errors with a set of historical net demand measurements consisting $N$ days or scenarios, and each scenario includes $T$ time steps. In this case, the system operator solves the following problem

in which the training / identification inputs include

• •

  $z_{i,t}$ is the net demand measurement over scenario $i$ and time period $t$;

• •

  $x_{i}$ is the baseline demand prediction features;

• •

  $\lambda_{i,t}$ is the demand response incentives over scenario $i$ and time period $t$;

• •

  $f(\cdot)$ is the MLP model used for baseline demand prediction;

• •

  $g(\cdot)$ is the agent demand response model.

and the problem has the following variables

• •

  $\theta$ is the DR agent model parameters;

• •

  $\beta$ is the MLP network weights used for the baseline demand model;

• •

  $\hat{D}_{i,t}$ is the baseline demand over scenario $i$ and time period $t$;

• •

  $y^{*}_{i,t}$ is the DR agent response decision over scenario $i$ and time period $t$.

The joint identification problem minimizes the mean squared loss between the predicted net demand $(\hat{D}_{i,t}+y^{*}_{i,t})$ and the observed net demand $z_{i,t}$ for $N$ training samples over $T$ horizon. Fig 1 visualizes the proposed end-to-end joint baseline demand estimation and agent model identification framework. Specifically, the top green colored block in Fig 1 is the baseline demand forecast module that maps the demand forecast features to baseline demand prediction; and the bottom green colored block is the demand response forecast, that takes DR signals and predict the DR agent response. The summation of these two intermediate predictions construct the net demand prediction $(\hat{D}_{i}+y^{*}_{i})$, which is compared against the true net demand data $z_{i}$ as the model loss function (3a).

We use the OptNet architecture [5] to represent the agent identification problem (1). The OptNet architecture leverages prior knowledge about the optimization and thereby, favouring sample efficiency. To obtain the demand response forecast, it solves the optimization problem in the forward loop, with the learnt parameters. Then it compares the difference between the optima of the optimization problem with the learnt parameters, and the optima of the true optimization problem, and use the difference as the loss function to update the unknown optimization parameters.

Our objective is to calculate $\frac{\partial y^{*}_{t}}{\partial\theta}$ and to backpropagate the training loss with respect to the agent model parameters $\theta$. To this end, we ignore the dependency of $\alpha$ over $D_{t}$ for now and treat it as a parameter instead of a function. In the next subsection on end-to-end training, we will discuss incorporate the dependency $\alpha$ over $D_{t}$ using derivative chain rules. We first converting the agent model into a set of equations using the Karush–Kuhn–Tucker (KKT) conditions. The KKT conditions of optimization problem (1) can be written as, for $t=1,...,T$,

where (5a) is the stationarity condition of the DR agent problem (1), (5b)–(5e) are the complementary slackness conditions associated with constraints (1b) and (1c). $y^{*},\underline{\mu}^{*},\overline{\mu}^{*},\underline{\nu}^{*},\overline{\nu}^{*}$ are the optimal primal and dual variables.

Next, we calculate the derivative of the optimal response $y_{t}^{*}$ with respect to model parameters $\theta=\{\alpha,\underline{P},\overline{P},\underline{E},\overline{E}\}$. We first take the total derivatives of the KKT conditions and summarize in compact matrix form, which is given in Eq (6), {strip}

where $\mathbbm{1}$ is an all-one vector and $\Gamma$ is a lower triangular matrix with appropriate dimension, and $\Lambda(\cdot)$ creates a diagonal matrix from a vector. The deviation process of Eq (6) is provided in the Appendix.

Our goal is to obtain the derivative of the optimal primal solution $y^{*}$ with respect to the optimization parameters $\theta=\{\alpha,\underline{P},\overline{P},\underline{E},\overline{E}\}$ (we can also compute the Jacobian of optimal dual solutions with respect to the optimization parameters using the same approach). We divide both sides of (6) with the parameter derivatives $d\alpha,d\underline{P},d\overline{P},d\underline{E},d\overline{E}$ respectively, and multiply both sides with $A^{-1}$:

Note that the derivative of a parameter with respect to itself (for example, $\frac{d\alpha}{d\alpha}$) is one while the derivative to different parameters (for example, $\frac{d\overline{P}}{d\alpha}$) are zeros. Therefore we obtain $\frac{\partial y^{*}}{\partial\alpha},\frac{\partial y^{*}}{\partial\overline{P}},\frac{\partial y^{*}}{\partial\underline{P}},\frac{\partial y^{*}}{\partial\overline{E}},\frac{\partial y^{*}}{\partial\underline{E}}$ each as the first element of (7a)–(7e).

Remark 1. An implicit assumption for the above derivation regarding the Jacobian terms $\frac{\partial y^{*}}{\partial\alpha},\frac{\partial y^{*}}{\partial\overline{P}},\frac{\partial y^{*}}{\partial\underline{P}},\frac{\partial y^{*}}{\partial\overline{E}},\frac{\partial y^{*}}{\partial\underline{E}}$ is that the optimal solution is continuously differentiable in a neighborhood of the optimization parameters. However, for a general convex optimization problem, the optimal solution sets are not necessarily continuously differentiable in the optimization parameters. Paper [27] provided sufficient conditions under which these Jacobian are well-defined. We consider the DR agent optimization model as a Quadratic Program (QP), which satisfies Assumptions 1-3 in [27], thus the Jacobians are well-defined. Beyond the quadratic form DR agent model, our framework can also be used in general convex DR agent models, where both the objective function and constraints are twice continuously differentiable in the optimization variable $y$, and continuously differentiable in the optimization parameters.

Remember that we consider the problem setting where the system operator only has access to historical demand response incentive signal $\lambda_{i}$, baseline demand feature $x_{i}$ and the net demand measurement $z_{i}$ (baseline demand + demand response). The system operator solves the joint learning problem as stated in (3), in order to jointly identify the baseline demand and agent DR model. To this end, we propose an end-to-end (E2E) training framework that combines both the gradient of the baseline demand forecast module and the demand response agent module. The model parameters are updated as follows during iteration $k$:

where $\eta_{1}$ and $\eta_{2}$ are learning rates for baseline and DR agent modules. $\frac{\partial L}{\partial D}$ and $\frac{\partial L}{\partial y^{*}}$ are gradients of the loss function with respect to the intermediate baseline demand and agent response predictions. $\frac{\partial D}{\partial\beta}$ is the gradient of $D$ with respect to the baseline demand forecast model weights. For the MLP network model used in this paper, we compute the gradient via backpropagation provided by PyTorch. The gradient update of the baseline module also includes $\frac{\partial y^{*}}{\partial\alpha}\frac{\partial\alpha}{\partial D}$ as we assume the DR agent’s discomfort coefficient $\alpha$ depends on the baseline demand. $\frac{\partial y^{*}}{\partial\alpha}$ is calculated in (7a), and $\frac{\partial\alpha}{\partial D}$ depends on the agent model. Gradient updates of the DR agent module is solely based on $\frac{\partial y^{*}}{\partial\theta}$ which we calculate according to (7). Notably, this training step has a separable structure that the demand forecast neural network module has no common parameters with the differentiable optimization module, hence their gradient will not impact each other.

Compared to classic two-step methods that first predict the baseline demand and then estimate the agent DR model and response amount, our end-to-end framework alleviates error accumulation since both modules are optimized simultaneously. Compared with a fully connected network overall features, our design embeds prior knowledge about the agent behavior (i.e., the utility optimization problem), which reduces the sample and computational complexity in magnitudes, and also providing better model explainability.

The training process is summarized as following for a given training data set $\{x_{i},\lambda_{i},z_{i}\}$

1. 1.

   Initialize unknown parameters for the baseline forecast module $\beta^{(1)}$ and the DR agent module $\theta^{(1)}$; set the iteration number to 1.

2. 2.

   At $k$th iteration

   1. (a)

      Forward the baseline demand prediction module (2) with feature input $x_{i}$ and parameters $\beta^{(k)}$, record the baseline prediction result $\hat{D}_{i}^{(k)}$.

   2. (b)

      Forward solve the DR problem (1) with price $\lambda_{i}$ and parameters $\theta^{(k)}$ using an optimization solver, record the primal results ${y_{i}^{*}}^{(k)}$ and the dual results ${\underline{\mu}_{i}^{*}}{(k)},{\overline{\mu}_{i}^{*}}{(k)},{\underline{\nu}_{i}^{*}}{(k)},{\overline{\nu}_{i}^{*}}{(k)}$.

   3. (c)

      Add up the two prediction as the net demand prediction $\hat{D}_{i}^{(k)}+{y_{i}^{*}}^{(k)}$, compare it against the observed net demand $z_{i}$ and obtain loss $L=||\hat{D}_{i}^{(k)}+{y_{i}^{*}}^{(k)}-z_{i}||_{2}^{2}$.

   4. (d)

      Update $\beta^{(k+1)}$ and $\theta^{(k+1)}$ with Eq (8a) and (8b). Proceed to the next iteration.

Remark 2. Warm start. In experiments, we observed that adding a warm start to the baseline demand forecast module using net measurements could improve the performance. In particular, the warm start is stated as

We note that the warm start trains the baseline demand prediction model using net demand measurements $z_{i}$ rather than historical baseline demand (which the system operator may not have access to). The optimized model parameter $\beta$ is then used as a warm starting point for the end-to-end training stage, in which $\beta$ is optimized together with DR agent model parameters $\theta$ following Eq (8a) and (8b).

Remark 3. Demand response dependency over baseline. We ends this section by discussing how the proposed framework handles the dependency of demand response over the baseline demand. In some demand response scenarios like load shedding and peak shaving, the demand response is dependent over the baseline demand. In other scenarios like behind-the-meter energy storage price response, the change of demand may not depend on the baseline demand. This difference is modeled as $\frac{\partial\alpha}{\partial D}$ that if the disutility parameter $\alpha$ is dependent/independent of the baseline demand estimation $D$. If there is no dependency, we can set this value to zero and then the baseline gradient update reduces to $\beta^{(k+1)}\leftarrow\beta^{(k)}+\eta_{1}\frac{\partial L}{\partial D}\frac{\partial D}{\partial\beta}$.

We generate DR data using real demand profiles as the baseline, with added response to price signals. For the baseline demand, we assume a demand response aggregator works directly with aggregated measurement which is viewed as one response agent. We use smart meter data from the Iowa Distribution Test Systems Dataset [28]. The system consists of 3 feeders, 240 nodes, and 1120 customers, all of which are equipped with smart meters. These smart meters measure hourly energy consumption (kWh) in the year 2017. The electricity consumption data is at the secondary distribution transformer level. In particular, we use the aggregated demand data from Feeder 1, Bus 1004 for the experiment. We also include the temperature data of Columbus from the Climate Data Online (CDO) portal in the National Climatic Data Center (NCDC) website [29] as baseline load forecast features. We use 200 days of data for training and 60 days for testing the demand response prediction results using the identified parameters. Table I provides the detailed statistics of the demand profile.

For demand response, we consider a demand-independent demand response agent model in this case study as follows

where $\alpha$ is the comfort deviation coefficient and $M$ represents the limit of the total demand change within the optimization period. In order to test the performance of different agent models, we randomly sample the ground truth $\alpha\sim\text{Uniform}[10,50]$ and $M\sim\text{Uniform}[1,10]$ at each test time. $\lambda$ is the time-of-use tariff, where we use the day ahead LMP data from NYISO [30]. We assume a day-ahead scheduled DR in which the optimization horizon is set to be $T=24$ (24 hours).

We use the above agent model (10) and the corresponding day ahead LMP [30] to generate synthetic demand response data, for a specific sampled combination ($\alpha$, $M$). The baseline demand data and the synthetic demand response data are combined to generate the net measurement data set. Fig 2 visualizes an example of the baseline demand and net demand for 1 day. Throughout our experiments, all MLP baseline demand forecast models are feed-forward neural networks consisting of three hidden layers with (200, 100, 100) hidden units per layer, and the demand response module uses OptNet architecture defined on optimization Eq. (10). We use Adam optimizer with a learning rate of $1e-3$ for the base load forecast module and $1e-1$ for the OptNet layer. For evaluations, we repeated for $N=10$ times with different choices of $(\alpha,M)$ to report the average performance and the standard deviations. Table II summarizes the mean absolute error (MAE) and standard deviations for the parameter identification results.

In addition to agent demand response model identification, we compare the baseline forecast performance over two settings: a-priori baseline prediction, in which we directly use the trained baseline forecast model to predict demand over a 24-hour horizon. A-priori baseline prediction is useful for resource scheduling and estimating feeder loads. The other setting is ex-post baseline estimation, in which we estimate the baseline by subtracting the demand response component (based on the identified response model parameters) from the net demand. The ex-post baseline estimation is useful to identify effective demand response and to settle payments.

Results in Table III show a relatively large a-priori baseline prediction error due to the challenge of performing day-ahead demand prediction with limited training data and features (note that we train the model over 200 days and test over the rest 60 days). But the ex-post baseline estimation results are substantially improved. As a base for comparison, the MAE and MAPE between net demand (baseline+response) and baseline are 0.44 kW / 4.77% for the aggregation case. Hence, our method provides a significant improvement over ex-post baseline estimation, where the baseline estimation error is reduced by about 90%. Finally, Fig 3 provides an example of the end-to-end training approach performance in demand response identification, a-priori baseline estimation and ex-post baseline estimation in 5 consecutive test days.

We conduct identification experiments using real world demand response experiment data from the low carbon London project [31], which included demand profile at 30-minute granularity during full year 2013 from 1,100 households receiving time-of-use (ToU) tariffs and 4,400 households receiving non-ToU flat rate (14.23 pence/kWh) for comparison, all located in the London area. The ToU tariff is set to 11.76 pence/kWh on default, 67.20 pence/kWh during high price periods, 3.99 pence/kWh during low price periods. High price periods were activated 69 times with a total duration of 394 hours, and low price periods were activated 93 times with a total duration of 830 hours. Customers are grouped into three categories based on income level: affluent, comfortable, and adversity, each category contains ToU and non-ToU customers. Fig. 4 shows an example of the average consumer demand from the ToU group and the non-ToU group.

We design the agent optimization model as

where $[x]^{+}=\max\{0,x\}$ is the positive value function, $D_{t}$ is the baseline demand, and $D_{\mathrm{min}}$ is the annual minimum demand that represents non-responsive demand, also called the “ghost demand”. Hence, we assume the disutility coefficient is inverse proportional to the responsive baseline demand and has different values for demand increments ($\alpha_{1}$) and demand reduction ($\alpha_{2}$). We do not include the daily total response constraints in this design because there is only at most one demand response event per day in the considered dataset.

The input features to the baseline demand forecast MLP model are 6 dimensional: temperature, humidity, demand of the previous day and calendar information including hour/day/month.

We compare our end-to-end model performance with a benchmark method to estimate $\alpha_{1}$ and $\alpha_{2}$ using control groups. The benchmark method follows a regression approach and features include non-ToU control group demand and temperature. The baseline prediction MAE of the benchmark method over the default tariff period is 0.009 kW. The MAPE is 5.8%, $\alpha_{1}$ and $\alpha_{2}$ in (11) is thus estimated using least square using the difference between the actual demand and the predicted baseline during the higher or low tariff period.

Table IV shows the estimation result of $\alpha_{1}$ and $\alpha_{2}$ as formulated in (11) for three aggregated demand groups. The result shows our method generates very similar results compared to the benchmark method, which has access to the historical baseline demand information using a control group. Note that we do not assume any prior knowledge about the baseline demand or use demand profile from the non-ToU control group in the end-to-end method. Both methods provide very similar identification results and all are consistent with the study result from the London Project [31]. E2E training generates higher $\alpha_{1}$ in the Comfortable and Adversity group, possibly due to different error biases in the baseline prediction compared to the benchmark method; note that the benchmark method still has a baseline prediction MAPE of 5.8%.

As observed in Section V-A, the agent DR model identification module demonstrates robustness against baseline estimation errors, i.e., even though the a-priori baseline load estimation contains large errors, the demand response identification, and ex-post baseline estimation are accurate. In this section, we conduct investigative studies to understand why the OptNet layer is robust to the baseline estimation error, and how robust the model is against different levels of prediction errors. To exam the effect of prediction errors, we consider a setting where the incentive signal $\lambda_{t}$ and noisy demand response measurements $\hat{y}^{*}$ are available. In particular, we inject a zero-mean Gaussian noise with different standard deviations to the ground-truth demand response optima $y^{*}$. By changing the standard deviations of noise, we aim to simulate the effect of varying prediction errors, i.e., a larger standard deviation corresponds to a prediction model with higher errors.

Table V shows the DR model identification results under different levels of Gaussian noise. Note that $\sigma=0$ is the setting of no Gaussian noise. In practice, this corresponds to the scenario where the true baseline demands of users are available. As the noise level (i.e., $\sigma$) increases, the parameter identification error increases. But in general, the results demonstrate that the implicit optimization layer is robust to erroneous input data and can recover the underlying optimization problem parameters even under high-magnitude noise.

In addition, we plot the residual error of the baseline demand forecast network trained with historical data in Section V-A. As shown in Fig 5, the residual error of the baseline demand forecast model roughly follows a normal distribution, where the mean is near zero (unbiased estimator), while the standard deviation equals $3.5$ (relatively high variance). The ablation result shows in Table V that, the agent model can still achieve accurate parameter estimation, as long as the residual error of baseline estimation is near zero. Robustness of the implicit differential layer is an attractive property since it guarantees the applicability of the proposed approach under incomplete information (i.e., no true baseline), erroneous and noisy prediction, and measurement. An interesting future direction would be to formally study the robustness of the OptNet for DR model identification.

In all the previous experiments, the demand response and baseline models are trained with data collected, assuming the agent models are time-variant. However, in practice, users’ preferences and flexibility for demand response may vary with respect to time. For instance, a user might be less willing to reduce the energy usage of his air conditioner on an extremely hot day than a mild day ($\alpha$ value is higher and $M$ value is smaller in the first case than the latter). This is challenging since the system operators don’t know when/how the agent DR model is changed.

We performed an experiment to mimic such effects whereby we generated data for training the demand response and baseline models from a mixture of two agent models. As shown in Figure 6, in the green area user model is $\alpha=20,M=3$ (more responsive) and in the red area, user model is $\alpha=50,M=0.5$ (less responsive). The identified result by our model is $\hat{\alpha}=28.57,\hat{M}=1.749$, which is close to the mean of the two models. It would be interesting to investigate different clustering or classification methods to first separate data from different user response modes and then run parameter identification for future works.