# Bayesian adversarial multi-node bandit for optimal smart grid protection
against cyber attacks

Jianyu Xu<EMAIL_ADDRESS>Bin Liu<EMAIL_ADDRESS>Huadong Mo
<EMAIL_ADDRESS>Daoyi Dong<EMAIL_ADDRESS>International Business
School Suzhou, Xi’an Jiaotong-Liverpool University, Suzhou, China  Department
of Management Science, University of Strathclyde, G1 1XQ, Glasgow, UK  School
of Engineering and Information Technology, University of New South Wales,
Canberra, ACT 2600, Australia

###### Abstract

The cybersecurity of smart grids has become one of key problems in developing
reliable modern power and energy systems. This paper introduces a non-
stationary adversarial cost with a variation constraint for smart grids and
enables us to investigate the problem of optimal smart grid protection against
cyber attacks in a relatively practical scenario. In particular, a Bayesian
multi-node bandit (MNB) model with adversarial costs is constructed and a new
regret function is defined for this model. An algorithm called Thompson-Hedge
algorithm is presented to solve the problem and the superior performance of
the proposed algorithm is proven in terms of the convergence rate of the
regret function. The applicability of the algorithm to real smart grid
scenarios is verified and the performance of the algorithm is also
demonstrated by numerical examples.

###### keywords:

Multi-node bandit; reinforcement learning; Bayesian updating; cyber attack;
smart grid.

††thanks: This work is partially supported by Fudan Scholar and Special
Research Grant, UNSW, and by the Australian Research Council’s Discovery
Projects funding scheme under Project DP190101566. Corresponding author
Huadong Mo. Tel. +61-2-62688683. Fax +61-2-62688683

, , ,

##  1 Introduction

The upgrade of traditional grids to smart grids has brought many benefits to
the overall management of power and energy systems, including higher
reliability, better efficiency, improved integration of renewable energy
resources, more flexible choice for stakeholders and lower operation cost
(Konstantelos et al., 2016; Pogaku et al., 2007; Yu et al., 2015). However,
the core technologies, $e.g.$, communication techniques and SCADA systems
(Abiri-Jahromi et al., 2020; Khalili et al., 2020; Rana, Xiang and Wang, 2018;
Todescato et al., 2020), which deliver advantages of smart grids, also open
the grids to vulnerabilities that already exist in the Information and
Communications Technology world. Now, those vulnerabilities pose threats to
smart grids, such as denial of service (DoS) attacks, false data injection,
replay attacks, privacy data theft and sabotage of critical infrastructure
(Gallo et al., 2020; Mo and Sansavini, 2017; Zhu and Martínez, 2013). In
addition, the failures in a smart grid caused by cyber attacks can easily
cascade to other highly dependent critical infrastructure sectors, such as
transportation systems, wastewater systems, health care systems and banking
systems, resulting in extensive physical damage and social and economic
disruption (Abiri-Jahromi et al., 2020; Che et al., 2018).

While governments, the private sector and academia are recognising the cyber
vulnerability of smart grids, the likelihood and impact of a cyber attack are
difficult to quantify. Furthermore, for a smart grid, there may be mandatory
standards and operation requirements from the grid stakeholders. Current risk
management strategies are generally qualitative or heuristic (Patsakis et al.,
2018). In these strategies, some assumptions, $e.g.$, constant reward with
respect to successful anti cyber attack (Rana, Li and Su, 2017; Smith and
Paté-Cornell, 2018), may be unrealistic for most smart grids.

This paper presents a probabilistic risk analysis framework to enhance the
smart grid cyber security. In particular, the dynamical and stochastic
characteristics of smart grids, such as uncertain demands, are taken into
account to investigate the effect of defending strategies on the real
operation cost. The Optimal Power Flow model (Zhang and Papachristodoulou,
2015) is applied to a 11-node radial smart grid originated from the Elia grid,
Belgium. Compared with the existing studies that focus on the inherent risk
(Mo and Sansavini, 2017; Zhang and Papachristodoulou, 2015), such as the
natural degradation and uncertain renewable energy resources for better
maintenance actions and power dispatch, this paper concerns about impact of
the external threat - cyber attacks on the operation cost for effective
deployment of cyber defense teams. In previous work, the cost of each attack
on a node is usually assumed to be a constant (Smith and Paté-Cornell, 2018).
Nevertheless, through investigating some practical scenarios, it has been
found that the costs are more likely to be decided by some adversarial factors
from the nature. Therefore, an adversarial cost sequence associated with each
node is assumed, and a widely used variation constraint is introduced on each
cost sequence. To cope with the objective of sequential decision strategy, the
problem is formulated using the reinforcement learning framework (Li et al.,
2020; Littman, 2015; Sutton and Barto, 2018). Specifically, the Bayesian prior
method (Smith and Paté-Cornell, 2018) is employed for the model parameters and
the problem is formulated as a Bayesian adversarial multi-node bandit (MNB)
model. In addition, a Bayesian minmax type regret function is constructed,
which is subject to the learning context.

Research on online algorithms for adversarial MNB problems started from Auer
et al. (2002), which is later extended to general adversarial reinforcement
learning models by (Even-Dar et al., 2009). Recently, Besbes et al. (2019)
implemented the classical method in Auer et al. (2002) for the adversarial
bandit problem assuming cost sequences with restricted variation and achieved
the state-of-the-art performance for such problems. An obvious gap when using
these methods in the problem of our work is that they do not adapt to the
Bayesian framework and cannot take advantages of the additional information
provided by the Bayesian assumptions. As a result, these algorithms may be
ineffective in Bayesian problems in terms of the convergence rate of the
regret. An alternative feasible method is proposed in Smith and Paté-Cornell
(2018), where model parameters are considered to be system states included in
the state space. Therefore, the problem is no longer Bayesian and can be
solved using existing algorithms for stochastic MNB. However, this method
suffers from a tremendous state space and is always computationally hard even
in problems with moderate sizes.

To cope with the abovementioned technical challenge, an online learning method
is developed in this study which integrates the Thompson sampling method
(Russo and Van Roy, 2014) with the classical Hedge algorithm (Auer et al.,
2002; Russo and Van Roy, 2014). Our algorithm takes advantages of both
algorithms. In particular, Thompson sampling method is used to cope with the
Bayesian framework and Hedge algorithm is applied to the adversarial costs. A
theoretical bound on the regret function is proved in the proposed algorithm,
which is superior to the typical regret bound by the state-of-the-art
algorithms for the same problem. The applicability and numerical performance
of the proposed algorithm is also illustrated through real data and simulation
studies.

The main contributions of our work are summarized as follows.

  * 1.

A relatively practical adversarial reward with restricted variation constraint
is proposed for the MNB model. The prior information of the model parameters
is incorporated through a Bayesian framework and a Bayesian adversarial MNB
model is formulated.

  * 2.

A Bayesian sup regret is defined as the criterion function for decision making
objective of the proposed problem. A new online algorithm, called Thompson-
Hedge algorithm, is developed to solve the problem. The convergence rate of
the Bayesian sup regret for the proposed algorithm is theoretically proven.

  * 3.

Real smart grid data are used to demonstrate the applicability of the proposed
method and numerical results show that the proposed algorithm can achieve the
state-of-the-art performance.

This paper is organized as follows. Section 2 formulates the problem as a
Bayesian reinforcement learning model and constructs the regret function. In
Section 3, an algorithm called Thompson-Hedge algorithm is developed for our
learning model and a theoretical upper bound on the regret function is
established. Section 4 verifies the feasibility of our model and algorithm
using a real data set. In Section 5, a comparative study is investigated on
the performance between our algorithm and a state-of-the-art algorithm.
Section 6 presents concluding remarks.

##  2 Problem formulation

In the electricity network with cyber attacks, an attacker launches a
coordinated attack using multiple attack vectors, because the smart grid
communication networks are physically distributed and highly heterogeneous.
The successful rate of such attack behavior is data-driven, which follows a
Poisson distribution, as demonstrated by the empirical study from U.S.
Department of Energy (Smith and Paté-Cornell, 2018). The network defender aims
to optimally allocate cyber defense teams among nodes in the network, via
probing one node per day. Such a defending action thwarts all attempted cyber
attacks to that node on that day, and also helps update his/her belief about
the uncertain successful rates of attack. The above interaction between the
attacker and the defender has the sequential decision-making nature and leads
itself to a Bayesian MNB model. This model employs proactively defense teams
that traditionally respond to cyber threats after they occur.

The considered attack scenario in a realistic smart grid is that the DoS
attacks block the demand response (DR) messages and dispatch commands
(Pillitteri and Brewer, 2014; Smith and Paté-Cornell, 2018). The DoS attacks
can be accomplished by flooding the communication channel, $e.g.$, the one
between the demand response automation server and customer systems or the one
between the control center and power plant router, with other messages and
commands, or by tampering with the communication channel (Amin et al., 2013;
Pillitteri and Brewer, 2014; Smith and Paté-Cornell, 2018). Above actions can
prevent legitimate DR messages and dispatch commands being received and
transmitted, $i.e.$, depriving authorized access or control to customer
systems and power plant, resulting in demand not being responded and power
plant not being controlled (Pillitteri and Brewer, 2014). Therefore, the
impact of successful DoS attacks on the optimal power flow model of the smart
grid is described by making the target node temporarily unavailable and
disconnected from the grid, which is illustrated in Section 4.1.

Consider a smart grid with $N$ nodes that suffers cyber attacks. Let
$\mathbb{N}\triangleq\left\\{1,\cdots,N\right\\}$ be the set of all nodes. At
each time $t=1,2,\ldots$, the operator’s action is to choose a node to probe.
If $i\in\mathbb{N}$ is probed, the operator observes a (random) number of
$K_{i,\,t}$ cyber attacks on node $i$.

The historical cyber incident data from US Department of Energy (Fig. 2 of
Smith and Paté-Cornell (2018)) has demonstrated that the arrival interval of
successful cyber attacks on per node can be well described by a truncated
Poisson distribution, and there is no record of multiple cyber attacks
(maximum is 3) on per node per day. Therefore, the physical meaning is that in
real applications, there should not be more than 3 cyber attacks arriving per
node in a defined time interval. In addition, the Palm-Khintchine theorem
(Heyman and Sobel, 2004; Smith and Paté-Cornell, 2018) justifies that the
aggregate arrivals from many sources (no need to be strict Poisson) approach a
truncated Poisson distribution in the limit time interval. It is assumed that
for all $i\in\mathbb{N}$, $\left\\{K_{i,\,t},\,t=1,2,\ldots\right\\}$ is a
sequence of independent and identically distributed random variable sequence
drawn from a Poisson distribution truncated at $m>0$ with rate $\lambda_{i}$,
$i.e.$, for all $k=1,2,\ldots$, and $t=1,2,\ldots$,

