Input Modeling and Uncertainty Quantification for Improving Volatile Residential Load Forecasting

TL;DR

This paper introduces a two-stage approach incorporating input modeling and uncertainty quantification to enhance residential load forecasting accuracy and reliability, addressing input errors often neglected by existing methods.

Contribution

It proposes a novel two-stage framework that models input errors and propagates their impact into load forecasts, improving accuracy and interval coverage.

Findings

Outperforms existing methods in predictive accuracy

Provides better interval coverage for load forecasts

Uses global sensitivity analysis for efficiency

Abstract

Load forecasting has long been recognized as an important building block for all utility operational planning efforts. Over the recent years, it has become ever more challenging to make accurate forecasts due to the proliferation of distributed energy resources, despite the abundance of existing load forecasting methods. In this paper, we identify one drawback suffered by most load forecasting methods: neglect to thoroughly address the impact of input errors on load forecasts. As a potential solution, we propose to incorporate input modeling and uncertainty quantification to improve load forecasting performance via a two-stage approach. The proposed two-stage approach has the following merits. (1) It provides input modeling and quantifies the impact of input errors, rather than neglecting or mitigating the impact, a prevalent practice of existing methods. (2) It propagates the impact of input errors into the ultimate point and interval predictions for the target customer's load to improve predictive performance. (3) A variance-based global sensitivity analysis method is further proposed for input-space dimensionality reduction in both stages to enhance the computational efficiency. Numerical experiments show that the proposed two-stage approach outperforms competing load forecasting methods in terms of both point predictive accuracy and coverage ability of the predictive intervals.

Figures

Tables

Table 1. TABLE I: Estimated total Sobol’ indices S ^ T j 2 subscript superscript ^ 𝑆 2 subscript 𝑇 𝑗 \widehat{S}^{2}_{T_{j}} ’s for the inputs θ 1 , 2 , θ 1 , 3 , … , θ 1 , 8 subscript 𝜃 1 2 subscript 𝜃 1 3 … subscript 𝜃 1 8 \theta_{1,2},\theta_{1,3},\ldots,\theta_{1,8} in the second stage of NNGP-NIGP. Predictions are performed for customer 1 on 3 days in the IEEE 8-bus test case.

${\hat{S}}_{T_j}^2$	$\theta$	Day 1	Day 2	Day 3
${\hat{S}}_{T_2}^2$	${\theta }_{1,2}$	$0.834$	$0.800$	$0.897$
${\hat{S}}_{T_3}^2$	${\theta }_{1,3}$	$0.015$	$0.021$	$0.002$
${\hat{S}}_{T_4}^2$	${\theta }_{1,4}$	$2.568\times {10}^{-10}$	$0.084$	$1.415\times {10}^{-9}$
${\hat{S}}_{T_5}^2$	${\theta }_{1,5}$	$3.435\times {10}^{-11}$	$1.348\times {10}^{-9}$	$1.558\times {10}^{-10}$
${\hat{S}}_{T_6}^2$	${\theta }_{1,6}$	$2.316\times {10}^{-16}$	$9.519\times {10}^{-4}$	$9.984\times {10}^{-5}$
${\hat{S}}_{T_7}^2$	${\theta }_{1,7}$	$7.247\times {10}^{-11}$	$1.460\times {10}^{-4}$	0.004
${\hat{S}}_{T_8}^2$	${\theta }_{1,8}$	$5.692\times {10}^{-11}$	$2.690\times {10}^{-6}$	$3.468\times {10}^{-10}$

Table 2. TABLE II: The average computational times (in seconds) of NNGP-NIGP with and without GSA for one day’s prediction in all test cases.

Scale	NNGP-NIGP without GSA	NNGP-NIGP with GSA
8-bus	$53.694$s	$2.469$s $\sim 4.228$s
14-bus	$75.387$s	$2.247$s $\sim 12.375$s
24-bus	$136.736$s	$3.104$s $\sim 43.583$s
123-bus	$241.381$s	$3.617$s $\sim 176.736$s

Equations

$${\mathbf{x}}_t={\left({\theta }_{i,1}^t,\dots ,{\theta }_{i,i-1}^t,{\theta }_{i,i+1}^t,{\theta }_{i,i+2}^t,\dots ,{\theta }_{i,N}^t\right)}^{\top },$$

$$0=-P_i+\sum _{k=1}^N|V_i||V_k|(G_{ik}\cos {\theta }_{ik}+B_{ik}\sin {\theta }_{ik}),$$

$$z_i^1(\mathbf{x})$$

$$\text{with}{a}_j^1(\mathbf{x})$$

$$K^1(\mathbf{x},{\mathbf{x}}^{\prime })=\mathbb{E}[z_i^1(\mathbf{x})z_i^1({\mathbf{x}}^{\prime })]={\sigma }_b^2+{\sigma }_w^2\mathbb{E}[a_i^1(\mathbf{x})a_i^1({\mathbf{x}}^{\prime })].$$

$${\mu }_1=(\mathbb{E}[{\theta }_{i,1}],\dots ,\mathbb{E}[{\theta }_{i,i-1}],\mathbb{E}[{\theta }_{i,i+1}],\dots ,\mathbb{E}[{\theta }_{i,N}]{)}^{\top }$$

$${\mathsf{V}}_1=(\mathbb{V}[{\theta }_{i,1}],\dots ,\mathbb{V}[{\theta }_{i,i-1}],\mathbb{V}[{\theta }_{i,i+1}],\dots ,\mathbb{V}[{\theta }_{i,N}]{)}^{\top }$$

$$y=\tilde{y}+{\epsilon }_y,$$

$$\mathbf{x}=\tilde{\mathbf{x}}+{\mathbf{\epsilon }}_x,$$

$$y=f(\mathbf{x}-{\mathbf{\epsilon }}_x)+{\epsilon }_y,$$

$$f(\mathbf{x}-{\mathbf{\epsilon }}_x)=f(\mathbf{x})-{\mathbf{\epsilon }}_x^{\top }\frac{\partial f(\mathbf{x})}{\partial \mathbf{x}}+\dots ,$$

$$y\approx f(\mathbf{x})-{\mathbf{\epsilon }}_x^{\top }\partial \bar{\mathbf{f}}+{\epsilon }_y.$$

