# Optimal Uncertainty-guided Neural Network Training

H M Dipu Kabira, Abbas Khosravia, Abdollah Kavousi-Fardb, Saeid Nahavandia,
Dipti Srinivasanc aInstitute for Intelligent Systems Research and Innovation
(IISRI), Deakin University, Australia.  bDepartment of Electrical and
Electronics Engineering, Shiraz University of Technology, Iran.  cDepartment
of Electrical & Computer Engineering, National University of Singapore,
Singapore.

###### Abstract

The neural network (NN)-based direct uncertainty quantification (UQ) methods
have achieved the state of the art performance since the first inauguration,
known as the lower-upper-bound estimation (LUBE) method. However, currently-
available cost functions for uncertainty guided NN training are not always
converging and all converged NNs are not generating optimized prediction
intervals (PIs). Moreover, several groups have proposed different quality
criteria for PIs. These raise a question about their relative effectiveness.
Most of the existing cost functions of uncertainty guided NN training are not
customizable and the convergence of training is uncertain. Therefore, in this
paper, we propose a highly customizable smooth cost function for developing
NNs to construct optimal PIs. The optimized average width of PIs, PI-failure
distances and the PI coverage probability (PICP) are computed for the test
dataset. The performance of the proposed method is examined for the wind power
generation and the electricity demand data. Results show that the proposed
method reduces variation in the quality of PIs, accelerates the training, and
improves convergence probability from 99.2% to 99.8%.

###### keywords:

Neural Network, Wind Power, Prediction Interval, Uncertainty Quantification,
Cost Function.

##  1 Introduction

All-natural quantities have some uncertainties. The quantity may slightly or
greatly vary for the same circumstances. The variance of the quantity at a
circumstance is the level of uncertainty for the circumstance. The level of
uncertainty is heteroscedastic. The predictability of the same quantity can be
different based on circumstances. Traditional point prediction systems predict
a value which is the most probable for the corresponding input combination.
The actual value may differ from the prediction slightly or greatly based on
circumstances [1]. For instance, electricity demand at off-peak and full-peak
hours is highly predictable but the demand on the transition times may largely
vary from one day to another day [2, 3]. The difference between the prediction
and the actual value may be caused by the modeling error, or the inherent
randomness of the system [4, 5]. The inclusion of some inputs such as the
current temperature and calendar information may reduce the uncertainty of
predictions. However, some portion of the uncertainty can be random and cannot
be predicted based on existing features. The uncertainty can also be
asymmetrically heteroscedastic and the point prediction with a certain error
possibility cannot provide adequate information to the user [6, 7].

Probabilistic forecast, such as an uncertainty bound is also popular in
decision-making [8, 9]. However, a single uncertainty bound is unable to
represent the level of uncertainty. Multiple uncertainty-bounds can be applied
to quantify uncertainties. PI is a recognized UQ method, applies the upper and
the lower bounds to quantify the level of uncertainty. Probabilistic forecasts
such as prediction intervals (PIs) with a certain coverage probability are
more appropriate for understanding the uncertain condition. Fig. 1 presents
the uncertainty captured by PIs with 99% coverage probability. The probability
density function changes from sample to sample [10, 11]. The width of the PI
varies from sample to sample based on the corresponding uncertainty. Decision-
makers get the most probable regions of targets from PIs, generated by NNs
even for an asymmetric and heteroscedastic system [12].

Figure 1:  Importance of the uncertainty quantification. The point prediction,
presented by the red line with a constant error possibility (the root mean
square error) cannot represent the heteroscedastic uncertainty. PIs,
represented by green lines becomes narrow in the less uncertain regions and
become wide in more uncertain regions. Therefore, PIs can represent
heteroscedastic uncertainty.

NN is commonly recognized as a black box to its end-user. It can provide
state-of-art performance with proper training. The designer of any NN chooses
NN-size, activation function and initial weights. Initial weights can also be
random or pre-defined. Proper training provides an optimal selection of
weights on interconnects and biases. The optimal NN output is a weighted sum
of inputs and functions of inputs. Weight optimization is performed through a
reward-based system. The reward is calculated based on the performance of NN
through a cost function in each cycle. Therefore, a cost function needs to be
designed considering the purpose, quality criteria, critical situations, and
the convergence of optimization [13].

This paper proposes an optimal PI construction technique considering different
aspects of recently proposed NN-based direct PI construction techniques. The
direct construction of PIs from the NN result in a sharp PI for any type of
probability distributions, such as skewed-Gaussian, log-normal and multimodal.
Several direct PI construction techniques have been proposed in the
literature. The relative performance of these techniques is questionable and
should be comprehensively investigated. Every new algorithm is claimed to be
the best with data analysis during its proposal. However, the result may vary
with different datasets. Therefore, we discuss the philosophy of developing
cost functions and provide a novel one. The proposed cost function combines
the important philosophies of recently proposed NN-based direct PI computation
methods. The paper presents a rigorous performance analysis for the wind power
generation and electricity demand data. The method is also applicable for the
UQ of several other datasets, such as electricity prices, and other renewable
generations. The improvements in the convergence for electricity demand, hydro
and solar generations are also analyzed. Moreover, weather, geographical
positions, and various human-made events can be considered as an input of the
NN. The effect of many events may not be expressed through mathematical
equations. NNs can find all of those hidden relations with reward-based
training.

The paper is organized with the following flow of information. Section 2
presents the basics of the uncertainty quantification and the advantage of the
NN-based direct uncertainty quantification. Section 3 presents all NN-based
direct PI construction methods. Section 4 presents the proposed PI
construction method. Section 5 reports the simulation results and performance
metrics. Section 6 is the concluding section.

##  2 Uncertainty and its Quantification

###  2.1 Increased Uncertainty in Power Grid

All real-world events consist of sub-events, among which many are random.
However, their combined effect can be interval predictable or even
deterministic. When tossing a coin for a single time, the probability of
getting the head or the tail is equal. However, the outcome of one hundred
tossings of a coin is quite predictable. There is a 97.9% chance of getting 40
to 60 heads and there is 72.88% chance of getting 45 to 55 heads. Therefore,
20% region of the output range contains about 97-98% of the probability
density. Therefore, that outcome is interval predictable. The outcome of ten
thousand tossings is more predictable. There is a 95.56% chance of getting
4900 to 5100 heads. Therefore, 2% region of the output range contains about
95-96% of the probability density and it is point predictable with 1-2% root
mean squared error (RMSE) [14].

