# Optimal Adaptive Prediction Intervals for Electricity Load Forecasting in
Distribution Systems via Reinforcement Learning

Yufan Zhang, Student Member, IEEE, Honglin Wen, Student Member, IEEE,
Qiuwei Wu, Senior Member, IEEE, and Qian Ai, Senior Member, IEEE Yufan Zhang,
Honglin Wen, and Qian Ai are with Department of Electrical Engineering,
Shanghai Jiao Tong University, Shanghai 200240, China.
Qiuwei Wu is with Tsinghua-Berkeley Shenzhen Institute,Tsinghua Shenzhen
International Graduate School, Tsinghua University, Shenzhen 518055, China.
Corresponding author: Qiuwei Wu (e-mail: qiuwu@sz.tsinghua.edu.cn).

###### Abstract

Prediction intervals (PIs) offer an effective tool for quantifying uncertainty
of loads in distribution systems. The traditional central PIs cannot adapt
well to skewed distributions, and their offline training fashion is vulnerable
to the unforeseen change in future load patterns. Therefore, we propose an
optimal PI estimation approach, which is online and adaptive to different data
distributions by adaptively determining symmetric or asymmetric probability
proportion pairs for quantiles of PIs’ bounds. It relies on the online
learning ability of reinforcement learning (RL) to integrate the two online
tasks, i.e., the adaptive selection of probability proportion pairs and
quantile predictions, both of which are modeled by neural networks. As such,
the quality of quantiles-formed PI can guide the selection process of optimal
probability proportion pairs, which forms a closed loop to improve PIs’
quality. Furthermore, to improve the learning efficiency of quantile
forecasts, a prioritized experience replay (PER) strategy is proposed for
online quantile regression processes. Case studies on both load and net load
demonstrate that the proposed method can better adapt to data distribution
compared with online central PIs method. Compared with offline-trained
methods, it obtains PIs with better quality and is more robust against concept
drift.

Keywords: Prediction interval; Online load forecasting; Reinforcement
learning; Uncertainty.

##  I Introduction

Load forecasting plays a pivotal role in distribution systems operation, such
as congestion management [1] and energy transaction [2]. However, it is widely
acknowledged that accurate short-term load forecast (STLF) in distribution
systems is challenging due to the uncertain characteristics of loads.
Moreover, with the increasing penetration of renewable energy resources (RESs)
in distribution systems, load measurements are blended with generation of RESs
(and thus referred to as net loads), which further complicates the forecasting
task. Although efforts have been devoted to accurate forecasts by focusing on
both features [3] and algorithms [4], point forecasts still cannot accommodate
the inherent uncertainty of loads. The recent decade has witnessed the trend
in load forecasting community, i.e., transition from point forecasts to
probabilistic forecasts; see recent review [5].

In general, probabilistic short-term load forecasting (PSTLF) is to
communicate the probability of load at future time, given information up to
current time [6]. It is usually developed in forms of quantile [7], prediction
interval (PI) [8], and probability density [9]. A PI is specified by two
bounds and a nominal coverage probability (NCP) $(1-\beta)\times 100\%$ ,
which indicates the probability that future observations fall in the interval.
It has been successfully employed in the decision-making of distribution
systems under uncertainty. For instance, PI delivers an uncertainty set for
robust optimization (RO), such that RO produces a decision reliable for any
realization in the set. In [10], multiple PIs with different NCP levels were
formulated, so that RO was performed for the multiple PIs-constrained dispatch
model, which reduced the conservatism of RO. Besides, PI offers information
about uncertainty for interval optimization. In [11], PI was used to describe
prediction uncertainty for the economic dispatch of virtual power plant, which
was then converted into a deterministic problem.

Traditionally, parametric approach with distributional assumption [12],
residual resampling [13], and quantile regression (QR) [14] are three popular
approaches to constructing PIs. The parametric approach relies on
distributional assumption, the shape parameters of which are estimated by
data-driven methods. The delta, Bayesian, and mean-variance methods have been
proposed [12], which assume that forecast errors have a Gaussian distribution.
As residual resampling and QR are distribution-free, they are also referred to
nonparametric approaches, whose successful applications in PI estimation can
be traced back to [14]. Resampling methods construct PI through sampling from
the empirical distribution of point forecast errors. Ref. [13] proposed a two-
stage bootstrap sampling method for PSTLF. The first stage used bootstrap to
characterize uncertainties from forecast models, and the second stage applied
bootstrap to resample from the regression errors. Although the resampling
methods are distribution-free, their performances are restricted by the point
forecast methods relying on. QR directly estimates quantiles for loads, and
uses two corresponding quantiles as the lower and upper bounds of PI. Ref.
[15] proposed a cost-oriented PI method, where the quantile proportion pair of
PI was selected to minimize operational cost. It is easy to combine QR with
machine learning methods, for instance, extreme learning machine (ELM) [16],
gradient boosting regression tree [17], and neural basis expansion model [8].
A pinball-loss guided QR for load responsiveness estimation was proposed in
[18], which was solved in a fully-distributed manner. In [19], an ensemble
method for combining QR models was proposed, which was designed as a linear
programming problem that minimized the pinball loss. Although QR method is
distribution-free and easy to use, it faces the problem, i.e., how to
determine the probability proportions for PI’s lower and upper bounds. A
typical approach is using the symmetric probability proportions. Such PI
formed by $q^{\beta/2}$ and $q^{1-\beta/2}$ is referred to as a central
prediction interval (CPI). Although these works have enriched tools for PSTLF
and offered complementary understanding, recent studies reveal that CPIs are
usually unnecessarily wide [20].

Therefore, increasing works have advocated the optimal prediction interval
(OPI) that aims to minimize PI width under the constraint of coverage. One of
the most popular approaches is to estimate the PI’s bounds directly, with the
objective accounting for both reliability and sharpness, which is therefore
free of quantile proportions selection. Ref. [21] developed a neural network
(NN) model to output the bounds. The NN was trained by minimizing the PI-based
objective function. However, for the NN-based method, the PI-based objective
function needs to be carefully designed, such that it is differentiable for NN
to learn. To solve this problem, [22] proposed a relaxed surrogate to the
original non-differentiable objective, making gradient descent a viable
solution method. However, the loss functions in these two works are not proper
metrics in the view of statistics, which translates into saying that the
quality of the forecasts issued by these two methods may be a concern in the
term of proper statistical scores. Ref. [23] proposed a hybrid method of OPI
based on ELM and particle swarm optimization. The multi-objective optimization
for the ELM’s output weights determination was developed. Although there is no
differentiable requirement for the PI-based objective by using particle swarm
optimization, the multi-objective optimization involving the soft constraints
is hard to obtain the Pareto optimal, which hinders accurate estimation of the
PI. Another approach is to optimize the selection of quantiles’ probability
proportions. For instance, [20] employed an ELM to generate nonlinear
features, and mixed-integer linear programming (MILP) was formulated to obtain
the OPI. The ELM does not involve backpropagations, which means the feature
learning stage is inefficient and leads to the obstacles in the later
estimation stage. Therefore, it remains an open issue to develop OPI.

Another limitation of the aforementioned works stem from offline-trained
algorithms. Their frameworks contain offline training and subsequent online
forecasting. However, offline-trained algorithms rely on the common assumption
that training and test data are drawn from the same distribution, also
referred to as independently identically distribution. In reality, concept
drift often occurs [24]. In other words, the distribution of test data differs
from that of training data. Usually, offline models cannot adapt well to the
changing distribution. Moreover, to retrain offline models, it usually
requires abundant recollected samples, and takes expensive computation.
Indeed, online learning is an appealing solution to tackle this issue [25]. In
this context, the model iteratively predicts a value, and then receives a cost
according to an adversarially chosen function. The parameters of the model are
updated in an online fashion, with only newly observed sample or a batch of
samples. An online ensemble learning approach for load forecasting was
proposed in [26]. Using the passive aggressive regression method, the convex
optimization problem was formulated for the new observation to update the
ensemble weights. However, as this approach is grounded in convex
optimization, it is not suitable for the parameter estimation of nonconvex
models. Using recurrent neural network (RNN) for online prediction, [27]
updated RNN’s parameters by a part of data batch sampled from the buffer.
However, since the samples stored in the buffer are assumed to have the equal
significance for model updating, their different contributions to improving
the forecast performance can hardly be distinguished.

Indeed, the OPI problem can be decoupled into two correlated sub-problems,
i.e., probability proportion selections and corresponding quantile forecasts
through QR models fitting. As such, the tasks of proportion selections and QR
fitting can be fulfilled separately, which frees the limitation imposed on QR
models in [20]. And the task of proportion selections is nontrivial and cannot
be searched in an exhaustive manner, as probability proportions are
continuous. Instead, it can be modeled as a decision-making problem that an
agent decides the proportions given some contexts, which can be further
naturally considered in an online setting. In this paper, the OPI problem at
each time-step is formulated as a bilevel program, where the upper level
performs the quantile proportion selection, and given the selected proportion,
the lower level predicts quantiles for PI construction. Reinforcement learning
(RL) is leveraged to solve the bilevel program in an online fashion, which
integrates upper and lower decision tasks in a closed-loop manner.
Specifically, the upper-level task is modeled by the decision making of the
agent, whereas the lower-level task is solved in the environment. RL has
strong theoretical background rooted in Markov Decision Process (MDP). And the
bounded reward value guarantees the online iterations converging to the
optimal value [28]. Recent years have witnessed its successful applications,
such as in the online energy management [29, 30]. Therefore, this paper
proposes an online OPI forecasting method that formulates PIs in an adaptive
manner. In general, it is a closed-loop feedback framework working toward the
quality improvement of estimated PIs. Concretely, the closed-loop framework is
established based on the value-based RL algorithm, where the agent is modeled
by a neural network. At each time-step, the agent receives an input state, and
adaptively determines the symmetric or asymmetric quantile proportion pairs
for PI as corresponding action. Based on the action, the environment estimates
PI’s bounds via the online QR models. Then, the parameters of online QR models
are updated through gradient descent optimization by a batch of data from the
buffer. In particular, an online quantile forecasting method with prioritized
experience replay (PER) [31] is proposed. Unlike treating samples in the
buffer equally, here the data is assigned with different probability of being
sampled to update the model parameters, which is derived by their importance
for improving the QR’s performance. After the quantile predictions, the
overall quality of the formed PI is used as the feedback to guide the agent in
return. Compared with the existing studies, the main contributions of the
paper are summarized as follows:

(1) The proposal of a new approach for OPI, which is adaptive to the
underlying data distribution, and uncovers the relationship between optimal
quantile proportion selection and quantile forecasting for OPI estimation;

(2) A RL-based online OPI prediction framework, which integrates the quantile
proportion selection and quantile estimations for PI’s lower and upper bounds
(LUB). As such, the two tasks work in a closed-loop manner toward PI quality
improvement. Also, under such a framework, a large spectrum of QR models are
applicable to the quantile forecasting task;

(3) An online quantile forecasting approach with the PER strategy, which
improves the data efficiency and QR models’ prediction quality by learning
more frequently on the prioritized samples.

The remainder of this paper is organized as follows. Section II introduces the
definition of the proposed problem and the overall framework. Section III
proposes the online OPI algorithm. Results are discussed and evaluated in
Section IV followed by the conclusions.

##  II Problem Statement

In this section, we formulate the OPI problem mathematically in subsection _A_
, and present the overall framework in subsection B.

###  II-A Problem Definition

Let $y_{t}$ denote the load at time $t$, which is a realization of the
corresponding random variable $Y_{t}$. The probabilistic load forecasting with
lead time $k$ translates to communicate the density function
$f_{t+k|t}(Y_{t+k}|\Omega_{t})$ and cumulative distribution function
$F_{t+k|t}(Y_{t+k}|\Omega_{t})$, where
$\Omega_{t}=\left\\{Y_{t+h-1},...,Y_{t}\right\\}$ is the information set up to
time $t$. Specifically, a quantile $q_{t+k}^{\alpha_{t+k}}$ is defined as
$F_{t+k|t}^{-1}(\alpha_{t+k})$ where the quantile’s probability proportion
$\alpha_{t+k}\in\left[0,1\right]$. Therefore, a PI with NCP $(1-\beta)\times
100\%$ can be developed as:

$\left[\hat{q}_{t+k}^{\alpha_{t+k}},\hat{q}_{t+k}^{\alpha_{t+k}+1-\beta}\right].$
$None$

Here the quantiles are the outputs of QR models