$$\displaystyle\mathbf{C}(\mathbf{x}_{i},\mathbf{x}_{j})=\sigma_{f}^{2}\exp\big{(}-\frac{1}{2}(\mathbf{x}_{i}-\mathbf{x}_{j})^{\top}\Lambda^{-1}(\mathbf{x}_{i}-\mathbf{x}_{j})\big{)},$$

$$\mathbb{E}[f({\mathbf{x}}_{\ast })|\mathbf{X},Y]$$

$$\mathbb{V}[f({\mathbf{x}}_{\ast })|\mathbf{X},Y]$$

$$\mathbb{E}[f({\mathbf{x}}_{\ast })]={\left([\mathbf{C}+{\sigma }_y^2{\mathbf{I}}_N+\mathrm{missing}diag(\mathbf{\Delta }\bar{\mathbf{f}}{\Sigma }_x{\mathbf{\Delta }\bar{\mathbf{f}}}^{\top }){]}^{-1}Y\right)}^{\top }\mathbf{q}.$$

$$\mathbf{q}=\int \mathbf{C}({\tilde{\mathbf{x}}}_{\ast },\mathbf{X})p({\tilde{\mathbf{x}}}_{\ast }|{\mathbf{x}}_{\ast },{\Sigma }_x)d{\tilde{\mathbf{x}}}_{\ast },$$

$$q_i$$

$$\mathbb{V}[f({\mathbf{x}}_{\ast })]$$

$${\mathbf{Q}}_{ij}$$

$$S_i=\frac{{\mathbb{V}}_{X_i}({\mathbb{E}}_{X_{\sim i}}[Y|X_i])}{\mathbb{V}(Y)},$$

$$S_{T_i}=\frac{{\mathbb{E}}_{X_{\sim i}}[{\mathbb{V}}_{X_i}(Y|X_{\sim i})]}{\mathbb{V}(Y)}=1-\frac{{\mathbb{V}}_{X_{\sim i}}({\mathbb{E}}_{X_i}[Y|X_{\sim i}])}{\mathbb{V}(Y)},$$

$$X=\sum _{i=1}^p{\beta }_i{\gamma }_i,$$

$${\hat{S}}_{T_i}=\frac{{(2n)}^{-1}\sum _{j=1}^n{\left(\hat{f}(\mathbf{A}{)}_j-\hat{f}({\mathbf{A}}_{\mathbf{B}}^{(i)}{)}_j\right)}^2}{\mathbb{V}(\hat{f}(\mathbf{C}))},$$

$$H_{\ell }=\sum _{i=24(\ell -1)+1}^{24\ell }{\hat{S}}_{T_i}^1,\ell =1,2,\dots ,N-1.$$

$$\text{MAPE}=T^{-1}\sum _{h=1}^T\left|\frac{P_h-{\hat{P}}_h}{P_h}\right|,$$

$$\text{CP}=n_c/24,$$

$$z_i^l(\mathbf{x})=b_i^l+\sum _{j=1}^{N_l}{W}_{ij}^l{a}_j^l(\mathbf{x}),\text{mbox}with{a}_j^l(\mathbf{x})=\varphi (z_j^{l-1}(\mathbf{x})).$$

$$K^l(\mathbf{x},{\mathbf{x}}^{\prime })$$

$$K^0(\mathbf{x},{\mathbf{x}}^{\prime })=\mathbb{E}[z_j^0(\mathbf{x})z_j^0({\mathbf{x}}^{\prime })]={\sigma }_b^2+{\sigma }_w^2\left(\frac{\mathbf{x}\cdot {\mathbf{x}}^{\prime \top }}{d_{in}}\right).$$

Full text

I Introduction

With the rapid growth of distributed energy resources (DERs), the high level of uncertainty in power grids has been steadily increasing, rendering responsive and reliable load forecasting more critical and challenging than ever. A plethora of methods have been proposed to meet the challenges in load forecasting, which generally fall into two categories: one-stage and two-stage methods. A typical one-stage method adopts a single prediction model for load forecasting, and various models have been considered in the literature. For example, in [1] and [2], standard multiple linear regression models are enhanced by different techniques to tackle large datasets for load forecasting. In [3] and [4], modified ARIMA models are used to capture customers’ load patterns. In [5] and [6], support vector regression algorithms are modified to achieve a high predictive accuracy especially when handling high-resolution distribution systems. In [7, 8, 9], different types of neural network (NN)-based methods (e.g., recurrent NN, wavelet NN, and fuzzy NN) are proposed for modeling the complex relationships between customers’ loads and input variables. Last but not least, Gaussian process (GP)-based methods that adopt various forms of covariance kernels and input variables are studied in [10, 11, 12] to obtain an informative interval prediction for the target customer’s load.

As the name suggests, two-stage methods, on the other hand, perform load forecasting in two stages, with each stage having a separate model for a different purpose. The second-stage model is dedicated to the target customer’s load prediction, while the first-stage model is typically set up to process input features (i.e., either select useful features or provide accurate input estimates) to carry out the second-stage load forecasting. Two-stage methods are known for their enhanced predictive performance thanks to the first-stage processing [13]. The interested reader is referred to [13, 14, 15, 16, 17] and references therein for a review of two-stage load forecasting methods.

Despite their successful applications reported, most existing methods have one major drawback: they fail to thoroughly address the impact of input errors when performing point and interval predictions for the target customer’s load. Input errors arise from the use of the estimates of relevant input features for load forecasting. The commonly used input features include future weather conditions (e.g., temperature and cloud coverage), future power states (e.g., voltages and phase angles), etc. As these future input values are not known at the time of prediction, they must be estimated, say, by domain experts and prediction models. The predicted input values are inevitably susceptible to errors, which can severely undermine the ultimate load forecasting performance achieved. However, a majority of existing methods neglect or underestimate the corresponding impact as if the input data used for prediction were free of errors.

There indeed exist some two-stage methods that tackle input errors to mitigate their impact on point predictions to some extent. The main idea is to apply noise filtering techniques (e.g., Kalman filter) to filter out input errors. However, this practice does not tackle the impact of input errors on interval predictions, neither does it provide any input uncertainty quantification. A better way is to appropriately quantify the impact of input errors and reflect it in both point and interval load predictions.