A grid is connected to millions of electrical appliances. Whether individual
equipment is consuming electricity or not is difficult to predict. However,
the total electricity consumption of a large grid is point predictable. A fair
coin has an equal probability of head or tail but an appliance may have a 10%
probability of consuming electricity and that probability varies depending on
time or weather or any other conditions. Therefore, the total power
consumption of a large grid is predictable and depends on major events, such
as weather, time, vacation, sports, etc.

Let us consider two events: the percentage of heads in tossing 1) five coins,
2) ten thousand coins. When the first event has a higher weight than the
second event on a quantity, the quantity becomes highly uncertain. The grid
also contains elements of different predictability. Wind power generation is
highly random. Installation of large scale wind power plants has made the
overall generation more unpredictable. The large-scale introduction of the
electric vehicle has made the overall consumption more unpredictable.

###  2.2 Uncertainty Quantification

All systems have some inherent randomness, known as the aleatory uncertainty.
The output value slightly or greatly varies for the same input combination.
The other type of uncertainty is epistemic or subjective which can be properly
captured by the precise modeling [15, 16]. This uncertainty can cause
significant prediction error when several secondary or tertiary effects are
overlooked during the modeling. The future value of the parameter or the
target ($t_{j}$) is represented as:

$t_{j}=y_{j}+\epsilon_{j}$ (1)

where $\epsilon_{j}$ is the zero expectation error signal and $j(=1,2,\cdots
n)$ is the sample number. Therefore, total uncertainty is represented as
follows:

$t_{j}-\hat{y}_{j}=[y_{j}-\hat{y}_{j}]+\epsilon_{j}$ (2)

where $\hat{y}_{j}$ is the true regression mean for the $j^{th}$ sample. When
the two terms in (2) are independent, the total variation associated with the
model outcome is represented as:

$\sigma_{j}^{2}=\sigma_{\hat{y}_{j}}^{2}+\sigma_{\hat{\epsilon}_{j}}^{2}$ (3)

The term $\sigma_{\hat{y}_{j}}^{2}$ represents the subjective uncertainty and
$\sigma_{\hat{\epsilon}_{j}}^{2}$ represents the inherent randomness.

Point prediction is the most widely used approach for predicting unknown
quantities. It provides a value corresponding to the lowest statistical error.
Statistical errors are measured by several criteria, such as the RMSE, mean
absolute percentage error (MAPE) or sum squared error (SSE). Therefore, the
user lacks information about the level of heteroscedastic uncertainty
($\sigma^{2}_{j}$). An interval forecast with a certain probability of
coverage expresses both uncertainties.

The PI for $j^{th}$ sample is traditionally represented as follows [17]:

$PI_{j,\ Traditional}^{(\alpha)}=[\mu_{j}-\lambda^{(\alpha)}\sigma_{j},\ \
\mu_{j}+\lambda^{(\alpha)}\sigma_{j}]$ (4)

Here, $\mu_{j}$ is the mean and $\lambda^{(\alpha)}$ is the target coverage
probability parameter. The value of $\lambda^{(\alpha)}$ is 1.15, 1.64 and
1.96 respectively for 75%, 90% and 95% theoretical coverage [18]. Therefore,
the width of the PI of $j^{th}$ sample is as follows:

$PIW_{j}=2\lambda^{(\alpha)}\sigma_{j}$ (5)

The distribution of uncertainties may not be purely Gaussian for any
combination of inputs. Some distributions can be log-normal, skewed Gaussian,
or multimodal. Therefore, PIs constructed through the Gaussian assumption
fails to maintain a narrow width and the required coverage simultaneously.

A relatively smarter approach to constructing PIs is considering the
cumulative distribution function. Conditional probability functions can
construct an interval of ($1-\alpha$) confidence level. The upper bound
($\overline{y}_{j}$) is a value greater than ($1-\alpha/2$) portion of the
probability density function. Therefore, $\overline{y}_{j}$ can be represented
as [19]:

$P_{j}(t_{j}<\overline{y}_{j})=1-\alpha/2$ (6)

where $P_{j}(condition)$ is the probability function for the target at
$j^{th}$ sample ($t_{j}$). The cumulative probability density function ($CP$)
can represent the relation as follows:

$CP(\overline{y}_{j})=1-\alpha/2$ (7)

Taking inverse:

$\overline{y}_{j}=CP^{-1}(1-\alpha/2)$ (8)

The lower bound ($\underline{y}_{j}$) and the upper bound ($\overline{y}_{j}$)
form the PI. Therefore, the PI (=$[\underline{y}_{j},\overline{y}_{j}]$)
becomes as follows:

$PI_{j,\ Conditional}^{(\alpha)}=[CP^{-1}(\alpha/2),CP^{-1}(1-\alpha/2)]$ (9)

The direct NN-based PI computation techniques construct an optimal PI
empirically without considering any theoretical probability distribution.

Figure 2: The advantage of the NN-based direct PI construction over
traditional and conditional PI construction techniques for a non-Gaussian
probability distribution. The NN-optimization technique finds an optimal PI
for any arbitrary probability distribution.

Fig. 2 presents the advantage of the NN-based direct PI construction over
traditional and conditional PI construction techniques for a non-Gaussian
probability distribution, known as asymmetric heteroscedastic uncertainties
[20]. Although, there exist multiple variables in our NN based PI construction
system, Fig. 2 presents a one-dimensional input range for an easier
visualization [14]. Traditional and conditional PIs follow Eq. (4) and Eq. (9)
respectively and the LUBE NN is directly trained to optimize the cost
function. Conditional PI is optimal compared to traditional PIs. The NN-based
direct PI moves upper and lower bounds slightly and evaluates the effect
during the training to achieve an optimal PI coverage probability (PICP).

##  3 Direct NN-based PI Construction Methods

Several direct NN-based methods have been proposed by different groups. The
current section discusses relevant methods where a cost function is proposed
or modified to overcome a limitation.

###  3.1 Initial LUBE Method

The lower upper bound estimation ($LUBE$) method is the very first method for
constructing PIs through direct training [21]. The LUBE method is based on the
following philosophies:

  * 1.

The PI should cover samples with an equal or higher probability compared to
the PI nominal coverage (PINC = $1-\alpha$).

  * 2.

When the required PI coverage probability (PICP) is achieved (PICP $\geq$
PINC), the quality of PI depends on the narrowness of the average width.