$\displaystyle
Pr\left(K_{i,\,t}=k\,|\,\lambda_{i}\right)=\left\\{\begin{array}[]{cc}\frac{\lambda_{i}^{k}}{k!}e^{-\lambda_{i}},&k<m\\\
\sum_{j=m}^{\infty}\frac{\lambda_{i}^{j}}{j!}e^{-\lambda_{i}},&k=m\\\
0,&\textup{otherwise}\end{array}\right..$

In particular, the average number of attacks on each node $i\in\mathbb{N}$,
denoted by $\mu_{i}$, is

$\mu_{i}=\sum_{k=1}^{m-1}k\cdot\frac{\lambda_{i}^{k}}{k!}e^{-\lambda_{i}}+m\sum_{k=m}^{\infty}\frac{\lambda_{i}^{k}}{k!}e^{-\lambda_{i}}.$

Let $\lambda=(\lambda_{1},\cdots,\lambda_{N})$ and assume a foresight belief
on $\lambda$ in terms of a prior distribution
$\mathcal{F}=\left(\mathcal{F}_{1},\cdots,\mathcal{F}_{N}\right)$ over
$\lambda$, where $\lambda_{i}\sim\mathcal{F}_{i}$ for all $i\in\mathbb{N}$. In
particular, it is assumed that each $\mathcal{F}_{i}$ is a gamma distribution
with parameters $\left(\alpha_{i},\,\beta_{i}\right)$ and all
$\mathcal{F}_{i}$’s are mutually independent with each other. Whenever $k$
cyber attacks are observed at node $i\in\mathbb{N}$, the parameters
$\left(\alpha_{i},\,\beta_{i}\right)$ are updated through Bayesian posterior
as
$\left(\alpha_{i},\,\beta_{i}\right)\rightarrow\left(\alpha_{i}+k,\,\beta_{i}+1\right)$.

At time $t=1,2,\ldots$, a cost $c_{i,\,t}$ is calculated from an optimization
process for all $i\in\mathbb{N}$. Without loss of generality, it is assumed
that all $c_{i,\,t}$’s are normalized such that $c_{i,\,t}\in[0,\,1/m]$ for
all $i\in\mathbb{N}$ and $t=1,2,\ldots$. The process uses some inputs ($e.g.$,
external factors) that are decided by an adversary (environment). Without
ambiguity, the cost vector $c_{t}=\left(c_{1,\,t},\cdots,c_{N,\,t}\right)$ is
assumed to be decided by the adversary at each time $t=1,2,\ldots$. If node
$i\in\mathbb{N}$ is probed, the operator avoids incurring a total cost of
$X(i,\,t)=K_{i,\,t}\cdot c_{i,\,t}$, $X(i,\,t)\in[0,\,1]$. Equivalently, a
reward of $K_{i,\,t}\cdot c_{i,\,t}$ is achieved by probing node $\pi_{t}$ at
time $t$. An admissible policy $\pi$ is a sequence of mappings
$\pi=\left\\{\pi_{1},\pi_{2},\cdots\right\\}$, where each $\pi_{t}$ is a
mapping from historical observations to the set of all probability
distributions on $\mathbb{N}$. To distinguish, let $\\{i_{1},i_{2},\cdots\\}$
represent the sequence of chosen nodes. The axioms of the problem are
formulated in the problem protocol below.

Algorithm 1 Problem protocol:

1:

2:The adversary decides an cost $c_{i,\,t}$ for all $i\in\mathbb{N}$.

3:The operator probes a node $i_{t}\in\mathbb{N}$ according to $\pi_{t}$ and
observes the number of cyber attacks on this node to be $k_{i_{t},\,t}$.

4:All $c_{i,\,t}$’s are revealed to the operator.

5:A reward $k_{i_{t},\,t}\cdot c_{i_{t},\,t}$ is achieved by the operator.

For notational convenience, let
$\theta=\left(\alpha_{1},\,\beta_{1},\cdots,\alpha_{N},\,\beta_{N}\right)$
denote the initial parameter vector and
$\vec{c}_{T}=\left\\{c_{1},\cdots,c_{T}\right\\}$ the cost sequence up to time
$T$. To proceed, the optimal policy and the regret function is firstly
formulated given the parameter vector $\theta$ and cost sequence
$\vec{c}_{t}$. The mean reward of node $i\in\mathbb{N}$ at time $t$ is
$\mathbf{E}\left(K_{i,\,t}\cdot c_{i,\,t}\right)=\mu_{i}\cdot c_{i,\,t}$.
Since the reward sequence from each node is not stationary due to $c_{i,\,t}$,
for all $t=1,2,\ldots$, the optimal node is defined as
$i_{t}^{*}\triangleq{\arg\max}_{i\in\mathbb{N}}\mu_{i}\cdot c_{i,\,t}$ and a
non-stationary optimal policy is supposed to choose node $i_{t}^{*}$ at time
$t=1,2,\ldots$. Thus, for any admissible policy $\pi$, given $\theta$ and
$\vec{c}_{T}$, the regret function up to time $T$ is defined as

$R^{\pi}\left(\lambda,\,\vec{c}_{T},\,T\right)\triangleq\sum_{t=1}^{T}\mu_{i_{t}^{*}}\cdot
c_{i_{t}^{*},\,t}-\mathbf{E}^{\pi}\left(\sum_{t=1}^{T}\mu_{\pi_{t}}\cdot
c_{\pi_{t},\,t}\,\Big{|}\,\lambda,\,\vec{c}_{T}\right),$ (1)

where $\mathbf{E}^{\pi}$ means the expectation taken with respect to the
(random) sequence $\left\\{i_{1},\cdots,i_{T}\right\\}$ generated by $\pi$.

A constraint is imposed upon the adversary by introducing a sequence of sets
$\left\\{\mathbb{C}_{t}\subset[0,\,1]^{t}\right\\}_{t\geq 2}$, that for all
time $t=1,2,\ldots$, $\vec{c}_{T}\in\mathbb{C}_{t}$. Then, given $\theta$, the
regret function with respect to the worst case up to time $T>0$ is

$\sup_{\vec{c}_{T}\in\mathbb{C}_{T}}R^{\pi}\left(\lambda,\,\vec{c}_{T},\,T\right).$
(2)

Note that
$\sup_{\vec{c}_{T}\in\mathbb{C}_{T}}R\left(\lambda,\,\vec{c}_{T},\,T\right)$
is supposed to measure the performance with the action sequence uniformly on
all possible cost sequences. Finally, incorporating the prior distribution
$\mathcal{F}$ on $\theta$, the Bayesian regret function with respect to the
worst case is defined as

$\displaystyle\mathcal{R}^{\pi}\left(\mathbb{C}_{T},\,T\right)\triangleq\mathbf{E}_{\lambda\sim\mathcal{F}}\left[\sup_{\vec{c}_{T}\in\mathbb{C}_{T}}R^{\pi}\left(\lambda,\,\vec{c}_{T},\,T\right)\right].$
(3)

To distinguish between the regret functions in (1), (2) and (3), they are
named as regret, sup regret and Bayesian sup regret, respectively, for
brevity.

By analyzing some real databases (such as the Elia Grid in Section 4), the
sequence of sets $\left\\{\mathbb{C}_{T}\subset[0,\,1]^{T}\right\\}_{T\geq 2}$
is formulated with constraint on the cost. It is concluded from the
statistical analysis of the real data that the temporal variations of the cost
sequence of each node $i\in\mathbb{N}$, $i.e.$,
$\left\\{|c_{i,\,t+1}-c_{i,\,t}|\right\\}_{t\geq 1}$, are stationary and
uniformly upper bounded by a value, which only relies on $i$ and is small
compared to the average value of $\left\\{c_{i,\,t}\right\\}_{t\geq 1}$.
Therefore, it is assumed that for each node the cumulative temporal variation
up to time $t\geq 1$ is upper bounded by a linear function of $t$. Based on
this assumption, a uniform upper bound on the cumulative temporal variation
for all nodes is introduced and
$\left\\{\mathbb{C}_{T}\subset[0,\,1]^{T}\right\\}_{T\geq 2}$ is constructed
as follows.

$\displaystyle\begin{split}\mathbb{C}_{T}\triangleq&\\{\left\\{c_{i,\,t}\right\\}_{t=1}^{T}\in[0,\,1]^{T},\,i\in\mathbb{N}:\>\\\
&\sum_{t=1}^{T}\max_{i\in\mathbb{N}}|c_{i,\,t+1}-c_{i,\,t}|\leq\mathcal{V}_{T}\\},\>\forall
T\geq 2,\end{split}$ (4)

where $\\{\mathcal{V}_{T}\\}_{T\geq 1}$ is a sequence of positive numbers.
Since $c_{i,\,t}\in[0,\,1/m]$, $\forall i\in\mathbb{N}$ and $t\geq 1$,
$\max_{i\in\mathbb{N}}|c_{i,\,t+1}-c_{i,\,t}|\ll 1/m$ holds for all $t\geq 1$.
Thus, $\mathcal{V}_{T}\leq O(T)$ and $\mathcal{V}_{T}/T\ll 1/m$. Meanwhile,
$\mathcal{V}_{T}$ is monotonically increasing in $T$. Based on these
considerations, the following assumption on $\mathcal{V}_{T}$ is proposed.

######  Assumption 1.

There exists $T_{0}\geq 1$, such that for all $T\geq T_{0}$,
$1/m\leq\mathcal{V}_{T}\leq T/m$.

The restriction of $T\geq T_{0}$ for some $T_{0}$ in Assumption 1 is necessary
to further present our results. Actually the condition
$1/m\leq\mathcal{V}_{T}$ cannot be directly concluded based on the previous
discussions. Since $\mathcal{V}_{T}$ is the accumulation of
$\max_{i\in\mathbb{N}}|c_{i,\,t+1}-c_{i,\,t}|$ and we suppose that
$\max_{i\in\mathbb{N}}|c_{i,\,t+1}-c_{i,\,t}|\ll 1/m$ . Hence, when $T$ is
small, $\mathcal{V}_{T}$ may not satisfy the condition
$1/m\leq\mathcal{V}_{T}$. It is worth noting that a justification of the
existence of $T_{0}$ is useful to guarantee the performance of our algorithm,
while it is not necessary to know the exact value of $T_{0}$. Prior knowledge
can be incorporated to ensure the existence of $T_{0}$. For example, if the
difference term $\max_{i\in\mathbb{N}}|c_{i,\,t+1}-c_{i,\,t}|$ fails to
rapidly become infinitely small (which is natural when the cost sequences are
non-stationary), then $\mathcal{V}_{T}$ will surely be larger than $1/m$ when
$T$ is large enough. Examples of securing an estimated value for $T_{0}$ can
be found in the smart grid application of Section 4, where $\mathcal{V}_{T}$
is estimated to grow linearly in $T$ and a threshold time for
$\mathcal{V}_{T}$ to exceed $1/m$ can be easily calculated.

##  3 Thompson-Hedge algorithm

In this section, an online learning algorithm is developed to optimize the
Bayesian sup regret for our problem. The algorithm has two layers. In the
outer layer, at each time $t\geq 1$, the algorithm uses Thompson sampling
(posterior sampling) method to sample a parameter vector
$\left(\lambda_{1},\cdots,\lambda_{N}\right)$ from the posterior distributions
of all $\lambda_{i}$’s based on the historical observations. Then, in the
inner layer, a so-called sub-algorithm is fed with the sampled parameters. The
sub-algorithm returns an action $\pi_{t}\in\mathbb{N}$. At the end of this
loop, the algorithm chooses node $\pi_{t}$, observes the reward from $\pi_{t}$
and updates the posterior distribution of $\lambda_{\pi_{t}}$. In the end of
this section, a theorem is provided to characterize the convergence rate of
the regret function.

Since Hedge algorithm is applied as the sub-algorithm in the inner layer, the
proposed algorithm is named as Thompson-Hedge algorithm. Before formally
presenting our algorithm, the Hedge algorithm is firstly introduced. Hedge
algorithm is a classical algorithm designed for the adversarial MNB problem
with full feedback (Freund and Schapire, 1999). It is a randomized algorithm
that maintains a weight $\omega_{t}(i)$ for each $i\in\mathbb{N}$. Hedge
algorithm then chooses node $i\in\mathbb{N}$ with a probability proportional
to $\omega_{t}(i)$ at time $t$. Subsequently, the weight for each node is
updated according to the observed reward of this node. The algorithm chooses a
size $\Delta$ and restarts updating the weights of each node every $\Delta$
times. The set of intervals between two successive restarting epochs
(including the first restarting time) is called a batch. Since Hedge algorithm
is designed for non-Bayesian adversarial MNB, in order to use it under the
Bayesian framework, one needs to modify the original Hedge algorithm and feed
it by the value of parameter $\lambda$. Denote Hedge($\lambda$) as the
modified Hedge algorithm used in our Bayesian adversarial MNB fed by
$\lambda$. The algorithm is summarized in Algorithm 1.

Algorithm 2 Algorithm 1: Hedge($\lambda$) Algorithm

1:Parameter vector $\lambda$,
$\varepsilon=\left(1-\sqrt{\ln{N}/2T}\right)^{-1}$, calculate $\mu_{i}$ using
$\lambda_{i}$.

2:for $b=1,2,\cdots,\left\lceil T/\Delta\right\rceil$ do

3: For all $i\in\mathbb{N}$, $\omega_{1}^{i}=1$.

4: while $t\leq\min\left\\{T,\,b\cdot\Delta\right\\}$ do

5: For each $i\in\mathbb{N}$, set
$p^{i}_{t}=\frac{\omega_{t}^{i}}{\sum_{j\in\mathbb{N}}\omega_{t}^{j}}.$

6: Choose a node $\pi_{t}$ according to the probability distribution
$\left\\{p^{i}_{t}\right\\}_{i\in\mathbb{N}}$ and receive a reward
$k_{\pi_{t},\,t}\cdot c_{\pi_{t},\,t}$.

7: For each node $i\in\mathbb{N}$, update
$\omega_{t+1}^{i}=\omega_{t}^{i}\varepsilon^{\mu_{i}\cdot c_{\pi_{t},\,t}}.$

8: Set $t\rightarrow t+1$.

9: end while

10:end for

Hedge($\lambda$) follows the paradigm of the classical Hedge algorithm. The
following lemma shows that if the value of parameter $\lambda$ fed to
Hedge($\lambda$) is the true parameter of the Bayesian adversarial MNB model,
then Hedge($\lambda$) retains a convergent sup regret with a known upper
bound.

######  Lemma 1.

If the input $\lambda$ in the Hedge algorithm is the true model parameter and
the batch size is chosen to be
$\Delta=\left\lceil\left(\ln{N}\right)^{1/3}\left({T}/{m\mathcal{V}_{T}}\right)^{2/3}\right\rceil$,
then for all $T\geq N$, the sup regret by the Hedge($\lambda$) algorithm is
upper bounded by

$\sup_{\vec{c}_{T}\in\mathbb{C}_{T}}R^{\pi}\left(\lambda,\,\vec{c}_{T},\,T\right)\leq\left(8+2\sqrt{2}\right)\left(m\mathcal{V}_{T}\ln{N}\right)^{1/3}\left(T\right)^{2/3}.$
(5)

Detailed proof of Lemma 1 is given in the Appendix. In the proof of our main
result of Theorem 1, the conclusion of Lemma 1 is used as an intermediate
benchmark result.

In the proposed Thompson-Hedge algorithm, parameter $\lambda$ is sampled from
its posterior distribution in each epoch. Then, the sampled parameter is fed
to the inner algorithm Hedge($\lambda$) and a node is chosen correspondingly.
In particular, the algorithm chooses a node according to the probability
weight of each node. Then, the algorithm observes the costs and updates the
posterior distribution of $\lambda$. Next, a sampled parameter $\hat{\lambda}$
is drawn from the posterior distribution of $\lambda$. Finally, the algorithm
updates the weight of each node for next epoch using both the observed costs
and $\hat{\lambda}$. $\textup{Gamma}(\alpha,\,\beta)$ is used to represent
Gamma distribution with parameters $(\alpha,\,\beta)$. The proposed Thompson-
Hedge algorithm is summarized in Algorithm 2.

Algorithm 3 Algorithm 2: Thompson-Hedge Algorithm

1:Parameter $\theta$, batch size $\Delta$,
$\varepsilon=\left(1-\sqrt{\ln{N}/2T}\right)^{-1}$, t=1.

2:for $b=1,2,\cdots,\left\lceil T/\Delta\right\rceil$ do

3: For all $i\in\mathbb{N}$, $\omega_{1}^{i}=1$.

4: while $t\leq\min\left\\{T,\,b\cdot\Delta\right\\}$ do

5: For each $i\in\mathbb{N}$, set
$p^{i}_{t}=\frac{\omega_{t}^{i}}{\sum_{j\in\mathbb{N}}\omega_{t}^{j}}.$

6: Choose a node $i_{t}$ according to the probability distribution
$\left\\{p^{i}_{t}\right\\}_{i\in\mathbb{N}}$ and receive a reward
$k_{i_{t},\,t}\cdot c_{i_{t},\,t}$.

7: For node $i_{t}$, update
$\left(\alpha_{i_{t}},\,\beta_{i_{t}}\right)\rightarrow\left(\alpha_{i_{t}}+k_{i_{t},\,t},\,\beta_{i_{t}}+1\right)$.

8: For each node $i\in\mathbb{N}$, sample
$\hat{\lambda}_{i}\sim\textup{Gamma}(\alpha_{i},\,\beta_{i})$, calculate
$\hat{\mu}_{i}$ using $\hat{\lambda}_{i}$, and update
$\omega_{t+1}^{i}=\omega_{t}^{i}\cdot\varepsilon^{\hat{\mu}_{i}\cdot
c_{i_{t},\,t}}.$

9: Set $t\rightarrow t+1$.

10: end while

11:end for

Before the main result (Theorem 1) is presented, the following lemma is
introduced that will be used for the proof of the main result.

######  Lemma 2.

[(Osband et al., 2013), Lemma 2] If $F$ is the true distribution of $\lambda$,
and $\hat{\lambda}$ is the sampled parameter in epoch $t$, then for any
$\sigma(\mathbb{H}_{t})$-measurable function $g$, it follows

$\mathbf{E}_{\lambda}\left[g(\lambda)\,|\,\mathbb{H}_{t}\right]=\mathbf{E}_{\hat{\lambda}}\left[g(\hat{\lambda})\,|\,\mathbb{H}_{t}\right].$

Lemma 2 shows a central result in Bayesian learning area. The sampled
parameter from the posterior distribution each time can be “considered” as the
true parameter in the sense that any deterministic function using it as an
argument independently has the same expectation. On the basis of Lemma 1 and
Lemma 2, the main theorem is presented as below, which establishes an upper
bound on the Bayesian sup regret for our Thompson-Hedge algorithm and thus
indicates its convergence rate.

######  Theorem 1.

Let $\pi$ be the Thompson-Hedge algorithm with batch size
$\Delta=\left\lceil\left(\ln{N}\right)^{1/3}\left({T}/{m\mathcal{V}_{T}}\right)^{2/3}\right\rceil$,
$\mathbb{C}_{T}$ be defined in (4), and $T_{0}$ be the same value as defined
in Assumption 1, then for all $T\geq T_{0}$, the Bayesian sup regret
$\mathcal{R}^{\pi}\left(\mathbb{C}_{T},\,T\right)$ can be upper bounded by

$\mathcal{R}^{\pi}\left(\mathbb{C}_{T},\,T\right)\leq\left(8+2\sqrt{2}\right)\left(m\mathcal{V}_{T}\ln{N}\right)^{1/3}\left(T\right)^{2/3}.$
(6)

The basic idea for the proof is sketched as follows. The Bayesian sup regret
by Thompson-Hedge algorithm is firstly decomposed into two parts: the regret
of the introduced Hedge($\lambda$) algorithm and the difference of performance
between Thompson-Hedge and Hedge algorithms. Since the sup regret by Hedge
algorithm is bounded by Lemma 1, Theorem 1 can then be proved by upper
bounding the performance difference between the two algorithms. To distinguish
between two algorithms in the work, $\pi^{TH}$ and $\pi^{H}$ are used for
Thompson-Hedge and Hedge algorithms correspondingly. Moreover, $R^{TH}$ and
$R^{H}$ are denoted as the regret function of $\pi^{TH}$ and $\pi^{H}$,
respectively.

a) Separation of the target regret $\mathcal{R}^{TH}$.