Motivated by the aforementioned discussion, we contribute to the load forecasting literature by proposing a two-stage approach to thoroughly address the impact of input errors through input modeling and uncertainty quantification. Specifically, in the first stage, we adopt a hybrid model of NN and GP to provide an accurate point prediction for future input feature values, together with an interval prediction to quantify the input uncertainty. In the second stage, the point and interval predictions obtained are seamlessly utilized by a GP model that can handle noisy inputs, which forms the ultimate point and interval predictions for the target customer’s future load with the impact of input errors taken into account. Figure 1 shows a schematic diagram of the proposed two-stage load forecasting approach aimed at improving load forecasts. To enhance the computational efficiency of the proposed two-stage approach, we further propose to perform input-space dimensionality reduction via a model-free global sensitivity analysis (GSA) method for input feature selection.

The proposed two-stage approach is tested on various IEEE distribution system test cases, namely, 8-bus, 14-bus, 24-bus and 123-bus systems, in comparison with representative one-stage and two-stage methods. The proposed two-stage approach outperforms the competing methods in terms of both point and interval predictions across all cases tested.

The remainder of the paper is structured as follows. Section II elaborates on the proposed two-stage approach. Section III presents the GSA method that aims at enhancing the computational efficiency. Section IV presents numerical experiments for comparing the proposed two-stage approach with competing forecasting methods. Section V concludes the paper.

II Load Forecasting with Input Modeling and Uncertainty Quantification

In this section, we elaborate on how the proposed two-stage approach can thoroughly address the issue of input errors by quantifying and propagating their impact into the ultimate load forecasts.

For ease of exposition, we first define the input features to be used in the second-stage model. Consider predicting customer $i$’s load at hour $t$ in a distribution grid comprising $N$ customers. We adopt the following set of input features:

[TABLE]

with $\theta_{i,k}$ denoting the difference in phase angles of customers $i$ and $k$. Such a choice was first adopted by the IGP framework proposed in [17], which is inspired by the well-known power flow equation:

[TABLE]

where $P_{i}$ denotes the load of bus $i$ in a distribution grid, $|V_{i}|$ is the voltage magnitude of bus $i$, and $G_{ik}$ and $B_{ik}$ are coefficients describing the relationship between buses $i$ and $k$, with $i$, $k=1,2,\ldots,N$. Equation (2) implies that a given target customer’s load is strongly related to the phase angle differences between other customers and herself.

It is known in the literature that there typically exists a strong coupling between load and phase angle, whereas the coupling between load and voltage magnitude is rather weak [18]. Hence, the use of phase angle information as inputs is arguably sufficient for load forecasting purposes. We note that other relevant input features such as weather conditions and geographic locations can be easily incorporated into the proposed two-stage approach as well. For simplicity, we keep the input features as given in (1) in this work.

II-A Input Modeling and Error Quantification in the First Stage

II-A1 Model Choice

Consider predicting target customer $i$’s hour-ahead load $P_{i}^{t+1}$ at hour $t$, the second-stage model in our two-stage approach requires an input vector $\mathbf{x}_{t+1}$ in the form of (1) for load prediction. Since the input values in $\mathbf{x}_{t+1}$ are unknown to us at hour $t$, the first-stage model is adopted to perform hour-ahead input prediction. To fully quantify the input prediction errors, the first-stage model must provide an interval prediction as well as a point prediction. GP becomes a natural choice to this end thanks to its inherent predictive uncertainty measure. However, GP is known for its incompetent extrapolation performance [19], hence unsuitable as the first-stage model for predicting future input values. In contrast, NNs are known for their ability to capture complex output patterns and their outstanding extrapolation performance. However, NNs do not come with a predictive uncertainty measure. We propose to use the neural network–Gaussian process (NNGP) model recently proposed in [20] to accomplish both tasks simultaneously. The underlying idea of NNGP is to place prior distributions on the weight and bias parameters of a deep neural network (DNN) with infinite width and derive its equivalent GP formulation. By using the equivalent GP for prediction, NNGP can combine the strengths of NN and GP, providing a highly accurate point prediction as well as an informative interval prediction for future input values.

Below we briefly explain the rationale behind NNGP. Consider a fully-connected DNN with $L$ hidden layers ($L\geq 2$), where the number of hidden units in each layer, $N_{u}$, are equal. Denote the pointwise nonlinearity function by $\phi$. Let $\mathbf{x}\in\mathbb{R}^{d_{in}}$ and $\mathbf{y}\in\mathbb{R}^{d_{out}}$ denote the input and output vectors of the network, respectively. Consider the $i$th hidden unit in the $l$th layer, and denote the post-nonlinearity and the post-affine transformation by $a_{i}^{l}$ and $z_{i}^{l}$, respectively. Hence, we have $\mathbf{y}=\left(z_{1}^{L}(\mathbf{x}),z_{2}^{L}(\mathbf{x}),\ldots,z_{d_{out}}^{L}(\mathbf{x})\right)^{\top}$. Denote the weight and bias parameters in the $l$th layer by $W_{ij}^{l}$ and $b_{i}^{l}$; we place normal prior distributions with zero means and variances $\sigma_{w}^{2}/N_{u}$ and $\sigma^{2}_{b}$ respectively on $W_{ij}^{l}$ and $b_{i}^{l}$.

To derive an equivalent GP for a generic DNN at hand, we first show how to do so for a one-layer NN with infinite width. Consider the output of the $i$th hidden unit of the first layer in the NN, we have

[TABLE]