The LUBE method minimizes the coverage width criterion (CWC). PICP, PINAW, and
CWC are defined as follows:

$PICP=\frac{1}{n}\sum_{j=1}^{n}c_{j}$ (10)

provided that,

$c_{j}=\begin{cases}1,t_{j}\in[\underline{y}_{j},\overline{y}_{j}]\\\
0,t_{j}\not\in[\underline{y}_{j},\overline{y}_{j}].\end{cases}$ (11)
$PINAW=\frac{1}{nR}\sum_{j=1}^{n}(\underline{y}_{j}-\overline{y}_{j}).$ (12)
$CWC=PINAW\\{1+\gamma^{(\alpha,PICP)}e^{\eta(PINC-PICP)}\\}$ (13)

provided that,

$\gamma^{(\alpha,PICP)}=\begin{cases}1,PICP<PINC\\\ 0,PICP\geq
PINC\end{cases}$

$\eta$ =50 is a hyperparameter used for penalizing low coverages. The LUBE
method also proposes the structure of the NN. Fig. 4 presents the structure of
a lower-upper-bound-based NN. In the LUBE NN, input and hidden layers are
shared. Therefore, the difference between the lower and upper bounds is the
weights of the connections between the last hidden layer to the output layer
and output layer biases.

###  3.2 $LUBE$ with Independent Width and Penalty Factors

The width factor (PINAW) exists in a multiplicative manner with the penalty
factor ($\gamma^{(\alpha,PICP)}$), as shown in (13) in the initial LUBE
method. Therefore, the optimization was controversial [22, 23]. One of the
concerns is the optimization at zero width (PINAW = 0), that happens
frequently with the initial LUBE method. Khosravi et. al. resolved most of
those controversies with the proposal of a new CWC equation of independent
width and penalty factor [24] as presented at Eq. (14).

$CWC=PINAW+\gamma^{(\alpha,PICP)}e^{\eta(PINC-PICP)}$ (14)

The function is modified based on the following philosophy:

  * 1.

A separate PINAW in an additive manner restricts optimization at PINAW = 0.

###  3.3 $LUBE$ with Continuous Cost Function

The NN training in LUBE method often fails to converge with the traditional
LUBE cost function. Therefore, the following equation is recently introduced
to achieve a continuous cost function at (PICP = PINC) [25]:

$CWC=PINAW+\gamma^{(\alpha,PICP)}$ (15)

provided that,

$\gamma^{(\alpha,PICP)}=\begin{cases}e^{\eta(PINC-PICP)}-1,PICP<PINC\\\
0,PICP\geq PINC\end{cases}$

where $\gamma^{(\alpha,PICP)}$ is the PI-failure penalty function. The value
of that function is called the penalty. $\gamma^{(\alpha,PICP)}$ is zero for
(PICP $\geq$ PINC). In contrast, its value exponentially increases with the
lowering of PICP when PICP $<$ PINC.

The concept of a continuous cost function is well known in the model-design
rules of circuit simulations [26, 27]. The philosophy behind the continuous
cost function is as follows:

  * 1.

A continuous model is the prerequisite for the convergence of the simulation.
Especially when the model is applied to optimize anything through an iterative
process.

  * 2.

A discontinuity on the model or it’s derivative may result in a very large
gradient ($\partial y/\partial w$) during an iteration, as the $\Delta y$
remains large for a very small $\Delta w$ at the point of discontinuity.
Therefore, the discontinuity may potentially result in the overflow while
nearby regions are mapped following the gradient.

As the most optimized point of the NN training is PICP = PINC, the step-sizes
are reduced near that point to achieve optimized weight with high precision
[28, 29]. That means, $\Delta x$ or a component of $\Delta x$ may become very
small while iterating near the most optimized point. Therefore, the cost
function needs to be smooth near the minima.

###  3.4 Can Wan’s Interval Score-based Cost Function

C Wan et al. [1] also propose a slightly different cost function for the NN-
based direct interval forecasting. They consider both width dependent and
coverage dependent components. The average coverage error (ACE) is the
coverage dependent component in C Wan’s cost function. ACE is defined as
follows:

$ACE=PICP-PINC$ (16)

where $PINC=1-\alpha$ is the target PI coverage, known as the PI nominal
coverage. PI-width of the $j^{th}$ sample is defined as follows:

$PIW_{j}=\overline{y}_{j}-\underline{y}_{j}$ (17)

C Wan et al. [30] also define a component named the interval score. The
interval score for $j^{th}$ sample is as follows:

$S_{j}=\begin{cases}-2\alpha\times
PIW_{j}-4(\underline{y}_{j}-t_{j}),t_{j}<\underline{y}_{j}\\\ -2\alpha\times
PIW_{j},\overline{y}\geq t_{j}\geq\underline{y}_{j}\\\ -2\alpha\times
PIW_{j}-4(t_{j}-\overline{y}_{j}),t_{j}>\overline{y}_{j}\end{cases}$ (18)

The average interval score is as follows:

$S_{AV}=\frac{1}{n}\sum_{j=1}^{n}S_{j}$ (19)

With the weighted summation of ACE and the interval score, the C Wan’s cost
function is as follows:

$\underset{NN}{\text{min}}\ \ \lambda|S_{AV}|+\gamma|ACE|.$ (20)

Both $\lambda$ and $\gamma$ are set to one to provide the equal weight towards
the PICP calibration and the PI sharpness.

According to our analysis, that cost function is based on the following
philosophies:

  * 1.

Both high and low PICPs are penalized with equal concentration to maintain a
gradient throughout the input domain.

  * 2.

Both of the width of the interval ($PIW_{j}$) and the failure distance
($\underline{y}_{j}-t_{j}$ or $t_{j}-\overline{y}_{j}$) are optimized.

  * 3.

The failure distance is considered with much higher priority ($2/\alpha$
times) compared to the width.

###  3.5 L. G. Mar´ıns’ Deviation from Mid Interval Consideration

L. G. Mar´ın et al. [31] consider the deviation from the mid interval and a
continuous cost function. Their cost function is as follows:

$\underset{NN}{\text{min}}\ \
\beta_{1}PINAW+\beta_{2}||e||^{2}+exp[-\eta\big{(}PICP-(1-\alpha)\big{)}]$
(21)

They also perform the normalization of $||e||^{2}$ and PINAW. They introduce
the deviation from mid interval as $||e||$ and named it as the error quantity
defined as follows:

$||e||=\sqrt{\sum_{j=1}^{n}\Big{|}t_{j}-\frac{\overline{y}_{j}+\underline{y}_{j}}{2}\Big{|}^{2}}.$
(22)

The function is proposed based on the following philosophy:

  * 1.

The most probable region of the target may not stay near the middle of the
interval in the direct PI construction method. Minimization of the deviation
of the target from the center can potentially shift the most probable region
near the center of PIs.

###  3.6 G. Zhang’s Deviation Information-based Criterion (DIC)

When the cost function only considers PICP and PINAW, the NN finds smart ways
to optimize the cost function. Though the PINAW becomes much smaller, the NN
often keeps a smaller width instead of covering the target in critical
situations. Therefore, target PICP is achieved with much smaller PINAW.
However, PI misses the target by a large distance in critical situations. The
situation is illustrated with a rough drawing in Fig. 3. Intervals covering
and missing the target are denoted by the green color and red color
respectively. In the critical situation, the NN aims to narrow down the PI
instead of covering it by increasing the width. However, the user of the
prediction algorithm needs indications of sudden rises or falls with higher
accuracy to manage critical situations.

G. Zhang et al. [32] tried to avoid the computation extensive exponential cost
function. They proposed deviation information-based criterion defined as:

$DIC=PINAW+\gamma^{(\alpha,PICP)}.pun$ (23)

provided that,

$pun=\sigma_{p}\sum_{j=1}^{N_{L}}(\underline{y}_{j}-t_{j})+\sigma_{p}\sum_{j=1}^{N_{U}}(t_{j}-\overline{y}_{j})$

and

$\gamma^{(\alpha,PICP)}=\begin{cases}1,PICP<PINC\\\ 0,PICP\geq
PINC\end{cases}$

The cost function is modified based on the following philosophy:

  * 1.

Avoiding computationally expensive exponential function.

  * 2.

Bringing the most probable region of the target near the middle of the PI.

Figure 3:  A rough diagram presenting NN based PIs when the optimization
considers only PICP and PINAW. PIs cover 80% to 90% targets but fail to
predict sharp changes. Successful and unsuccessful PIs are represented by
green and red lines respectively with upper and lower bound marks.

##  4 Proposed Method

In this paper, a smooth and customizable cost function is proposed for the
uncertainty guided NN training. Uncertainty guided NNs predict the upper bound
and the lower bound. Bounds are computed with NNs without any assumption on
the distribution.

Studying motives of all proposed algorithms, the following key criteria of a
good cost function are concluded:

  1. 1.

PICP, PINAW and PI normalized average failure distance (PINAFD) are important
parameters of an ideal cost function. The consideration of the deviation from
the mid-interval also results in a slightly lower failure distance but the PI
becomes much wider.

  2. 2.

Cost function and its derivatives need to be continuous for a better
convergence profile of the NN-training.

  3. 3.

Often a trained NN fails to maintain $PICP\geq PINC$ with the test dataset due
to the slight variation between the training and test datasets. Therefore, a
small coverage margin ($\delta$) is required during the training.

  4. 4.

The entire input domain needs to have single minima at $PICP=PINC+\delta$.
Except at the minima, the input domain needs to have non-zero gradients
directing the optimization towards the minima.

  5. 5.

A simpler or less computation extensive cost function is preferred. However,
the simplicity of the function is not related to the quality of PIs.

  6. 6.

Different users may have different preferences towards width, failure
distance, and coverage penalty.

PIs may fail to cover the target during a sudden fluctuation but the NN should
try to bring the nearby bound of the PI close to the target. When the nearby
bound of the PI is close to the uncovered target, the user can manage the
situation with a lower difficulty. The proposed NN training considers the
minimization of the prediction interval normalized failure distance (PINAFD).
The expression of PINAFD is as follows:

$PINAFD=\frac{\sum_{j=1}^{n}(1-c_{j})\times
min(|t_{j}-\overline{y}_{j}|,|\underline{y}_{j}-t_{j}|)}{R\times\Big{\\{}\sum_{j=1}^{n}(1-c_{j})\Big{\\}}+\epsilon}.$
(24)

where $c_{j}$ contains the same meaning as of Eq. (10). The minimum distance
of the target from bounds is considered to be the failure distance for the
corresponding sample when the target is not bounded. The total failure
distance is divided by the total number of missing samples
($\sum_{j=1}^{n}(1-c_{j})$) and normalized by the range ($R$) to achieve
normalized average value. $\epsilon$ is a small value to avoid an undefined
PINAFD value for 100% PICP. During an iteration of the training, NN may cover
all samples, resulting in $\sum_{j=1}^{n}(1-c_{j})=0$. In such a situation,
the value of PINAFD becomes zero by zero. To avoid that undefined value, a
small value, $\epsilon=1e-10$ is added to the denominator to make zero PINAFD
for 100% PICP.

We formulate the proposed optimization parameter by adding a weighted PINAFD.
The proposed optimization parameter is as follows:

$\underset{NN}{\text{min}}\ \ PINAW+\rho.PINAFD+\beta\times(1-\alpha+\delta-
PICP)^{2}$ (25)

where $\rho$ is the failure distance resistance parameter, $\beta$ is the PI
coverage penalty factor, and $\delta$ is the coverage margin. Usually, $\rho$
is set to one to provide an equal concentration towards the width and the
failure distance. The coverage penalty $\beta$ needs to be high enough to
provide higher concentration towards PICP than a slight shrink of average
width or failure distance. Therefore, the value of $\beta$ is usually set to
more than 200. The test PICP can be slightly lower than the train PICP due to
a slight variation among datasets. The initial LUBE method strictly maintains
the PICP higher than the nominal PICP. However, in a linear or a polynomial
penalty, the test PICP can be slightly lower than the nominal. Therefore, a
slight margin of $\delta=\alpha/50$ is kept. The variation of PICP is related
to the nominal PICP. When $\alpha$ is small, the sample density near the edge
of PIs is lower and the variation in PICP is also lower. With $\alpha=5\%$, a
slight margin of $\delta=0.1\%$ is enough for obtaining a test PICP of 95%
most of the time. Therefore, $\delta=\alpha/50$ is kept for increasing the
margin with increasing $\alpha$.

The proposed optimization parameter, presented in Eq. (25) is named as the
Coverage Width Failure Distance Criteria (CWFDC). The proposed cost function
is considering the coverage (PICP), the normalized average width (PINAW), and
the failure distance (FD). Besides these considerations, it is smooth and
considers a small margin to withstand a slight PICP variation.