In the first step, the Bayesian sup regret $\mathcal{R}^{TH}$ is separated
into a combination of two terms. In Step b) and Step c), the upper bound of
these two terms are obtained separately and therefore, the upper bound of
$\mathcal{R}^{TH}$ can be readily obtained. The separation of
$\mathcal{R}^{TH}$ is given as below.

$\displaystyle\mathcal{R}^{TH}\left(\mathbb{C}_{T},\,T\right)=\mathbf{E}_{\lambda\sim\mathcal{F}}\Bigg{[}\sup_{\vec{c}_{T}\in\mathbb{C}_{T}}\Bigg{\\{}\sum_{t=1}^{T}\mu_{i_{t}^{*}}\cdot
c_{i_{t}^{*},\,t}$
$\displaystyle-\mathbf{E}^{TH}\left(\sum_{t=1}^{T}\mu_{i_{t}}\cdot
c_{i_{t},\,t}\,\Big{|}\,\lambda,\,\vec{c}_{T}\right)\Bigg{\\}}\Bigg{]}$
$\displaystyle\leq\mathbf{E}_{\lambda\sim\mathcal{F}}\Bigg{[}\sup_{\vec{c}_{T}\in\mathbb{C}_{T}}\Bigg{\\{}\sum_{t=1}^{T}\mu_{i_{t}^{*}}\cdot
c_{i_{t}^{*},\,t}$
$\displaystyle-\mathbf{E}^{H}\left(\sum_{t=1}^{T}\mu_{i_{t}i}\cdot
c_{i_{t},\,t}\,\Big{|}\,\lambda,\,\vec{c}_{T}\right)\Bigg{\\}}\Bigg{]}$
$\displaystyle+\mathbf{E}_{\lambda\sim\mathcal{F}}\Bigg{[}\sup_{\vec{c}_{T}\in\mathbb{C}_{T}}\Bigg{\\{}\mathbf{E}^{H}\left(\sum_{t=1}^{T}\mu_{i_{t}}\cdot
c_{i_{t},\,t}\,\Big{|}\,\lambda,\,\vec{c}_{T}\right)$
$\displaystyle-\mathbf{E}^{TH}\left(\sum_{t=1}^{T}\mu_{i_{t}}\cdot
c_{i_{t},\,t}\,\Big{|}\,\lambda,\,\vec{c}_{T}\right)\Bigg{\\}}\Bigg{]}.$ (7)

For notational convenience, denote the above two terms in the two square
brackets by $\Lambda_{1}$ and $\Lambda_{2}$, respectively. The relation in (3)
can be rewritten as

$\mathcal{R}^{TH}=\Lambda_{1}+\Lambda_{2}$

Note that $\Lambda_{1}$ is closely associated with the sup regret of $\pi^{H}$
and can be bounded based on the result in Lemma 1. $\Lambda_{2}$ is the
difference between the return of $\pi^{TH}$ and $\pi^{H}$ that will be bounded
in the following step.

b) Upper bound of $\Lambda_{1}$

$\Lambda_{1}$ is rewritten as

$\displaystyle\Lambda_{1}=\mathbf{E}_{\lambda\sim\mathcal{F}}\left[\sup_{\vec{c}_{T}\in\mathbb{C}_{T}}R^{H}\left(\lambda,\,\vec{c}_{T},\,T\right)\right].$

Note that the upper bound in Lemma 1 holds for any true model parameter
$\lambda$. Thus the upper bound still holds after taking expectation on
$\lambda$ on both sides of (5). Therefore, it holds that

$\displaystyle\Lambda_{1}\leq\left(8+2\sqrt{2}\right)\left(m\mathcal{V}_{T}\ln{N}\right)^{1/3}\left(T\right)^{2/3}.$

c) Upper bound $\Lambda_{2}$

The following clarification is made for notational convenience. $\lambda_{t}$
is used for the sampled parameter by Thompson-Hedge algorithm at time $t$ and
$\lambda$ for the true model parameter which is the input of Hedge algorithm.
To bound the Bayesian sup regret by Thompson-Hedge algorithm, the difference
between the Bayesian sup regret functions by Thompson-Hedge and Hedge
algorithms, respectively, is firstly bounded, namely,

$\displaystyle\begin{split}\mathbf{E}_{\lambda\sim\mathcal{F}}\Bigg{[}\sup_{\vec{c}_{T}\in\mathbb{C}_{T}}\Bigg{\\{}&\mathbf{E}^{TH}\left(\sum_{t=1}^{T}\mu_{i_{t}}\cdot
c_{i_{t},\,t}\,\Big{|}\,\lambda,\,\vec{c}_{T}\right)\\\
&-\mathbf{E}^{H}\left(\sum_{t=1}^{T}\mu_{i_{t}}\cdot
c_{i_{t},\,t}\,\Big{|}\,\lambda,\,\vec{c}_{T}\right)\Bigg{\\}}\Bigg{]},\end{split}$
(8)

where $\mathbf{E}^{TH}$ and $\mathbf{E}^{H}$ have the same meaning as
$\mathbf{E}^{\pi}$ given $\pi$, for Thompson-Hedge algorithm and Hedge
algorithm, respectively. Denote the observation history before time $t\geq 2$
as
$\mathbb{H}_{t}=\left\\{i_{1},\,k_{i_{1},\,1},\,c_{1},\cdots,i_{t-1},\,k_{i_{t-1},\,t-1},\,c_{t-1}\right\\}$.
Conditioned on the observation history, the probability weight function
$p_{t}(\cdot)$ in both algorithms is functions of $\lambda$ and $\lambda(t)$
correspondingly. Moreover, let $p_{t}(\cdot)$ and $\tilde{p}_{t}^{i}(\cdot)$
be the probability functions of Thompson-Hedge and Hedge algorithms,
respectively. For any fixed $\vec{c}_{T}\in\mathbb{C}_{T}$, the following
relation holds,

$\displaystyle\mathbf{E}^{TH}\left(\sum_{t=1}^{T}\mu_{i_{t}}\cdot
c_{i_{t},\,t}\,\Big{|}\,\lambda,\,\vec{c}_{T}\right)$
$\displaystyle-\mathbf{E}^{H}\left(\sum_{t=1}^{T}\mu_{i_{t}}\cdot
c_{i_{t},\,t}\,\Big{|}\,\lambda,\,\vec{c}_{T}\right)$ $\displaystyle=$
$\displaystyle\mathbf{E}^{TH}\left(\sum_{t=1}^{T}\mu_{i_{t}}\cdot
c_{i_{t},\,t}\,\Big{|}\,\lambda,\,\vec{c}_{T},\,\mathbb{H}_{T}\right)$
$\displaystyle-\mathbf{E}^{H}\left(\sum_{t=1}^{T}\mu_{i_{t}}\cdot
c_{i_{t},\,t}\,\Big{|}\,\lambda,\,\vec{c}_{T},\,\mathbb{H}_{T}\right)$
$\displaystyle=$
$\displaystyle\sum_{t=1}^{T}\mathbf{E}_{\lambda_{t}}\left[p_{t}^{i}(\lambda_{t})-\tilde{p}_{t}^{i}(\lambda)|\mathbb{H}_{t}\right]\mu_{i}\cdot
c_{i,\,t}.$ (9)

At any time $t\geq 1$, note that $\lambda_{t}$ is the sampled parameter from
the same posterior distribution as the true $\lambda$. Meanwhile,
$p_{t}^{i}(\cdot)$ and $\tilde{p}_{t}^{i}(\cdot)$ are the same deterministic
function based on $\mathbb{H}_{t}$. According to Lemma 2, it follows

$\mathbf{E}_{\lambda_{t},\,\lambda}\left[p_{t}^{i}(\lambda_{t})-\tilde{p}_{t}^{i}(\lambda)|\mathbb{H}_{t}\right]\mu_{i}\cdot
c_{i,\,t}=0.$

Therefore, for any fixed $\vec{c}_{T}\in\mathbb{C}_{T}$, it holds that

$\displaystyle\mathbf{E}_{\lambda\sim\mathcal{F}}\Bigg{[}\Bigg{\\{}\mathbf{E}^{TH}\left(\sum_{t=1}^{T}\mu_{i_{t}}\cdot
c_{i_{t},\,t}\,\Big{|}\,\lambda,\,\vec{c}_{T}\right)$
$\displaystyle-\mathbf{E}^{H}\left(\sum_{t=1}^{T}\mu_{i_{t}}\cdot
c_{i_{t},\,t}\,\Big{|}\,\lambda,\,\vec{c}_{T}\right)\Bigg{\\}}\Bigg{]}$
$\displaystyle=\mathbf{E}\Bigg{(}\mathbf{E}_{\lambda\sim\mathcal{F}}\Bigg{[}\Bigg{\\{}\mathbf{E}^{TH}\left(\sum_{t=1}^{T}\mu_{i_{t}}\cdot
c_{i_{t},\,t}\,\Big{|}\,\lambda,\,\vec{c}_{T},\,\mathbb{H}_{T}\right)$
$\displaystyle-\mathbf{E}^{H}\left(\sum_{t=1}^{T}\mu_{i_{t}}\cdot
c_{i_{t},\,t}\,\Big{|}\,\lambda,\,\vec{c}_{T},\,\mathbb{H}_{T}\right)\Bigg{\\}}\Bigg{]}\Bigg{)}$
$\displaystyle=\mathbf{E}\left(\mathbf{E}_{\lambda_{t},\,\lambda}\left[p_{t}^{i}(\lambda_{t})-\tilde{p}_{t}^{i}(\lambda)|\mathbb{H}_{t}\right]\mu_{i}\cdot
c_{i,\,t}\right)$ $\displaystyle=0.$