where the superscript “0” denotes the input layer of the NN, and $x_{k}$ denotes the $k$th component of the input vector $\mathbf{x}$. As $z_{i}^{1}(\mathbf{x})$ is the sum of independent and identically distributed terms, the standard Central Limit Theorem indicates that $z_{i}^{1}(\mathbf{x})$ becomes normally distributed as $N_{u}\rightarrow\infty$. Hence, given any input vectors $\mathbf{x}^{1}$, $\mathbf{x}^{2}$, …, $\mathbf{x}^{n}$, the vector $\mathbf{z}_{i}^{1}=\left(z_{i}^{1}(\mathbf{x}^{1}),z_{i}^{1}(\mathbf{x}^{2}),\ldots,z_{i}^{1}(\mathbf{x}^{n})\right)^{\top}$ follows a multivariate normal distribution. That is, $\mathbf{z}_{i}^{1}\sim\mathcal{N}(\mu^{1},K^{1})$, where $K^{1}$ denotes the $n\times n$ covariance matrix for the first layer, with its $(g,h)$th component given by $K^{1}(g,h)=\mathrm{Cov}(z_{i}^{1}(\mathbf{x}^{g}),z_{i}^{1}(\mathbf{x}^{h}))$, $g,h=1,2,\ldots,n$. Since the parameters $W_{ij}^{1}$’s and $b_{i}^{1}$’s have zero means, we have $\mu^{1}(\mathbf{x})=\mathbb{E}[z_{i}^{1}(\mathbf{x})]=0$ for any $\mathbf{x}\in\mathbb{R}^{d_{in}}$. In particular, given any two inputs $\mathbf{x}$, $\mathbf{x}^{\prime}\in\mathbb{R}^{d_{in}}$, it follows that

[TABLE]

We can see from the above that each component of the output vector of a one-layer NN with infinite width is equivalent to a GP with a covariance function specified by (3).

The aforementioned result can be generalized to NNs with more than one layer in a recursive manner by utilizing the relationship between $a_{j}^{l}$ and $z_{j}^{l-1}$, i.e. $a_{j}^{l}(\mathbf{x})=\phi(z_{j}^{l-1}(\mathbf{x}))$. The resulting covariance function in the $L$th layer, $K^{L}(\mathbf{x},\mathbf{x}^{\prime})$, hence serves as the covariance kernel of the equivalent GP for a DNN with $L$ hidden layers. For the sake of brevity, we refer the interested reader to Appendix A for a detailed derivation.

II-A2 Input Prediction and Error

Quantification

In the first stage, we predict the hour-ahead input vector component by component via NNGP. For predicting each component, we utilize not only their own values observed in the past but also the observed values of other components.

Let us consider predicting the input $\theta_{i,j}$ at hour $h^{*}$, for $j=1,2,\ldots,i-1,i+1,\ldots,N$ (see the definition given in (1)). We construct a training set for NNGP as $\{D_{j}^{h^{*}-24n_{t_{1}}},D_{j}^{h^{*}-24n_{t_{1}}+1},\ldots,D_{j}^{h^{*}-1}\}$, where $D_{j}^{h}$ denotes an input-output pair $(X_{j}^{h},Y_{j}^{h})$ for $h=h^{*}-24n_{t_{1}},h^{*}-24n_{t_{1}}+1,\ldots,h^{*}-1$, and $n_{t_{1}}$ denotes the number of days from which the training points are sampled. For our purpose, the output layer $Y_{j}^{h}$ is in fact $\theta_{i,j}^{h}$, namely, the phase angle difference between customer $i$ and customer $j$ at hour $h$. Taking into account the correlations between different customers (recall (1), (2) and discussions therein), we set $X_{j}^{h}=(\Theta_{i,j}^{h\top},\Theta_{i,j^{1}}^{h\top},\ldots,\Theta_{i,j^{n_{j}}}^{h\top})^{\top}$, where $\Theta_{i,j}^{h}=(\theta_{i,j}^{h-24n_{in}},\ldots,\theta_{i,j}^{h-1})^{\top}$, with $n_{in}$ denoting the number of days from which the data in the input layer are sampled, and the indices $j^{1},j^{2},\ldots,j^{n_{j}}$ are selected from $\{1,2,\ldots,N\}$ excluding $i$ and $j$.

Upon constructing the training set, NNGP can be trained using the covariance kernel obtained following the steps as given in Section II-A1. The predictive mean and variance of each component in the hour-ahead input vector can be obtained via NNGP in the same manner as using a GP model. Hence, the predictive mean of the input vector

[TABLE]

and predictive variance of the input vector

[TABLE]

are readily available, and can then be used in the second-stage model for making ultimate load forecasts.

The choice of the parameters for NNGP (i.e., $n_{t_{1}}$, $n_{in}$, $\sigma_{b}$, $\sigma_{w}$, and the number of hidden layers $L$) greatly impacts the predictive accuracy achieved. Appropriate parameter values can be obtained by a cross-validation procedure, and we refer the interested reader to [20] for details.

II-B Input Error Propagation and Uncertainty Quantification in the Second-Stage Load Forecasting

II-B1 Model Choice

In the second stage, we aim to appropriately propagate the impact of input errors assessed in the first stage into the ultimate load forecasting results, reflecting it in both point and interval predictions. The noisy input Gaussian process (NIGP) model proposed in [21] is arguably an adequate choice as the second-stage model. Compared to standard GP models, NIGP can propagate the impact of input errors into both point and interval predictions, providing a thorough solution to the issue of input errors in load forecasting.

Below we briefly describe the rationale behind NIGP. Assume that the observed output $y\in\mathbb{R}$ is a noisy measurement of an actual output $\tilde{y}$, that is,

[TABLE]

where the output noise is assumed to be normally distributed, i.e., $\varepsilon_{y}\sim\mathcal{N}(0,\sigma^{2}_{y})$, with $\sigma^{2}_{y}$ denoting the output noise variance. Assume that the observed input vector $\mathbf{x}\in\mathbb{R}^{d}$ can be modeled as the actual input vector $\tilde{\mathbf{x}}$ corrupted by some input noise:

[TABLE]

where the input noise vector is assumed to be normally distributed, i.e., $\bm{\varepsilon}_{x}\sim\mathcal{N}(0,\Sigma_{x})$, with $\Sigma_{x}$ denoting the $d\times d$ input noise covariance matrix. Based on (6) and (7), we can establish the following relationship between the observed input-output pair $(\mathbf{x},y)$ as follows:

[TABLE]

where $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ denotes the target input-output relationship to estimate. Suppose that $f(\cdot)$ can be modeled as a GP with some covariance kernel. The first-order Taylor expansion of $f(\mathbf{x}-\bm{\varepsilon}_{x})$ yields

[TABLE]