##  5 Result Evaluation

The wind power generation and the electricity demand samples from August 2012
to August 2019 are downloaded from the UK-grid website [33]. Four recent
samples and the time in the hour ($Time_{Day}=hour+minutes/60$) on a
corresponding day is provided as the input to the NN. Fig. 4 presents the
structure of the NN with input-output combinations. We apply the simulated
annealing technique for NN training. NNs of different sizes and
initializations are trained with different cost functions to evaluate the
result. Different steps of the result evaluation are as follows:

  1. 1.

Finding the optimal NN-sizes for different cost functions, different PICP, and
different datasets.

  2. 2.

Evaluation of PINAW, PICP, and PINAFD for NNs with optimal sizes.

The NN size optimization is vital to avoid overfitting or underfitting [34].
The optimal neuron size of any NN-based prediction system depends on both the
data and the cost function. Therefore, size of NNs is optimized for each cost
function at first. Then, NNs of optimized sizes are trained to construct PIs
and to evaluate its performance for the test set.

###  5.1 The LUBE method

Single hidden layer NNs with different neuron numbers are trained to find an
optimal NN size. The NN is trained with four random initializations and one
with the lowest CWC is selected. Fig. 5(a) presents the lowest CWC values for
different NN sizes and three different PINC values for wind power generation
data. The number of neurons is varied between 5 and 15. The optimal NN size is
found to be 9, 8, and 7 for PINC = 95%, 90%, and 80% respectively. Similarly,
the optimum size of NNs is found to be 8, 8 and 6 for PINC = 95%, 90%, and 80%
respectively for the electricity demand data.

Khosravi et al. [24] previously observed that the training of the LUBE method
fails to converge once in twelve cases. Therefore, the training converges at
roughly 91.7% situations. 5.2% converged training provide a much wider or too
narrow PIs. Therefore, only acceptable PIs are considered. The performance of
the LUBE method is presented as the first segment of Table 1. Five hundred NNs
are trained for each of the wind power and the electricity demand data and NNs
providing logical PIs are considered. The LUBE method generates high-quality
PIs on average. However, some PIs exhibit correct PICP on cross-validation
data but provide slightly lower PICP compared to PINC for the test data. This
happens due to the slight variation between the test dataset and the cross-
validation dataset.

Figure 4: The structure of the NN with input-output combinations. Four recent
samples and the time is applied to quantify the uncertainty on the next
sample. Table 1: PI Optimization Performance for the 5-minutes ahead forecast
on the Wind Power Generation and Electricity Demand Data of the UK grid.

$LUBE$ Method
---
Data | $1-\alpha$ | $N_{Trial}$ | $\mu_{PINAW}$ | $\mu_{PICP}$ | $\sigma_{PICP}$ | $\mu_{PINAFD}$ | $\mu_{CWC}$ | $\mu_{CWFDC}$ | $\tilde{N}_{iter}^{(PICP(1\%))}$ | $\tilde{N}_{iter}^{(PINAW(1.5))}$
Wind Power | 0.95 | 500 | 5.42 | 95.14 | 1.42 | 1.23 | 5.42 | 8.25 | 91 | 126
0.90 | 500 | 4.61 | 90.12 | 2.61 | 2.08 | 4.62 | 13.09 | 113 | 131
0.80 | 500 | 3.23 | 80.15 | 2.32 | 2.74 | 3.24 | 28.47 | 119 | 137
E. De-mand | 0.95 | 500 | 2.07 | 95.06 | 0.70 | 0.51 | 2.09 | 3.11 | 78 | 104
0.90 | 500 | 2.01 | 90.09 | 1.26 | 0.74 | 2.05 | 4.08 | 92 | 119
0.80 | 500 | 1.90 | 80.14 | 1.53 | 0.77 | 1.96 | 8.27 | 99 | 127
Direct Interval Forecasting by Can Wan
Data | $1-\alpha$ | $N_{Trial}$ | $\mu_{PINAW}$ | $\mu_{PICP}$ | $\sigma_{PICP}$ | $\mu_{PINAFD}$ | $\mu_{CWC}$ | $\mu_{CWFDC}$ | $\tilde{N}_{iter}^{(PICP(1\%))}$ | $\tilde{N}_{iter}^{(PINAW(1.5))}$
Wind Power | 0.95 | 500 | 7.74 | 95.11 | 0.33 | 0.41 | 7.74 | 8.25 | 47 | 106
0.90 | 500 | 6.87 | 90.31 | 0.52 | 0.39 | 6.87 | 19.36 | 51 | 105
0.80 | 500 | 5.06 | 80.24 | 0.97 | 0.24 | 5.06 | 30.9 | 56 | 110
E. De-mand | 0.95 | 500 | 2.59 | 95.09 | 0.31 | 0.29 | 2.51 | 3.19 | 41 | 88
0.90 | 500 | 2.41 | 90.04 | 0.44 | 0.22 | 2.44 | 3.64 | 46 | 97
0.80 | 500 | 2.32 | 80.01 | 0.61 | 0.18 | 2.36 | 7.36 | 50 | 102
Mid-interval Deviation Consideration by L. G. Mar´ın
Data | $1-\alpha$ | $N_{Trial}$ | $\mu_{PINAW}$ | $\mu_{PICP}$ | $\sigma_{PICP}$ | $\mu_{PINAFD}$ | $\mu_{CWC}$ | $\mu_{CWFDC}$ | $\tilde{N}_{iter}^{(PICP(1\%))}$ | $\tilde{N}_{iter}^{(PINAW(1.5))}$
Wind Power | 0.95 | 500 | 9.02 | 95.21 | 1.42 | 1.09 | 9.02 | 9.21 | 68 | 98
0.90 | 500 | 7.36 | 90.22 | 1.51 | 1.56 | 7.36 | 9.32 | 64 | 103
0.80 | 500 | 6.51 | 80.30 | 2.03 | 1.83 | 6.51 | 18.34 | 64 | 107
E. De-mand | 0.95 | 500 | 3.03 | 95.16 | 0.72 | 0.39 | 2.89 | 3.33 | 42 | 81
0.90 | 500 | 2.78 | 90.19 | 1.30 | 0.34 | 2.65 | 3.49 | 45 | 94
0.80 | 500 | 2.72 | 80.25 | 1.68 | 0.33 | 2.40 | 3.81 | 48 | 97
Deviation Information-based Criterion by G. Zhang
Data | $1-\alpha$ | $N_{Trial}$ | $\mu_{PINAW}$ | $\mu_{PICP}$ | $\sigma_{PICP}$ | $\mu_{PINAFD}$ | $\mu_{CWC}$ | $\mu_{CWFDC}$ | $\tilde{N}_{iter}^{(PICP(1\%))}$ | $\tilde{N}_{iter}^{(PINAW(1.5))}$
Wind Power | 0.95 | 500 | 5.23 | 94.15 | 3.42 | 0.89 | 5.25 | 8.62 | 110 | 178
0.90 | 500 | 4.82 | 90.11 | 4.61 | 1.02 | 5.83 | 13.94 | 102 | 159
0.80 | 500 | 4.14 | 79.53 | 6.89 | 1.50 | 8.12 | 22.5 | 106 | 162
E. De-mand | 0.95 | 500 | 2.10 | 95.06 | 2.09 | 0.28 | 2.20 | 2.71 | 78 | 152
0.90 | 500 | 2.05 | 90.09 | 2.67 | 0.22 | 2.09 | 2.99 | 83 | 155
0.80 | 500 | 1.93 | 80.01 | 3.49 | 0.20 | 2.01 | 4.43 | 89 | 174
Proposed Method ($\rho$ =1, $\beta$ =1000, $\delta$ = $\alpha$/50)
Data | $1-\alpha$ | $N_{Trial}$ | $\mu_{PINAW}$ | $\mu_{PICP}$ | $\sigma_{PICP}$ | $\mu_{PINAFD}$ | $\mu_{CWC}$ | $\mu_{CWFDC}$ | $\tilde{N}_{iter}^{(PICP(1\%))}$ | $\tilde{N}_{iter}^{(PINAW(1.5))}$
Wind Power | 0.95 | 1000 | 7.40 | 95.09 | 0.11 | 0.53 | 7.40 | 8.03 | 45 | 92
0.90 | 1000 | 6.01 | 90.18 | 0.32 | 0.58 | 6.01 | 9.99 | 41 | 97
0.80 | 1000 | 4.31 | 80.41 | 0.38 | 0.67 | 4.31 | 10.08 | 42 | 88
E. De-mand | 0.95 | 1000 | 2.34 | 95.10 | 0.14 | 0.31 | 2.36 | 2.67 | 43 | 93
0.90 | 1000 | 2.29 | 90.21 | 0.41 | 0.23 | 2.35 | 3.02 | 40 | 90
0.80 | 1000 | 2.27 | 80.40 | 0.59 | 0.20 | 2.33 | 3.48 | 28 | 84