$\begin{split}&\hat{q}_{t+k}^{\alpha_{t+k}}=f^{\alpha_{t+k}}(\bm{s}_{t};\hat{\bm{W}}^{\alpha_{t+k}}_{t})\\\
&\hat{q}_{t+k}^{\alpha_{t+k}+1-\beta}=f^{\alpha_{t+k}+1-\beta}(\bm{s}_{t};\hat{\bm{W}}^{\alpha_{t+k}+1-\beta}_{t}),\end{split}$
(2)

where the vector $\bm{s}_{t}$ is the realization of $\Omega_{t}$ which is
composed of the $h$ previous values, and $\hat{\bm{W}}^{\alpha_{t+k}}_{t}$,
$\hat{\bm{W}}^{\alpha_{t+k}+1-\beta}_{t}$ are QR models’ parameters for the
estimations of LUB’s quantiles. The estimated parameters of QR models can be
obtained by solving the optimization problem of minimizing a pinball loss

$\begin{split}&\hat{\bm{W}}_{t}^{\alpha_{t+k}}=\mathop{\arg\min}_{\bm{W}_{t}^{\alpha_{t+k}}}\mathbb{E}_{F_{t+k|t}}L^{\alpha_{t+k}}(Y_{t+k},q_{t+k}^{\alpha_{t+k}})\\\
&\qquad\qquad
s.t.q_{t+k}^{\alpha_{t+k}}=f^{\alpha_{t+k}}(\bm{s}_{t};\bm{W}^{\alpha_{t+k}}_{t})\\\
&\hat{\bm{W}}_{t}^{\alpha_{t+k}+1-\beta}=\mathop{\arg\min}_{\bm{W}_{t}^{\alpha_{t+k}+1-\beta}}\mathbb{E}_{F_{t+k|t}}L^{\alpha_{t+k}+1-\beta}(Y_{t+k},q_{t+k}^{\alpha_{t+k}+1-\beta})\\\
&\qquad\qquad
s.t.q_{t+k}^{\alpha_{t+k}+1-\beta}=f^{\alpha_{t+k}+1-\beta}(\bm{s}_{t};\bm{W}^{\alpha_{t+k}+1-\beta}_{t}),\end{split}$
(3)

where $L^{\alpha}$, $\alpha\in\\{\alpha_{t+k},\alpha_{t+k}+1-\beta\\}$ is the
pinball loss, which is defined as

$\begin{split}L^{\alpha}(y_{t+k},q_{t+k}^{\alpha})=(1-\alpha)(q_{t+k}^{\alpha}-y_{t+k})I_{\left\\{q_{t+k}^{\alpha}-y_{t+k}\geq
0\right\\}}\\\
+\alpha(y_{t+k}-q_{t+k}^{\alpha})I_{\left\\{q_{t+k}^{\alpha}-y_{t+k}\textless
0\right\\}},\end{split}$ (4)

where $I_{\left\\{.\right\\}}$ is the indicator function.

For symmetric distributions like Gaussian distribution, $\alpha_{t+k}$ is
often determined as $\beta/2$, for any time indiscriminately. However, for
asymmetric distributions such as log-normal distribution, Beta distribution,
such CPIs cannot adapt well to the real distribution. Indeed, conditioned on a
certain data distribution $F_{t+k|t}(Y_{t+k}|\Omega_{t})$, there exists an OPI
$\left[\hat{q}_{t+k}^{\alpha_{t+k}^{*}},\hat{q}_{t+k}^{\alpha_{t+k}^{*}+1-\beta}\right]$
with optimal proportions $\alpha_{t+k}^{*},\alpha_{t+k}^{*}+1-\beta$, whose
width is smaller than that of CPI with the same NCP. The illustration of the
differences between CPI and OPI under a skewed probability distribution is
shown in Fig. 1. The shaded parts have the same area under the density curve.
Since the probability density function is with positive skewness, the width of
OPI is smaller than that of CPI.

Figure 1: Illustration of CPI and OPI under a skewed probability density
function.

In this paper, instead of treating the quantile proportion as the
predetermined constant, the problem of OPIs formulation involves selection of
quantile proportions and corresponding quantile estimations. That is, we wish
to find the optimal proportion pair at time $t$, such that the Winkler score
$w(y_{t+k},\hat{q}_{t+k}^{\alpha_{t+k}^{*}},\hat{q}_{t+k}^{\alpha_{t+k}^{*}+1-\beta})$
that measures the overall quality of PI is minimized. The selected proportion
pair at time $t$ is then used to specify corresponding QR models. Therefore,
the OPI problem at time $t$ can be formulated as a bilevel programming
problem,

$\begin{split}&\mathop{\min}_{\bm{\theta}_{t}}w(y_{t+k},\hat{q}_{t+k}^{\alpha_{t+k}^{*}},\hat{q}_{t+k}^{\alpha_{t+k}^{*}+1-\beta})\\\
&s.t.\\\ &\alpha_{t+k}^{*}=f(\bm{s}_{t};\bm{\theta}_{t})\\\
&\hat{q}_{t+k}^{\alpha_{t+k}^{*}}=f^{\alpha_{t+k}^{*}}(\bm{s}_{t};\hat{\bm{W}}^{\alpha_{t+k}^{*}}_{t})\\\
&\hat{q}_{t+k}^{\alpha_{t+k}^{*}+1-\beta}=f^{\alpha_{t+k}^{*}+1-\beta}(\bm{s}_{t};\hat{\bm{W}}^{\alpha_{t+k}^{*}+1-\beta}_{t})\\\
&\hat{\bm{W}}_{t}^{\alpha_{t+k}^{*}}=\mathop{\arg\min}_{\bm{W}_{t}^{\alpha_{t+k}^{*}}}\mathbb{E}_{F_{t+k|t}}L^{\alpha_{t+k}^{*}}(Y_{t+k},q_{t+k}^{\alpha_{t+k}^{*}})\\\
&\qquad\qquad
s.t.q_{t+k}^{\alpha_{t+k}^{*}}=f^{\alpha_{t+k}^{*}}(\bm{s}_{t};\bm{W}^{\alpha_{t+k}^{*}}_{t})\\\
&\hat{\bm{W}}_{t}^{\alpha_{t+k}^{*}+1-\beta}=\mathop{\arg\min}_{\bm{W}_{t}^{\alpha_{t+k}^{*}+1-\beta}}\mathbb{E}_{F_{t+k|t}}L^{\alpha_{t+k}^{*}+1-\beta}(Y_{t+k},q_{t+k}^{\alpha_{t+k}^{*}+1-\beta})\\\
&\qquad\qquad
s.t.q_{t+k}^{\alpha_{t+k}^{*}+1-\beta}=f^{\alpha_{t+k}^{*}+1-\beta}(\bm{s}_{t};\bm{W}^{\alpha_{t+k}^{*}+1-\beta}_{t}).\end{split}$
(5)

As each two adjacent states are auto-correlated, the state at time $t$
transits to the state at time $t+1$ according to the state transition function
$S^{M}(.)$, i.e., $\bm{s}_{t+1}=S^{M}(\bm{s}_{t})$. It involves a similar
problem that minimizes the Winkler score at the time $t+1$.

In (5), the upper level performs the probability proportion selection by
learning a mapping $f(\bm{s}_{t};\bm{\theta}_{t})$ with parameters
$\bm{\theta}_{t}$. The lower-level problem performs the quantile forecasting.
In general, the two problems in two levels are interdependent. Usually, the QR
model is highly nonconvex and the bilevel problem in (5) is difficult to
solve. However, in some special case where the lower-level QR model is linear,
the bilevel problem can be transformed to a single level mathematical
programming problem. For the general problem, we propose to solve it in an
iterative manner via RL framework. Also, as random variable’s cumulative
distribution function $F_{t+k|t}$ changes with on-going input data, the upper-
and lower-level optimization problem changes from time to time. Therefore, the
selection of optimal proportion should be updated with time. Besides, the
estimation of QR models should also be updated with time, due to the
nonstationary characteristics of the time series. Therefore, at time $t$, we
optimize over the probability proportion and corresponding quantiles.

###  II-B The Proposed Framework

In this paper, RL serves as the basic closed-loop framework for integrating
the aforenamed tasks. Specifically, to issue the quantile forecasts for load
at time $t$ , one uses previous values $\left\\{y_{t-h},...,y_{t-1}\right\\}$
as input features and estimates $\hat{q}_{t-1+k}^{\alpha_{t-1+k}^{*}}$,
$\hat{q}_{t-1+k}^{\alpha_{t-1+k}^{*}+1-\beta}$ via parameters
$\hat{\bm{W}}^{\alpha_{t-1+k}^{*}}_{t-1}$,
$\hat{\bm{W}}^{\alpha_{t-1+k}^{*}+1-\beta}_{t-1}$. Then after $y_{t}$ is
revealed, the parameter of the corresponding QR models
$f^{\alpha_{t-1+k}^{*}}(.)$, $f^{\alpha_{t-1+k}^{*}+1-\beta}(.)$ are updated
by gradient descent at the time $t$ as:
$\hat{\bm{W}}^{\alpha_{t-1+k}^{*}}_{t}$,
$\hat{\bm{W}}^{\alpha_{t-1+k}^{*}+1-\beta}_{t}$, while the parameters of the
remaining QR models at the time $t$ remain the same as that of their
parameters at the time $t-1$. Then, after $\alpha_{t+k}^{*}$ at the time $t$
is revealed, quantiles $\hat{q}_{t+k}^{\alpha_{t+k}^{*}}$,
$\hat{q}_{t+k}^{\alpha_{t+k}^{*}+1-\beta}$ are estimated by
$f^{\alpha_{t+k}^{*}}(.;\hat{\bm{W}}^{\alpha_{t+k}^{*}}_{t})$,
$f^{\alpha_{t+k}^{*}+1-\beta}(.;\hat{\bm{W}}^{\alpha_{t+k}^{*}+1-\beta}_{t})$.
For simplicity of notations, $\left\\{y_{t-h},...,y_{t-1}\right\\}$ is denoted
as $\bm{x}_{t}$. The whole framework is illustrated in Fig. 2.

Figure 2: Illustration of proposed RL-based OPI framework.

As shown, an agent, modeled by NN, is responsible for selecting the optimal
quantile proportion $\alpha_{t-1+k}^{*}$ for lower bound’s quantile under the
current state/context $\bm{x}_{t}$. As the probability proportion is
continuous in the range of $(0,\beta)$ , the optimal proportion cannot be
found by exhaustive enumeration, which leads to the adoption of RL method. And
corresponding QR models with different quantile proportions are in the
environment for predicting the LUB’s quantiles for PIs. Specifically, the
states of the agent are also the input features fed into QR models. And the
quality of the constructed PIs is used as the feedback reward to link the two
tasks. Since we focus on online estimation, the lead time $k$ equals one in
the following discussion. We detail the RL formulation of the online OPIs
construction problem. The key elements of which are defined as follows:

State: For state at the time $t$, the QR’s input feature $\bm{x}_{t}$ is also
used as the state for the RL’s agent.

State transition function: To issue the quantile forecasts for load at time
$t$, $\bm{x}_{t}=\\{y_{t-h},...,y_{t-1}\\}$ is the state of the agent as well
as input features of the QR model. Likewise, to issue the quantile forecasts
for load at time $t+1$, the state of the agent is
$\bm{x}_{t+1}=\\{y_{t-h+1},...,y_{t}\\}$. The state transition from
$\bm{x}_{t}$ to $\bm{x}_{t+1}$ in this problem is a mapping of continuous
variables and independent of the action. Since the state transition of high-
dimensional continuous states is difficult to model, it is hard to solve the
formulated DRL problem via a model-based RL method that involves learning a
dynamic model. Therefore, we decide to model the abovementioned RL problem in
a model-free approach, rather than a traditional model-based RL approach.

Action: To choose optimal probability proportion from the continuous range of
$(0,\beta)$, the probability proportion is transformed into discrete
variables. Concretely, the discrete action space $\mathcal{A}$ is a set of
probability proportions with the size of $|\mathcal{A}|$, which is equal to
the number of quantile proportions. With the increasing number of actions, the
range of $(0,\beta)$ is modeled more precisely for the probability proportion
of lower bound’s quantile. Given the certain NCP, the action space is:

$\left\\{\frac{i\cdot\beta}{|\mathcal{A}|+1}\right\\}^{|\mathcal{A}|}_{i=1},$
$None$

where $|\mathcal{A}|=2^{n}-1,n=1,2,...$ is the number of actions. So, at the
time $t$, the optimal $\alpha^{*}_{t}$ is chosen from the set defined in (6).
When $n=1$, the problem returns back to CPI. The illustration of the action
space with the change of action numbers is shown in Fig. 3. The action space
with the larger dimension contains all the actions in the previous action
space with the smaller dimension. Although the number of actions has impact on
the agent’s policy, as the difference between adjacent quantile proportions
under relatively large action size is small, the impact of the action space’s
size on the selection policy is marginal.