where $\partial{f(\cdot)}/\partial{\mathbf{x}}$ denotes the partial derivative process of $f(\cdot)$, which is also a GP [22]. Given that the distribution of the product of two normal random vectors $\bm{\varepsilon}_{x}$ and $\partial{f(\mathbf{x})}/\partial{\mathbf{x}}$ is analytically intractable, we approximate $\partial{f(\mathbf{x})}/\partial{\mathbf{x}}$ by the gradient of the posterior mean of $f(\mathbf{x})$, and denote it by $\bm{\partial}\mathbf{\bar{f}}$. It then follows from (8) and (9) that

[TABLE]

We see from (10) that the observed output $y$ is approximately normally distributed, i.e., $P(y|f)=\mathcal{N}(f,\sigma_{y}^{2}+\bm{\partial}\mathbf{\bar{f}}^{\top}\Sigma_{x}\bm{\partial}\mathbf{\bar{f}})$.

Now suppose that the training set for NIGP comprises $N$ observed input-output pairs, $\{(\mathbf{x}_{i},y_{i})\}_{i=1}^{N}$. Denote the $N\times d$ input matrix by $\mathbf{X}$ and the $N$-dimensional output vector by $Y$. Let $\mathbf{C}$ denote the $N\times N$ covariance matrix obtained by evaluating the covariance kernel at $\mathbf{X}$, and denote the gradient of the posterior means at the $N$ training inputs by $\bm{\Delta}\mathbf{\bar{f}}$, which is an $N\times d$ matrix. When the popular squared exponential covariance kernel is used, the covariance between any two inputs $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$ is given by

[TABLE]

where $\sigma_{f}^{2}$ represents the process variance, and $\Lambda$ denotes the $d\times d$ diagonal matrix with its main diagonal elements being the length-scale hyper-parameters.

The expressions for the predictive mean and variance of NIGP can be derived when the test input is deterministic and stochastic. Denote the observed test input by $\mathbf{x}_{*}$ and the actual test input by $\tilde{\mathbf{x}}_{*}$. For the deterministic case, we have $\mathbf{x}_{*}=\tilde{\mathbf{x}}_{*}$. Hence,

[TABLE]

where $\mathop{\mathrm{missing}}{diag}(\mathbf{\Delta\bar{f}}\Sigma_{x}\mathbf{\Delta\bar{f}}^{\top})$ denotes the diagonal of the matrix $\mathbf{\Delta\bar{f}}\Sigma_{x}\mathbf{\Delta\bar{f}}^{\top}$, and $\mathbf{I}_{N}$ denotes the $N\times N$ identity matrix. Notice that $\mathbf{\Delta\bar{f}}\Sigma_{x}\mathbf{\Delta\bar{f}}^{\top}$ is a correction term added to the covariance matrix of a standard GP model to reflect the uncertainty in the inputs $\mathbf{X}$.

In the stochastic case, assume that the observed test input $\mathbf{x}_{*}$ is multivariate normally distributed, i.e., $\mathbf{x}_{*}\sim\mathcal{N}(\tilde{\mathbf{x}}_{*},\Sigma_{x})$. It follows that the predictive mean of $f(\mathbf{x}_{*})$ can be given as

[TABLE]

Compared to the predictive mean given in (12), $\mathbf{q}$ is the counterpart of $\mathbf{C}(\tilde{\mathbf{x}}_{*},\mathbf{X})$ in the stochastic case, which is given by

[TABLE]

where $p(\tilde{\mathbf{x}}_{*}|\mathbf{x}_{*},\Sigma_{x})$ is the posterior of $\tilde{\mathbf{x}}_{*}$. When the squared exponential covariance kernel in (11) is adopted, $\mathbf{q}$ is an $N\times 1$ vector whose $i$th component is given by

[TABLE]

where $\mathbf{I}_{d}$ denotes the $d\times d$ identity matrix. Utilizing the law of total variance, the predictive variance of $f(\mathbf{x}_{*})$ follows as

[TABLE]

where $\mathop{\mathrm{missing}}{tr}(\bm{M})$ denotes the trace of matrix $\bm{M}$, and $\mathbf{Q}$ is an $N\times N$ matrix whose $(i,j)$th component is given by

[TABLE]

and $\bm{\alpha}=(\mathbf{C}+\sigma^{2}_{y}\mathbf{I}_{N}+\mathop{\mathrm{missing}}{diag}(\mathbf{\Delta\bar{f}}\Sigma_{x}\mathbf{\Delta\bar{f}}^{\top}))^{-1}Y.$ We refer the interested reader to [21, 19] for more details on the derivation of the predictive mean and variance functions of NIGP. A MATLAB package is available for practical implementation of NIGP [23].

II-B2 Input Error Propagation and Quantification in Load Forecasting

In the second stage, NIGP can be conveniently implemented for making ultimate load forecasts. As NIGP is essentially a variant of GP, its training and prediction can be performed in the same manner as using a standard GP. Consider predicting the target customer $i$’s hour-ahead load $P_{i}^{h^{*}}$ at hour $h^{*}-1$. We first estimate the hyper-parameters $\sigma^{2}_{f}$, $\Lambda$ and $\sigma^{2}_{y}$ via maximum likelihood estimation with a set of training input-output pairs $\{(\mathbf{x}_{h},P_{i}^{h})\}_{h=1}^{n_{t_{2}}}$, where $n_{t_{2}}$ denotes the number of days from which the training points are sampled, and $\mathbf{x}_{h}$ is in the form of (1). With the estimates of the hyper-parameters, the load forecast $\widehat{P}_{i}^{h_{*}}$ can be made with the predictive mean and variance of the input vector, $\mu_{1}$ and $\mathsf{V}_{1}$, respectively given by (4) and (5). Specifically, set $\Sigma_{1}=\mathop{\mathrm{missing}}{diag}(\mathsf{V}_{1})$, a diagonal matrix whose main diagonal entries are given by the components of $\mathsf{V}_{1}$. Then the predictive mean and variance for the target customer’s hour-ahead load $P_{i}^{h^{*}}$ can be obtained by replacing $\tilde{\mathbf{x}}_{*}$ and $\Sigma_{x}$ in (13) and (II-B1) respectively by $\mu_{1}$ and $\Sigma_{1}$. The predictive mean can be used as a point prediction for $P_{i}^{h^{*}}$, and the prediction variance can be utilized to construct an interval prediction for $P_{i}^{h^{*}}$.

III Improving Computational Efficiency via Global Sensitivity Analysis