$\mu_{CWC}=\mu_{PINAW}+\gamma^{(\alpha,PICP)}(\mu_{PICP})$;
$\mu_{CWFDC}=\mu_{PINAW}+\mu_{PINAFD}+1000*(1-\alpha+\delta-\mu_{PICP})^{2}$
$\tilde{N}_{iter}^{(PICP(1\%))}$ is median of the number of iterations to
reach $|1-\alpha+\delta-PICP|<1\%$
$\tilde{N}_{iter}^{(PINAW(1.5))}$ is median of the number of iterations to
reach $PINAW<1.5\times PINAW_{Opt}$

Figure 5: NN-size optimization for the uncertainty quantification of the wind
power generation of the UK grid. (a) The LUBE cost function. Optimized NN-
sizes are 9 for $\alpha$=5%, 8 for $\alpha$=10% and 7 for $\alpha$=20%. (b) C.
Wan’s cost function. Optimized NN-sizes are 10 for $\alpha$=5%, 9 for
$\alpha$=10% and 9 for $\alpha$=20%. (c) L. G. Marn’s cost function. Optimized
NN-sizes are 11 for $\alpha$=5%, 10 for $\alpha$=10% and 10 for $\alpha$=20%.
(d) G Zhang’s cost function. Optimized NN-sizes are 10 for $\alpha$=5%, 9 for
$\alpha$=10% and 9 for $\alpha$=20%. (e) the proposed cost function. Optimized
NN-size is 10 for $\alpha$=5%, 10%, and 20%.

###  5.2 C. Wan’s method

Four neurons with different structures and random initialization are trained
and the one corresponding to the lowest $\gamma|ACE|+\lambda|S_{AV}|$ value is
selected. Fig. 5(b) presents those lowest cost function values. The number of
neurons is varied from 5 to 15. The optimal NN size is found to be 10, 9, and
9 for PINC = 95%, 90%, and 80% respectively for the wind power data.
Similarly, NN size is found to be 11, 8, and 7 for PINC = 95%, 90%, and 80%
respectively for the electricity demand data.

Five hundred NNs of optimal size are trained for both of the wind power and
the electricity demand data. The reported result considers NNs which provide
logical PIs. With the C. Wan’s cost function, roughly 99% NN training
converges and all converged NNs provide a logical PI. The performance of that
method is presented as the second segment of Table 1. C. Wan’s method provides
PIs of slightly higher width but reduces the average failure distance.

###  5.3 L. G. Marn’s method

L. G. Marn’s method considers the deviation from the mid-interval. This time,
NNs with the lowest $PINAW+||e||^{2}/n$ values are considered. Fig. 5(c)
presents the lowest cost function values. Following the same process, optimal
NN size is found to be 11, 10, and 10 for PINC = 95%, 90%, and 80%
respectively for the wind power data. NN size is found to be 11, 9, and 8 for
PINC = 95%, 90%, and 80% respectively for the electricity demand data.

As this cost function and its derivatives are continuous, the convergence is
higher (99.2%). The performance of the method is presented as the third
segment in Table 1. As a slightly higher PICP compared to PINC is also
penalized in this cost function, the average PICP becomes slightly greater
compared to other methods. As Marn’s method does not optimize the failure
distance, the failure distance is greater compared to C. Wan’s method. In the
contrary, consideration of the deviation from the mid-interval results in a
lower failure distance compared to the LUBE method. Also, optimization of the
deviation from the mid-interval brings the most probable regions of targets
close to the middle of PIs. That results in slightly wider PIs.