Note that the relation above holds for any $\vec{c}_{T}\in\mathbb{C}_{T}$,
which leads to

$\displaystyle\begin{split}\mathbf{E}_{\lambda\sim\mathcal{F}}\Bigg{[}\sup_{\vec{c}_{T}\in\mathbb{C}_{T}}\Bigg{\\{}&\mathbf{E}^{TH}\left(\sum_{t=1}^{T}\mu_{i_{t}}\cdot
c_{i_{t},\,t}\,\Big{|}\,\lambda,\,\vec{c}_{T}\right)\\\
&-\mathbf{E}^{H}\left(\sum_{t=1}^{T}\mu_{i_{t}}\cdot
c_{i_{t},\,t}\,\Big{|}\,\lambda,\,\vec{c}_{T}\right)\Bigg{\\}}\Bigg{]}\equiv
0.\end{split}$ (10)

d) Upper bound of the regret $\mathcal{R}^{TH}$.

Finally, the upper bound of the Bayesian sup regret of Thompson-Hedge
algorithm is obtained by combining the results in a), b) and c). In
particular, it holds that

$\displaystyle\mathcal{R}^{TH}\left(\mathbb{C}_{T},\,T\right)$
$\displaystyle\leq\Lambda_{1}+\Lambda_{2}=\Lambda_{1}+0$
$\displaystyle\leq\left(8+2\sqrt{2}\right)\left(m\mathcal{V}_{T}\ln{N}\right)^{1/3}\left(T\right)^{2/3}.$

Therefore, the proof is concluded.

######  Remark 1.

Relevant research that considers a similar problem can be found in Besbes et
al. (2019). In Besbes et al. (2019), the classical EXP3 type algorithm was
used and an upper bound of the order
$O\left(\left(m\mathcal{V}_{T}N\ln{N}\right)^{1/3}\left(T\right)^{2/3}\right)$
was obtained for the sup regret. Since the upper bound holds uniformly on the
parameter space, the same upper bound also holds for the Bayesian sup regret.
If the constant $\left(8+2\sqrt{2}\right)$ in (6) is neglected, our bound
outperforms the bound in Besbes et al. (2019) by a term of $O(N^{1/3})$, which
implies that the performance of our Thompson-Hedge algorithm is less sensitive
to the number of nodes $N$. It indicates that when considering problems with
large scales, our proposed algorithm is supposed to retain a Bayesian sup
regret that converges relatively faster. Meanwhile, Besbes et al. (2019)
constructs a lower bound of the order
$O\left(\mathcal{V}_{T}^{1/3}T^{2/3}\right)$ on the regret by any algorithm.
Similarly, the lower bound also holds uniformly on all the parameters and
adapts to our problem. Therefore, the lower bound shows that our algorithm
achieves the order optimality.

##  4 Applicability in smart grids

The model is formulated following the learning context in Section 2 and an
online learning algorithm is developed in Section 3. This section presents an
application case to demonstrate the applicability of the proposed model and
method in practical smart grids. In particular, Assumption 1 plays a central
role in our model. A real data set is used to verify that Assumption 1 may
hold in reality, and thus the proposed method is practical. Section 4.1
introduces the procedure for calculating the reward $c_{i,\,t}$ in the smart
grid. In Section 4.2, linear regression method is used to show a linear growth
rate of the critical quantity $\mathcal{V}_{T}$ in Assumption 1, and thus it
is concluded that Assumption 1 holds for the selected data set.

###  4.1 Operation cost of the smart grid

To facilitate reading, meanings of the variables used in calculating the
operation cost of the smart grid, are displayed in different categories as
follows. For the operational variables, $\vartheta_{t}$ is the operation state
of the smart grid at time $t$, $P_{i,j}^{t}$ [kW] is the output of type $j$
power source at Node $i$, $\bar{P}_{i,j}$ is the rate power of type $j$ power
source at Node $i$, $L_{i}^{t}$ is the power load at Node $i$,
$\bar{L}_{i}^{t}$ is the load shedding, $i.e.$, power demand not supplied, at
Node $i$, $\varphi$ is the state of charge of the energy storage system (ESS)
and $\bar{C}_{ESS}$ is the power capacity of the ESS. For the configuration
parameters of the smart grid, $B_{i,i^{{}^{\prime}}}$ [$1/\Omega$] is the
susceptance of the pairs of Nodes $(i,i^{{}^{\prime}})$, $\delta_{i}$ is the
voltage angle at Node $i$, $\Delta\delta$ is the voltage angle difference
$(\delta_{i}-\delta_{i^{{}^{\prime}}})$ of two nearby nodes,
$A_{i,i^{{}^{\prime}}}$ [A] is the ampacity of the pairs of Nodes
$(i,i^{{}^{\prime}})$ and $V$ [kW] is the nominal voltage of the smart grid.
For the cost and price coefficients, $Co^{\vartheta_{t}}$ is the operation
cost of the smart grid subject to operation state $\vartheta_{t}$, $Cs_{j}$ is
the variable operation cost for power source $j$, $Cf_{i,i^{{}^{\prime}}}$ is
the variable operation cost for feeder $(i,i^{{}^{\prime}})$, $Cp$ is the
penalty cost for power demand not supplied, and $Ep^{\vartheta_{t}}$ is the
energy price associated with operation state $\vartheta_{t}$. $T_{s}$ is the
duration of $\vartheta_{t}$, $\mathbb{P}$ denotes a subset of power sources,
$\mathbb{N}$ denotes the set of all nodes, and $\mathbb{F}$ denotes the set of
node pairs with transmission line between them. $s_{i,i^{{}^{\prime}}}=1$
indicates that power cannot be transmitted between Node $i$ and Node
$i^{{}^{\prime}}$ due to the successful cyber attack on Node $i$.

The linear Direct Current (DC) power flow model (Mo and Sansavini, 2019;
Sahraei-Ardakani and Hedman, 2016) is introduced to interprete the physical
meaning of $c_{i,\,t}$ in the smart grid. At time $t$, operation state of the
distributed generation system (DGS) is denoted by the following vector:

$\vartheta_{t}=\left[P_{i,\,j}^{t}\,L_{i}^{t}\right],$ (11)