In this section, we aim to improve the computational efficiency of the proposed two-stage load forecasting approach. Due to the use of spatial information (i.e. phase angle difference) for load forecasting, the input-space dimensionality can increase dramatically with the number of customers in the power grid, resulting in a high computational cost. To reduce the input-space dimensionality, one can consider retaining only a few most important input dimensions for load forecasting. Although many feature selection techniques (e.g., filter methods, wrapper methods, and embedded methods, see [24, 25]) are available, most of them are parametric approaches, hence pose restrictive assumptions on the underlying input-output relationship. In this work, we adopt the variance-based GSA method, which is model-free. GSA can be applied to both stages of our two-stage approach to improve the computational efficiency achieved.

III-A *Variance-Based GSA for Systems with Functional Inputs *

GSA focuses on quantifying how sensitive the output is to each individual input feature and their interactions, and identifies those inputs that contribute the most to the output variability. GSA methods fall into two categories: regression-based and variance-based. The main idea of the variance-based GSA methods is to decompose the variance of the output as a sum of contributions of each input feature. One of the most widely used variance-based GSA methods is the Sobol’ index method, on which we provide a brief overview below.

Consider the following model: $X\mapsto Y=f(X)$, where $f(\cdot)$ denotes the underlying input-output function, $X=(X_{1},X_{2}\ldots,X_{m})^{\top}$ denotes the $m\times 1$ vector of input features, and $Y$ represents the output. The first-order Sobol’ index of $X_{i}$ is defined as

[TABLE]

which quantifies the impact of $X_{i}$ on $Y$ only, excluding any interactions with other inputs. The total Sobol’ index of $X_{i}$ is defined as

[TABLE]

where $X_{\sim i}$ denotes the collection of all inputs except $X_{i}$. The total Sobol’ index measures the total impact of $X_{i}$ on $Y$, including higher-order impacts through interactions with other inputs. By definition, a Sobol’ index takes a value in $[0,1]$. The larger an index value is, the greater impact the associated input has on the output.

A majority of existing methods proposed for Sobol’ index estimation focus on tackling spatial inputs, see, e.g., [26, 27, 28, 29]. In our problem context, however, methods capable of handling functional inputs (i.e., inputs that vary with time) are required. GSA methods that tackle systems with functional inputs are relatively scarce. In [30], the authors proposed an effective Sobol’ index estimation method for GSA of dynamic systems, which is suitable to be applied to the first- and second-stage analyses to improve the computational efficiency achieved.

To estimate the Sobol’ indices, we first decompose the functional inputs through a simultaneous principle component analysis as follows:

[TABLE]

where $\bm{\gamma}=\{\gamma_{1},\gamma_{2},\ldots,\gamma_{p}\}$ denotes the collection of basis functions such as B-splines and wavelets, and $\bm{\beta}=\{\beta_{1},\beta_{2}\ldots,\beta_{p}\}$ denotes the set of corresponding coefficients, with $p$ denoting the number of basis functions used. Then we approximate the joint distribution of the coefficients in $\bm{\beta}$ using a Gaussian mixture model. Subsequently, we can construct Monte Carlo (MC) samples of the functional inputs based on (15) by sampling from the approximated distribution of $\bm{\beta}$. Given a sample of input features, we can obtain the corresponding outputs $Y$ using an approximated input-output function $\widehat{f}$. The first-order and total Sobol’ indices can then be estimated via MC sampling. Taking into account the potential interactions between the components of the input vector, we adopt the total Sobol’ index as the importance measure and use the following estimator proposed in [31]:

[TABLE]

where $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{C}$ respectively denote three independently generated $n\times m$ matrices of input samples, each row of which gives a randomly sampled input vector. In (16), $\widehat{f}(\mathbf{A})$ evaluates the approximated function $\widehat{f}$ at the input vector given by each row of matrix $\mathbf{A}$, and the subscript $j$ of $\widehat{f}(\mathbf{A})_{j}$ denotes the $j$th row of the resulting matrix $\widehat{f}(\mathbf{A})$. $\mathbf{A}_{\mathbf{B}}^{(i)}$ denotes the matrix whose $i$th column is taken from matrix $\mathbf{B}$ while the other $m-1$ columns remain from matrix $\mathbf{A}$.

III-B *Speeding Up Input Prediction and Uncertainty

Quantification*

The purpose of performing GSA in the first stage is to select a set of customer indices $\{j^{1},j^{2},\ldots,j^{n_{j}}\}$ from $\{1,2,\ldots,N\}$ (see Section II-A2) to help with input prediction. Consider predicting the value of input $\theta_{i,j}$ at each hour on a particular day, namely, $\theta_{i,j}^{h^{*}},\ldots,\theta_{i,j}^{h^{*}+23}$. Prior to prediction, we first construct a training set with full-dimensional inputs for the GSA purpose. Specifically, we set $n_{j}=N-2$ and $n_{in}=1$ in Section II-A2. Hence, each input dimension corresponds to the phase angle difference of the target customer and another customer at one of the previous 24 hours and the input-space dimensionality is $24(N-1)$. Recall that in the first stage, the input-output relationship is approximated by NNGP, hence the first-stage GSA is carried out based on the NNGP model to obtain the total index estimates, $\widehat{S}^{1}_{T_{i}}$, $i=1,2,\ldots,24(N-1)$. Since each index $i$ corresponds to a particular hour and a particular customer, we then quantify the impact of each component of the input vector $(\theta_{i,1},\ldots,\theta_{i,i-1},\theta_{i,i+1},\ldots,\theta_{i,N})^{\top}$ by the following summary statistic:

[TABLE]

Here $H_{\ell}$ measures the impact of each input feature on the output in the first-stage input prediction, and the larger the value, the greater the impact. We sort the $H_{\ell}$’s in a non-increasing order and select the first $n_{j}$ corresponding components (excluding $\Theta_{i,j}^{h}$) to be included in $X_{j}^{h}$ for input prediction as detailed in Section II-A2. The value of $n_{j}$ can be determined together with other parameters of NNGP by running a cross-validation procedure.

III-C Speeding Up Load Forecasting