Remark 1: In the general solution of the model in (5), the optimal solution of
the upper-level problem in (5) can be continuous probability proportion.
However, we choose to approximate the continuous proportion by discrete action
space out of two practical concerns:

Firstly, as we train QR models separately for each quantile proportion in the
action space, with finite discrete probability proportions, the number of QR
models is also finite. So, it is only required to optimize the parameters of
limited number of QR models. In contrast, if probability proportions are
modeled in a continuous action space, the number of QR models will be
infinite, which considerably increases the computational burden. Secondly,
since the adjacent probability proportions in the action space are close to
each other, the discrete probability proportions are adequate for representing
the continuous probability proportions in a range.

Figure 3: Illustration of action space.

Reward: The Winkler score is a comprehensive evaluation index of PI quality,
which measures both reliability and sharpness simultaneously. The smaller the
Winkler score, the better the quality of the constructed PI. It is defined as:

$w_{t}=\left\\{\begin{aligned}
&\hat{q}_{t}^{\alpha_{t}^{*}+1-\beta}-\hat{q}_{t}^{\alpha_{t}^{*}},\qquad\quad\quad\qquad\hat{q}_{t}^{\alpha_{t}^{*}}\leq
y_{t}\leq\hat{q}_{t}^{\alpha_{t}^{*}+1-\beta}\\\
&\hat{q}_{t}^{\alpha_{t}^{*}+1-\beta}-\hat{q}_{t}^{\alpha_{t}^{*}}+2(\hat{q}_{t}^{\alpha_{t}^{*}}-y_{t})/\beta,\quad\quad\quad\quad
y_{t}<\hat{q}_{t}^{\alpha_{t}^{*}}\\\
&\hat{q}_{t}^{\alpha_{t}^{*}+1-\beta}-\hat{q}_{t}^{\alpha_{t}^{*}}+2(y_{t}-\hat{q}_{t}^{\alpha_{t}^{*}+1-\beta})/\beta,\hat{q}_{t}^{\alpha_{t}^{*}+1-\beta}<y_{t}\end{aligned}\right.$
$None$

Therefore, the negative value of Winkler score for the constructed PI at time-
step $t$ is used as the reward:

$r_{t}=-w_{t}.$ $None$

Agent: The optimal probability proportion selection problem is solved by
agent, which gradually learns how to choose the optimal probability proportion
for the lower bound quantile to construct good quality PI.

Environment: The LUB’ quantiles forecasting problem is in the environment,
which is consisted of $2\cdot|\mathcal{A}|$ quantile predictors which are QR
models. The online forecasting process is illustrated in the next Section and
each predictor has its own model parameters and buffer.

##  III Online Optimal Construction of PI

In this Section, a general form of online quantile forecasting with PER
strategy is proposed in subsection A. Based on it, the proposed online OPI
method is presented in subsection B.

###  III-A Online Quantile Forecasting with PER Strategy

Contrary to generating the best quantile predictor on the entire training set
at once, online QR updates the predictor for future load at each step.
Therefore, in the case of the $\alpha_{th}$ quantile online forecasting, the
QR model $f^{\alpha}(.;\bm{W}_{t}^{\alpha})$ with parameters
$\bm{W}_{t}^{\alpha}$ is updated in a recursive and adaptive manner for data
in a sequential order. Here, we update the model using gradient descent
recursively by a small batch of data. Specifically, the new data revealed is
stored into the predictor’s buffer $D_{t}^{\alpha}$, and a small batch of
feature and label pairs, which are referred as experiences, are sampled from
the buffer to update the model parameters $\bm{W}_{t}^{\alpha}$. It has
several advantages over the standard online stochastic gradient descent.
First, the data efficiency is improved, because the historical data pairs can
be potentially used in many updates. Second, since the parameter updating is
based on a small batch of data not a single data, the parameters learning
process can be more stable.

Here, instead of sampling data from the buffer uniformly at random, the online
quantile forecasting with PER strategy is proposed, where the data are sampled
with priority. The priority of sampling experience $e_{j}=(\bm{x}_{j},y_{j})$
is defined as:

$p^{\alpha}_{j}=L^{\alpha}(y_{j},\hat{q}_{j}^{\alpha}),$ $None$

where
$\hat{q}_{j}^{\alpha}=f^{\alpha}(\bm{x}_{j};\hat{\bm{W}}_{t-1}^{\alpha})$ is
the estimated quantile and $y_{j}$ is the realization of the load.
$L^{\alpha}(y_{j},\hat{q}_{j}^{\alpha})$ is the pinball loss of the $e_{j}$,
and $p_{j}^{\alpha}$ is the priority.

And for the new data, its priority is assigned with the maximum priority in
the predictor’s buffer, thus favoring the new data during sampling later.

The probability of the experience $e_{j}$ being sampled is defined based on
the priority as follows:

$P_{j}=(p_{j}^{\alpha})^{\sigma}\big{/}\sum_{k}(p_{k}^{\alpha})^{\sigma},$
$None$

where $\sigma$ controls the extent that the prioritization is used. When
$\sigma=0$, it corresponds to the case of uniform sampling. Therefore, by
assigning data with probability proportional to the priority, the experience
with large priority is more likely to be drawn frequently for parameter
updating. Therefore, with PER strategy, the sample with larger pinball loss
and the new sample are more likely to be sampled for parameters updating. That
is, the samples that are harder to learn and that bring new information are
learned more frequently.

With PER strategy, the loss function $\ell(\hat{\bm{W}}^{\alpha}_{t-1})$ of a
small batch of experiences with the size $B$ can be written as follows:

$\ell(\hat{\bm{W}}^{\alpha}_{t-1})=\frac{1}{B}\sum_{j}w_{j}(p_{j}^{\alpha})^{2},$
$None$ $w_{j}=(N\cdot P_{j})^{-\rho}\big{/}w^{max},$ $None$

where $w_{j}$ is the normalized importance-sampling (IS) weight, which
corrects the bias introduced by prioritized sampling [31]. Through the
correction, (12) shows that the experience with larger probability being
sampled is assigned with the smaller weight. $N$ is the number of experiences
in the predictor’s buffer, and $w^{max}$ is the maximum IS weight. The
parameter $\rho$ controls the extent of importance-sampling correction, which
fully compensates for the priority sampling when $\rho=1$. When $\rho=0$,
there is no correction.

Then, the predictor’s parameters update is as follows:

$\hat{\bm{W}}^{\alpha}_{t}=\hat{\bm{W}}^{\alpha}_{t-1}-\eta_{\alpha}\cdot\nabla_{\hat{\bm{W}}^{\alpha}_{t-1}}\ell(\hat{\bm{W}}^{\alpha}_{t-1}).$
$None$

Here, $\eta_{\alpha}$ is the learning rate of the QR model. The final online
quantile forecasting method with PER strategy is outlined in Algorithm 1, and
the output of Algorithm 1 is the predicted quantile at each time-step. The
additional remarks are presented below.

Algorithm 1 Online Quantile Forecasting with PER Strategy

1:Batch size $B$, learning rate $\eta_{\alpha}$, parameters $\sigma$,$\rho$,
the quantile’s probability proportion $\alpha$

2:Initialize the quantile predictor $f^{\alpha}(.;\hat{\bm{W}}_{0}^{\alpha})$
with random weights $\hat{\bm{W}}_{0}^{\alpha}$

3:for $t$=1,2,… do

4: Given the input feature $\bm{x}_{t}$, predict the quantile
$\hat{q}_{t}^{\alpha}=f^{\alpha}(\bm{x}_{t};\hat{\bm{W}}_{t-1}^{\alpha})$

5: Reveal true value $y_{t}$

6: Store the experience $e_{t}=(\bm{x}_{t},y_{t})$ in the predictor’s buffer
$D_{t}^{\alpha}$ with maximal priority
$p_{t}^{\alpha}=max_{i<t}p_{i}^{\alpha}$

7: if the size of predictor’s buffer $|D_{t}^{\alpha}|\geq B$ then

8: for $j=1,...,B$ do

9: Sample experience $e_{j}$ from $D_{t}^{\alpha}$ according to $P_{j}$
calculated by (10)

10: Compute IS weight by (12)

11: Update priority by (9)

12: end for

13: Update the model parameters by (13)

14: end if

15: Update the input feature $\bm{x}_{t+1}$ at the next time-step based on the
true load value $y_{t}$

16:end for

Line 3 – Line 5: The input feature $\bm{x}_{t}$ and target variable $y_{t}$
are the new data at the time-step $t$, which are stored into the predictor’s
buffer $D_{t}^{\alpha}$ at each time-step.

Line 6 – Line 13: By drawing samples from the predictor’s buffer
$\left\\{e_{j}\right\\}_{j=1}^{B}\sim D_{t}^{\alpha}$ according to the
probability, the predictor’s parameters are updated on the small batch of
experiences.

Line 14: After parameter updating, the feature at the next time-step is
updated based on the true load value $y_{t}$ at time $t$ using the moving-
split window method. Then, the predictor estimates the next time-step quantile
based on the updated feature and model parameters.

###  III-B The Proposed Online OPI Method

The proposed online OPI method is shown in Algorithm 2. The value-based RL
method dueling deep Q-networks (DQN) is used. We note that the specific RL
method is not the main focus of the work, and many other RL methods with
discrete action space can also be applied. Dueling DQN [32] divides the
estimation of state-action value function $Q_{\pi}(\bm{x}_{t},\bm{\alpha})$
under the policy $\pi$ into two streams: the respective estimations of state
value $V_{\pi}(\bm{x}_{t})$ and state-dependent action advantages
$A_{\pi}(\bm{x}_{t},\bm{\alpha})$. Here, $\bm{\alpha}$ is the action vector
formed by the discrete proportions in the set defined in (6). Therefore, the
state-action value for selecting any proportion in $\bm{\alpha}$ at time $t$
can be expressed as:

$Q_{\pi}(\bm{x}_{t},\alpha_{t})=V_{\pi}(\bm{x}_{t})+(A_{\pi}(\bm{x}_{t},\alpha_{t})-\frac{1}{|\mathcal{A}|}\sum_{\alpha_{t}^{\prime}\in\mathcal{A}}A_{\pi}(\bm{x}_{t},\alpha_{t}^{\prime})),$
$None$

where $\alpha_{t}$ is the selected proportion.

Based on the state-action value (14), to balance the exploration and
exploitation, the action is chosen in a $\epsilon$-greedy manner.

Since the forecast errors measured by Winkler score at different time (for
instance $t$ and $t+v$, where $v$ is a positive integer) are autocorrelated
[33], the reward at time $t$ is correlated with the reward at time $t+v$. So
the agent needs to be farsighted to consider the future reward at the current
time-step, which matches the delayed reward feature of RL. Therefore, it is
natural to model the state-action value at time $t$ as the expected discounted
reward function:

$Q_{\pi}(\bm{x}_{t},\alpha_{t})=E_{\pi}[r_{t}+\gamma
r_{t+1}+\gamma^{2}r_{t+2}+...|S_{t}=\bm{x}_{t},A_{t}=\alpha_{t}],$ $None$

where $\gamma$ is the discount factor.

Using the Bellman equation, the current state-action value is decomposed into
immediate reward plus discounted state-action value of successor state.
Therefore, the agent of dueling DQN performs the minimization of temporal
difference (TD) error:

$L=\mathbb{E}[(r_{t}+\gamma\mathop{max}\limits_{\alpha_{t+1}}Q_{\pi^{\prime}}(\bm{x}_{t+1},\alpha_{t+1};\bm{\theta}^{-})-Q_{\pi}(\bm{x}_{t},\alpha_{t};\bm{\theta}))^{2}],$
$None$

where $Q_{\pi^{\prime}}(\bm{x}_{t+1},\alpha_{t+1};\bm{\theta}^{-})$ is
estimated by the target network, and
$Q_{\pi}(\bm{x}_{t},\alpha_{t};\bm{\theta})$ is estimated by the local
network, and $\bm{\theta}^{-},\bm{\theta}$ are the parameters of target and
local networks respectively. Based on (16), the local network is updated by:

$\begin{split}&\frac{\partial
L}{\partial\bm{\theta}}=\eta_{Q}\cdot\mathbb{E}[(r_{t}+\gamma\mathop{max}\limits_{\alpha_{t+1}}Q_{\pi^{\prime}}(\bm{x}_{t+1},\alpha_{t+1};\bm{\theta}^{-})-\\\
&Q_{\pi}(\bm{x}_{t},\alpha_{t};\bm{\theta}))\frac{\partial
Q_{\pi}(\bm{x}_{t},\alpha_{t};\bm{\theta})}{\partial\bm{\theta}}],\end{split}$
(17)

where $\eta_{Q}$ is the agent’s learning rate.

Then, the target network’s parameters are updated by:

$\bm{\theta}\cdot\tau+(1-\tau)\cdot\bm{\theta}^{-}\rightarrow\bm{\theta}^{-},$
$None$

where $\tau$ is the soft update parameter.

Algorithm 2 Online OPI Method

1:Batch size $B$, learning rates $\eta_{\alpha}$,$\eta_{Q}$, parameters
$\sigma$,$\rho$,$\tau$,$\gamma$, NCP $1-\beta$

2:Initialize $2\cdot|\mathcal{A}|$ QR models’ and RL agent’s parameters with
random weights

3:for $t$=1,2,… do

4: Given the input state $\bm{x}_{t}$, with the probability $\epsilon$, the
agent selects a random action; otherwise the agent selects
$\alpha_{t}^{*}=\mathop{\arg\max}Q_{\pi}(\bm{x}_{t},\bm{\alpha})$

5:// Update parameters for online quantile predictors

6: Execute action in the environment: Select the QR models for LUB’s quantiles
with probability proportions $\alpha_{t}^{*}$ and $\alpha_{t}^{*}+1-\beta$,
and predict the corresponding quantiles:
$\hat{q}_{t}^{\alpha_{t}^{*}}=f^{\alpha_{t}^{*}}(\bm{x}_{t};\hat{\bm{W}}_{t-1}^{\alpha_{t}^{*}})$,$\hat{q}_{t}^{\alpha_{t}^{*}+1-\beta}=f^{\alpha_{t}^{*}+1-\beta}(\bm{x}_{t};\hat{\bm{W}}_{t-1}^{\alpha_{t}^{*}+1-\beta})$

7: Reveal true value $y_{t}$

8: Respectively store the experience $e_{t}=(\bm{x}_{t},y_{t})$ in the
selected predictors’ buffers
$D_{t}^{\alpha_{t}^{*}}$,$D_{t}^{\alpha_{t}^{*}+1-\beta}$ with maximal
priority $p_{t}^{\alpha}=\mathop{max}_{i<t}p_{i}^{\alpha_{t}^{*}}$,
$p_{t}^{\alpha_{t}^{*}+1-\beta}=\mathop{max}_{i<t}p_{i}^{\alpha_{t}^{*}+1-\beta}$

9: Execute the Line 6 – Line 13 of Algorithm 1 to update QR models’ parameters
for the selected LUB’s predictors, respectively.

10:// Update parameters for RL’s agent

11: Calculate the reward in (8) and update the input feature $\bm{x}_{t+1}$
for the next time-step based on the true load value $y_{t}$

12: Store the transition $(\bm{x}_{t},r_{t},\alpha_{t}^{*},\bm{x}_{t+1})$ in
the dueling DQN’s buffer $D_{t}^{Q}$

13: Sample random minibatch of transitions from $D_{t}^{Q}$, and update the
agent’s network using (17)

14:end for

##  IV Case Study

This section testifies the effectiveness of the proposed method based on real-
world datasets. And the aim is to show the proposed online method (1) compared
with online CPI method, adjusts adaptively to either symmetric or unsymmetric
data distributions, (2) by learning from newly arriving data, achieves better
performance than offline-trained PI estimation methods, and (3) compared with
offline-trained CPI and OPI methods, is robust against the concept drift.
Also, the effectiveness of the PER strategy is demonstrated in subsection E,
and convergence performance is analyzed in subsection F.

###  IV-A Implementation Details

An hourly ahead forecasting is investigated. Smart meter data from Low Carbon
London trail [34] is used to verify the proposed method. In the trail, the
customers receive the Time-of-Use (ToU) tariff are in the ToU group, while
those receive the flat tariff are in the non-ToU group. To validate the
proposed method is effective to data with different underlying distributions
and to PIs with different NCP, the two cases are considered. And PIs with
95$\%$ and 90$\%$ NCP are constructed for the two cases respectively. The
first case is based on the aggregated hourly load in the year of 2013 of 251
customers in the ToU group. The second case considers the net load, which is
the aggregated hourly load of 251 customers in the ToU group minus the hourly
wind output. The probability distributions of data in the two cases are
different. Specifically, it is usually believed that the probability
distribution of the load is approximated with the Normal distribution.
Therefore, the distribution of the data in the first case has relatively small
skewness. And the probability distribution of the wind is believed to be
approximated with the Weibull distribution, which is positive skew with large
skewness. Therefore, the distribution of the net load in the second case is
positively skewed and has relatively larger skewness.

Table I: Summary of Quantile Predictors’ Parameters Item | Value
---|---
Batch size | 128
No. of neurons in each hidden layer | 128
No. of hidden layers | 1
No. of neurons in input layer | 168
No. of neurons in output layer | 1
Optimizer | Adam
Learning rate | 1e-3
Table II: Summary of Dueling DQN’s Agent’s Parameters Item | Value
---|---
Batch size | 128
No. of hidden layers | 2
No. of neurons in the first hidden layer | 512
No. of neurons in the second hidden layer | 256
No. of neurons in input layer | 168
Optimizer | Adam
Learning rate | 1e-4

Since the specific forecasting model is not the main focus of the work,
quantile multi-layer perception (QMLP) is used for the 2$|\mathcal{A}|$
quantile predictors. There is no restriction to the type of machine learning-
based predictor in the proposed framework. As long as the predictor can
incrementally update its model when new data is collected, it can be used in
the proposed setting. For example, NN can update parameters from online
experiences fairly well, so many other types of NN-based models can be used
here. The historical load 168 hours ago is used as the input for both
predictors’ feature and agent’s state. The quantile predictors’ parameters are
summarized in Table I, along with the parameters of Dueling DQN’s agent
summarized in Table II, where the number of neurons in output layer is equal
to the number of actions which will be discussed in the later subsection.
Winkler, average coverage deviation (the reliability/calibration score), and
sharpness scores are used as the three criteria [20] for comparison. The
program is implemented on the laptop with Intel®CoreTM i5-10210U 1.6 GHz CPU,
and 8.00 GM RAM111Codes will be available at
https://github.com/YufanChang/Optimal-Adaptive-Prediciton-Interval after
publication..

###  IV-B Comparison with Online CPI Forecasting Methods

Under the cases of load and net load, the proposed online closed-loop OPI
forecasting with PER is compared with online CPI forecasting with PER. The
number of quantile proportions changes from 3 to 63. The results are shown in
Fig. 4, where the results of the load are in Fig. 4 (a), and the results of
the net load are in Fig. 4 (b). The Winkler score is used as the overall
quality criteria of PIs, which is calculated as the average values over the
entire forecasting steps.

Figure 4: Comparison of Winkler scores of the proposed method and online CPI
in two cases.

In terms of the Winkler score, the proposed method and online CPI forecasting
with PER show interesting results under the cases of load and net load. With
the change of the number of nominal proportions, in the case of the load, the
Winkler scores of the proposed method fluctuate around the Winkler score of
online CPI forecasting with PER; in the case of the net load, the Winkler
scores of the proposed method are lower than that of online CPI forecasting
with PER, which indicates the better performance. This is because the data in
the first case has small skewness, the optimal lower and upper quantiles are
approximated with that of the CPI. Therefore, the proposed method and online
CPI forecasting with PER have the similar performances. However, for the
second case, with the penetration of the wind, the distribution of the net
load becomes more complex and has larger skewness. The optimal lower and upper
quantiles are very different from that of the CPI. Thus, online CPI
forecasting with PER cannot adapt to the distribution of the net load, which
results in the poor performance. Also, we notice that there are some
fluctuations in results of Fig. 4 under different sizes of action space. It is
brought by the exploration mechanism of the RL algorithm and uncertainty
caused by quantile crossing, a common concern in QR models. With the increase
penetration of the wind, the skewness of the net load distribution becomes
more evident. Thus, the online CPI forecasting with PER becomes more
unsuitable. To prove this point, the noises sampled from the Beta distribution
$4\cdot Beta(\alpha=2,\beta=5)$ are added to the net load to increase the
skewness. The relative Winkler scores, which is the Winkler score of online
CPI forecasting with PER minus the scores of the proposed method, are shown in
Fig. 5. It shows that the improvements of the Winkler scores become more
obvious when added with the noise. So the results demonstrate the
effectiveness of the proposed method which adaptively selects optimal
quantiles conditioned on the data distribution.

Figure 5: Relative Winkler score of the net load under different levels of
noise.

###  IV-C Comparison with Offline-trained PI Forecasting Methods on Test Set

To further evaluate the forecast performance of the proposed approach,
offline-trained PI forecasting comparison candidates belonging to the
categories of CPI and OPI forecasting approaches are investigated.
Specifically, the mixed integer programming model (MLMIP) [20], belonging to
adaptive OPI methods, is used. And the CPI methods include the light quantile
gradient boosting regression tree (Light QGBRT), QMLP guided by the pinball
loss, and a naive benchmark that is an extension of persistence model in
probabilistic setting, which constructs PIs for each hour based on symmetric
quantiles derived from historical statistical data [8]. Here, the network
parameters of QMLP method are the same as that of the quantile predictors of
the proposed method. To make fair comparison, although the proposed online
method is performed in a stream fashion and there is no need of the dataset
division, the data is partitioned into the training and test sets of size
70$\%$ and 30$\%$ for the comparison. The numbers of quantile proportions of
the proposed method are chosen as 7 and 3, respectively for the two cases. For
the proposed method and comparison candidates, the results on the test set
after the learning on training set are presented in Fig. 6.

Figure 6: Quality criteria of PIs of the proposed method and offline
comparison candidates on test set in the two cases.

Due to the constraint of restricting the empirical coverage probability no
less than the NCP, MLMIP method always has good calibration performance.
However, the Winkler score of MLMIP method is not satisfying, especially for
the PI forecasting task in the second case. It is inferred that although MLMIP
method is independent of distribution assumption, and demonstrates its good
performance for the load level in MW, the relatively simple forecasting model,
ELM that the approach relies on, may dampen its ability to capture the complex
data distribution that the net load has, which is more variable and uncertain.
Also, if like [20], multiple MLMIP models are trained for each month
separately, the performance of it is expected to be improved. Since the
comparisons in this subsection focus on the performance of one particular
model for the whole dataset, this case is not investigated.

In terms of Winkler score, in the two cases, the proposed method has the
lowest Winkler score on the test set. Compared with the naive benchmark, it
improves the performance by 33$\%$ in the first case, and by 22$\%$ in the
second case. Among the offline methods, QMLP performs the best. However, the
QMLP-based method sacrifices the sharpness for the high reliability.
Specifically, in the first case, although the CPI is approximated with the
OPI, the proposed method’s merit, i.e., consecutive learning from the new
data, enables it to better capture the future trend of the load. Therefore,
the proposed method has the smaller Winkler score with sharper PIs than the
QMLP-based method. In the second case, the improvement of the Winkler score
compared with the best comparison candidate is mainly due to the CPI
forecasting method’s incapability of adapting to data with large skewness.

Fig. 7 illustrates the 168 hours PIs on test set predicted by the proposed
method, superimposed with the CPIs predicted by the QMLP in the two cases. In
Fig. 7 (a), the proposed method achieves the lower upper bounds and higher
lower bounds than QMLP. When the load displays large variability, the width of
PI is relatively large to capture the uncertainty. During the other periods,
the upper and lower bounds of the proposed method are very close to the true
load value. The good performance can be attributed to its adaptive adjustment
to the stream of data. In Fig. 7 (b), the net load displays the larger
uncertainty than the load in Fig. 7 (a). And the negative value of the net
load means the power flows from the demand side to the grid. The width of PI
produced by QMLP is relatively larger than the proposed method. And the
proposed method adapts well to the sudden decrease of the net load, which is
partly caused by the sudden increase of the wind power.

Figure 7: 168 hours PIs forecast on test set obtained via the proposed method
and QMLP in two cases.

The time for updating the model parameters of the proposed method and offline-
trained comparison candidates is shown in Table III. Since the updating time
for the load and net load cases is similar, only the time in the net load case
is displayed. It shows that the proposed method takes the smallest time for
updating the model parameters. Among the offline-trained methods, since Light
QGBRT is a light model designed to be computationally efficient, it takes the
smallest time. The results in Table III are reasonable. The model updating of
the proposed method only relies on a small batch of data. In contrast, the
model updating of the offline methods involves the large numbers of
iterations, and the model can only be updated until enough samples are
collected.

Table III: Comparison of time for updating the model parameters The proposed method | Light QGBRT | QMLP | MLMIP
---|---|---|---
96.8ms | 930ms | 21.5s | 1min 44s

###  IV-D Performance Comparison under the Concept Drift

The concept drift means the statistical properties of the target variable,
which the model is trying to predict, changes in unforeseen ways. To prove
that the proposed online method is robust against the concept drift, the
methods are examined on the scenarios where the training and test sets have
different data distributions. Concretely, for the first case, the training set
is replaced with the synchronous aggregated hourly load in the year of 2013 of
150 customers in the non-ToU group. For the second case, the training set is
replaced with the same aggregated load in the non-ToU group minus the hourly
wind output. And the test sets of the two cases remain the same as those in
the subsection C. Therefore, for the offline-trained comparison candidates,
the results presented in the subsection C are regarded as the benchmarks to
show the results under no concept drift. For the proposed method, the
performances after online learning on the same test sets under no concept
drift are the benchmarks.

The Winkler scores on the test set are shown in Fig. 8. In the two cases, the
proposed method has the similar Winkler scores with the benchmarks, and has
the best performance among the comparison candidates. Thanks to the PER
strategy, the data following the new distribution are more likely to be
sampled frequently for the QR models to learn. Therefore, the method can
efficiently adapt to the changes in load patterns. For the offline-trained
comparison candidates, the concept drift deteriorates the performances.
Especially in the second case where the data distribution is more complex, the
offline-trained methods have larger Winkler score under the concept drift. For
example, the Winkler scores of MLMIP, Light QGBRT, and Naive method are nearly
twice as large as the benchmarks. Therefore, the results solidly verify the
better performance of the proposed method.

Figure 8: Winkler score of PIs of the proposed method and offline comparison
candidates on test set in the two cases under concept drift.

###  IV-E Comparison with Online OPI Forecasting Without PER Strategy

To demonstrate the effectiveness of PER strategy, the proposed method is
compared with its counterpart without PER strategy, where a small batch of
data are randomly sampled from buffer to update the QR model recursively. The
results over the entire forecasting steps are shown in Fig. 9, where the
results of the load are in the left column, and the results of the net load
are in the right column.

Figure 9: Quality criteria of PIs of the proposed method and its counterpart
without PER strategy in the two cases with the change of quantile proportions.

It is evident that in the two cases, the counterpart without PER strategy has
the most satisfactory sharpness performance. However, the calibration of this
approach undergoes the significant discrepancy from the NCP. Since the
sharpness and calibration only display the one aspect of the PI quality, the
overall metric such as Winkler score is needed for the performance evaluation.
In terms of the performance evaluated by the Winkler score, the proposed
method outperforms the online counterpart without PER strategy. The result
demonstrates the effectiveness of the PER technique. Therefore, the ways that
the experiences are sampled to update predictors’ parameters can affect the
performance of QR models for online quantile forecasting. Compared with
updating the predictors’ parameters by random experiences in the buffer, PER
enables the learning to be more effective, because the samples hard to predict
and the new sample are learned more frequently by the QR models.

###  IV-F Convergence Analysis

To demonstrate that the proposed method has good convergence performance, we
show PIs in the first 168 hours of QR models’ learning processes in Fig. 10,
and the moving average rewards gained by agent in Fig. 11. It is shown that
the upper and lower quantiles produced by QR models cannot capture the trends
in load or net load profiles at the very beginning, but QR models learn
quickly to issue accurate predictions in a few steps and obtain PIs with good
quality. For agent’s policy learning, the achieved average rewards increase
during the learning process, and converge soon to the maximum value. The
learning processes of QR models and RL’s agent can affect each other. However,
since QR models and RL’s agent work toward the same goal, namely improving
PI’s quality, the dynamic learning processes of QR models and RL’s agent are
stable and converge fast.

Figure 10: PIs in the first 168 hours of QR models’ learning processes.
Figure 11: RL’s agent’s moving average rewards during the learning process.

##  V Conclusions

In this paper, we propose a new approach to determine the optimal adaptive
prediction intervals for load forecasting in distribution networks via
reinforcement learning. It provides optimal prediction intervals adaptive to
different data distributions. The approach uses the RL method to unify two
correlated sub-problems, i.e., probability proportion selection and
corresponding quantile forecasting, which allows it to be implemented in an
online closed-loop fashion. Case studies on both load and net load forecasting
have shown that the proposed method can produce good quality PIs. Compared
with the online CPI method, it selects the optimal quantiles’ probability
proportions with good flexibility, and thus adapts well to the skewness of
data distributions. Also, by learning from the load incrementally, it obtains
lower Winkler score and is more robust against concept drift compared with
offline-trained methods. The case studies also reveal that the PER strategy
improves the online learning efficiency, and helps to achieve better
performance.

In this work, we approximate the probability proportions in a discrete action
space. Therefore, it is still required to consider representing probability
proportions in a continuous action space in the future. Besides, it is worthy
to note that the proposed closed-loop framework is a general way that links
two separate tasks working toward the same goal. Its potential application is
not restricted to forecasting problems in power systems; instead, it may be
valuable to investigate the applications of the closed-loop framework on other
problems consisting of two separate tasks or stages.

## Acknowledgement

This work was supported by the National Natural Science Foundation of China
(U1866206).

## References

  * [1] W. Liu, Q. Wu, F. Wen, and J. Østergaard, “Day-ahead congestion management in distribution systems through household demand response and distribution congestion prices,” _IEEE Transactions on Smart Grid_ , vol. 5, no. 6, pp. 2739–2747, 2014.
  * [2] E. Vega-Fuentes, J. Yang, C. Lou, and N. K. Meena, “Transaction-oriented dynamic power flow tracing for distribution networks—definition and implementation in gis environment,” _IEEE Transactions on Smart Grid_ , vol. 12, no. 2, pp. 1303–1313, 2021.
  * [3] T. Hong, _Short term electric load forecasting_. North Carolina State University, 2010.
  * [4] J. W. Taylor and P. E. McSharry, “Short-term load forecasting methods: An evaluation based on european data,” _IEEE Transactions on Power Systems_ , vol. 22, no. 4, pp. 2213–2219, 2007.
  * [5] T. Hong, P. Pinson, Y. Wang, R. Weron, D. Yang, and H. Zareipour, “Energy forecasting: A review and outlook,” _IEEE Open Access Journal of Power and Energy_ , vol. 7, pp. 376–388, 2020.
  * [6] T. Gneiting and M. Katzfuss, “Probabilistic forecasting,” _Annual Review of Statistics and Its Application_ , vol. 1, pp. 125–151, 2014.
  * [7] B. Liu, J. Nowotarski, T. Hong, and R. Weron, “Probabilistic load forecasting via quantile regression averaging on sister forecasts,” _IEEE Transactions on Smart Grid_ , vol. 8, no. 2, pp. 730–737, 2015.
  * [8] H. Wen, J. Gu, J. Ma, L. Yuan, and Z. Jin, “Probabilistic load forecasting via neural basis expansion model based prediction intervals,” _IEEE Transactions on Smart Grid_ , vol. 12, no. 4, pp. 3648–3660, 2021.
  * [9] T. Li, Y. Wang, and N. Zhang, “Combining probability density forecasts for power electrical loads,” _IEEE Transactions on Smart Grid_ , vol. 11, no. 2, pp. 1679–1690, 2019.
  * [10] H. Qiu, W. Gu, Y. Xu, Z. Wu, S. Zhou, and J. Wang, “Interval-partitioned uncertainty constrained robust dispatch for ac/dc hybrid microgrids with uncontrollable renewable generators,” _IEEE Transactions on Smart Grid_ , vol. 10, no. 4, pp. 4603–4614, 2018.
  * [11] C. Huang, D. Yue, J. Xie, Y. Li, and K. Wang, “Economic dispatch of power systems with virtual power plant based interval optimization method,” _CSEE journal of power and energy systems_ , vol. 2, no. 1, pp. 74–80, 2016\.
  * [12] A. Khosravi, S. Nahavandi, D. Creighton, and A. F. Atiya, “Comprehensive review of neural network-based prediction intervals and new advances,” _IEEE Transactions on neural networks_ , vol. 22, no. 9, pp. 1341–1356, 2011\.
  * [13] J. Zhang, Y. Wang, M. Sun, and N. Zhang, “Two-stage bootstrap sampling for probabilistic load forecasting,” _IEEE Transactions on Engineering Management_ , pp. 1–9, 2020.
  * [14] T. Heskes, W. Wiegerinck, and H. Kappen, “Practical confidence and prediction intervals for prediction tasks,” _PROGRESS IN NEURAL PROCESSING_ , pp. 128–135, 1997.
  * [15] C. Zhao, C. Wan, and Y. Song, “Cost-oriented prediction intervals: On bridging the gap between forecasting and decision,” _IEEE Transactions on Power Systems_ , 2021.
  * [16] C. Wan, Z. Xu, P. Pinson, Z. Y. Dong, and K. P. Wong, “Probabilistic forecasting of wind power generation using extreme learning machine,” _IEEE Transactions on Power Systems_ , vol. 29, no. 3, pp. 1033–1044, 2013\.
  * [17] M. Landry, T. P. Erlinger, D. Patschke, and C. Varrichio, “Probabilistic gradient boosting machines for gefcom2014 wind forecasting,” _International Journal of Forecasting_ , vol. 32, no. 3, pp. 1061–1066, 2016\.
  * [18] Y. Zhang, Q. Ai, and Z. Li, “Admm-based distributed response quantity estimation: a probabilistic perspective,” _IET Generation, Transmission & Distribution_, vol. 14, no. 26, pp. 6594–6602, 2020.
  * [19] Y. Wang, N. Zhang, Y. Tan, T. Hong, D. S. Kirschen, and C. Kang, “Combining probabilistic load forecasts,” _IEEE Transactions on Smart Grid_ , vol. 10, no. 4, pp. 3664–3674, 2019.
  * [20] C. Zhao, C. Wan, Y. Song, and Z. Cao, “Optimal nonparametric prediction intervals of electricity load,” _IEEE Transactions on Power Systems_ , vol. 35, no. 3, pp. 2467–2470, 2020.
  * [21] A. Khosravi, S. Nahavandi, D. Creighton, and A. F. Atiya, “Lower upper bound estimation method for construction of neural network-based prediction intervals,” _IEEE transactions on neural networks_ , vol. 22, no. 3, pp. 337–346, 2010.
  * [22] T. Pearce, A. Brintrup, M. Zaki, and A. Neely, “High-quality prediction intervals for deep learning: A distribution-free, ensembled approach,” in _International conference on machine learning_. PMLR, 2018, pp. 4075–4084.
  * [23] C. Wan, Z. Xu, P. Pinson, Z. Y. Dong, and K. P. Wong, “Optimal prediction intervals of wind power generation,” _IEEE Transactions on Power Systems_ , vol. 29, no. 3, pp. 1166–1174, 2013.
  * [24] K. Miyaguchi and H. Kajino, “Cogra: Concept-drift-aware stochastic gradient descent for time-series forecasting,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 33, no. 01, 2019, pp. 4594–4601.
  * [25] E. Hazan, “Introduction to online convex optimization,” _arXiv preprint arXiv:1909.05207_ , 2019.
  * [26] L. Von Krannichfeldt, Y. Wang, and G. Hug, “Online ensemble learning for load forecasting,” _IEEE Transactions on Power Systems_ , vol. 36, no. 1, pp. 545–548, 2020.
  * [27] M. N. Fekri, H. Patel, K. Grolinger, and V. Sharma, “Deep learning for load forecasting with smart meter data: Online adaptive recurrent neural network,” _Applied Energy_ , vol. 282, p. 116177, 2021.
  * [28] T. Jaakkola, M. I. Jordan, and S. P. Singh, “On the convergence of stochastic iterative dynamic programming algorithms,” _Neural computation_ , vol. 6, no. 6, pp. 1185–1201, 1994.
  * [29] E. Mocanu, D. C. Mocanu, P. H. Nguyen, A. Liotta, M. E. Webber, M. Gibescu, and J. G. Slootweg, “On-line building energy optimization using deep reinforcement learning,” _IEEE transactions on smart grid_ , vol. 10, no. 4, pp. 3698–3708, 2018.
  * [30] T. Liu, B. Wang, and C. Yang, “Online markov chain-based energy management for a hybrid tracked vehicle with speedy q-learning,” _Energy_ , vol. 160, pp. 544–555, 2018.
  * [31] T. Schaul, J. Quan, I. Antonoglou, and D. Silver, “Prioritized experience replay,” _arXiv preprint arXiv:1511.05952_ , 2015.
  * [32] Z. Wang, T. Schaul, M. Hessel, H. Hasselt, M. Lanctot, and N. Freitas, “Dueling network architectures for deep reinforcement learning,” in _International conference on machine learning_. PMLR, 2016, pp. 1995–2003.
  * [33] F.-K. Sun, C. Lang, and D. Boning, “Adjusting for autocorrelated errors in neural networks for time series,” in _Advances in Neural Information Processing Systems_ , vol. 34. Curran Associates, Inc., 2021, pp. 29 806–29 819.
  * [34] J. R. Schofield, R. Carmichael, S. Tindemans, M. Bilton, M. Woolf, G. Strbac _et al._ , “Low carbon london project: Data from the dynamic time-of-use electricity pricing trial, 2013,” _uK Data Service, SN_ , 2015.


# Role of spin-orbit coupling effects in rare-earth metallic tetra-borides : a
first principle study

Ismail Sk  Bajkul Milani Mahavidyalaya, Purba Medinipur-721655,West Bengal,
India  Kazi Nazrul University, Asansol-713340,West Bengal, India  Nandan
Pakhira  Kazi Nazrul University, Asansol-713340,West Bengal, India

###### Abstract

We have investigated the electronic structure of rare-earth tetraborides,
$\textrm{RB}_{4}$, using first-principle electronic structure methods (DFT)
implemented in Quantum Espresso (QE). In this article we have studied heather-
to neglected strong spin-orbit coupling (SOC) effects present in these systems
on the electronic structure of these system in the non-magnetic ground state.
The calculations were done under GGA and GGA+SO approximations using ultrasoft
pseudopotentials and fully relativistic ultrasoft pseudopotentials (for SOC
case). Perdew-Burke-Ernzerhof generalized gradient approximation (PBE-GGA)
exchange-correlation functionals within the linearized plane-wave (LAPW)
method as implemented in QE were used. The projected density of states
consists of 3 distinct spectral peaks well below the Fermi energy and
separated from the continuum density of states around the Fermi energy. The
discrete peaks arises due to rare-earth $s$-orbital, rare-earth $p$ \+ B $p$
and B $p$-orbitals while the continuum arises due to hybridized B $p$, rare-
earth $d$ orbitals. Upon inclusion of SOC the peak arising due to rare-earth
$p$-orbitals gets split into two peaks corresponding to $j=0.5$ and $j=1.5$
configurations. In case of $\textrm{LaB}_{4}$, in the presence of SOC, spin-
split $4f$ orbitals contributes to density of states at the Fermi level while
the density of states at the Fermi level largely remains unaffected for all
other materials under consideration.

##  I Introduction

The strong Coulomb correlations present in $3d$ and $4d$ transition metal
compounds as well as in $4f$ lanthanides and $5f$ actinides are key to
understanding nonvel and exotic properties. The rare earth lanthanides except
Pm are good conductors of heat and electricity. Pm is radioactive with very
short life and its occurance in nature is extremely rare. The rare-earth
metallic tetra-borides exhibit various valency such as di, tri and tetravalent
state [1]. Cerium (Ce) and Terbium (Tb) primarily show tetravalent state where
as the other metallic tetraborides mostly show trivalent state [1]. Recently
intermediated valance state of Yb between $\textrm{Yb}^{2+}$ and
$\textrm{Y}b^{3+}$ is experimentally observerd and Kondo interaction is
significant in this system [2].

Recent observation of fractional magnetic plateau in $\textrm{TmB}_{4}$ and
$\textrm{ErB}_{4}$ have created lot of interest in these class of materials.
Stable magnetization plateau occuring at 1/2 fraction (of saturation
magnetization) and fractional plateaus at $1/7,1/8,\cdots$ etc [3,4].
fractions are similar to the plateaus observed in the Hall resistivity of two
dimensional degenerate electron gases subject to a perpendicular magnetic
field.

It is interesting to mention that the position of the rare-earth atoms as
shown in Fig. 1 forms a two dimensional Archimedian Shastry-Sutherland lattice
(SSL) [5]. SSL consisting of localized spin-1/2 is an example of geometrically
frustrated system with huge spin degeneracy and the observation of
magnetization plateaus is often attributed to this degenearacy. Insulating
$\textrm{SrCu}_{2}(\textrm{BO}_{3})$ [6] is a well studied system which can be
effectively mapped onto a nearest-neighbour SSL. However in metallic rare-
earth tetraborides “localized” spins interacts only through long range RKKY
[7] type of interactions. Hence the mapping of interacting fermionic model
onto an effective spin-1/2 models on SSL with nearest neighbour interaction is
highly non-trivial [5,8]. Correlated and frustrated systems are of great
academic interest as well as they have many potential technologicl
applications like memory device, spintronics, quantum computation etc. [9].
The very first step towards understanding the intriguing thermodynamic and
transport properties in these complex systems is to study their electronic
band structure. In an earlier work [10] electronic structure of
$\textrm{RB}_{4}$ (except $\textrm{TmB}_{4}$) have been studied using first
principle method. However strong atomic spin-orbit coupling effects present in
rare-earth atoms have been neglected. Inclusion of SOC for certain systems [R=
Yb, Pr, Gd, Tb, Dy] in the mangetic state have been considered and also
considered of $4f$ for $\textrm{RB}_{4}$ but there is no systematic study of
such effects in the non-magnetic (paramagnetic) state [11]. In the present
work we make a detailed study of SOC effects on the electronic structure of
rare earth tetra-borides. In particular we have chosen systems (R=La, Ce, Nd,
Sm) with relatively low SOC effects as well as systems (R=Ho,Er,Tm,Lu) with
relatively high SOC effects.

The organization of the rest of the paper is as follows. In Sec. II we discuss
crystal structure of the system. In Sec. III we elaborate the computational
details for band structure. In Sec. IV we discuss the results for various
systems and finally in Sec. V we conclude.

##  II Crystal structure

$\textrm{RB}_{4}$ crystallizes in the tetragonal symmetrey with space group
P4/mbm [12,13]. Fig. 1 summarizes crystal structure of $\textrm{RB}_{4}$ from
different perspectives. Fig. 1 (a) displays the full tetragonal structure
which consists of alternate layers of rare-earth (R) and B ions stacked along
$c$-axis. Fig. 1(b) shows the top view of the crystal structure. There are two
distinct types of B atoms - (i) planar and (ii) octahedral. Boron atoms form
octahedra as well as 7-atom rings in the $a-b$ plane [14]. Ring forming planar
B atoms (shown in blue) which are not part of octahedra also forms dimers and
these dimers are aranged in a regular pattern. In Fig. 1(c) we show one unit
cell formed by four such B octahedra. In Fig. 1(d) we show SSL formed by the B
atoms. From Fig. 1(b) it is clear that out of the 4 B atoms two are nearer
than other two. The exchange interaction between the two near B atoms mimics
the nearest neighbour interaction (J) and the interaction between the distant
B atoms mimics interaction along alternate diagonals. It is important to
mention that B atoms play a crucial role in the electronic structure of these
systems as they are in the $sp$-hybridized state.

Figure 1: (Color online) Tetragonal crystal structure of $RB_{4}$. Fig. 1(a)
represents the full structure consisting of different layers of rare-earth, R,
ions (red) and B (green) stacked along $c$-axis. (b) Top view of the B
sublattice (in the $a-b$ plane) comprising of 7 atom ring and a square formed
by the position of the R atoms. (c) Side view of the B sublattice (along
$c$-axis) showing two different types of B ; one forming dimer (shown in blue)
and the other part of the B octahedra (shown in green). (d) Shastry-Sutherland
lattice in two dimension.

##  III Computational Details

First-principle calculations were performed using density functional theory
(DFT)[15,16] as implemented in the open source package Quantum Espresso [17]
under the Burai [18] framework. The calculations are done within GGA and
GGA+SO approximation. We have used Ultra soft pseudo-potentials [19] Marzari-
Vanderbit smaering [20] for structural optimization and total energy
calculation of the system. Further, Perdew-Burke-Ernzerhof Generalized
Gradient Approximation (PBE-GGA) exchange-correlation functional within the
linearized augmented plane wave (LAPW) method is employed [21,22]. For the
case with SOC effect full relativistic Ultra soft pseudo-potentials were used.
The total Hamiltonian of the Kohn-sham DFT calculations with Spin-Orbit
coupling can be written as [23]

$\hat{H}=\hat{T}+\hat{V}_{ext}+\hat{V}_{es}+\hat{V}_{xc}+\hat{H}^{SOC}=\hat{T}+\hat{V}_{a}+\hat{H}^{SOC},$
(1)

where, $\hat{T}$,$\hat{V}_{ext}$,$\hat{V}_{es}$,$\hat{V}_{xc}$ and
$\hat{H}^{SOC}$ are the kinetic energy operator, external potential operator,
electrostatic or hartree potential operator, exchange-correlation potential
operator and spin-orbit coupling operator respectively. $\hat{V}_{a}$ is the
applied field or Kohn-Sham potential operator. The Hamiltonian $\hat{H}^{SOC}$
for relativistic limit in terms of momentum and spin operator can be expressed
as [23]

$\hat{H}^{SOC}=\frac{i}{4c^{2}}({\nabla}\hat{V}_{a}\times\hat{p}).\hat{s}$ (2)

For the central field approximation the Hamiltonian $\hat{H}^{SOC}$ [23] can
be written as

$\hat{H}^{SOC}=\zeta\hat{l}.\hat{s}$ (3)

where $\hat{l}$ is the angular momentum and
$\zeta=\frac{1}{2m^{2}c^{2}r}\frac{d\hat{V}_{a}}{dr}$, where $c$ is the speed
of light.

The lattice information were taken from the materials research project site
[24]. The lattice constants, the kinetic energy cut-off (Ecutwfc) and charge
density cut-off (Ecutrho) values used for $\textrm{RB}_{4}$ are mentioned in
Table 1. All the calculations were performed on three dimensional crystals
consisting of primitive tetragonal lattice with 20 atoms. The energy
conservation was achieved using $8^{3}$-points in the full Brillouin zone for
sampling. Energy convergence criteria of $10^{-6}$ Ry were used for self-
consistent calculations. The band structure is plotted along the path
involving high symmetry points. The high symmetry points for tetragonal
lattice system in the first Brillouin zone are $\Gamma$=$(0,0,0)$, $X$=
($\frac{\pi}{a}$,0,0), $M$=($\frac{\pi}{a}$,$\frac{\pi}{a}$,0),
$Z$=(0,0,$\frac{\pi}{c}$), $R$=($\frac{\pi}{a}$,0,$\frac{\pi}{c}$),
$A$=($\frac{\pi}{a}$,$\frac{\pi}{a}$,$\frac{\pi}{c}$). Calculated band
structures were plotted along the high symmetry directions $\Gamma-
X-M-\Gamma$, $Z-R-A-Z$, $X-R$, $M-A$.

##  IV Results and Discussion

In this study we have considered 4 canonical systems $\textrm{LaB}_{4}$,
$\textrm{CeB}_{4}$, $\textrm{NdB}_{4}$ and $\textrm{SmB}_{4}$ with relatively
small spin-orbit coupling strength and 4 canonical systems $\textrm{HoB}_{4}$,
$\textrm{ErB}_{4}$, $\textrm{TmB}_{4}$ and $\textrm{LuB}_{4}$ with much larger
spin-orbit coupling effect. In Table 1. we have summarized the lattice
constants for systems under consideration. In Table 2. we have summarized the
atomic spin-orbit coupling energy [25](in units of $\textrm{cm}^{-1})$ of
various rare-earth atoms under consideration.

Materials | a($\AA$) | c($\AA$) | Ecutwfc(Ry) | Ecutrho(Ry)
---|---|---|---|---
$\textrm{LaB}_{4}$ | $7.31066$ | $4.18269$ | $25$ | $225$
$\textrm{CeB}_{4}$ | $7.17377$ | $4.07463$ | $40$ | $340$
$\textrm{NdB}_{4}$ | $7.23842$ | $4.11996$ | $38$ | $342$
$\textrm{SmB}_{4}$ | $7.18656$ | $4.08152$ | $35$ | $315$
$\textrm{HoB}_{4}$ | $7.08619$ | $4.00815$ | $42$ | $340$
$\textrm{ErB}_{4}$ | $7.06973$ | $3.99708$ | $37$ | $332$
$\textrm{TmB}_{4}$ | $7.05321$ | $3.98405$ | $38$ | $340$
$\textrm{LuB}_{4}$ | $7.02687$ | $3.96821$ | $42$ | $378$
Table 1: Lattice constants and the parameters used for the calculations. Elements | SOC Energy($\textrm{cm}^{-1}$) | Elements | SOC Energy($\textrm{cm}^{-1}$)
---|---|---|---
La | $5.6\times 10^{3}$ | Ho | $8.1\times 10^{3}$
Ce | $5.8\times 10^{3}$ | Er | $8.4\times 10^{3}$
Nd | $6.3\times 10^{3}$ | Tm | $8.7\times 10^{3}$
Sm | $6.8\times 10^{3}$ | Lu | $9.3\times 10^{3}$
Table 2: Tabulation for spin-orbit energy Materials | Fermi energy $E_{F}$(eV)
---|---
| Without SOC | With SOC
$\textrm{LaB}_{4}$ | $12.447$ | $12.523$
$\textrm{CeB}_{4}$ | $13.103$ | $13.103$
$\textrm{NdB}_{4}$ | $12.290$ | $12.303$
$\textrm{SmB}_{4}$ | $12.285$ | $12.301$
Table 3: Fermi Energy due to effect of low spin-orbit interaction materials.

It is important to mention that the choice of kinetic energy cut-off and the
number of $k$-points chosen over the irreducible Brillouin zone are extremely
crucial in determining crystal structure and band structure calculation. We
have calculated the total energy as a function of the plane wave kinetic
energy cut-off as well as the number of $k$-points over irrudicible Brillouin
zone. In Fig. 2(a), (b) we show the convergence of the total energy as a
function of kinetic energy cutoff and number of $k$-points for one canonical
system $\textrm{LaB}_{4}$ with low spin-orbit coupling strength. In Fig. 2(c),
(d) we have shown the same for $\textrm{LuB}_{4}$, a material with much larger
SOC strength. It is clear that the kinetic energy cut-off in the range of
20-50 Ry and 30-50 Ry are deemed to be sufficient for convergence of total
energy in these two systems, respectively. Also, we have found that $4\times
4\times 4$ $k$-mesh (defined over irreducible Brillouin zone) is sufficient
for relative stability of tetragonal structure. For the entire calculation we
have chosen a $k$-mesh of size $8\times 8\times 8$.

Figure 2: Total Energy as a function of kinetic energy cut-off and $k$-points
along irreducible edges. Convergence of self-consistent field for (a,b)
$\textrm{LaB}_{4}$ and (c,d) $\textrm{LuB}_{4}$, respectively.

###  IV.1 Systems with low SOC effect

Taking the optimized crystal structure, we have calculated the electronic band
structures and projected density of states (PDOS) with and without spin-orbit
coupling effects under generalized gradient approximations (GGA) and GGA+SO,
respectively. In Table 3. we compare Fermi energy for systems with and with
out SOC. The Fermi energy for $\textrm{LaB}_{4}$ changes significantly but for
other systems change is only at the second decimal place. The main reason is
that except for $\textrm{LaB}_{4}$ (with SOC) the pseudo-potentials in the
non-magnetic state does not involve highly localised $4f$ orbitals and SOC
strongly affects $4f$ orbitals and its effect on other orbitals are only
secondary through hibridization with $4f$ orbitals. In Fig. 3(a) we have shown
the band structure for $\textrm{LaB}_{4}$ with and without SOC effect. The
Fermi level is set to zero for both the cases. As can be clearly observed from
Fig. 3(a), (b) except at discrete symmetry points $\Gamma$, $Z$ and $R$ there
is no significant SOC effect especially near the Fermi energy. However SOC
lifts degeneracy at special symmetry points. Also, it can be observed that
along the path $R-A$ bands are very flat and there is wide gap (of about 4 eV)
between the top and bottom bands in this region. Flat bands correspond to non-
dispersive localized bands arising mainly from deep core level state. In Fig.
4(a), (b) we have shown projected DOS from various orbitals at a given site in
the absence of SOC. At the Fermi level the contribution is predominantly from
B $2p$ and La $3d$. Discrete spectral peaks at -32 eV, -15 eV etc. arises due
to deep core level states like B $1s$, La $1s$, $2p$. In Fig. 4(c) we show
combined PDOS from all atoms as well as total DOS. When we switch on SOC the B
$2p$ state gets split into two peaks corresponding to $j=0.5$ and $j=1.5$.
Also La $3p$ state gets split into two peaks to $j=0.5$ and $j=1.5$. In the
presence of SOC there is contribution of $4f$ state (split into $j=2.5$ and
$j=3.5$) at the Fermi energy. This is an unique feature in the case of
$\textrm{LaB}_{4}$ and is absent in all other systems we have considered in
this study. PDOS corresponding to $4f$ is spread over wide range of energy
from -10 eV to 7 eV but the total spectral weight is much smaller than B $2p$
and La $3d$ contributions. Just above the Fermi level, in the range 0 to 7.5
eV, PDOS arises due to strong hybridization between La $3d$ orbitals and B
$1s$, $2p$ orbitals.

Figure 3: (a) (Color online) Electronic band structure for $\textrm{LaB}_{4}$.
(a) represents band structure without (black) and with SOC (red). (b)
represents same band structure in the narrower energy window about the Fermi
level (set to zero). Figure 4:  (c) PDOS of $\textrm{LaB}_{4}$ in the absence
of SOC. (a) and (b) shows contribution from different orbitals from B and La,
respectively. Top left panel shows combined contribution from all orbitals.
Dicrete peaks at -32 eV, -17 eV arises mainly due to La $1s$ and La $2p$
orbitals while the peak at -15 eV due to B $2p$. (c) and (d) shows spin-split
contribution from varios B and La orbitals in the presence of SOC. The peak at
-17 eV gets split into two peaks with $j=0.5$ and $j=1.5$. Inset of (d) shows
contribution of spin-split $4f$ orbitals about Fermi level. Top right panel
shows the combined contribution from all orbitals.

In Fig. 5 we summarize the band structure and projected density of states of
$\textrm{CeB}_{4}$ with and without SOC effects. $Ce$ is the first atom in the
lanthanides series which contains $4f$ orbital. As can be clearly seen from
Fig. 5(a) and 5(b) in the presence of SOC, otherwise degenerate bands split at
$\Gamma$ and R points but the bands remain degenrate at Z point. As in the
case of $\textrm{LaB}_{4}$ there exists non-dispersive flat bands along
$R-A-Z$ directions and there is a gap of around 4.5 eV between the top and the
bottom bands. In Fig. 5(c) we show the PDOS arising from various atomic
orbitals in the absence of SOC effects. The distinct spectral peaks appearing
at -14 eV and at -17 eV are due to B $2p$ and Ce $3p$ orbitals, respectively.
The extremely narrow spectral peak at -34 eV arises due to deep core level Ce
$1s$ state and B $2p$. The continuum density of states in the energy window
-10 eV to 8 eV arises due to hybridized B $2p$ and Ce $3p$, Ce $4p$, Ce $5d$
orbitals. In Fig. 5(d) we show the effect of SOC on PDOS for various atoms. As
in the case of $\textrm{LaB}_{4}$ there is appearance of extremely narrow peak
at -19 eV due to splitting of spin-degeneate B $2p$ and Ce $3p$ orbitals into
$j=0.5$ and $j=1.5$ manifolds. The DOS in the energy window -10 eV to 8 eV
remains largely unaffected as in the case of $\textrm{LaB}_{4}$ and there is
no additional contribution due to spin split Ce $4f$ orbitals.

Figure 5: (a) (Color online) Combined Band structure of $\textrm{CeB}_{4}$
with out (black) and with SOC (red), (b) same band structure in the narrower
energy window about the Fermi level (set to zero). Degenerate band splits at
$\Gamma$ and R points. Fig. (c) and (d) The PDOS of $\textrm{CeB}_{4}$. The
distinct spectral peaks appearing at -14 eV and at -17 eV are due to B $2p$
and Ce $3p$ orbitals, respectivelyIn the energy window -10 eV to 8 eV
represents continuum DOS due to hybridized B $2p$ and Ce $3p$, Ce $4p$, Ce
$5d$ orbitals.

The band structure and projected density of states of $\textrm{NdB}_{4}$ with
and without SOC effects have been summarised in Fig. 6. As shown in Fig. 6(a)
and 6(b), the spin degenerate bands splits in various regions due to SOC
effects present in these systems. Band splitting is more explicit along the
direction $\Gamma-Z-R$ and $A-Z$. Very few bands cros the Fermi level and far
from Fermi level most of the bands are much less dispersed and nearly flat. In
Fig. 6(c) we show the PDOS from various atoms without SOC effects. As in the
earlier cases the distinct spectral peak at -15 eV arises due to B $2p$
orbitals and the spectral peak at -19.5 eV arises due to B $2p$ and Nd $3p$
orbitals, respectively. The extremely narrow spectral peak at -38 eV arises
mainly due to non-dispersive deep core-level Nd $1s$ orbital. However B $2p$
orbitals have also contribution towards the peak at -38 eV. The continuum
density of states in the energy range between -10.5 eV to 7.5 eV arises due to
bybridized B $2p$ and Nd $3p$, Nd $4p$, Nd $5d$ orbitals. Finally, in Fig.
6(d) we show the PDOS in the presence of SOC effect. The continuum DOS in the
range -10.5 eV to 7.5 eV remains largely unaffected. However the peak at -19.5
eV gets split into two peaks at -18 eV and -21 eV. This arises due to
otherwise degenrate B $2p$ and Nd $3p$ orbitals splitting into $j=0.5$ and
$j=1.5$ manifolds due to SOC effects.

Figure 6: (a) (color online) Combined Band structure of $\textrm{NdB}_{4}$
with out (black) and with (red) SOC, (b) same band structure in the narrower
energy window about Fermi level (set to zero). Fig. (c) and (d) The PDOS of
$\textrm{NdB}_{4}$. Very small narrow spectral peak at -38 eV is found due to
non-dispersive deep core-level Nd $1s$ orbital. Other peaks at -38 eV arises
for B $2p$ orbitals. The continuum DOS in the energy range between -10.5 eV to
7.5 eV arises due to bybridized B $2p$ and Nd $3p$, Nd $4p$, Nd $5d$ orbitals.

In Fig. 7 repregents the band structure and density of states of
$\textrm{SmB}_{4}$ with and without SOC effects. It is interesting to mention
that $\textrm{SmB}_{4}$ is metallic whereas $\textrm{SmB}_{6}$ is a Kondo
insulator where Sm shows mixed valency $\textrm{Sm}^{+2}$ and
$\textrm{Sm}^{+3}$ at the ratio 3:7. In Fig. 7(a) and 7(b) we show electronic
band structure. Splitting of energy bands in the $\Gamma-Z-R$ direction is
much more prominent due to much larger SOC effects. Energy bands along $R-A-Z$
continues to remain flat. In Fig. 7(c) we show PDOS. The discrete peak arising
due to Sm $1s$ shifts further down to -41 eV. The spectral peak at -20.5 eV
and -14.5 eV arises due to B $2p$, Sm $3p$ and B $1s$, $2p$ orbitals,
respectively. The origin of continuum states in the range -10.5 eV to 8 eV is
same as earlier. When we switch on the SOC the spectral peak arising due to
$p$-orbitals of B and Sm gets split into $j=0.5$ and $j=1.5$ states and the
corresponding spectral peaks apears at -23 eV and -19 eV, respectively.

Figure 7: (a) (color online) Combined Band structure of $\textrm{SmB}_{4}$
with out (black) and with (red) SOC, (b) same band structure in the narrower
energy window about the Fermi level (set to zero). Fig. (c) and (d) the PDOS
of $\textrm{SmB}_{4}$. The isolated spectral peak at -20.5 eV and -14.5 eV
corresponds to B $2p$, Sm $3p$ and B $1s$, $2p$ orbitals, respectively.

###  IV.2 Systems with large SOC Effect

In the previous section we have considered SOC effects on 4 canonical systems
with relatively low SOC effect. In this section we consider SOC effects on 4
canonical systems with relatively large SOC effects. In Table. IV we have
summarized the Fermi energy with and without SOC effects.

In Fig. 8 we have shown the band structure and projected DOS for
$\textrm{HoB}_{4}$ with and without SOC effects. SOC effects on the splitting
of energy bands are prominent for wide range of energies. Bands far away from
Fermi energy are also affected due to strong SOC effects. Degeneracy lifting
effect along $\Gamma-Z-R$ are now quite explicit. Fermi level crossing bands
along $Z-R$ are also affected. However flat bands along $R-A-Z$ are not
affected by SOC. In Fig. 8(c) we show PDOS due to various atoms as earlier.
Continuum DOS in the range -10.5 eV to 7.5 eV arizes due to strong
hybridization between $5d$ orbitals of Ho and $2p$ orbitals of B atoms. The
spectral peak at -14 eV arises due to $2p$ orbitals of B atoms while the peak
at -24 eV arises due to Ho $3p$ orbitals. Extremely narrow and isolated peak
at -48 eV arises due to deep core level $1s$ orbital of Ho atom. In Fig. 8(d)
we show the effect of SOC on PDOS. There is enhancement of PDOS around Fermi
level. The spectral peak at -24 eV gets split into two peaks at -27 eV and -22
eV which arises due to Ho orbitals with $j=0.5$ and $j=1.5$, respectively.

Fig. 9 indicates the band structure and projected density of states of
$\textrm{ErB}_{4}$ in the presence and absence of SOC effects. As seen in Fig.
9(a) and 9(b), spin-split bands are quite visible in the energy range -4 eV to
-6 eV along $\Gamma-Z-R$ direction. Band splitting effects near the Fermi
level also starts showing up. Projected density of states as shown in Fig.
9(c) follows similar trend as in the case of other tetra-borides. The
continuum density of states in the range -10 eV to 7 eV arises from the
hybridized B $2p$ and Er $5d$ orbitals. The spectral peak due to Er $1s$ is
now at -51 eV. While the Er $3p$\+ B $2p$ spectral peak is at -25 eV. The
smaller peak arising due to B $2p$ is at -15 eV. As shown in Fig. 9(d),
inclusion of SOC effect causes splitting of the -25 eV peak into $j=0.5$ and
$j=1.5$ states situated at -29 eV and -23.5 eV, respectively.

Materials | Fermi energy $E_{F}$(eV)
---|---
| Without SOC | With SOC
$\textrm{HoB}_{4}$ | $12.271$ | $12.279$
$\textrm{ErB}_{4}$ | $12.280$ | $12.286$
$\textrm{TmB}_{4}$ | $12.274$ | $12.278$
$\textrm{LuB}_{4}$ | $12.251$ | $12.251$
Table 4: Fermi Energy due to effect of large spin-orbit interaction materials.
Figure 8: (a) (Color online) Combined Band structure of $\textrm{HoB}_{4}$
without (black) and with (red) SOC, (b) same band structure in the narrower
energy window about the Fermi level. The spectral peak at -14 eV corresponds
to $2p$ orbitals of B atoms while the peak at -24 eV arises due to Ho $3p$
orbitals. Extremely narrow and isolated peak is found at -48 eV due to deep
core level $1s$ orbital of Ho atom. Fig. (c) and (d) The PDOS of
$\textrm{HoB}_{4}$. Showing the individual contributions from each orbital
atom with and without SOC effect. The Fermi energy is set to zero assigned by
the dotted red line.

Figure 9: (a) (Color online) Combined Band structure of $\textrm{ErB}_{4}$
with (red) and without SOC (black), (b) same band structure in the narrower
energy window about the Fermi level (set to zero). The spin-split bands
appearing in the energy range -4 eV to -6 eV along $\Gamma-Z-R$ direction. (c)
and (d) PDOS of $\textrm{ErB}_{4}$. The spectral peak due to Er $1s$ is now at
-51 eV. While the Er $3p$\+ B $2p$ spectral peak is at -25 eV . At -15 eV due
to B $2p$ orbital. (d) Inclusion of SOC effect causes splitting of the -25 eV
peak into $j=0.5$ and $j=1.5$ states situated at -29 eV and -23.5 eV,
respectively

.

In Fig. 10 we summarize the electronic band structure and projected density of
states of $\textrm{TmB}_{4}$. In an earlier study band structure for
$\textrm{TmB}_{4}$ in the anti-ferromagnetic state was reported. So the
present study is relevant in the paramagnetic state of this system. As shown
in Fig. 10(a) and 10(b) energy bands far from the Fermi level are strongly
affected due to SOC. Energy bands in the energy range -4 eV to -6 eV show
significant splitting especially along $\Gamma-Z-R$ direction. Similar
features are also observable for energy bands in the window 1 eV to 2 eV. Some
of the Fermi level crossing bands show degeneracy lifting effects near Fermi
level. The spectral features are similar to the other tetraborides. The peaks
arising due to Tm $1s$, $3p$ and B $2p$ are at -52.5 eV and -26 eV
respectively. The continuum DOS in the energy range -11 eV to 7 eV arises due
to hybridized Tm $4p$, $5d$ orbitals with B $2p$ orbitals. Inclusion of SOC,
as shown in Fig. 10(d), causes splitting of the -26 eV spectral peaks into a
$j=0.5$ peak at -30 eV and a $j=1.5$ peak at -24 eV.

Finally, in Fig. 11 we show our results for $\textrm{LuB}_{4}$. Incidentally
Lu is the last member of the lanthanide series with completely filled $4f$
orbitals. As in the case of $\textrm{TmB}_{4}$ there is strong SOC effects on
the energy bands in the energy window -6 eV to -4 eV as well as in the window
1 eV to 2 eV. SOC effects on the Fermi level crossing bands near Fermi energy
are less compared to $\textrm{TmB}_{4}$. These features are well summarised in
Fig. 11(a) and 11(b). In Fig. 11(c) and 11(d) we show the projected DOS in the
absence and presence of SOC effects, respectively. The spectral peak at -56 eV
is due to Lu $1s$ orbital while the peak at -15 eV is due to B $1s$ and $2p$.
There is a strong peak at -27 eV arising due to Lu $3p$ orbital. The height of
this peak is much more than the other two discrete peaks. The continuum of
density of states around Fermi level arises due to hybridized B $2p$ and Lu
$5d$ orbitals. In the presence of SOC the peak at -26 eV gets split onto two
peaks at -32 eV and -25 eV with $j=0.5$ and $j=1.5$, respectively.

Figure 10: (a) (color online) Combined Band structure of $\textrm{TmB}_{4}$
without (black) and with (red) SOC, (b) same band structure in the narrower
energy window about Fermi level (set to zero). Fig. (c) and (d) The PDOS of
$\textrm{TmB}_{4}$. Showing the individual contributions from each orbital
atom with and without SOC effect. The continuum DOS in the energy range -11 eV
to 7 eV arises due to hybridized Tm $4p$, $5d$ orbitals with B $2p$ orbitals.
With SOC, the p-orbital split into a $j=0.5$ peak at -30 eV and a $j=1.5$ peak
at -24 eV.  Figure 11: (a) (color online) Combined Band structure of
$\textrm{LuB}_{4}$ without (black) and with (red) SOC, (b) same band structure
in the narrower energy window about Fermi level (indicates dotted red line).
Fig. (c) and (d) The PDOS of $\textrm{LuB}_{4}$. In the presence of SOC the
peak at -26 eV gets split onto two peaks at -32 eV and -25 eV with $j=0.5$ and
$j=1.5$, respectively like other system.

##  V Conclusion

We have investigated the electronic structure of $\textrm{RB}_{4}$ with non-
magnetic ground state. The electronic band structure shows splitting due to
interaction between spin and angular momentum. The bands splitting has been
interpreted with the help of PDOS. It has also been observed that the two new
branches for p-orbital appearing due to SOC effect. In case of
$\textrm{LaB}_{4}$ with SOC, the contribution of $4f$ orbitals to the DOS
about the Fermi level has been observed. This work is partially supported by
WB- DSTBT research grant no. STBT- 11012(26)/31/2019-ST SEC. One of us (NP)
would like to acknowledge hospitality of IIT, Kharagpore. One of us (IS) would
like to thank Bajkul Milani Mahavidyalaya (College) authority for giving me an
opportunity to pursue research as a Ph. D. scholar.

## References

  * (1) Goryachev, Y.M., Kovenskaya, B.A., Dudnik, E.M. et al. Electronic structure of tetraborides of the rare earth elements. J Struct Chem $\bf 16$, 951–954 (1975).
  * (2) Felix Frontini et al. Intermediate valence state in $\textrm{YbB}_{4}$ revealed by resonant x-ray emission spectroscopy. J. Phys.: Condens. Matter $\bf 34$ 34560 (2022)
  * (3) ai Swaroop Sunku, Tai Kong, Toshimitsu Ito, Paul C.Canfield, B. Sriram Shastry, Pinaki Sengupta, and Christos Panagopoulos, Phys. Rev. B $\bf 93$, 174408 (2016)
  * (4) Nandan Pakhira, Jyoti Krishna, S. Nandy, T. Maitra, and A Taraphder. Electronic structure of metallic tetra-boride TmB 4 : An LDA+DMFT study.arXiv:1807.05388
  * (5) B. S. Shastry and B. Sutherland :Exact ground state of a quantum mechanical antiferromagnet. Physica B $\bf 108$ 1069(1981)
  * (6) Miyahara, S. Ueda, K. Theory of the orthogonal dimer Heisenberg spin model for $\textrm{SrCu}_{2}({BO}_{3})_{2}$. J. Physics: Condens. Matter $\bf 15$, R327 (2003).
  * (7) Lawrence Baylor Robinson, Lloyd N. Ferguson, Jr., and Frederick Milstein. Indirect Exchange Coupling of Magnetic Moments in Rare-Earth Metals. Phys. Rev. B $\bf 3$, 1025-1033 (1971)
  * (8) T. Mtasumurai, D.Okuyama, T. Mouri, and Y. Murakami J. Phys. Soc. Jpn $\bf 80$, 074701 (2011)
  * (9) Pradhan, Subhasree and A. Taraphder. “First-Principles Study of The Metallic Rare-Earth Tetraboride $\textrm{TmB}_{4}$.” Materials Today: Proceedings $\bf 4$ (2017): 5532-5536.
  * (10) Z. P. Yin and W. E. Pickett, Rare-earth–boron bonding and 4f state trends in $\textrm{RB}_{4}$ tetraborides.Physical Review B $\bf 77$, 035135 2008
  * (11) Choi, Hongchul and Laref, Amel and Shim, Ji and Kwon, S. and Min, B.(2009). Electronic structures and magnetic properties of $\textrm{RB}_{4}$ (R=Yb,Pr,Gd,Tb,Dy). Journal of Applied Physics. $\bf 105$. 07E107 - 07E107. 10.1063/1.3058707.
  * (12) Olsen,J. Waśkowska, A. Gerward, L. Vaitheeswaran, Ganapathy Venkatakrishnan, Kanchana Svane, A. Shitsevalova, Natalya Fillipov, V.. (2011). HoB4 at high pressure and low temperature: An experimental and theoretical study. High Pressure Research. $\bf 31$. 3-6. 10.1080/08957959.2010.513683.
  * (13) Z. Fisk and M. B. Maple, Multiple Phase Transition in Rare Earth Tetraborides at Low Temperature. Solid State Commun. $\bf 39$, 1189 1981.
  * (14) WILL, G. and SCHFER, W.. ”The crystal structures of $\textrm{ErB}_{4}$ and $\textrm{DyB}_{4}$ studied by neutron diffraction” Zeitschrift fr Kristallographie - Crystalline Materials $\bf 144$, 217-225 (1976)
  * (15) P. Hohenberg and W. Kohn Phys. Rev. $\bf 136$, B864 (1964).
  * (16) W. Kohn and L.J. Sham, Phys. Rev. $\bf 140$, A1133 (1965).
  * (17) Paolo Giannozzi et al. QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. J. Phys.: Condens. Matter $\bf 21$ 395502 (2009)
  * (18) https://nisihara.wixsite.com/burai/
  * (19) Vanderbilt David. Soft self-consistent pseudopotentials in a generalized eigenvalue formalism. Phys. Rev. B $\bf 41$, 7892 –7897 (1990)
  * (20) N. Marzari, D. Vanderbilt, A. De Vita, and M. C. Payne, “Thermal contraction and disordering of the Al(110) surface,” Phys. Rev. Lett. $\bf 82$, 3296–3299 (1999).
  * (21) Perdew, J. P., Zunger, Alex. Self-interaction correction to density-functional approximations for many-electron systems. Phys. Rev. B $\bf 23$, 5048 (1981)
  * (22) Baroni, Stefano and de Gironcoli, Stefano and Dal Corso, Andrea and Giannozzi, Paolo. Phonons and related crystal properties from density-functional perturbation theory. Rev. Mod. Phys. $\bf 73$, 515-520(2001)
  * (23) Islam, M.R., Islam, M.S., Ferdous, N. et al. Spin–orbit coupling effects on the electronic structure of two-dimensional silicon carbide. J Comput Electron $\bf 18$, 407–414 (2019).
  * (24) A. Jain, S.P. Ong, G. Hautier, et al. The Materials Project: A materials genome approach to accelerating materials innovation APL Materials, 2013, $\bf 1$(1), 011002.
  * (25) Martin WC. Table of Spin-Orbit Energies for p-Electrons in Neutral Atomic (core)np Configurations. J Res Natl Bur Stand A Phys Chem. 1971; $\bf 75A$(2):109-111. doi:10.6028/jres.075A.010