where the power sources considered in this work consists of natural gas plant,
biomass plant, wind farm, Photovoltaic farm and ESS. These data are all
sampled from the historical database provided by the Elia Grid,
Belgium111(URL: https://www.elia.be/en/grid-data).

Before calculating $c_{i,\,t}$, the first objective is to achieve the minimal
cost in the presence of load shedding, denoted by $Co^{\vartheta_{t}}$, by
solving the following linear optimization problem

$\displaystyle\min Co^{\vartheta_{t}}=$
$\displaystyle\sum_{i\in\mathbb{N}}\sum_{j\in\mathbb{P}}\left(Cs_{j}-Ep^{\vartheta_{t}}\right)P_{i,\,j}^{t}$
$\displaystyle+\sum_{(i,\,i^{{}^{\prime}})\in\mathbb{F}}Cf_{i,\,i^{{}^{\prime}}}\left|B_{i,\,i^{{}^{\prime}}}\left(\delta_{i}-\delta_{i^{{}^{\prime}}}\right)\right|$
$\displaystyle+\left(Cp+Ep^{\vartheta_{t}}\right)\sum_{i=1}^{N}\bar{L}_{i}^{t}$
(12)

subject to

$L_{i}^{t}-\bar{L}_{i}^{t}-\sum_{j\in\mathbb{P}}P_{i,\,j}^{t}-\sum_{i=1}^{N}B_{i,\,i^{{}^{\prime}}}(\delta_{i}-\delta_{i^{{}^{\prime}}})=0,$
(13) $0\leq P_{i,\,j}\leq\bar{P}_{i,\,j}^{t},\>j\neq ESS,$ (14)
$\left|P_{i,\,j}^{t}\cdot
T_{s}\right|\leq\min\left\\{(1-\varphi)\cdot\bar{C}_{ESS},\,\bar{P}_{i}^{j}\cdot
T_{s}\right\\},\>j=ESS,$ (15)
$\left|B_{i,\,i^{{}^{\prime}}}(\delta_{i}-\delta_{i^{{}^{\prime}}})\right|\leq(1-s_{i,\,i^{{}^{\prime}}})VA_{i,\,i^{{}^{\prime}}},$
(16)

where constraint (13) requires that the power generated and consumed is
balancing at each node of DGS, constraint (14) requires that the power
generated should not be larger than the rated power, constraint (15) indicates
that the charging or discharging of ESS should not be larger than the
remaining capacity or nominal rate, and constraint (16) indicates that the
power flow between two nodes should not be larger than the capacity of the
transmission line. Eq. (4.1) is the objective function of a typical linear DC
Optimal Power Flow model, which aims to minimize the total operation cost
including the cost of generating power, the cost of running feeders and the
penalty cost of demand not supplied, when satisfying physical constraints
(13)-(16) of DGS.

The linear DC optimal power flow model is configurated in the Matlab and can
be solved using the well-known Simplex method, where the values of
configuration parameters have been given in Mo and Sansavini (2019) and the
computation time of operation cost for each practical operation state is
around 0.006 second in Gurobi. This computation time contributes most to the
total simulation time and therefore our proposed algorithm can be implemented
in real time.

The physical meaning of $c_{i,\,t}$ in the smart grid can be defined as the
difference between the operation cost of the DGS without cyber attacks and the
operation cost of the DGS given the Node $i^{*}$ is temporarily unavailable
caused by a successful cyber attack. Therefore, when probing Node $i$ at time
$t$, the reward function $c_{i,\,i}$ is defined as

$\begin{split}c_{i,\,t}\triangleq&\Big{|}Co^{\vartheta_{t}}\left(P_{i,j}^{t},\,L_{i}^{t},\,\Delta\delta;\,s_{i,\,i^{{}^{\prime}}}=0\right)\\\
&-Co^{\vartheta_{t}}\left(P_{i,j}^{t},\,L_{i}^{t},\,\Delta\delta,\,s_{i,\,i^{{}^{\prime}}}\,;\,s_{i,\,i^{{}^{\prime}}}=1\right)\Big{|},\end{split}$
(17)

where
$Co^{\vartheta_{t}}\left(P_{i,j}^{t},\,L_{i}^{t},\,\Delta\delta;\,s_{i,\,i^{{}^{\prime}}}=0\right)$
is the operation cost of the DGS without cyber attacks, and
$Co^{\vartheta_{t}}\left(P_{i,j}^{t},\,L_{i}^{t},\,\Delta\delta;\,s_{i,\,i^{{}^{\prime}}}=1\right)$
is the operation cost of the DGS given Node $i$ is unavailable caused by the
successful cyber attack on Node $i$. That is to say, the reward $c_{i,t}$ can
be calculated via solving the linear optimization problem, defined by (4.1)
through (16) (with/without cyber attacks), where the input $\vartheta_{t}$ is
drawn from the dataset of Elia Grid, Belgium.

###  4.2 Numerical analysis of cost sequences

In this section, numerical analysis is presented based on a real data set to
verify that Assumption 1 holds in reality. Note that Assumption 1 implies a
linear or sub-linear upper bound in terms of $T$ on $\mathcal{V}_{T}$.
Therefore, if $\mathcal{V}_{T}$ has a linear or sub-linear growth rate in $T$,
then Assumption 1 is supposed to hold by choosing a proper value for $m$.
Thus, linear regression is performed on the sequence of $\mathcal{V}_{T}$
against time $T$. The realistic grid data from the Elia is used, which
provides data from the Belgian electricity market system. The underlying
electricity network (a subgrid of Elia Grid) is shown in Figure 1. In the DGS,
the dataset of Elia Grid, Belgium is recorded every 15 minutes, which means
that the optimal power flow model will be performed and generate one attack
cost $c_{i,t}$ every 15 minutes. The DGS under cyber attacks is investigated
over one week, which indicates that a specific dataset with 672 successive
observations of the attack cost is selected for illustration. $m$ is chosen to
be $4$ and the data has been normalized such that $c_{i,\,t}\in[0,\,1/4]$. The
scatter plot is presented in Figure 2.

Figure 1: The $11$-node radial DGS originated from the IEEE 13 node test
feeder. Figure 2: Data scatter plot and linear prediction of
$\mathcal{V}_{T}$.

From Figure 2, a significant linear relation between $\mathcal{V}_{T}$ and $T$
can be observed. For further justification, the numerical results of the
linear regression are presented in Table 1.

Table 1: Regression results. | Coefficients | Standard error | t Stat | P-value
---|---|---|---|---
Intercept | 1.560 | 0.04711 | 32.04 | $<0.01$
$T$ | 0.02648 | 0.0001212 | 218.5 | 0

According to the results in Table 1, the coefficient for $T$ is about $0.03$,
which means that Assumption 1 holds with $T_{0}$ chosen to be $9$. As a
significance test for the regression model, results of ANOVA for the
regression model is shown in Table 2.

Table 2: ANOVA results. | df | SS | MS | F | Significance
---|---|---|---|---|---
Regression | 1 | 17655 | 17655 | 47734 |  $\leq$0.01
Residual | 669 | 247.4 | 0.3699 |  |
Total | 670 | 17902.4 |  |  |

From the results in Table 2, the significance level is below $0.01$, which
means that the linear relation between $\mathcal{V}_{T}$ and $T$ for this data
set is significant and the model is justified.

##  5 Numerical results

In this section, the performance of our Thompson-Hedge algorithm is
illustrated through simulation studies. In Section 5.1, a comparative study is
conducted between the performances of our Thompson-Hedge algorithm and the
R.EXP3 algorithm which was recently proposed in Besbes et al. (2019).
Subsequently, in Section 5.2, sensitivity analysis is conducted to investigate
the influence of constraint $\mathcal{V}_{T}$ on the cost sequence.

###  5.1 Comparative study

This section compares the performance of the proposed Thompson-Hedge algorithm
with the R.EXP3 algorithm that was also designed for the same adversarial
problem with constrained variation on cost sequence. Both algorithms use batch
methods from the original EXP3 algorithm, $i.e.$, the time horizon is divided
into small batches and the algorithm “restarts” at the beginning of each
batch. Meanwhile, the R.EXP3 algorithm is also a randomized algorithm, retains
a weight function for each node and updates the weights each time. Different
from our Thompson-Hedge algorithm, in the R.EXP3 algorithm, if node
$i\in\mathbb{N}$ is chosen at any time, the total reward $K_{i,\,t}\cdot
c_{i,\,t}$ is considered as a whole adversarial reward. Since the attack
number $k_{i^{{}^{\prime}},\,t}$ is unknown if $i\neq i^{{}^{\prime}}$, the
R.EXP3 is supposed to ignore $c_{i^{{}^{\prime}},\,t}$ and treat the problem
as a typical MNB model with bandit feedback.

The simulation is executed under $N=10$ and $N=20$. To initialize, the
parameter
$\theta=\left(\alpha_{1},\,\beta_{1},\cdots,\alpha_{N},\,\beta_{N}\right)$ is
set for all $i\in\mathbb{N}$, and two large numbers for $Q$ and $L$. The
preset parameter $\theta$ is used to generate
$\lambda=(\lambda_{1},\cdots,\lambda_{N})$ for $Q$ times in total. Under each
generated $\lambda$, $\lambda$ is used to generate the sequence
$\\{K_{i,\,t}\\}_{t=1}^{T}$, and the adversarial cost sequence
$\\{c_{i,\,t}\\}_{t=1}^{T}$ is artificially generated. Then, both algorithms
are run based on $\\{K_{i,\,t}\\}_{t=1}^{T}$ and $\\{c_{i,\,t}\\}_{t=1}^{T}$,
$\forall i\in\mathbb{N}$. Finally, the Bayesian sup regret in the two
algorithms are calculated and compared. The simulation process is summarized
using pseudo codes in Algorithm 3. Experience from the previous work Smith and
Paté-Cornell (2018) and the real database from our numerical example are
integrated to determine the parameter vector $\theta$. In Smith and Paté-
Cornell (2018), the prior distribution is set to be Gamma$(2,\,2)$; in our
numerical example. It is estimated that $\lambda\approx 0.25$, which may
correspond to the prior distribution Gamma$(1,\,4)$. Therefore, two values for
$\theta$ are selected to implement the simulation:
$\theta_{1}:\>\left(\alpha_{1},\,\beta_{1}\right)=\cdots=\left(\alpha_{N},\,\beta_{N}\right)=(2,\,2)$
and
$\theta_{2}:\>\left(\alpha_{1},\,\beta_{1}\right)=\cdots=\left(\alpha_{N},\,\beta_{N}\right)=(1,\,4)$.
According to the choices of $\theta$, the truncation parameter $m$ is set as
$3$, which ensures that the probability $Pr(K_{i,\,t}>m)$ is small.

Algorithm 4 Algorithm 3: Simulation for comparison

1:Preset parameter
$\theta=(\alpha_{1},\,\beta_{1},\cdots,\alpha_{N},\,\beta_{N})$; numbers of
simulation trials $Q$, $L$; time horizon $T$.

2:for $q=1:Q$ do

3: For all $i\in\mathbb{N}$, generate $\lambda_{i}$ from
$Gamma(\alpha_{i},\,\beta_{i})$.

4: for $l=1:L$ do

5: for $t=1:T$ do

6: for $i\in\mathbb{N}$ do

7: Generate $K_{i,\,t}$ from $Poisson(\lambda_{i})$.

8: Generate $c_{i,\,t}$ artificially.

9: end for

10: end for

11: Run Thompson-Hedge algorithm and R.EXP3 algorithm independently based on
the sequences $\\{K_{i,\,t}\\}_{t=1}^{T}$ and $\\{c_{i,\,t}\\}_{t=1}^{T}$,
$\forall i\in\mathbb{N}$.

12: Calculate the regret function
$R_{l}^{TH}\left(\lambda,\,\vec{c}_{T},\,T\right)$ for Thompson-Hedge
algorithm and $R_{l}^{R3}\left(\lambda,\,\vec{c}_{T},\,T\right)$ for R.EXP3
algorithm.

13: end for

14: Calculate the sup regret $\tilde{R}_{q}^{TH}=\max_{l=1:L}R_{l}^{TH}$ for
Thompson-Hedge algorithm and $\tilde{R}_{q}^{R3}=\max_{l=1:L}R_{l}^{R3}$ for
R.EXP3 algorithm.

15:end for

16:Calculate the Bayesian sup regret
$\mathcal{R}^{TH}\left(\mathbb{C}_{T},\,T\right)=1/Q\sum_{q=1}^{Q}\tilde{R}_{q}^{TH},$
for Thompson-Hedge algorithm and,
$\mathcal{R}^{R3}\left(\mathbb{C}_{T},\,T\right)=1/Q\sum_{q=1}^{Q}\tilde{R}_{q}^{R3},$
for R.EXP3 algorithm.

Note that in the simulation, the empirical estimation of the sup regret and
Bayesian sup regret functions is used in the two algorithms. Therefore, it is
necessary to set the values of $Q$ and $L$ large enough to make the estimation
with good precision. In the simulation, it is set as $Q=L=10^{5}$. To simulate
the real situations, the procedures of sketching the cost sequence are given
as below:

Step 1: For all $i\in\mathbb{N}$, generate $c_{i,\,1}$ uniformly on
$(0,\,1/m)$ independently.

Step 2: For $t=2,\cdots,T$, generate $c_{i,\,t}$ independently from the
uniform distribution on the overlapping interval between
$(c_{i,\,t-1}-1/100m,\,c_{i,\,t-1}+1/100m)$ and $(0,\,1/m)$.

By Step 1, simulate the random initial value of each cost sequence. Then, by
Step 2, set $|c_{i,\,t}-c_{i,\,t-1}|\leq 1/50m$. It can be verified that any
cost sequence generated by Step 1 and Step 2 satisfies Assumption 1 with
$T_{0}=50$. The simulation results under
$\left(\alpha_{1},\,\beta_{1}\right)=\cdots=\left(\alpha_{N},\,\beta_{N}\right)=(2,\,2)$
are given in Figures 3 and 4.

Figure 3: Comparison of regrets ($\theta=\theta_{1}$, $N=10$). Figure 4:
Comparison of regrets ($\theta=\theta_{1}$, $N=20$).

Using the same process, simulation is implemented for
$\left(\alpha_{1},\,\beta_{1}\right)=\cdots=\left(\alpha_{N},\,\beta_{N}\right)=(1,\,4)$.
The results are given in Figures 5 and 6.

Figure 5: Comparison of regrets ($\theta=\theta_{2}$, $N=10$). Figure 6:
Comparison of regrets ($\theta=\theta_{2}$, $N=20$).

In Figures 3-6, our Thompson-Hedge algorithm outperforms the existing R.EXP3
algorithm in terms of Bayesian sup regret. In particular, Thompson-Hedge
algorithm is less sensitive to the problem scale $N$ by comparing the regret
curves under $N=10$ and $N=20$, which is consistent with our discussion in
Remark 1. Our algorithm has advantages over a typical algorithm designed for
MNB with bandit feedback, such as the EXP3 or R.EXP3 algorithm. The usual
convergence order of the regret function by a typical algorithm is
$O(N\log{N})^{1/3}$ for MNB with bandit feedback in terms of the problem scale
$N$. Our Thompson-Hedge algorithm feeds a set of sampled parameters to Hedge
algorithm. In the proof of Theorem 1, it is shown that using the sampled
parameters is “as good as” using the true parameters under the Bayesian
framework. Moreover, if the true parameters are known, then Hedge algorithm
can solve the problem, which has a convergence order of $O(\log{N})$.
Therefore, our Thompson-Hedge algorithm is supposed to retain a convergence
rate of $O(\log{N})$ for the Bayesian sup regret in the special case.

Meanwhile, upper bound in Theorem 1 is also drawn in Figures 3-6. It is worth
noting that the upper bound in Theorem 1 holds uniformly on all possible
values of $\theta$. Hence, it is concluded that when $\theta$ varies, the
deviation of the regret from the upper bound also varies. In particular, as
shown in the figures, the difference between the upper bound and the regret is
comparatively large under $\theta=\theta_{1}$, $N=20$ while small or moderate
under other three parameters. However, it can be observed from the figures
that our upper bound is obviously sharper than that for the R.EXP3 algorithm
in the sense that our upper bound is below the regret under R.EXP3 at some
time points in all the four figures.

###  5.2 Sensitivity to the variation

This section presents the sensitivity analysis of our Thompson-Hedge algorithm
to the variation constraint $\mathcal{V}_{T}$. In reality, according to the
specific environment and workload under which the network functions, the
adversary may generate the cost sequence subject to various rules. In our
model, the variation constraint $\mathcal{V}_{T}$ is used to characterize the
cost sequence. Thus, numerical results for the proposed Thompson-Hedge
algorithm are presented under four levels of variation. In particular, the
same steps given in Section 5.1 are employed, but $4$ levels are selected for
the variation scale in Step 2: $1/20m$, $1/50m$, $1/100m$ and $1/200m$. Other
parameters are set to be identical as those in Section 5.1. The simulation
results are illustrated under both $\theta_{1}$ and $\theta_{2}$ as shown in
Figures 7 and 8.

Figure 7: Regrets under different variations ($\theta=\theta_{1}$, $N=20$).
Figure 8: Regrets under different variations ($\theta=\theta_{2}$, $N=20$).

The results in Figures 7 and 8 illustrate the different performances of our
Thompson-Hedge algorithm under different variation scales. The Bayesian regret
monotonically increases with the variation. This is because larger variation
leads to more uncertainties of the cost sequence and makes the experience from
historical observations unreliable. Thus, as an adaptive online learning
algorithm, the Bayesian sup regret in our Thompson-Hedge algorithm may
increase and the performance may have fluctuations. However, it should be
noted that the incremental difference between the regret curves under $1/100m$
and $1/200m$ are not so significant compared to that between the curves under
$1/50m$ and $1/20m$. This implies that when the variation is within
$\mathcal{V}_{T}\leq T/50m$, the regret function is less sensitive to the
variation than the situation that $\mathcal{V}_{T}$ is around $1/10m$. In the
numerical example, the real data set (from Elia Grid, Belgium) is examined and
it appears that the corresponding $\mathcal{V}_{T}$ is around $T/100m$.
Therefore, it can be concluded that the sensitivity of our Thompson-Hedge
algorithm to the variation may not be high in practice.

##  6 Concluding remarks

This study investigates the sequential control problem in modern energy and
power systems with smart grids. The existing work is extended by introducing
an adversarial cost sequence with a variation constraint. A Bayesian MNB model
is constructed to cope with the problem and an online learning algorithm named
Thompson-Hedge algorithm is proposed to retain a converging regret function.
In addition, it is proved that the convergence rate of the regret in the
proposed algorithm is superior to an existing algorithm that can be used for
this problem.

In developing our algorithm, a basic idea is to firstly obtain sampled
parameters from the posterior distribution when the reinforcement learning
model is a partly parameterized Bayesian model. Subsequently, the sampled
parameters are used instead of the true parameters in the model. It is worth
pointing out that although specific models are used in this work, the paradigm
of our algorithm framework can also be used for other classical models.

From a conceptual perspective, future research will explore the applications
of the proposed algorithms to more situations in reality. For example, the
proposed algorithm can be used to mitigate the impact of failure in tracking
goods, or obtain suppliers to exchange inventory information due to successful
cyber attacks. Other potential applications can be found in the domains of
multi-robot systems (Liu et al., 2017), Internet of Things (IOT) (Perera et
al., 2015), $e.g.$, traffic management, water distribution, waste management,
online manufacturing and smart supply chain. It is also worth investigating
the impact of other types of cyber attacks such as the botnet command and
control, data exfiltration, data tampering, data destruction or even physical,
destruction via alternation of critical software/data (Pillitteri and Brewer,
2014) on the smart grid by properly modeling their influences on the Optimal
Power Flow model according to their attack mechanisms. For example, the replay
attacks maliciously repeat DR messages or dispatch commands, and the false
data injection attacks purposely alter the integrity of a smart grid by
compromising a subset of transmitted data packages and sending out inaccurate
DR messages or dispatch commands.

##  Appendix Appendix A Proof of Lemma 1

In this section, detailed proof of Lemma 1 is provided, which is used in the
proof of Theorem 1. Before reaching the proof of Lemma 1, a preliminary result
is presented in the following lemma Slivkins et al. (2019).

######  Lemma 3.

[Slivkins et al. (2019), Theorem 5.16] If all per-time costs are in $[0,\,1]$
and $\varepsilon$ is chosen to be $\left(1-\sqrt{\ln{N}/2T}\right)^{-1}$, then
the Hedge algorithm satisfies

$\displaystyle\max_{i\in\mathbb{N}}\left\\{\sum_{t=\tau_{b-1}}^{\tau_{b}-1}K_{i,\,t}\cdot
c_{i,\,t}\right\\}$
$\displaystyle-\mathbf{E}^{H}\left(\sum_{t=\tau_{b-1}}^{\tau_{b}-1}K_{i_{t},\,t}\cdot
c_{i_{t},\,t}\,\Big{|}\,\lambda,\,\vec{c}_{T}\right)$ $\displaystyle\leq
2\sqrt{2}\cdot\sqrt{\Delta\ln{N}}.$

Note that the original theorem in Slivkins et al. (2019) concerns a reward-
based problem and the parameter $\varepsilon$ is chosen to be
$\sqrt{\ln{N}/2T}$. Actually, by using the transformation of $cost=1-reward$
between the per-time cost and reward, the result can be transferred in the
original theorem into the conclusion in Lemma 3. Detailed proof of Lemma 1 is
given as follows.

Proof of Lemma 1. The proof is divided into four steps. In the first step, the
best single action policy is used to decompose the regret function into two
parts. In the second and third steps, upper bounds on the two parts from the
first step are constructed, respectively. Finally, the upper bound on the
original regret function is presented. Denote the sub-sequence of times when
Hedge algorithm restarts as
$1=\tau_{0}\leq\tau_{1}\leq\cdots\leq\tau_{\left\lceil T/\Delta\right\rceil}$,
which implies that $\tau_{b+1}-\tau_{b}=\Delta$ for all
$b=1,\cdots,\left\lceil T/\Delta\right\rceil-1$ and $\tau_{\left\lceil
T/\Delta\right\rceil}-\tau_{\left\lceil T/\Delta\right\rceil-1}\leq\Delta$.
For all $b=1,\cdots,\left\lceil T/\Delta\right\rceil$, note that the set of
times $\left\\{\tau_{b-1},\cdots,\tau_{b}-1\right\\}$ exactly forms the $b$th
batch.

Step 1 Regret decomposition.
For each $b=1,\cdots,\left\lceil T/\Delta\right\rceil$, the best single node
$I_{b}^{*}$ in $b$th batch is defined as
$I_{b}^{*}\triangleq{\arg\max}_{i\in\mathbb{N}}\left\\{\sum_{t=\tau_{b-1}}^{\tau_{b}-1}\mu_{i}\cdot
c_{i,\,t}\right\\}$. A specific sequence
$\\{c_{i,\,t}\\}_{t=1}^{T}\in\mathbb{C}_{T}$ is fixed and the regret function
within $b$th batch is decomposed as

$\displaystyle\begin{split}&\sup_{\vec{c}_{T}\in\mathbb{C}_{T}}\sum_{t=\tau_{b-1}}^{\tau_{b}-1}\mu_{i_{t}^{*}}\cdot
c_{i_{t}^{*},\,t}-\mathbf{E}^{H}\left(\sum_{t=\tau_{b-1}}^{\tau_{b}-1}\mu_{i_{t}}\cdot
c_{i_{t},\,t}\,\Big{|}\,\lambda,\,\vec{c}_{T}\right)\\\
&\leq\sup_{\vec{c}_{T}\in\mathbb{C}_{T}}\left\\{\sum_{t=\tau_{b-1}}^{\tau_{b}-1}\mu_{i_{t}^{*}}\cdot
c_{i_{t}^{*},\,t}-\sum_{t=\tau_{b-1}}^{\tau_{b}-1}\mu_{I_{b}^{*}}\cdot
c_{I_{b}^{*},\,t}\right\\}\\\
&+\sup_{\vec{c}_{T}\in\mathbb{C}_{T}}\Bigg{\\{}\sum_{t=\tau_{b-1}}^{\tau_{b}-1}\mu_{I_{b}^{*}}\cdot
c_{I_{b}^{*},\,t}\\\
&-\mathbf{E}^{H}\left(\sum_{t=\tau_{b-1}}^{\tau_{b}-1}\mu_{i_{t}}\cdot
c_{i_{t},\,t}\,\Big{|}\,\lambda,\,\vec{c}_{T}\right)\Bigg{\\}}.\end{split}$
(A.18)

The two terms in the right hand side of (A.18) are defined as

$R_{b,\,1}\triangleq\sup_{\vec{c}_{T}\in\mathbb{C}_{T}}\left\\{\sum_{t=\tau_{b-1}}^{\tau_{b}-1}\mu_{i_{t}^{*}}\cdot
c_{i_{t}^{*},\,t}-\sum_{t=\tau_{b-1}}^{\tau_{b}-1}\mu_{I_{b}^{*}}\cdot
c_{I_{b}^{*},\,t}\right\\},$

and

$\displaystyle R_{b,\,2}\triangleq$
$\displaystyle\sup_{\vec{c}_{T}\in\mathbb{C}_{T}}\Bigg{\\{}\sum_{t=\tau_{b-1}}^{\tau_{b}-1}\mu_{I_{b}^{*}}\cdot
c_{I_{b}^{*},\,t}$
$\displaystyle-\mathbf{E}^{H}\left(\sum_{t=\tau_{b-1}}^{\tau_{b}-1}\mu_{i_{t}}\cdot
c_{i_{t},\,t}\,\Big{|}\,\lambda,\,\vec{c}_{T}\right)\Bigg{\\}}.$

Step 2 Upper bounding $R_{b,\,1}$.
Let
$\mathcal{V}_{T}^{b}\triangleq\sum_{t=\tau_{b-1}}^{\tau_{b}-2}\left|c_{i,\,t+1}-c_{i,\,t}\right|$
be the total variation of $c_{i,\,t}$ in $b$th batch. For all
$\vec{c}_{T}\in\mathbb{C}_{T}$, the following relation holds

$\max_{\tau_{b-1}\leq t\leq\tau_{b}-1}\left\\{\mu_{i_{t}^{*}}\cdot
c_{i_{t}^{*},\,t}-\mu_{I_{b}^{*}}\cdot c_{I_{b}^{*},\,t}\right\\}\leq
2m\mathcal{V}_{T}^{b}.$ (A.19)

By contradiction, if (A.19) does not hold, then there is at least one time
$\tau_{b-1}\leq t_{0}\leq\tau_{b}-1$ such that $\mu_{i_{t}^{*}}\cdot
c_{i_{t}^{*},\,t_{0}}-\mu_{I_{b}^{*}}\cdot c_{I_{b}^{*},\,t_{0}}\geq
2m\mathcal{V}_{T}^{b}$. Let $i_{t_{0}}^{*}=i_{0}$. Since $\mu_{i}\leq m$,
$\forall i\in\mathbb{N}$, for all $\tau_{b-1}\leq t\leq\tau_{b}-1$, it follows

$\displaystyle\mu_{i_{0}}\cdot c_{i_{0},\,t}$
$\displaystyle>\mu_{i_{0}}\left(c_{i_{0},\,t_{0}}-\mathcal{V}_{T}^{b}\right)>\mu_{i_{0}}\cdot
c_{i_{0},\,t_{0}}-m\mathcal{V}_{T}^{b}$ $\displaystyle>\mu_{I_{b}^{*}}\cdot
c_{I_{b}^{*},\,t_{0}}+m\mathcal{V}_{T}^{b}>\mu_{I_{b}^{*}}\left(c_{I_{b}^{*},\,t_{0}}+\mathcal{V}_{T}^{b}\right)$
$\displaystyle>\mu_{I_{b}^{*}}\cdot c_{I_{b}^{*},\,t}.$

However, this contradicts with the fact that $I_{b}^{*}$ is the optimal single
node of $b$th batch. Thus, inequality (A.19) holds. Therefore, it can be
obtained that

$\displaystyle R_{b,\,1}$
$\displaystyle=\sup_{\vec{c}_{T}\in\mathbb{C}_{T}}\sum_{t=\tau_{b-1}}^{\tau_{b}-1}\left(\mu_{i_{t}^{*}}\cdot
c_{i_{t}^{*},\,t}-\mu_{I_{b}^{*}}\cdot c_{I_{b}^{*},\,t}\right)$
$\displaystyle\leq\Delta\sup_{\vec{c}_{T}\in\mathbb{C}_{T}}\max_{\tau_{b-1}\leq
t\leq\tau_{b}-1}\left\\{\mu_{i_{t}^{*}}\cdot
c_{i_{t}^{*},\,t}-\mu_{I_{b}^{*}}\cdot c_{I_{b}^{*},\,t}\right\\}$
$\displaystyle\leq 2\Delta m\mathcal{V}_{T}^{b}.$

Step 3 Upper bounding $R_{b,\,2}$.
Note that

$\displaystyle\sum_{t=\tau_{b-1}}^{\tau_{b}-1}\mu_{I_{b}^{*}}\cdot
c_{I_{b}^{*},\,t}$
$\displaystyle=\sum_{t=\tau_{b-1}}^{\tau_{b}-1}\mathbf{E}\left(K_{I_{b}^{*},\,t}\cdot
c_{I_{b}^{*},\,t}\right)$
$\displaystyle=\max_{i\in\mathbb{N}}\left\\{\sum_{t=\tau_{b-1}}^{\tau_{b}-1}\mathbf{E}\left(K_{i,\,t}\cdot
c_{i,\,t}\right)\right\\}$
$\displaystyle\leq\mathbf{E}\left(\max_{i\in\mathbb{N}}\left\\{\sum_{t=\tau_{b-1}}^{\tau_{b}-1}K_{i,\,t}\cdot
c_{i,\,t}\right\\}\right).$

Therefore, it holds that

$\displaystyle R_{b,\,2}$
$\displaystyle\leq\sup_{\vec{c}_{T}\in\mathbb{C}_{T}}\Bigg{\\{}\mathbf{E}\Bigg{[}\left(\max_{i\in\mathbb{N}}\left\\{\sum_{t=\tau_{b-1}}^{\tau_{b}-1}K_{i,\,t}\cdot
c_{i,\,t}\right\\}\right)$
$\displaystyle-\mathbf{E}^{H}\left(\sum_{t=\tau_{b-1}}^{\tau_{b}-1}\mu_{i_{t}}\cdot
c_{i_{t},\,t}\,\Big{|}\,\lambda,\,\vec{c}_{T}\right)\Bigg{]}\Bigg{\\}}$
$\displaystyle\leq\sup_{\vec{c}_{T}\in\mathbb{C}_{T}}\Bigg{\\{}\mathbf{E}\Bigg{[}\left(\max_{i\in\mathbb{N}}\left\\{\sum_{t=\tau_{b-1}}^{\tau_{b}-1}K_{i,\,t}\cdot
c_{i,\,t}\right\\}\right)$
$\displaystyle-\mathbf{E}^{H}\left(\sum_{t=\tau_{b-1}}^{\tau_{b}-1}k_{i_{t},\,t}\cdot
c_{i_{t},\,t}\,\Big{|}\,\lambda,\,\vec{c}_{T}\right)\Bigg{]}\Bigg{\\}}.$

According to Lemma 1 and noting that the relation in Lemma 1 holds for
arbitrary $\vec{c}_{T}$ satisfying $K_{i,\,t}\cdot c_{i,\,t}\leq 1$, $\forall
i,\,t$, it follows

$\displaystyle R_{b,\,2}\leq 2\sqrt{2}\cdot\sqrt{\Delta\ln{N}}.$

Step 4 Upper bounding the sup regret.
In the final step, discussions in the previous two steps and the upper bound
on the sup regret by Hedge algorithm are summarized. Note that

$\displaystyle\sup_{\vec{c}_{T}\in\mathbb{C}_{T}}R^{H}\left(\lambda,\,\vec{c}_{T},\,T\right)\leq\sum_{b=1}^{\left\lceil
T/\Delta\right\rceil}\left(R_{b,\,1}+R_{b,\,2}\right)$
$\displaystyle\leq\sum_{b=1}^{\left\lceil T/\Delta\right\rceil}\left(2\Delta
m\mathcal{V}_{T}^{b}+2\sqrt{2}\cdot\sqrt{\Delta\ln{N}}\right)$
$\displaystyle\leq 2\Delta
m\mathcal{V}_{T}+\left(\frac{T}{\Delta}+1\right)2\sqrt{2}\cdot\sqrt{\Delta\ln{N}},$
(A.20)

where the third inequality follows from the relation $\sum_{b=1}^{\left\lceil
T/\Delta\right\rceil}\mathcal{V}_{T}^{b}=\mathcal{V}_{T}$.

Embedding
$\Delta=\left\lceil\left(\ln{N}\right)^{1/3}\left({T}/{m\mathcal{V}_{T}}\right)^{2/3}\right\rceil$
into (Appendix A), it follows

$\displaystyle\sup_{\vec{c}_{T}\in\mathbb{C}_{T}}R^{\pi}\left(\lambda,\,\vec{c}_{T},\,T\right)$
$\displaystyle\leq
2\left[\left(\ln{N}\right)^{1/3}\left(\frac{T}{m\mathcal{V}_{T}}\right)^{2/3}+1\right]m\mathcal{V}_{T}$
$\displaystyle+2\sqrt{2}\cdot
T\cdot\sqrt{\ln{N}\left(\ln{N}\right)^{-1/3}\left(\frac{T}{m\mathcal{V}_{T}}\right)^{-2/3}}$
$\displaystyle+2\sqrt{2}\sqrt{\ln{N}\left[\left(\ln{N}\right)^{1/3}\left(\frac{T}{m\mathcal{V}_{T}}\right)^{2/3}+1\right]}.$

Since $N>2$, Assumption 1 leads to
$\left(\ln{N}\right)^{1/3}\left({T}/{m\mathcal{V}_{T}}\right)^{2/3}>1$ and
$\left[\ln{N}/T(m\mathcal{V}_{T})^{2}\right]<1$. It follows

$\displaystyle\sup_{\vec{c}_{T}\in\mathbb{C}_{T}}R^{\pi}\left(\lambda,\,\vec{c}_{T},\,T\right)$
$\displaystyle\leq\left(4+2\sqrt{2}\right)\left(m\mathcal{V}_{T}\ln{N}\right)^{1/3}\left(T\right)^{2/3}$
$\displaystyle+4\left(\ln{N}\right)^{2/3}\left(\frac{T}{m\mathcal{V}_{T}}\right)^{1/3}$
$\displaystyle\leq\left(4+2\sqrt{2}\right)\left(m\mathcal{V}_{T}\ln{N}\right)^{1/3}\left(T\right)^{2/3}$
$\displaystyle+4\left(m\mathcal{V}_{T}\ln{N}\right)^{1/3}\left(T\right)^{2/3}\left[\frac{\ln{N}}{T\left(m\mathcal{V}_{T}\right)^{2}}\right]^{1/3}$
$\displaystyle\leq\left(4+2\sqrt{2}\right)\left(m\mathcal{V}_{T}\ln{N}\right)^{1/3}\left(T\right)^{2/3}$
$\displaystyle+4\left(m\mathcal{V}_{T}\ln{N}\right)^{1/3}\left(T\right)^{2/3}$
$\displaystyle\leq\left(8+2\sqrt{2}\right)\left(m\mathcal{V}_{T}\ln{N}\right)^{1/3}\left(T\right)^{2/3},$

which concludes the proof.

## References

  * (1)
  * Abiri-Jahromi et al. (2020) Abiri-Jahromi, A., Kemmeugne, A., Kundur, D. and Haddadi, A. (2020), ‘Cyber-physical attacks targeting communication-assisted protection schemes’, IEEE Transactions on Power Systems 35(1), 440–450.
  * Amin et al. (2013) Amin, S., Schwartz, G. A. and Sastry, S. S. (2013), ‘Security of interdependent and identical networked control systems’, Automatica 49(1), 186–192.
  * Auer et al. (2002) Auer, P., Cesa-Bianchi, N., Freund, Y. and Schapire, R. E. (2002), ‘The nonstochastic multiarmed bandit problem’, SIAM journal on computing 32(1), 48–77.
  * Besbes et al. (2019) Besbes, O., Gur, Y. and Zeevi, A. (2019), ‘Optimal exploration–exploitation in a multi-armed bandit problem with non-stationary rewards’, Stochastic Systems 9(4), 319–337.
  * Che et al. (2018) Che, L., Liu, X., Shuai, Z., Li, Z. and Wen, Y. (2018), ‘Cyber cascades screening considering the impacts of false data injection attacks’, IEEE Transactions on Power Systems 33(6), 6545–6556.
  * Even-Dar et al. (2009) Even-Dar, E., Kakade, S. M. and Mansour, Y. (2009), ‘Online markov decision processes’, Mathematics of Operations Research 34(3), 726–736.
  * Freund and Schapire (1999) Freund, Y. and Schapire, R. E. (1999), ‘Adaptive game playing using multiplicative weights’, Games and Economic Behavior 29(1-2), 79–103.
  * Gallo et al. (2020) Gallo, A. J., Turan, M. S., Boem, F., Parisini, T. and Ferrari-Trecate, G. (2020), ‘A distributed cyber-attack detection scheme with application to dc microgrids’, IEEE Transactions on Automatic Control DOI: 10.1109/TAC.2020.2982577.
  * Heyman and Sobel (2004) Heyman, D. and Sobel, M. (2004), ‘Superposition of renewal processes’, Stochastic Models in Operations Research: Stochastic Processes and Operating Characteristics p. 158.
  * Khalili et al. (2020) Khalili, M., Zhang, X., Cao, Y., Polycarpou, M. M. and Parisini, T. (2020), ‘Distributed fault-tolerant control of multiagent systems: An adaptive learning approach’, IEEE Transactions on Neural Networks and Learning Systems 31(2), 420–432.
  * Konstantelos et al. (2016) Konstantelos, I., Giannelos, S. and Strbac, G. (2016), ‘Strategic valuation of smart grid technology options in distribution networks’, IEEE Transactions on Power Systems 32(2), 1293–1303.
  * Li et al. (2020) Li, J.-A., Dong, D., Wei, Z., Liu, Y., Pan, Y., Nori, F. and Zhang, X. (2020), ‘Quantum reinforcement learning during human decision-making’, Nature Human Behaviour 4(3), 294–307.
  * Littman (2015) Littman, M. L. (2015), ‘Reinforcement learning improves behaviour from evaluative feedback’, Nature 521(7553), 445–451.
  * Liu et al. (2017) Liu, Z., Wang, L., Wang, J., Dong, D. and Hu, X. (2017), ‘Distributed sampled-data control of nonholonomic multi-robot systems with proximity networks’, Automatica 77, 170–179.
  * Mo and Sansavini (2017) Mo, H. and Sansavini, G. (2017), ‘Dynamic defense resource allocation for minimizing unsupplied demand in cyber-physical systems against uncertain attacks’, IEEE Transactions on Reliability 66(4), 1253–1265.
  * Mo and Sansavini (2019) Mo, H. and Sansavini, G. (2019), ‘Impact of aging and performance degradation on the operational costs of distributed generation systems’, Renewable energy 143, 426–439.
  * Osband et al. (2013) Osband, I., Russo, D. and Van Roy, B. (2013), (more) efficient reinforcement learning via posterior sampling, in ‘Advances in Neural Information Processing Systems’, pp. 3003–3011.
  * Patsakis et al. (2018) Patsakis, G., Rajan, D., Aravena, I., Rios, J. and Oren, S. (2018), ‘Optimal black start allocation for power system restoration’, IEEE Transactions on Power Systems 33(6), 6766–6776.
  * Perera et al. (2015) Perera, C., Liu, C. H. and Jayawardena, S. (2015), ‘The emerging internet of things marketplace from an industrial perspective: A survey’, IEEE Transactions on Emerging Topics in Computing 3(4), 585–598.
  * Pillitteri and Brewer (2014) Pillitteri, V. Y. and Brewer, T. L. (2014), Guidelines for smart grid cybersecurity, Technical report.
  * Pogaku et al. (2007) Pogaku, N., Prodanovic, M. and Green, T. C. (2007), ‘Modeling, analysis and testing of autonomous operation of an inverter-based microgrid’, IEEE Transactions on Power Electronics 22(2), 613–625.
  * Rana, Li and Su (2017) Rana, M. M., Li, L. and Su, S. W. (2017), ‘Cyber attack protection and control of microgrids’, IEEE/CAA Journal of Automatica Sinica 5(2), 602–609.
  * Rana, Xiang and Wang (2018) Rana, M. M., Xiang, W. and Wang, E. (2018), ‘Smart grid state estimation and stabilisation’, International Journal of Electrical Power & Energy Systems 102, 152–159.
  * Russo and Van Roy (2014) Russo, D. and Van Roy, B. (2014), ‘Learning to optimize via posterior sampling’, Mathematics of Operations Research 39(4), 1221–1243.
  * Sahraei-Ardakani and Hedman (2016) Sahraei-Ardakani, M. and Hedman, K. W. (2016), ‘Computationally efficient adjustment of facts set points in dc optimal power flow with shift factor structure’, IEEE Transactions on Power Systems 32(3), 1733–1740.
  * Slivkins et al. (2019) Slivkins, A. et al. (2019), ‘Introduction to multi-armed bandits’, Foundations and Trends® in Machine Learning 12(1-2), 1–286.
  * Smith and Paté-Cornell (2018) Smith, M. D. and Paté-Cornell, M. E. (2018), ‘Cyber risk analysis for a smart grid: how smart is smart enough? a multiarmed bandit approach to cyber security investment’, IEEE Transactions on Engineering Management 65(3), 434–447.
  * Sutton and Barto (2018) Sutton, R. S. and Barto, A. G. (2018), Reinforcement learning: An introduction, MIT press.
  * Todescato et al. (2020) Todescato, M., Bof, N., Cavraro, G., Carli, R. and Schenato, L. (2020), ‘Partition-based multi-agent optimization in the presence of lossy and asynchronous communication’, Automatica 111, 108648.
  * Yu et al. (2015) Yu, K., Ai, Q., Wang, S., Ni, J. and Lv, T. (2015), ‘Analysis and optimization of droop controller for microgrid system based on small-signal dynamic model’, IEEE Transactions on Smart Grid 7(2), 695–705.
  * Zhang and Papachristodoulou (2015) Zhang, X. and Papachristodoulou, A. (2015), ‘A real-time control framework for smart power networks: Design methodology and stability’, Automatica 58, 43–50.
  * Zhu and Martínez (2013) Zhu, M. and Martínez, S. (2013), ‘On distributed constrained formation control in operator–vehicle adversarial networks’, Automatica 49(12), 3571–3582.


# Large peak-to-valley ratio of negative-differential-conductance
in graphene p-n junctions

V. Hung Nguyen1,2111E-mail<EMAIL_ADDRESS>A. Bournel1 and P.
Dollfus1 1Institut d’Electronique Fondamentale, UMR8622, CNRS, Univ. Paris
Sud, 91405 Orsay, France
2Center for Computational Physics, Institute of Physics, VAST, P.O. Box 429 Bo
Ho, Hanoi 10000, Vietnam

###### Abstract

We investigate the transport characteristics of monolayer graphene p-n
junctions by means of the non-equilibrium Green’s function technique. It is
shown that thanks to the high interband tunneling of chiral fermions and to a
finite bandgap opening when the inversion symmetry of graphene plane is
broken, a strong negative-differential-conductance behavior with peak-to-
valley ratio as large as a few tens can be achieved even at room temperature.
The dependence of this behavior on the device parameters such as the Fermi
energy, the barrier height, and the transition length is then discussed.

##  I Introduction

Since the discovery of the Esaki tunnel diode esak58 , the effect of negative-
differential-conductance (NDC) has given rise to an intense research activity
from fundamental aspects of transport to possible applications including
oscillator, frequency multiplier, memory, fast switching, etc. mizu95 .
Motivated by the recent development of graphene nanoelectronics neto09 ;
schw10 ; lin010 , this effect has been also investigated and discussed in some
nanostructures based on monolayer graphene vndo08 ; chau09 , bilayer graphene
hung09 , and graphene nanoribbons ren009 ; lian10 ; vndo10 . However, these
studies have shown that due to their zero-bandgap, the peak-to-valley ratio
(PVR) of NDC in 2D-graphene structures is relatively small, while, though it
can be improved significantly in narrow graphene nanoribbons, the results
obtained are severely affected by increasing the ribbon’s width and/or,
especially, by the edge defects lian10 ; vndo10 .

A bandgap opening may be a key-step to improve the operation of graphene
devices. Beyond the technique which consists in patterning a graphene sheet
into nanoribbons mhan07 , some recent works have suggested that a bandgap can
open when the inversion symmetry of the graphene plane is broken, e. g., by
the interaction with the substrate zhou07 ; vita08 ; ende10 , the patterned
hydrogen adsorption balo10 , the adsorption of water molecules yava10 , or by
the controlled structural modification echt08 . In particular, the experiment
reported in zhou07  has shown that graphene epitaxially grown on SiC can
exhibit a bandgap up to $260$ meV. Though relatively small compared to that in
conventional semiconductors, it is about ten times greater than the thermal
energy at room temperature, which has stimulated further investigations of
graphene-based transistors lian07 ; kedz08 ; mich10 . In this work, our aim is
to look for possibilities of achieving a strong NDC behavior in monolayer
graphene p-n junctions. This expectation comes from the fact that the
appearance of a finite bandgap can result in a low valley current, while, the
high interband tunneling of chiral fermions in graphene may lead to a high
current peak. Therefore, a large PVR is expected to be achieved in these
junctions.

##  II Model and calculation

Figure 1: (Color online) Schematic of p-n junction based on monolayer graphene
(a) and its potential profile along the transport direction (b). The n-doped
(left) and the p-doped (right) regions are used to form a potential barrier in
the device. The transition between two regions is characterized by the length
$L$.

Graphene has a honeycomb lattice structure with a unit cell consisting of two
carbon atoms - normally referred to as $A$ and $B$ atoms. To describe the
charge states in the system, a simple nearest neighbor tight binding model can
be conveniently used neto09 , with $a_{c}=0.142$ nm as carbon-carbon distance,
$t=2.7$ eV as hopping energy between nearest neighbor sites, and
$\varepsilon_{A}=-\varepsilon_{B}=\Delta$ as the onsite energies in the two
sublattices. When $\Delta=0$, this model results in a band structure with zero
gap, i.e., the conduction and valence bands meet at the K and K’ points in the
Brillouin zone. By making $\Delta\neq 0$, the inversion symmetry is broken and
therefore a finite gap opens in graphene’s band-structure. The energy
dispersion close to the K-point simply writes

$E\left({\vec{k}}\right)=\pm\sqrt{\hbar^{2}v_{F}^{2}\left({k_{x}^{2}+k_{y}^{2}}\right)+\Delta^{2}}$
(1)

where $v_{F}=3a_{c}t/2\hbar\approx 10^{6}$ m/s is the Fermi velocity,
$\vec{k}=\left({k_{x},k_{y}}\right)$ is the 2D-momentum, and the sign +/-
stands for the conduction/valence band, respectively. From eq. (1), the
bandgap is determined as $E_{g}=2\Delta$. To describe the excited states
around such K-point, one can use the following effective (massive Dirac-like)
Hamiltonian:

$H=-i\hbar
v_{F}\left({\sigma_{x}\partial_{x}+\sigma_{y}\partial_{y}}\right)+\Delta\sigma_{z}+U$
(2)

where $U$ is the external potential energy and $\sigma_{x,y,z}$ are the Pauli
matrices. The Hamiltonian (2) is now used to study the transport
characteristics of the p-n junction schematized in Fig. 1. The p-doped and
n-doped graphene regions can be generated by electrostatic doping huar07 ;
zhan08  or by chemical doping farm09 ; bren10 . The back-gate is used to
control the carrier density in the graphene layer by applying a voltage
$V_{sub}$. The junction is characterized by the potential barrier $U_{0}$ and
the transition length $L$. This length is the size of the region across which
the charge density changes monotonically from n-type to p-type. Though
expected to be short farm09 , it is finite and is considered in this work as a
parameter.

In principle, the presence of defects and impurities is unavoidable and may
result in a substantial amount of disorder in the graphene. However, such
disorder affects strongly the electronic characteristics only when the
graphene sheets are narrow cres08 . Thus, we assume that it can be negligible
in this work where we consider the transport in wide systems vozm05 . Our
study addresses the ballistic transport through the junction and the non-
equilibrium Green’s function (NEGF) technique (see the formulation in ref.
vndo08 ) is then used to solve eq. (2) and to determine the transport
quantities.

Figure 2: (Color online) Local density of states (a) and transmission
coefficient (b) in a graphene p-n junction. The transition length in (a)
$L=20$ nm; in (b): $L=5$ nm (dotted-dashed), 10 nm (solid) and 20 nm (dashed
lines). All results are for $\tilde{\Delta}=65$ meV
($\tilde{\Delta}=\sqrt{E_{y}^{2}+\Delta^{2}}$ and $E_{y}=\hbar v_{F}k_{y}$).

##  III Results and discussion

Figure 3: (Color online) (a) $I-V$ characteristics of graphene p-n junction
with different energy bandgaps. Other parameters are: $E_{F}=0.26$ eV,
$U_{0}=0.52$ eV, and $L=10$ nm. (b) and (c) illustrate the band diagrams at
low bias and in the current valley, respectively.

Before investigating the behavior of the electrical current through these
junctions, we plot a map of the local density of states and the transmission
coefficient as a function of energy in Fig. 2(a) and Fig. 2(b), respectively.
The results show clearly three important transport regions: (i) thermionic
when $E<U_{N}-\tilde{\Delta}$ or $E>U_{P}+\tilde{\Delta}$, (ii) interband
tunneling when $U_{N}+\tilde{\Delta}<E<U_{P}-\tilde{\Delta}$, and (iii)
transmission gap when $\left|{E-U_{\alpha}}\right|\leq\tilde{\Delta}$, where
$U_{\alpha}\equiv U_{N,P}$ denote the potential energies in the n-doped and
p-doped regions, respectively. The appearance of transmission gap is
essentially due to the fact that the longitudinal momentum $k_{x}$ defined
from eq. (1) is imaginary in such energy region, i.e., the carrier states are
evanescent and therefore decay rapidly when going to one of the two junction
sides. Note that the same results are obtained for
$\tilde{\Delta}={\rm{const}}$ even with different $E_{y}$ and/or different
$\Delta$. When $\Delta=0$, the transmission gap is just a function of $E_{y}$
and disappears for normal incident particles. By rising $\Delta$, this gap
increases and gets its minimum value $E_{g}=2\Delta$ when $E_{y}=0$. Moreover,
due to the appearance of evanescent states around the neutral points in the
transition region, Fig. 2(b) shows clearly that the larger the transition
length, the lower the interband tunneling. However, it is worth noting that
the interband tunneling of chiral fermions observed here is still very high in
comparison with typical values of the order of $10^{-7}$ observed in
conventional Esaki tunnel diodes carl97 . Therefore, a high current peak and
then a large PVR are expected to be observed in the considered graphene
junctions.

Now, in Fig. 3(a), we display the $I-V$ characteristics obtained for different
energy bandgaps. Note that, throughout the work, the current is computed at
room temperature. It is shown that the NDC behavior appears more clearly and
its PVR increases when increasing the bandgap though the current peak reduces.
The form of the $I-V$ curve can be explained well by looking at the diagrams
in Fig. 3(b,c): at low bias (3(b)), the contribution of interband tunneling
processes makes the current rising; when further rising the bias (3(c)), the
interband tunneling is suppressed due to the transmission gaps and therefore
the current decreases; when the bias is high enough, the contribution of
charge carriers in the thermionic region makes the current rising very
rapidly. Because the transmission gap increases, the current, especially in
the valley region, decreases when increasing the bandgap. This results in a
large PVR of NDC, for instance, about 46 in the case of $E_{g}=260$ meV shown
in Fig. 3(a).

Figure 4: (Color online) $I-V$ characteristics of graphene p-n junction with
different Fermi energies (a) and different potential barriers (b). Other
parameters are: $E_{g}=260$ meV, $E_{F}=U_{N}+U_{0}/2$ in (b), $U_{0}=0.52$ eV
in (a), and $L=10$ nm, where $U_{N}$ is the potential energy in the n-doped
graphene. Figure 5: (Color online) $I-V$ characteristics of graphene p-n
junction with different transition lengths. The inset shows the evolution of
peak-to-valley ratio with respect to $L$. Other structure parameters are:
$E_{g}=260$ meV, $E_{F}=0.26$ eV, and $U_{0}=0.52$ eV.

To evaluate the possibilities of having a strong NDC behavior in the
junctions, we consider the behavior of the $I-V$ characteristics with
different Fermi energies in Fig. 4(a) and different barrier heights in Fig.
4(b). Fig. 4(a) shows that the current valley shifts to the low bias region
while the current peak decreases when $E_{F}$ goes from the value of
$U_{m}=U_{N}+U_{0}/2$ to the top of the valence band in the p-doped region or
to the bottom of conduction band in the n-doped one. Simultaneously, the
current is a symmetric function of $E_{F}$ around $U_{m}$. This result leads
to the fact that the strongest NDC is achieved when $E_{F}=U_{m}$ and the
rectification behavior around the zero bias can be observed when $E_{F}$ goes
away from such value, e.g., see the cases of $E_{F}=0.13$ or 0.39 eV in Fig.
4(a). The latter behavior is essentially due to the role of the transmission
gap (as illustrated in the diagrams in Figs. 3(b,c)), which induces a strong
reduction of current in the low-positive bias region. Besides, when decreasing
the barrier height, we find that the transmission gap around $U_{P}$ moves
downward in energy, and as a consequence, the interband tunneling and the
current peak are reduced. Indeed, this point is clearly illustrated in Fig.
4(b), where $E_{F}$ is chosen to be $U_{N}+U_{0}/2$ to achieve the strongest
NDC. Moreover, it is shown that the current valley shifts to the low bias and
its value is reduced, which finally results in an increase of the PVR when
decreasing $U_{0}$. For instance, the PVR is $\sim 30$, $46$, and $79$ for
$U_{0}=0.6$, $0.52$, and $0.44$ eV, respectively. To obtain the best device
operation (either a high current peak or a large PVR), the study suggests that
the barrier height about two times greater than the energy bandgap should be
used.

Finally, the $I-V$ characteristics in the junctions with different transition
lengths are displayed in Fig. 5. It is shown that because of the suppression
of interband tunneling (as in Fig. 2(b)), the current peak is reduced when
increasing $L$. Therefore, as seen in the inset of Fig. 5, the PVR decays as a
function of the length $L$. However, the PVR as large as about 10 is still
obtained even for $L$ up to 20 nm. This suggests the possibility of achieving
a very large PVR with a short transition length, which can be realized by
controlling the device geometry, e.g., by appropriately reducing the gate
dielectric thickness in the case of the electrostatic doping zhan08 , or by
using the chemical doping to generate the junction as mentioned in ref. farm09
. As expected above the large PVR, e.g., $\sim$ 123 obtained for $L=5$ nm in
this work, is much higher than for conventional Esaki tunnel diodes where its
highest reported value is just about 16 (see in ref. oehm10  and references
therein).

##  IV Conclusion

Using the NEGF technique, we have investigated the transport characteristics
of monolayer graphene p-n junctions. Even at room temperature a negative-
differential-conductance with peak-to-valley ratio as large as a few tens can
be observed in these junctions thanks to a finite bandgap and to the high
interband tunneling of chiral fermions. The dependence of this behavior on the
device parameters was analyzed. It is shown that the strong negative-
differential-conductance behavior is achieved when the Fermi level is in the
center of the potential barrier, the barrier height is about two times greater
than the bandgap, and the transition length is short. We hope that obtained
results can be helpful for designing efficient room-temperature negative-
differential-conductance devices based on finite-bandgap graphene.

## Acknowledgments

This work was partially supported by the French ANR through the projects
NANOSIM-GRAPHENE (ANR-09-NANO-016) and MIGRAQUEL (ANR-10-BLAN-0304). One of
the authors (VHN) acknowledges the Vietnam’s National Foundation for Science
and Technology Development (NAFOSTED) for financial support under the Project
No. 103.02.76.09.

## References

  * (1) L. Esaki, Phys. Rev. 109, 603 (1958).
  * (2) H. Mizuta and T. Tanoue, _The Physics and Applications of Resonant Tunneling Diodes_ (Cambridge University Press, Cambridge, 1995).
  * (3) A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
  * (4) F. Schwierz, Nat. Nanotech. 5, 487 (2010).
  * (5) Y.-M. Lin, C. Dimitrakopoulos, K. A.Jenkins, D. B. Farmer, H.-Y. Chiu, A. Grill, and Ph. Avouris, Science 327, 662 (2010); L. Liao, Y.-C. Lin, M. Bao, R. Cheng, J. Bai, Y. Liu, Y. Qu, K. L. Wang, Y. Huang, and X. Duan, Nature 467, 305 (2010).
  * (6) V. Nam Do, V. Hung Nguyen, P. Dollfus, and A. Bournel, J. Appl. Phys. 104, 063708 (2008).
  * (7) H. Chau Nguyen and V. Lien Nguyen, J. Phys.: Condens. Matter 21, 045305 (2009).
  * (8) V. Hung Nguyen, A. Bournel, V. Lien Nguyen, and P. Dollfus, Appl. Phys. Lett. 95, 232115 (2009).
  * (9) H. Ren, Q.-X. Li, Y. Luo, and J. Yang, Appl. Phys. Lett. 94, 173110 (2009).
  * (10) G. Liang, S. B. Khalid, and K.-T. Lam, J. Phys. D: Appl. Phys. 43, 215101 (2010).
  * (11) V. Nam Do and P. Dollfus, J. Appl. Phys. 107, 063705 (2010).
  * (12) M. Y. Han, B. Ozyilmaz, Y. Zhang, and P. Kim, Phys. Rev. Lett. 98, 206805 (2007).
  * (13) S. Y. Zhou, G.-H. Gweon, A. V. Fedorov, P. N. First, W. A. de Heer, D.-H. Lee, F. Guinea, A. H. Castro Neto, and A. Lanzara, Nat. Mater. 6, 770 (2007).
  * (14) L. Vitali, C. Riedl, R. Ohmann, I. Brihuega, U. Starke, and K. Kern, Surf. Sci. 602, L127 (2008).
  * (15) C. Enderlein, Y. S. Kim, A. Bostwick, E. Rotenberg, and K. Horn, New J. Phys. 12, 033014 (2010)
  * (16) R. Balog, B. Jørgensen, L. Nilsson, M. Andersen, E. Rienks, M. Bianchi, M. Fanetti, E. Lægsgaard, A. Baraldi, S. Lizzit, Z. Sljivancanin, F. Besenbacher, B. Hammer, T. G. Pedersen, P. Hofmann, and L. Hornekær, Nat. Mater. 9, 315 (2010).
  * (17) F. Yavari, C. Kritzinger, C. Gaire, L. Song, H. Gullapalli, T. Borca-Tasciuc, P. M. Ajayan, and N. Koratkar, Small 6, 2535 (2010).
  * (18) T. J. Echtermeyer, M. C. Lemme, M. Baus, B. N. Szafranek, A. K. Geim, and H. Kurz, IEEE Electron Device Lett. 29, 952 (2008).
  * (19) G. Liang, N. Neophytou, M. S. Lundstrom, and D. E. Nikonov, J. Appl. Phys. 102, 054307 (2007).
  * (20) J. Kedzierski, P.-L. Hsu, P. Healey, P. W. Wyatt, C. L. Keast, M. Sprinkle, C. Berger, and W. A. de Heer, IEEE Trans. Electron Devices 55, 2078 (2008).
  * (21) P. Michetti, M. Cheli, and G. Iannaccone, Appl. Phys. Lett. 96, 133508 (2010).
  * (22) B. Huard, J.A. Sulpizio, N. Stander, K. Todd, B. Yang, and D. Goldhaber-Gordon, Phys. Rev. Lett. 98, 236803 (2007).
  * (23) L. M. Zhang and M. M. Fogler, Phys. Rev. Lett. 100, 116804 (2008); T. Low, S. Hong, J. Appenzeller, S. Datta, and M. S. Lundstrom, IEEE Trans. Electron Devices 56, 1292 (2009).
  * (24) D. B. Farmer, Y.-M. Lin, A. Afzali-Ardakani, and P. Avouris, Appl. Phys. Lett. 94, 213106 (2009).
  * (25) K. Brenner and R. Murali, Appl. Phys. Lett. 96, 063104 (2010).
  * (26) A. Cresti, N. Nemec, B. Biel, G. Niebler, F. Triozon, G. Cuniberti, and S. Roche, Nano Res. 1, 361 (2008).
  * (27) In wide graphene sheets, the most remarkable change in the electronic structure induced by the disorder appears at the K-point where it can give rise to localized states (see in M. A. H. Vozmediano _et al._ , Phys. Rev. B 72, 155121 (2005)).
  * (28) A. D. Carlo, P. Lugli, and P. Vogl, Solid State Commun. 101, 921 (1997).
  * (29) M. Oehme, M. Šarlija, D. Hähnel, M. Kaschel, J. Werner, E. Kasper, and J. Schulze, IEEE Trans. Electron Devices 57, 2857 (2010).