Figure 6:  Typical representative PICP values over iterations for the proposed
method and two existing methods. Only the first 100 iterations are presented.
These values vary from simulation to simulation. Figure 7:  PIs of 80%, 95%,
and 99% PICP with targets for (a) the electricity demand and (b) wind power
generation of UK grid on 17/7/2019. A 5-minutes ahead UQ is performed. Two
samples are drawn in each hour for better visualization. Intervals are
computed through an optimally trained NN. The structure of the NN is presented
in Fig. 4. Four previous samples and the time is used to compute the interval.
Figure 8:  PIs of 80%, 95%, and 99% PICP with targets for (a) the electricity
demand and (b) wind power generation of UK grid on 17/7/2019. A 30-minutes
ahead UQ is performed. Intervals are wider compared to the 5-minutes ahead UQ.

###  5.4 G. Zhang’s method

The same simulations are performed for NNs trained using G. Zhang’s cost
function. The lowest $DIC$ values are considered for drawing Fig. 5(d).
Optimal NN size is found to be 10, 9, and 9 for PINC = 95%, 90%, and 80%
respectively for the wind power data. Optimal NN size is 11, 9, and 8 for PINC
= 95%, 90%, and 80% respectively and for the electricity demand data.

Due to the discontinuity of the cost function at PICP=PINC, the convergence of
that simulation is low (96.2%). The performance of that method is presented as
the fourth segment in Table 1. This function also provides a good PICP on
average and a high failure distance. As the low PICP is poorly penalized and
the sum of failure distance and the average width may potentially result in a
wrong gradient direction, the variation in PICP is higher compared to other
methods.

###  5.5 The proposed method

NNs are trained using the proposed cost function as a part of the performance
evaluation. The lowest $PINAW+PINAFD$ values among four NNs are considered for
drawing Fig. 5(e). Optimal NN size is found to be 10 for all PINC values (95%,
90%, and 80%) for the wind power generation. Optimal NN size is found to be
12, 10, and 9 respectively for 95%, 90%, and 80% PINC for the electricity
demand data.

The failure distance resistance parameter ($\rho$) is set to one to provide an
equal concentration towards width and the failure distance. The PI coverage
penalty factor ($\beta$) is set to one thousand to reduce the PICP variation
among different NNs. According to observations, the PICP variation is higher
for a higher $\alpha$. Therefore, the coverage margin ($\delta$) is set to
($\alpha/50$). One thousand NNs of optimal size are trained and NNs providing
logical PIs are considered. 99.8% of NN training is converged and provides a
logical PI. The probable reason for 0.02% convergence failure can be PICP = 0
for the first thousand iterations. Although the NN training process changes
weights in each iteration, the PICP remains zero and the change in PINAW +
PINAFD is negligible compared to the large failure penalty at PICP = 0. Thus,
the gradient becomes small and causing slow convergence.

Although the simulation does not converge in 0.2% situations, the convergence
of the proposed system is much better compared to any existing method. The
simulation does not converge in 0.8% to 9% situations in currently available
methods. Fig. 6 presents PICP values over iterations. The value of PICP
reaches a feasible range (PINC - $\alpha$ $<$ PICP $<$ 1) within the first
thirty iterations most of the time in the proposed method. Usually, more than
fifty iterations are required to get such a feasible range in the initial LUBE
method. Other continuous cost functions also require more than thirty
iterations to get such a feasible range. Moreover, the PICP iterates near 100%
when a greater PICP (PICP $>$ PINC + $\delta$) is not penalized. Variation in
the finalized PICP is higher in LUBE and other continuous cost functions due
to discontinuity at the gradient near PICP = PINC. As the proposed cost
function is continuous along with its high order derivatives, the convergence
is faster and possesses a low variance.

The performance of this method is presented as the last segment of Table 1.
PICP variation becomes much lower with the proposed method when $\beta$ is set
to 1000. The mean value of PICP remains very close to PINC + $\delta$ and the
PICP variation is close to $\delta$. As a result, about 85% of NNs maintains
PICP $>$ PINC. On other systems, less than 70% of NNs maintains PICP $>$ PINC.
The proposed system has much lower PINAW + PINAFD as it is the performance
criterion during PICP = PINC + $\delta$. However, the user may choose any
desired failure distance resistance parameter ($\rho$) according to their
preference. Table 1 also present the convergence of the neural network
training. Both PICP and PINAW reaches near to their optimum value with fewer
iterations with the proposed cost function. A few high values can change the
average greatly, therefore we consider median values to compare convergences.

Fig. 7 presents PIs with targets for the 5-minutes-ahead prediction with the
proposed NN training. Fig. 7(a) presents PIs of three different PICPs with
targets for the electricity demand data. Fig. 7(b) presents PIs of three
different PICPs with targets for the wind power generation data. PIs and
corresponding targets are drawn on the same plot to visualize the PI of the
heteroscedastic system.

Fig. 8 presents PIs with targets for the 30-minutes-ahead prediction with the
proposed NN training. Fig. 8(a) presents PIs of three different PICPs with
targets for the electricity demand data. Fig. 8(b) presents PIs of three
different PICPs with targets for the wind power generation data.

##  6 Conclusion

Uncertainty is inescapable but uncertainty aware decisions bring higher
sustainability and profitability. The NN based LUBE PI construction method
achieved state of the art performance in quantifying asymmetrically
heteroscedastic uncertainty in terms of narrow width and required PICP.
However, NNs need to be re-trained for new types of signals and the non-smooth
LUBE cost function often fails to achieve an efficient uncertainty guided NN.
Several improvements to the LUBE method have proposed by different researchers
for optimal NN training. Each one of them has some limitations in terms of
convergence, understandability, smoothness, parameter insufficiency, or
customizability. Therefore, a smooth optimization function is proposed. A low
failure distance results in a low non-coverage penalty. Moreover, the user may
prefer to minimize the failure distance. Therefore, the cost function
considers coverage, width, and failure distance criteria to train NNs for the
construction of PIs with higher consistency. Researchers may bring 100%
convergence of the training in the future.