In the seccond stage, we aim to identify and select a few components in the input vector $(\theta_{i,1},\ldots,\theta_{i,i-1},\theta_{i,i+1},\ldots,\theta_{i,N})^{\top}$ to be used in the NIGP model for predicting the target customer $i$’s hour-ahead load. We carry out GSA based on the NIGP model trained in the second stage following the steps detailed in Section III-A. We obtain the total Sobol’ index estimate, $\widehat{S}_{T_{j}}^{2}$, for each input $\theta_{i,j}$. The higher $\widehat{S}_{T_{j}}^{2}$ is, the more important $\theta_{i,j}$ is considered for load prediction of the target customer $i$. Upon identifying the most important components in the input vector, only those components will be retained as inputs for future load prediction. Hence, a higher computational efficiency can be achieved.

IV Numerical Experiments

In this section, we numerically evaluate the predictive performance of the proposed two-stage load forecasting approach against two benchmark methods, i.e., a representative two-stage model — IGP proposed in [17] and a representative one-stage model — SARIMA studied in [3]. We refer to the proposed two-stage approach as NNGP-NIGP. As mentioned in Section II, IGP adopts the same input features as NNGP-NIGP. Nonetheless, the first stage of IGP only provides a point prediction for the hour-ahead input vector via $k$-means clustering, and fails to quantify and propagate the impact of input errors into ultimate load forecasts. Without a mechanism to properly tackle the uncertain inputs, IGP constructs interval predictions for the target customer’s load via a heuristic approach. As will be shown later on, the comparison of NNGP-NIGP and IGP clearly manifests the substantial improvement in the predictive performance brought about by proper input-error propagation and quantification. The SARIMA model, on the other hand, modifies the traditional seasonal ARIMA model with a multi-model partitioning algorithm, and uses the target customer’s data for forecasting. This one-stage model serves as a baseline for comparison in the numerical experiments.

IV-A Experimental Setup

Experiments are conducted on four IEEE test cases, respectively, 8-bus, 14-bus, 123-bus test cases as well as a 24-bus test case which is constructed to simulate the phase unbalance problem based on the 8-bus test case. To simulate highly uncertain load behaviors caused by DERs in real-life power systems, historical load profiles drawn from real-world data sources are used for simulations. Taking into account the uncertain renewable generation behaviors of DERs, we first pre-process the hourly PV generation data over a year drawn from Renewable.ninja [32], and then subtract the pre-processed data from the load data of each bus. To obtain phase angle values, we perform power flow analysis to generate the states of the power system hourly over a one-year time frame using the MATLAB Power System Simulation Package (MATPOWER) [33, 34], based on the processed load data.

The predictive performance of each approach is evaluated by the accuracy of the point prediction and the coverage ability of the predictive interval obtained. To evaluate point predictive accuracy, we use the mean absolute percentage error (MAPE), which is defined below:

[TABLE]

where $P_{h}$ and $\widehat{P}_{h}$ respectively denote the load actually observed and the predicted load at hour $h$, and $T$ denotes the number of hourly predictions made. In our experiments, $T=24$ is used to calculate MAPE for every prediction day. Since the SARIMA model does not provide an interval prediction, we only compare the coverage abilities of the predictive intervals given by NNGP-NIGP and IGP. The coverage abilities are assessed by the coverage probability (CP):

[TABLE]

where $n_{c}$ denotes the number of observed load data points covered by the predictive interval on a prediction day. Predictions are performed over 30 days in each season for all customers in each test case considered.

To set appropriate parameter values for the first-stage NNGP model and GSA, a cross-validation procedure was adopted. For the 8-bus, 14-bus, and 24-bus test cases, the parameter setting adopted is $n_{t}=60$, $n_{in}=3$, $\sigma_{b}=1$, $\sigma_{w}=1$, and $n_{j}=1$ (the number of neighboring customers whose information is included for the first-stage input prediction; see Section III-B for details). For the 123-bus test case, the parameter setting used is $n_{t}=90$, $n_{in}=3$, $\sigma_{b}=1$, $\sigma_{w}=1$, and $n_{j}=2$. Regarding the second-stage NIGP model, the only parameters that require estimation are those in the covariance kernel, which are learned via maximum likelihood estimation. The numerical experiments are performed on a laptop with 6th generation Intel®  CoreTM i7 processor and 8.0GB DDR4 memory.

IV-B Comparison of Point and Interval Predictions

The MAPEs obtained by NNGP-NIGP, IGP, and SARIMA in all test cases are summarized in Fig. 2. We see that across all test cases, NNGP-NIGP outperforms IGP and SARIMA in terms of the point estimation accuracy. This is because NNGP-NIGP can properly incorporate the impact of input errors in the ultimate point prediction for the target customer’s load, while the competing methods fail to do so.

The CPs obtained by NNGP-NIGP and IGP in all test cases are summarized in Fig. 3. We see that NNGP-NIGP has a higher coverage probability than IGP. To closely examine the ability of NNGP-NIGP in appropriately propagating the impact of input errors into the predictive uncertainty associated with the target customer’s load forecasts, we show in Fig. 4 the predictive intervals obtained by NNGP-NIGP and IGP for one arbitrarily chosen customer on a typical prediction day in all test cases. Notice that in Fig. 4, the predictive interval of NNGP-NIGP can cover more observations than IGP due to suitable input error propagation. Moreover, the predictive interval given by NNGP-NIGP is much narrower than that given by IGP. In fact, the predictive interval of IGP is found sometimes wider than actually needed, which may result in unnecessary waste of resources if used for planning and operation purposes.

IV-C *Applying GSA to Improve Computational Efficiency *

We demonstrate that the proposed GSA method can effectively identify those important input features for load forecasting and improve the computational efficiency of NNGP-NIGP. Table I shows the estimated total Sobol’ indices obtained in the second stage of NNGP-NIGP for load prediction of customer 1 on three different days in the 8-bus test case. Since the estimated Sobol’ index value corresponding to the input $\theta_{1,2}$ is significantly higher than those for the other inputs, $\theta_{1,2}$ is identified as the most important input for load prediction of customer 1. This conclusion is consistent with that of IGP given in [17], where a heuristic feature selection method relying on the automatic relevance determination kernel is used. As presenting the GSA results for the other test cases (especially for the 123-bus test case) would require considerable space, we mention without showing details that the GSA results are highly consistent with those obtained by IGP in other test cases as well. Table II shows the average computational times consumed by NNGP-NIGP with and without GSA in all test cases. We see that the computational time of NNGP-NIGP is significantly reduced when only those important input features are used for load forecasting.

V Conclusion

This work represents one of the first efforts to thoroughly address the impact of input errors on both point and interval predictions for the target customer’s load. We identify and verify the opportunity of improving load forecasting performance by incorporating suitable input modeling and uncertainty quantification in a two-stage approach. The first stage delivers point and interval estimates for future input feature values, which are to be used in the second-stage model. The second stage propagates the impact of input errors into the ultimate point and interval predictions for the target customer’s load. We select a representative two-stage load forecasting model, IGP, as the main benchmark model, which also performs input prediction, nevertheless via a much less sophisticated approach. The proposed two-stage approach is found to outperform IGP and a representative one-stage model, SARIMA, with respect to both point and interval predictions in various test cases. In addition, a model-free functional GSA method is considered to reduce the input-space dimensionality for an enhanced computational efficiency. The numerical results show that the GSA method can significantly improve the computational efficiency of the proposed two-stage approach.

Appendix A Derivation of the equivalent GP for NNGP

In Section II-A1, we have derived the equivalence between a one-layer neural network with infinite width and a GP. The covariance function of the corresponding GP is given in (3). In this section, we proceed to generalize the derivation to an NN with more than one layer (i.e., a DNN) by induction, and the key lies in deriving the covariance kernel of the equivalent GP that corresponds to a DNN. This step is particularly important for load forecasting, as it ensures that a DNN can be regarded as a GP, which can be conveniently used for obtaining point and interval estimates for future input values in a standard manner [22]. Consider the $l$th layer $(l\geq 2)$ of the DNN, we have

[TABLE]

It follows that $z_{i}^{l}\sim\mathcal{N}(0,K^{l})$, where $K^{l}$ is the covariance matrix in the $l$th layer, whose entry is given by the covariance between the $z_{i}^{l}$ values corresponding to a pair of inputs. For any two inputs $\mathbf{x}$ and $\mathbf{x}^{\prime}$, we have

[TABLE]

Notice that the joint distribution of $z_{i}^{l-1}(\mathbf{x})$ and $z_{i}^{l-1}(\mathbf{x}^{\prime})$ is bivariate normal with zero mean and a $2\times 2$ covariance matrix whose distinct entries are given by $K^{l-1}(\mathbf{x},\mathbf{x}^{\prime})$, $K^{l-1}(\mathbf{x},\mathbf{x})$, and $K^{l-1}(\mathbf{x}^{\prime},\mathbf{x}^{\prime})$. To show the recursive relationship between $K^{l}$ and $K^{l-1}$, we rewrite (17) as

[TABLE]

where $F_{\phi}$ is a deterministic function that depends only on the nonlinearity function $\phi$. For the base case $K^{0}$, suppose that $W_{ij}^{0}\sim\mathcal{N}(0,\sigma_{w}^{2}/d_{in})$ and $b_{j}^{0}\sim\mathcal{N}(0,\sigma_{b}^{2})$, we have

[TABLE]

The covariance matrix for the $l$th layer can thus be obtained recursively via (18) and (19). The resulting covariance matrix $K^{L}$ for the $L$th layer is the ultimate covariance kernel adopted by the NNGP model.

To efficiently calculate the covariance kernel corresponding to the GP model equivalent to the DNN, an algorithm which utilizes the recursive relationship between the covariance matrix in adjacent layers (see (18)) is available.

Constructing the covariance matrix $K^{L}$ for the equivalent GP involves computing the Gaussian integral in (17) for all pairs of training-training and training-test points, recursively for all layers. For some forms of nonlinearity functions used by the NNGP model, the integration can be performed analytically; but for most forms of nonlinearity functions, numerical integrations are required. Computing integrals independently for each pair of data points and each layer can be computationally expensive, and a more efficient algorithm to recursively calculate the covariance kernel $K^{l}$ in each layer $l$ is available in [20]. We briefly summarize the algorithm as follows:

1. Generate: pre-activations $u=[−u_{max},\ldots,u_{max}]$ consisting of $n_{g}$ elements evenly spaced between $−u_{max}$ and $u_{max}$; variances $s=[0,s_{max}]$ with $n_{v}$ evenly spaced elements, where $s_{max}<u^{2}_{max}$; and correlations $c=(-1,\ldots,1)$ with $n_{c}$ linearly spaced elements.

2. Populate a matrix $F$ containing a lookup table for the function $F_{\phi}$ given in (18). This involves numerically approximating a Gaussian integral, in terms of the marginal variances $s$ and correlations $c$. The entries in the matrix $F$ are computed as

[TABLE]

where $U_{ab}=[u_{a},u_{b}]^{\top}$ with $u_{a}\in u$ and $u_{b}\in u$, and

[TABLE]

3. For every pair of data points $\mathbf{x}$ and $\mathbf{x}^{\prime}$ in layer $l$, compute $K^{l}(\mathbf{x},\mathbf{x}^{\prime})$ using (18). Approximate the function $F_{\phi}(K^{l-1}(\mathbf{x},\mathbf{x}^{\prime}),K^{l-1}(\mathbf{x},\mathbf{x}),K^{l-1}(\mathbf{x}^{\prime},\mathbf{x}^{\prime}))$ by bilinear interpolation into the matrix $F$ from Step 2, where we interpolate into $s_{i}$ using the value of $K^{l-1}(\mathbf{x},\mathbf{x})$, and interpolate into $c_{j}$ using $K^{l-1}(\mathbf{x},\mathbf{x}^{\prime})/K^{l-1}(\mathbf{x},\mathbf{x})$. Remember that $K^{l-1}(\mathbf{x},\mathbf{x})=K^{l-1}(\mathbf{x}^{\prime},\mathbf{x}^{\prime})$, due to data pre-processing to guarantee constant norm.

4. Repeat the previous step recursively for all layers.

Following the algorithm above, the covariance of any two arbitrary inputs in the $L$th layer can be calculated efficiently. The resulting covariance matrix in the $L$th layer is then used for prediction for the corresponding GP. We refer the interested reader to [20] for more details on this algorithm. A Python package is available for practical implementation of NNGP [35].