## References

  * [1] H. D. Kabir, A. Khosravi, M. A. Hosen, S. Nahavandi, Neural network-based uncertainty quantification: A survey of methodologies and applications, IEEE Access 6 (2018) 36218–36234.
  * [2] N. A. Shrivastava, A. Khosravi, B. K. Panigrahi, Prediction interval estimation of electricity prices using pso-tuned support vector machines, IEEE Transactions on Industrial Informatics 11 (2) (2015) 322–331.
  * [3] I. Koprinska, M. Rana, V. G. Agelidis, Correlation and instance based feature selection for electricity load forecasting, Knowledge-Based Systems 82 (2015) 29–40.
  * [4] Y.-Q. Zhang, X. Wan, Statistical fuzzy interval neural networks for currency exchange rate time series prediction, Applied Soft Computing 7 (4) (2007) 1149–1156.
  * [5] M. P. Clements, Evaluating the bank of england density forecasts of inflation, The Economic Journal 114 (498) (2004) 844–866.
  * [6] A. Kendall, Y. Gal, What uncertainties do we need in bayesian deep learning for computer vision?, in: Advances in Neural Information Processing Systems, 2017, pp. 5580–5590.
  * [7] Y. Gal, Uncertainty in deep learning, University of Cambridge.
  * [8] S. Roy, S. B. Roy, I. N. Kar, A new design methodology of adaptive sliding mode control for a class of nonlinear systems with state dependent uncertainty bound, in: 2018 15th International Workshop on Variable Structure Systems (VSS), IEEE, 2018, pp. 414–419.
  * [9] N. Laptev, J. Yosinski, L. E. Li, S. Smyl, Time-series extreme event forecasting with neural networks at uber, in: International Conference on Machine Learning, no. 34, 2017, pp. 1–5.
  * [10] H. Quan, D. Srinivasan, A. Khosravi, Uncertainty handling using neural network-based prediction intervals for electrical load forecasting, Energy 73 (2014) 916–925.
  * [11] H. Quan, D. Srinivasan, A. Khosravi, Integration of renewable generation uncertainties into stochastic unit commitment considering reserve and risk: A comparative study, Energy 103 (2016) 735–745.
  * [12] R. Ak, Y.-F. Li, V. Vitelli, E. Zio, Adequacy assessment of a wind-integrated system using neural network-based interval predictions of wind power generation and load, International Journal of Electrical Power & Energy Systems 95 (2018) 213–226.
  * [13] A. Altınten, F. Ketevanlioğlu, S. Erdoğan, H. Hapoğlu, M. Alpbaz, Self-tuning pid control of jacketed batch polystyrene reactor using genetic algorithm, Chemical Engineering Journal 138 (1) (2008) 490–497.
  * [14] H. D. Kabir, A. Khosravi, S. Nahavandi, Partial adversarial training for neural network-based uncertainty quantification, IEEE Transactions on Emerging Topics in Computational Intelligence.
  * [15] B. M. Ayyub, G. J. Klir, Uncertainty modeling and analysis in engineering and the sciences, Chapman & Hall/CRC Press, 2006.
  * [16] M. Alam, B. Grimm, J. P. Parmigiani, Effect of incident angle on crack propagation at interfaces, Engineering Fracture Mechanics 162 (2016) 155–163.
  * [17] R. W. Johnson, An introduction to the bootstrap, Teaching Statistics 23 (2) (2001) 49–54.
  * [18] B. R. Kirkwood, J. A. Sterne, Essential medical statistics, John Wiley & Sons, 2010\.
  * [19] P. Pinson, G. Kariniotakis, Conditional prediction intervals of wind power generation, IEEE Transactions on Power Systems 25 (4) (2010) 1845–1856.
  * [20] M. A. Carnero, A. Pérez, E. Ruiz, Identification of asymmetric conditional heteroscedasticity in the presence of outliers, SERIEs 7 (1) (2016) 179–201.
  * [21] A. Khosravi, S. Nahavandi, D. Creighton, Construction of optimal prediction intervals for load forecasting problems, Power Systems, IEEE Transactions on 25 (3) (2010) 1496–1503.
  * [22] P. Pinson, J. Tastu, Discussion of “prediction intervals for short-term wind farm generation forecasts” and “combined nonparametric prediction intervals for wind power generation”, IEEE Transactions on Sustainable Energy 5 (3) (2014) 1019–1020.
  * [23] C. Wan, Z. Xu, J. Østergaard, Z. Y. Dong, K. P. Wong, Discussion of “combined nonparametric prediction intervals for wind power generation”, IEEE Transactions on Sustainable Energy 5 (3) (2014) 1021–1021.
  * [24] A. Khosravi, S. Nahavandi, Closure to the discussion of prediction intervals for short-term wind farm generation forecasts and combined nonparametric prediction intervals for wind power generation and the discussion of combined nonparametric prediction intervals for wind power generation , Sustainable Energy, IEEE Transactions on 5 (3) (2014) 1022–1023.
  * [25] H. M. D. Kabir, A. Khosravi, A. Hosen, S. Nahavandi, Partial adversarial training for prediction interval, 2018 International Joint Conference on Neural Networks (IJCNN) (2018) 1–6.
  * [26] G. J. Coram, How to (and how not to) write a compact model in verilog-a, in: Behavioral Modeling and Simulation Conference, 2004. BMAS 2004. Proceedings of the 2004 IEEE International, IEEE, 2004, pp. 97–106.
  * [27] H. D. Kabir, Z. Ahmed, R. Kariyadan, L. Zhang, M. Chan, Modeling of fringe current for semiconductor-extended organic tfts, in: Electron Devices and Solid-State Circuits (EDSSC), 2016 IEEE International Conference on, IEEE, 2016, pp. 177–180.
  * [28] G. D. Magoulas, M. N. Vrahatis, G. S. Androulakis, Effective backpropagation training with variable stepsize, Neural networks 10 (1) (1997) 69–82.
  * [29] G. D. Magoulas, M. N. Vrahatis, G. S. Androulakis, Improving the convergence of the backpropagation algorithm using learning rate adaptation methods, Neural Computation 11 (7) (1999) 1769–1796.
  * [30] C. Wan, Z. Xu, P. Pinson, Direct interval forecasting of wind power, IEEE Transactions on Power Systems 28 (4) (2013) 4877–4878.
  * [31] L. G. Marín, F. Valencia, D. Sáez, Prediction interval based on type-2 fuzzy systems for wind power generation and loads in microgrid control design, in: Fuzzy Systems (FUZZ-IEEE), 2016 IEEE International Conference on, IEEE, 2016, pp. 328–335.
  * [32] G. Zhang, Y. Wu, K. P. Wong, Z. Xu, Z. Y. Dong, H. H.-C. Iu, An advanced approach for construction of optimal wind power prediction intervals, IEEE Transactions on Power Systems 30 (5) (2015) 2706–2715.
  * [33] Gridwatch, Uk national grid status (Oct. 2017).
URL http://www.gridwatch.templar.co.uk/

  * [34] D. M. Blei, A. Y. Ng, M. I. Jordan, Latent dirichlet allocation, Journal of machine Learning research 3 (Jan) (2003) 993–1022.