††thanks: These authors contributed equally††thanks: These authors contributed
equally

# Discovery of an electronic crystal in a cuprate Mott insulator

Mingu Kang  Department of Physics, Massachusetts Institute of Technology;
Cambridge, MA 02139, USA  Center for Complex Phase of Materials, Max Planck
POSTECH Korea research Initiative; Pohang 790-884, Republic of Korea  Charles
Zhang  Department of Physics, University of Toronto, Toronto, Ontario M5S1A7,
Canada  Enrico Schierle  Helmholtz-Zentrum Berlin für Materialien und Energie;
Albert-Einstein Strabe 15, D-12489 Berlin, Germany  Stephen McCoy  Material
Science and Engineering, University of California, Riverside; Riverside,
California 92521, USA  Department of Electrical Engineering, University of
California, Riverside, Riverside; California 92521, USA  Jiarui Li  Department
of Physics, Massachusetts Institute of Technology; Cambridge, MA 02139, USA
Ronny Sutarto  Canadian Light Source; Saskatoon, Saskatchewan S7N 2V3, Canada
Feizhou He  Canadian Light Source; Saskatoon, Saskatchewan S7N 2V3, Canada
Andreas Suter  Laboratory for Muon-Spin Spectroscopy, Paul Scherrer Institute,
CH-5232 Villigen, Switzerland  Thomas Prokscha  Laboratory for Muon-Spin
Spectroscopy, Paul Scherrer Institute, CH-5232 Villigen, Switzerland  Zaher
Salman  Laboratory for Muon-Spin Spectroscopy, Paul Scherrer Institute,
CH-5232 Villigen, Switzerland  Eugen Weschke  Helmholtz-Zentrum Berlin für
Materialien und Energie; Albert-Einstein Strabe 15, D-12489 Berlin, Germany
Shane Cybart  Material Science and Engineering, University of California,
Riverside; Riverside, California 92521, USA  Department of Electrical
Engineering, University of California, Riverside, Riverside; California 92521,
USA  John Y. T. Wei  Department of Physics, University of Toronto, Toronto,
Ontario M5S1A7, Canada  Riccardo Comin<EMAIL_ADDRESS>Department of Physics,
Massachusetts Institute of Technology; Cambridge, MA 02139, USA

###### Abstract

Copper oxide high temperature superconductors universally exhibit multiple
forms of electronically ordered phases that break the native translational
symmetry of the CuO2 planes. The interplay between these orders and the
superconducting ground state, as well as how they arise through doping a Mott
insulator, is essential to decode the mechanisms of high-temperature
superconductivity. Over the years, various forms of electronic liquid crystal
phases – including charge/spin stripes and incommensurate charge-density-waves
(CDWs) – were found to emerge out of a correlated metallic ground state in
underdoped cuprates. Early theoretical studies also predicted the emergence of
a Coulomb-frustrated ’charge crystal’ phase in the very lightly-doped,
insulating limit of the CuO2 planes. This charge crystal phase was observed
only recently at the surface of Bi2Sr2CaCu2O8+y, raising new questions on
whether such a state may also exist in the bulk and thus be a common
instability of cuprates. Here, we use resonant X-ray scattering, electron
transport, and muon spin rotation measurements to fully resolve the electronic
and magnetic ground state and search for signatures of charge order in very
lightly hole-doped cuprates from the 123 family RBa2Cu3O7-δ (RBCO; R: Y or
rare earth). X-ray scattering data from RBCO films reveal a breaking of
translational symmetry more pervasive than was previously known, extending
down to the insulating and magnetically ordered Mott limit. The ordering
vector of this charge crystal state is linearly connected to the charge-
density-waves of underdoped RBCO, suggesting that the former phase is a
precursor to the latter as hole doping is increased. Most importantly, the
coexistence of charge and spin order in RBCO suggests that this electronic
symmetry-breaking state is common to the CuO2 planes in the very lightly-doped
regime. These findings bridge the gap between the Mott insulating state and
the underdoped metallic state, and underscore the prominent role of Coulomb-
frustrated electronic phase separation among all cuprates.

##  I  I. Introduction

Following the discovery of high-temperature superconductivity [1], early
theoretical studies proposed that multiple forms of electronic symmetry-
breaking exist in lightly-doped Mott insulators [2, 3, 4, 5]. These
theoretical predictions were confirmed by the discovery of charge and spin
stripes in La1.48Nd0.4Sr0.12CuO4 [6]. During the last decade,multiple forms of
these symmetry-breaking states have been discovered in all cuprate high-
temperature superconductors [7, 8, 9, 10, 11, 12, 13, 14, 15, 16]. In many
cases, these emergent phases are intertwined with superconductivity [17, 18,
19] in ways that have been investigated from multiple perspectives [20, 21,
22, 23, 24, 25, 26, 27].

The driving mechanisms behind electronic order in the CuO2 planes have been
the subject of longstanding questions that are essential for understanding of
quantum matter and high-temperature superconductivity in cuprates [28].
However, a consensus on the origin of electronic order has not formed fully,
partly because of key phenomenological differences between the cuprate
families of YBa2Cu3Oy (123) Bi2Sr2CaCu2O8+y (2212), and La2CuO4 (214). In the
latter, charge order cohabits with static spin order to form electronic
stripes characterized by Coulomb-frustrated real-space phase separation and a
spatial periodicity directly correlated with the hole density. However, this
’real-space Coulomb frustration’ scenario has not found wide applicability in
other cuprate families, because of the competing hypothesis of a ’momentum-
space instability’ involving scattering between itinerant states mediated by
strong antiferromagnetic fluctuations.

Figure 1:  Characterization of Pr-substituted YBCO and its phase diagram. (a)
Crystal structure of Y1-xPrxBCO containing a CuO2-Y,Pr-CuO2 bilayer. (b,c)
Phase diagram of YBCO7-δ (b) and Y1-xPrxBCO (c). Doping dependent data points
for superconductivity in YBCO7-δ, superconductivity in Y1-xPrxBCO, and charge
order in YBCO7-δ are adapted from Refs. [9, 29, 30]. (d) X-ray absorption
spectra of Y1-xPrxBCO films. Two main peaks located at 931.3 eV and 930 eV
corresponds to the Cu-$L_{3}$ and Pr-$M_{5}$ absorption edges respectively.
(e) Resistance of Y1-xPrxBCO films as a function of temperature (note the use
of two different linear scales on the $y$-axis above and below 250 $\Omega$).
(f) Asymmetry spectrum of muon spin rotation data from PrBCO films at 75 K.
(g) Temperature-dependent internal field and depolarization rate. Each data
point is obtained by fitting the asymmetry spectra taken at different
temperatures. Grey-shaded bars indicate the onset of magnetic order at
$T_{N}\approx$ 250 K. A magnetic volume fraction close to 100% has been
deduced from weak transverse field experiments (${B}_{ext}=75G$) measurements
(not shown), indicating that the film is fully magnetic.

A clear discriminant between these two frameworks is the existence of an
electronic crystal phase in the very lightly hole-doped and insulating CuO2,
before the onset of superconductivity. This state has been predicted under the
real-space scenario [31] but it is fundamentally incompatible with the
momentum-space scenario which requires a correlated metallic state with low-
energy quasiparticles and a Fermi surface. A recent microscopy study of Bi-
based cuprates has found periodic spatial modulations in the tunneling
conductance of an insulating Bi-2212 surface, providing fresh evidence for an
electronic crystal state [32]. This study demonstrates that the phenomenology
of electronic orders in underdoped Bi-based cuprates needs to be understood on
the basis of the real-space Coulomb frustration picture. If confirmed in the
family of Y-based cuprates, these findings would have a significant impact on
our description of electronic orders in the CuO2 planes across all cuprates.

To address this question, we turn our attention to the 123 family of
$R$Ba2Cu3O7-δ compounds (RBCO; R: Y or rare earth). Among cuprate high-
temperature superconductors, the RBCO family is one of the most widely
studied, owing to its low degree of electronic disorder that is enabled by
sourcing the carrier doping far away from the CuO2 planes. This possibility is
uniquely enabled by the crystal structure of RBCO, with the CuO chain layer
providing a charge reservoir for the CuO2 planes, which form a bilayer around
the spacer rare earth layer [Fig. 1(a)]. To sidestep challenges in the control
of oxygen stoichiometry in the chain layer at very low hole doping levels, we
have taken an alternative route of using Pr-substitution at the Y site in
nearly fully-oxygenated ($\delta\ll 1$) samples of (Y,Pr)Ba2Cu3O7
(Y1-xPrxBCO7-δ), taking cues from previous studies of doping control and
suppression of superconductivity in these systems [see Fig. 1(c) for the phase
diagram of Y1-xPrxBCO] [29, 33, 34]. In Y1-xPrxBCO, the concentration of
mobile carriers in the CuO2 planes is controlled via rare earth substitution
(Pr-for-Y) using the known ability of low-energy Pr $4f$ orbitals to
effectively sequester holes from Cu-${3d}_{x^{2}-y^{2}}$-O${}_{{2p}_{\sigma}}$
states [35, 36]. The fully-substituted, end compound PrBCO is known to be
insulating and magnetic [Fig. 1(c)], therefore providing an ideal platform to
study electronic orders down to the Mott limit.

Using resonant soft X-ray scattering (RXS) at Cu-$L_{3}$ (2$p\rightarrow$
3$d$) and Pr-$M_{5}$ (3$d\rightarrow$ 4$f$) edges, we reveal an electronically
modulated state across the extended phase diagram of Y1-xPrxBCO, from the
underdoped superconducting regime ($x$ = 0.275 and 0.5) to the Mott limit ($x$
= 1). In PrBCO, the observation of a Cu-site density wave in the absence of a
Fermi surface and its coexistence with static Neel order crucially underscores
the need for Coulomb frustration to support charge order in RBCO.

##  II II. Charge order
in magnetic insulating PrBa2Cu3O7-δ

We first present the basic characterization of Y1-xPrxBCO films with Pr
content $x$ = 0.275, 0.5, and 1. X-ray absorption spectra are shown in Figure
1(d). Two prominent absorption peaks are observed at 931.3 eV and 930.0 eV,
corresponding to the Cu-$L_{3}$ and Pr-$M_{5}$ resonances, respectively. The
intensity of the Pr-$M_{5}$ peak grows in accordance with higher Pr content
$x$. Transport properties of the Y1-xPrxBCO films are also displayed in Fig.
1(e). The $T_{c}$ of pristine YBCO at this oxygen content is $\approx$ 90 K in
accordance with prior reports [9, 37, 30]. With Pr substitution, $T_{c}$ is
reduced to 55 K at $x$ = 0.275 [Fig. 1(e), $T_{c}$ is estimated from the
midpoint of the transition], and is further suppressed to 12 K for $x$ = 0.5.
In the latter sample, an upturn in resistivity is manifested at low
temperatures, signaling the proximity to a metal-insulator transition.
Finally, $x$ = 1 PrBCO behaves as a strong insulator with its resistivity
diverging at low temperature. We also performed zero field muon spin rotation
(ZF-$\mu$SR) measurements to characterize the magnetic properties of the end
compound PrBCO. As shown in Fig. 1(f), the ZF-$\mu$SR time spectrum at 75 K
exhibits time oscillations that are indicative of the presence of static
magnetic order in the PrBCO film. By fitting the ZF-$\mu$SR spectra using a
standard function (see Appendix for detail), we extract the internal field and
depolarization rate as a function of temperature in Fig. 1(g). The internal
field shows a sharp onset of static spin order at $T_{N}\approx$ 250 K, where
the decay rate also peaks. The latter shows a second rise below 25 K due to
the additional ordering of Pr magnetic moments. Overall, the electronic and
magnetic properties of the Y1-xPrxBCO thin films studied here are
quantitatively consistent with previous reports [34, 29, 36, 38], confirming
that Pr-substitution effectively underdopes the CuO2 planes all the way to the
magnetically ordered, insulating limit [Fig. 1(c)].

Figure 2:  Charge order in PrBCO. Wide momentum-range RXS scans on PrBCO along
the Cu-O bond direction ($Q_{x}$) obtained at the Cu-$L_{3}$ edge. The double
peak structure is marked by red solid lines. The grey-dashed lines are used to
model the fluorescence background, while the red- and blue-shaded areas are
fitted Gaussian peaks. Figure 3:  Resonant profile of charge order and
hybridization between Pr-4$f$ and CuO2 orbitals. (a) Photon energy-dependence
of the X-ray scattering data (momentum scans) for positive $Q_{x}$. The curves
are vertically shifted for clarity. Two curves with data points marked with
open circles highlight the RXS scans at the peak of the Pr- and Cu-resonances.
(b) Intensity map representing the resonant behavior of charge order after
subtraction of the fluorescence background. The red-dashed and blue-solid
curves on the right show the resonant profile at
$Q_{\mathrm{CDW}}^{\mathrm{Cu}}$ and $Q_{\mathrm{CDW}}^{\mathrm{Pr}}$
respectively. (c,d) Schematics of the two density-waves with different Cu-3$d$
and Pr-4$f$ orbital contents. (e) Schematic of the Zhang-Rice singlet state
formed upon hybridization of planar Cu-3$d_{x2-y2}$ orbital and O-2$p$ hole
states. (f) Schematic of the Fehrenbacher-Rice state formed upon hybridization
of Pr-4$f_{z(x2-y2)}$ orbital and surrounding O-2$p$ hole states.

Figure 2 shows a RXS scan taken in PrBCO over a broad range of in-plane
momentum along the Cu-O bond direction ($Q_{x}$). On top of a smooth
background due to fluorescence and the tail of a specular reflection near
$Q_{x}$ = 0, a two-peak structure appears symmetrically at positive and
negative momentum transfers between $|Q_{x}|\approx$ 0.2 and 0.4 r.l.u..
Gaussian fits of the RXS profile reveal that the diffraction signatures are
composed of two distinct peaks centered at 0.36 and 0.28 r.l.u.. The full-
width-at-half-maximum in both cases is estimated to be $\approx$ 0.06-0.08
r.l.u., reflecting the presence of short-range charge correlations with a
spatial coherence length $\xi\approx$ 4-5 unit cells or $\approx$ 15-20 Å.
Both diffraction peaks originate from ordering of the charge degrees of
freedom as confirmed by the photon polarization analysis (see Fig. S1 in
Supplemental Material) [7, 39]. Additional azimuthal angle-dependent
experiments indicate that the double-peak structure appears only along the
Cu-O bond direction (see Fig. S2 in Supplemental Material) [39]. The
correlation length, polarization dependence, and ordering axis of these
diffraction signatures are all consistent with the characteristics of charge
order in other cuprates. However, the coexistence of Cu-$3d$ charge order with
static Cu-$3d$ spin order departs from all previous studies, signaling the
emergence of a new state.

To understand these initial observations, we performed RXS experiments as a
function of photon energy across the Cu-$L_{3}$ (931.3 eV) and Pr-$M_{5}$
(930.0 eV) edges. As shown in Fig. 3(a), the two main diffraction peaks
display a surprising evolution with photon energy. When X-rays are tuned more
than 2 eV away from the Pr and Cu resonances, the peaks vanish (see the
profiles at 928.0 and 933.0 eV for examples), indicating that the diffraction
signatures arise from charge modulations within the CuO2-Pr-CuO2 bilayer. Most
importantly, and as can be seen in the intensity color map of Fig. 3(b), the
first diffraction peak, hereafter termed
$Q_{\mathrm{CDW}}^{\mathrm{Pr}}\approx$ 0.28 r.l.u., resonates at the Pr-edge,
while the second peak $Q_{\mathrm{CDW}}^{\mathrm{Cu}}\approx$ 0.36 r.l.u. is
most prominent at the Cu-edge (931.3 eV).

Figure 4:  Charge order as a function of doping. (a) Doping-dependent RXS
spectra from Y1-xPrxBCO films ($x$ = 0.275, 0.5, 1). RXS spectra from nearly
optimally doped films of HoBCO and EuBCO (our ’control’ samples) are also
displayed. (b) Evolution of the charge order wavevectors in the RBCO family.
Open and filled blue diamonds represent $Q_{CDW}$ of YBCO7-δ along the $a$\-
and $b$-crystallographic axes, respectively (adopted from Ref. [9]). The
relationship between $Q_{CDW}$ and $p$ is obtained from the fit to the YBCO7-δ
series (grey-dashed line), from which we evaluate the hole-doping levels of
Y1-xPrxBCO.

##  III  III. Doping dependence of the diffraction data

X-ray scattering data indicate the presence of two electronic modulations. The
$Q_{\mathrm{CDW}}^{\mathrm{Pr}}$ peak can be ascribed to a spatial modulation
of the Pr-4$f$ orbital occupation [see schematics in Fig. 3(d)]. We interpret
this signal as arising from Pr-4$f$-O${}_{{2p}_{\pi}}$ hybridized valence
states [Fig. 3(f)] which have been proposed within the Fehrenbacher-Rice model
[40].

On the other hand, the $Q_{\mathrm{CDW}}^{\mathrm{Cu}}$ peak corresponds to a
periodic modulation of the occupation of Cu-${3d}_{x^{2}-y^{2}}$ orbitals
[Fig. 3(c,e)] and thus is closely related to the CDW order observed in YBCO.
This connection is validated by comparing the CDW signatures of PrBCO with
those we measured in ’control’ samples from the RBCO family. Namely, in
superconducting EuBCO and HoBCO films with near-optimal oxygen content, the
RXS profiles (See Supplemental Materials Fig. S4) resonate at the same energy
but have a very different ordering vector
$Q_{\mathrm{CDW}}^{\mathrm{Cu}}\approx$ 0.3. For samples at intermediate
doping within the (Y,Pr)BCO series, we find
$Q_{\mathrm{CDW}}^{\mathrm{Cu}}\approx$ 0.32 r.l.u. and 0.33 r.l.u. for $x$ =
0.275 and $x$ = 0.5, respectively. This progressive increase of the ordering
vector as doping is reduced can be understood on the basis of the well-
established dependence of the ordering vector on hole doping, which has been
previously used to ’gauge’ the doping level in oxide heterostructures [37].
Here, we use this approach to infer the carrier filling of the CuO2 planes of
PrBCO from the relationship between $Q_{\mathrm{CDW}}^{\mathrm{Cu}}$ and hole
doping ($p$) in the YBCO7-δ series [Fig. 4(b); data from Ref. [9]]. Using
$Q_{\mathrm{CDW}}^{\mathrm{Cu}}$ = 0.32 and 0.33 r.l.u. for $x$ = 0.275 and
0.5 and $Q_{\mathrm{CDW}}^{\mathrm{Cu}}$ = 0.36 r.l.u. for $x$ = 1 yields hole
concentrations $p$ = 0.11, 0.085, and 0.01, respectively. This latter estimate
indicates that the end compound PrBCO is in the very lightly hole-doped Mott
limit, in excellent accordance with the observed insulating and magnetic
ground state (with high Neel temperature $T_{N}$).

Figure 5:  Evolution of charge order in the Y1-xPrxBCO series. (a,b)
Temperature-dependent momentum scans from $x$ = 0.275 and $x$ = 0.5 Y1-xPrxBCO
films, measured at Cu-$L_{3}$ edge. (c,d) Temperature-dependent momentum scans
from PrBCO films, measured at Cu-$L_{3}$ (c) and Pr-$M_{5}$ edge (d). (e-g),
Resonant profile of the charge order peaks in $x$ = 0.275 (e), $x$ = 0.5 (f),
and $x$ = 1 (g). In (g), the resonant behavior of
$Q_{\mathrm{CDW}}^{\mathrm{Cu}}$ and $Q_{\mathrm{CDW}}^{\mathrm{Pr}}$ is
displayed by the red and blue lines. Overlaid filled grey curves represent the
X-ray absorption spectra of each sample. (h-j) Doping-dependence of wave
vectors (h), cutoff temperatures (i), and widths (j) of the charge order peaks
in Y1-xPrxBCO. YBCO ($x$ = 0, $\delta\approx$ 0) data are adopted from
previous RXS studies (Ref. [9, 37]). The data points plotted with grey square
and grey circle correspond to bulk YBCO (Ref. [9]) and YBCO thin film (Ref.
[37], 50 nm thick on La2/3Ca1/3MnO3), respectively. Rose shades in (h,i) are
guides-to-the-eye.

##  IV  IV. Photon energy and temperature dependence of the diffraction data

The resonant profiles (X-ray scattering intensity vs. photon energy) for
Y1-xPrxBCO samples at intermediate doping levels are shown in Fig. 5(e) and
5(f). The presence of a single resonance (at the Cu-edge) in these samples, in
combination with the linear dependence of the ordering vector on doping, lends
further support to the close connection between the charge crystal
$Q_{\mathrm{CDW}}^{\mathrm{Cu}}$ peak observed in PrBCO and the charge-
density-wave signatures observed at higher hole doping.

For the temperature dependence, we first focus on PrBCO, where the diffraction
signatures at the Cu and Pr edge display a differing temperature dependence.
As shown in Fig. 5(c,d), the resonant scattering intensity
$I_{\mathrm{CDW}}^{\mathrm{Pr}}$ clearly decreases with increasing
temperature, while $I_{\mathrm{CDW}}^{\mathrm{Cu}}$ evolves more gradually all
the way up to 300 K (see Fig. S3 for details on the quantitative analysis of
the temperature dependence [39]). This difference suggests that, while
coexisting, these two modulated states are independent as it can be expected
from their association to different degrees of freedom [20, 21, 22, 23, 24,
25, 26, 27]. We also note that the very gradual evolution of
$I_{\mathrm{CDW}}^{\mathrm{Cu}}$ in PrBCO is reminiscent of the temperature
dependence of charge density fluctuations seen in [41] and therefore could
suggest a partly fluctuating nature of the observed spatial modulations. The
temperature dependence for the intermediate doping levels shows an evolution
of $I_{\mathrm{CDW}}^{\mathrm{Cu}}$ consistent with a very gradual onset of
the CDW signal around 150-200 K, as observed in bulk YBCO in previous studies
[9].

To summarize the trends for the cutoff temperature and the correlation
lengths, while the former tends to increase with Pr-doping [Fig. 5(i)], the
latter remains similarly short-ranged for all doping levels [Fig. 5(j)].

##  V  V. Charge order in the magnetic insulating limit

In Figure 6, we summarize our results in the extended phase diagram of RBCO.
The latter includes long-range Cu (Neel) spin order (from our $\mu$SR data and
Ref. [29, 33]), low-temperature Pr spin order (from Ref. [29, 33, 42]),
superconductivity [from our transport data in Fig. 1(e) and Ref. [29, 33]],
pseudogap [from resistivity anomalies, Fig. 1(e)], and charge order (from our
RXS data, Figs. 2-5).

The present RBCO phase diagram demonstrates the coexistence between static Cu
charge and spin orders in the very lightly-doped limit. Combined with the
recent observation of electronic modulation on a very lightly-doped insulating
surface of Bi2Sr2CaCu2O8+y [32], our results support the existence of a common
’charge crystal’ phase emerging in the strongly-correlated ground state of
lightly-doped CuO 2 planes. This observation is highly significant for our
understanding of electronic symmetry-breaking in cuprates, as it reveals that
electronic liquid or crystal phases pervade the phase diagram of the CuO2
planes all the way into the Mott limit.

Figure 6:  Extended phase diagram of RBCO. Green region denotes long-range Cu
(Neel) spin order (from our $\mu$SR data and Ref. [29, 33]). Dashed red line
marks the onset of low-temperature Pr spin order (from Ref. [29, 33, 42]).
Blue region denotes superconductivity (from transport data in Fig. 1(e) and
Ref. [29, 33]). Cyan region denotes charge order (from RXS data, Figs. 2-5).
Magenta line on the left edge denotes the “charge crystal” phase. Grey squares
mark the pseudogap onset temperature (from non-linear deviations in
resistivity, see Fig. S5 in supplementary material). The top horizontal axis
(hole concentration $p$) uses the $p$ values as extrapolated from the doping-
dependence of the CDW wavevector $Q_{\mathrm{CDW}}$ (see text and Fig. 4). The
${T}_{\mathrm{CDW}}$ data point at $x=0$ and $p=0.15$ is from Ref. [9] and is
consistent with our results on HoBCO and EuBCO (see Supplementary Figure 4).

Regarding the origin of charge order in cuprates, two main viewpoints have
been adopted over the years [20, 21, 22, 23, 24, 25, 26, 27]: the real space
Coulomb-frustration scenario; and the momentum-space instability scenario.
These competing hypotheses have been invoked to successfully explain different
aspects of the phenomenology of charge order in cuprates. In this context, the
present results offer important new insights, as the observation of a
spatially modulated electronic phase in the insulating limit is fundamentally
incompatible with the momentum-space instability scenario due to the absence
of low-energy itinerant states in a gapped system. On the other hand, the
real-space Coulomb-frustration scenario has the merit that it does not require
the presence of a Fermi surface and is compatible with the insulating ground
state of PrBCO.

##  VI  Acknowledgements

This work was supported by the Air Force Office of Scientific Research Young
Investigator Program under grant FA9550-19-1-0063 and National Research
Foundation of Korea, Ministry of Science, Grant No. 2016K1A4A4A01922028. Work
at the University of Toronto was supported by the Natural Sciences and
Engineering Research Council (NSERC) and the Canada Foundation for Innovation
(CFI). Part of the research described in this paper was performed at the
Canadian Light Source, a national research facility of the University of
Saskatchewan, which is supported by CFI, NSERC, the National Research Council
(NRC), the Canadian Institutes of Health Research (CIHR), the Government of
Saskatchewan, and the University of Saskatchewan. M.K. acknowledges a Samsung
Scholarship from the Samsung Foundation of Culture

##  VII  Appendix: Materials and methods

###  VII.1  1\. Pulsed Laser Deposition of YPrBCO thin films

(Y,Pr)BCO thin films used in this study were epitaxially grown along the
$c$-axis on (LaAlO3)0.3(Sr2TaAlO6)0.7 (LSAT) substrates using pulsed laser-
ablated deposition (PLD). LSAT was chosen as substrate material for this study
because its in-plane lattice constant (3.868 Å) is well matched with the in-
plane lattice constants of (Y,Pr)BCO films [43]. The PLD growth was done using
a 248 nm KrF excimer laser with a fluence of 2 J/cm2. Films were grown using a
substrate temperature of 800 °C and 200 mTorr O2 chamber pressure. Further
details of the PLD growth process are given in Ref. [44]. All the films grown
had nominal thickness of 40 nm, which was set using target-specific growth
rates. The growth rates were determined using atomic force microscopy, by
measuring the step heights of test films after they were chemically etched to
create a step edge. The films are expected to be $\approx$ 50 % twinned,
following the results of other similarly grown YBCO films (see Ref. [45]).

###  VII.2  2\. Plasma assisted reactive sputter deposition of HoBCO and EuBCO
films

HoBCO and EuBCO films were deposited on LSAT (100), polished single crystal
substrates. The substrates were heated from 825 °C during deposition. A DC
magnetron sputtering source was used with 50 mm HoBCO and EuBCO stoichiometric
targets that had been verified by XRD. The chamber during growth was 330 mTorr
with an Argon to Oxygen ratio of 1.5 sccm. Once the chamber settings were
stable the DC magnetron sputter gun was ignited and ramped up to 55 watts over
the course of 10 minutes. Typical growth rates were around 1 nm/min with a
final film thickness of 55 nm. After deposition the chamber was filled with
500 Torr of ultra-high purity oxygen and allowed to stand at deposition
temperatures for 10 minutes. After, power was cut from the heater and the
samples were allowed to cool to room temperature in 500 Torr of oxygen. The
critical temperatures were 92-94K with current density of 3MA/cm2 at 77 K.
Electrical transport measurements for HoBaCuO Josephson devices are reported
here [46].

###  VII.3  2\. X-ray diffraction and transport characterization of YPrBCO
films

The as-grown films were characterized using X-ray diffraction (XRD) and
resistance vs. temperature ($R$ vs. $T$) measurements. XRD measurements were
done using the $\theta$-2$\theta$ method and show the expected substrate and
film peaks. The $R$ vs. $T$ measurements were done in a He-4 cryostat using an
AC resistance bridge and silver contacts in the Van der Pauw configuration.
Superconducting transitions were observed in the $x$ = 0.275 and $x$ = 0.5
films, and insulating behavior in the $x$ = 1 film. The critical temperatures
($T_{c}$) and $R$ vs. $T$ behaviors are similar to other YPrBCO films of
comparable Pr and oxygen doping (see Refs. [29, 47, 48, 49]). Similarly grown
YBCO films were shown to be epitaxially bonded to the substrate using
transmission electron microscopy (TEM) in a previous study [43].

###  VII.4  3\. Resonant X-ray scattering experiments

Resonant X-ray scattering experiments RXS experiments were performed in two
different beamlines, the UE46-PGM1 beamline of BESSY II at Helmholtz Zentrum
Berlin ($x$ = 0.275 and 1) and the REIXS (10ID-2) beamline of Canadian Light
Source ($x$ = 0.5). All RXS measurements were conducted in ultra-high vacuum
diffractometers with the pressure better than 10-9 torr. Samples were aligned
$insitu$ using (0 0 2) and (1 0 2) Bragg diffractions. Momentum-space scans
were obtained by rocking the sample angle at a fixed detector position. X-ray
absorption spectra (XAS) were measured by monitoring total electron yields. We
used out-of-scattering-plane ($\sigma$) polarizations for both RXS and XAS
measurements unless specified.

###  VII.5  4\. Muon spin spectroscopy experiments

The Low Energy $\mu$SR experiments where performed at the $\mu$E4 beamline
[50] which is part of the Swiss Muon Source, Paul Scherrer Institut, Villigen,
Switzerland. The $\mu$SR data have been analyzed with the help of MUSRFIT
[51]. Details about the $\mu$SR technique can be found in Ref. [52]. All the
measurements were conducted in ultra-high vacuum at a pressure smaller than
$10^{-8}$ mbar and zero magnetic field conditions. The decay positron
asymmetry spectra were analyzed with following function

$A(t)=A_{\rm osc}\,\cos(\gamma_{\mu}B_{\rm int}t+\phi)\cdot e^{-\lambda_{\rm
T}t}+A_{1}\,e^{-(\lambda_{\rm L}t)^{\beta}},$

where $A_{\rm osc}$ is the coherent zero field oscillation amplitude,
$\gamma_{\mu}$ is the muon gyromagnetic ratio, $B_{\rm int}$ the internal
magnetic field at the muon site, and $\phi$ is the phase between the initial
muon spin direction and the decay positron detector. $\lambda_{\rm T,L}$ are
the transverse and longitudinal depolarization rates. $\beta\simeq 1$ except
close to the phase transition where it is close to 2.

## References

  * Bednorz and MÃ¼ller [1986] J. G. Bednorz and K. A. MÃ¼ller, Possible highT c superconductivity in the Ba-La-Cu-O system, Zeitschrift fÃ¼r Physik B Condensed Matter 64, 189 (1986).
  * Zaanen and Gunnarsson [1989] J. Zaanen and O. Gunnarsson, Charged magnetic domain lines and the magnetism of high-Tc oxides, Phys. Rev. B 40, 7391 (1989).
  * Machida [1989] K. Machida, Magnetism in La2CuO4 based compounds, Physica C: Superconductivity 158, 192 (1989).
  * Poilblanc and Rice [1989] D. Poilblanc and T. M. Rice, Charged solitons in the Hartree-Fock approximation to the large- \textit{U} Hubbard model, Physical Review B 39, 9749 (1989).
  * Emery  _et al._ [1990] V. J. Emery, S. A. Kivelson, and H. Q. Lin, Phase separation in the t-J model, Physical Review Letters 64, 475 (1990).
  * Tranquada  _et al._ [1995] J. M. Tranquada, B. J. Sternlieb, J. D. Axe, Y. Nakamura, and S. Uchida, Evidence for stripe correlations of spins and holes in copper oxide superconductors, Nature 375, 561 (1995).
  * Ghiringhelli  _et al._ [2012] G. Ghiringhelli, M. L. Tacon, M. Minola, C. Mazzoli, N. B. Brookes, G. M. D. Luca, A. Frano, D. G. Hawthorn, F. He, T. Loew, M. M. Sala, D. C. Peets, M. Salluzzo, E. Schierle, R. Sutarto, G. a. Sawatzky, E. Weschke, B. Keimer, and L. Braicovich, Long-Range Incommensurate Charge Fluctuations in (Y,Nd)Ba2Cu3O6+x, Science 337, 821 (2012).
  * Chang  _et al._ [2012] J. Chang, E. Blackburn, A. T. Holmes, N. B. Christensen, J. Larsen, J. Mesot, R. Liang, D. a. Bonn, W. N. Hardy, A. Watenphul, M. V. Zimmermann, E. M. Forgan, and S. M. Hayden, Direct observation of competition between superconductivity and charge density wave order in YBa2Cu3O6.67, Nature Physics 8, 871 (2012), arXiv:1206.4333 .
  * Blanco-Canosa  _et al._ [2014] S. Blanco-Canosa, A. Frano, E. Schierle, J. Porras, T. Loew, M. Minola, M. Bluschke, E. Weschke, B. Keimer, and M. L. Tacon, Resonant X-ray Scattering Study of Charge Density Wave Correlations in YBa2Cu3O6+x, Physical Review B 90, 054513 (2014), arXiv:1406.1595 .
  * Comin  _et al._ [2014] R. Comin, A. Frano, M. M. Yee, Y. Yoshida, H. Eisaki, E. Schierle, E. Weschke, R. Sutarto, F. He, A. Soumyanarayanan, Y. He, M. Le Tacon, I. S. Elfimov, J. E. Hoffman, G. A. Sawatzky, B. Keimer, and A. Damascelli, Charge Order Driven by Fermi-Arc Instability in Bi2Sr2-xLaxCuO6+d, Science 343, 390 (2014).
  * da Silva Neto  _et al._ [2014] E. H. da Silva Neto, P. Aynajian, A. Frano, R. Comin, E. Schierle, E. Weschke, A. Gyenis, W. Jinsheng, J. Schneeloch, Z. Xu, S. Ono, G. Gu, M. Le Tacon, and A. Yazdani, Ubiquitous interplay between charge ordering and high temperature superconductivity in cuprates, Science (New York, N.Y.) 343, 393 (2014).
  * Tabis  _et al._ [2014] W. Tabis, Y. Li, M. L. Tacon, L. Braicovich, A. Kreyssig, M. Minola, G. Dellea, E. Weschke, M. J. Veit, M. Ramazanoglu, A. I. Goldman, T. Schmitt, G. Ghiringhelli, N. Barišić, M. K. Chan, C. J. Dorow, G. Yu, X. Zhao, B. Keimer, and M. Greven, Charge order and its connection with Fermi-liquid charge transport in a pristine high-Tc cuprate, Nature Communications 5, 5875 (2014), arXiv:1404.7658 .
  * da silva Neto  _et al._ [2015] E. H. da silva Neto, R. Comin, F. He, R. Sutarto, Y. Jiang, R. L. Greene, G. A. Sawatzky, and A. Damascelli, Charge ordering in electron doped superconductor Nd2-xCexCuO4, Science 347, 282 (2015).
  * da Silva Neto  _et al._ [2016] E. H. da Silva Neto, B. Yu, M. Minola, R. Sutarto, E. Schierle, F. Boschini, M. Zonno, M. Bluschke, J. Higgins, Y. Li, G. Yu, E. Weschke, F. He, M. L. Tacon, R. L. Greene, M. Greven, G. a. Sawatzky, B. Keimer, and A. Damascelli, Doping dependent charge order correlations in electron-doped cuprates, Science Advances 2, e1600782 (2016), arXiv:1607.06094 .
  * Kang  _et al._ [2019] M. Kang, J. Pelliciari, A. Frano, N. Breznay, E. Schierle, E. Weschke, R. Sutarto, F. He, P. Shafer, E. Arenholz, M. Chen, K. Zhang, A. Ruiz, Z. Hao, S. Lewin, J. Analytis, Y. Krockenberger, H. Yamamoto, T. Das, and R. Comin, Evolution of charge order topology across a magnetic phase transition in cuprate superconductors, Nature Physics 15, 335 (2019).
  * Comin [2016] R. Comin, Resonant X-ray scattering studies of charge order in cuprates, Annu. Rev. Condens. Matter Phys. 7, 369 (2016), arXiv:1406.3533 .
  * Kivelson  _et al._ [2003] S. A. Kivelson, I. P. Bindloss, E. Fradkin, V. Oganesyan, J. M. Tranquada, A. Kapitulnik, and C. Howald, How to detect fluctuating stripes in the high-temperature superconductors, Reviews of Modern Physics 75, 1201 (2003).
  * Tranquada [2015] J. M. Tranquada, Exploring intertwined orders in cuprate superconductors, Physica B: Condensed Matter Special Issue on Electronic Crystals (ECRYS-2014), 460, 136 (2015).
  * Fradkin  _et al._ [2015] E. Fradkin, S. A. Kivelson, and J. M. Tranquada, Colloquium : Theory of intertwined orders in high temperature superconductors, Rev. Mod. Phys. 87, 457 (2015).
  * Kivelson  _et al._ [1998a] S. Kivelson, E. Fradkin, and V. Emery, Electronic liquid-crystal phases of a doped Mott insulator, Nature 393, 550 (1998a), arXiv:9707327v1 [arXiv:cond-mat] .
  * Vojta and Sachdev [1999] M. Vojta and S. Sachdev, Charge Order, Superconductivity, and a Global Phase Diagram of Doped Antiferromagnets, Phys. Rev. Lett. 83, 3916 (1999).
  * Sachdev and La Placa [2013] S. Sachdev and R. La Placa, Bond order in two-dimensional metals with antiferromagnetic exchange interactions, Physical Review Letters 111, 027202 (2013), arXiv:1303.2114 .
  * Davis and Lee [2013] J. C. S. Davis and D.-H. Lee, Concepts relating magnetic interactions, intertwined electronic orders and strongly correlated superconductivity, Proceedings of the National Academy of Sciences 110, 10.1073/pnas.1316512110 (2013), arXiv:1309.2719 .
  * Lee [2014] P. A. Lee, Amperean pairing and the pseudogap phase of cuprate superconductors, Physical Review X 4, 031017 (2014), arXiv:1401.0519 .
  * Wang and Chubukov [2014] Y. Wang and A. Chubukov, Charge-density-wave order with momentum (2Q,0) and (0,2Q) within the spin-fermion model: Continuous and discrete symmetry breaking, preemptive composite order, and relation to pseudogap in hole-doped cuprates, Phys. Rev. B 90, 35149 (2014).
  * Dalla Torre  _et al._ [2015] E. G. Dalla Torre, Y. He, D. Benjamin, and E. Demler, Exploring quasiparticles in high-Tc cuprates through photoemission, tunneling, and x-ray scattering experiments, New Journal of Physics 17, 22001 (2015), arXiv:1312.0616 .
  * Duong and Das [2017] L. Q. Duong and T. Das, Correlation between Fermi arc and charge order resulting from the momentum-dependent self-energy correction in cuprates, Physical Review B 96, 125154 (2017).
  * Keimer  _et al._ [2015] B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, From quantum matter to high-temperature superconductivity in copper oxides, Nature 518, 179 (2015).
  * Radousky [1992] H. B. Radousky, A review of the superconducting and normal state properties of Y1-xPrxBa2Cu3O7, Journal of Materials Research 7, 1917 (1992).
  * Liang  _et al._ [2006] R. Liang, D. A. Bonn, and W. N. Hardy, Evaluation of CuO2 plane hole doping in YBa2 Cu3 O6+x single crystals, Physical Review B 73, 180505(r) (2006).
  * Kivelson  _et al._ [1998b] S. A. Kivelson, E. Fradkin, and V. J. Emery, Electronic liquid-crystal phases of a doped Mott insulator, Nature 393, 550 (1998b).
  * Zhao  _et al._ [2019] H. Zhao, Z. Ren, B. Rachmilowitz, J. Schneeloch, R. Zhong, G. Gu, Z. Wang, and I. Zeljkovic, Charge-stripe crystal phase in an insulating cuprate, Nature Materials 18, 103 (2019).
  * Akhavan [2002] M. Akhavan, The question of Pr in HTSC, Physica B 321, 265 (2002).
  * Soderholm  _et al._ [1987] L. Soderholm, K. Zhang, D. G. Hinks, M. A. Beno, J. D. Jorgensen, C. U. Segre, and I. K. Schuller, Incorporation of Pr in YBa2Cu3O7-d: electronic effects on superconductivity, Nature 328, 604 (1987).
  * Neumeier  _et al._ [1989] J. J. Neumeier, T. BjÃ¸rnholm, M. B. Maple, and I. K. Schuller, Hole filling and pair breaking by pr ions in y1-xprxba2cu3o6.95±0.02, Physical Review Letters 63, 2516 (1989), publisher: American Physical Society.
  * Sun  _et al._ [1994] A. G. Sun, L. M. Paulius, D. A. Gajewski, M. B. Maple, and R. C. Dynes, Electron tunneling and transport in the high-$t_{c}$ superconductor y1-xprxba2cu3o7-δ, Physical Review B 50, 3266 (1994), publisher: American Physical Society.
  * Frano  _et al._ [2016] A. Frano, S. Blanco-Canosa, E. Schierle, Y. Lu, M. Wu, M. Bluschke, M. Minola, G. Christiani, H. U. Habermeier, G. Logvenov, Y. Wang, P. A. van Aken, E. Benckiser, E. Weschke, M. Le Tacon, and B. Keimer, Long-range charge-density-wave proximity effect at cuprate/manganate interfaces, Nature Materials 15, 831 (2016).
  * Wojek  _et al._ [2012] B. M. Wojek, E. Morenzoni, D. G. Eshchenko, A. Suter, T. Prokscha, H. Keller, E. Koller, O. Fischer, V. K. Malik, C. Bernhard, and M. Döbeli, Magnetism, superconductivity, and coupling in cuprate heterostructures probed by low-energy muon-spin rotation, Phys. Rev. B 85, 024505 (2012).
  * [39] See Supplemental Material at [DOI] for Methods and Supplemental Figures, .
  * Fehrenbacher and Rice [1993] R. Fehrenbacher and T. M. Rice, Unusual electronic structure of PrBa2Cu3O7, Physical Review Letters 70, 3471 (1993).
  * Arpaia  _et al._ [2019] R. Arpaia, S. Caprara, R. Fumagalli, G. D. Vecchi, Y. Y. Peng, E. Andersson, D. Betto, G. M. D. Luca, N. B. Brookes, F. Lombardi, M. Salluzzo, L. Braicovich, C. D. Castro, M. Grilli, and G. Ghiringhelli, Dynamical charge density fluctuations pervading the phase diagram of a Cu-based high-Tc superconductor, Science 365, 906 (2019).
  * Hill [2000] J. P. Hill, X-ray-scattering study of copper magnetism in nonsuperconducting PrBa2Cu3O6.92, Physical Review B 61, 1251 (2000).
  * Zhang  _et al._ [2013] H. Zhang, N. Gauquelin, G. A. Botton, and J. Y. Wei, Attenuation of superconductivity in manganite/cuprate heterostructures by epitaxially-induced CuO intergrowths, Applied Physics Letters 103, 052606 (2013).
  * Zhang  _et al._ [2018] H. Zhang, N. Gauquelin, C. McMahon, D. G. Hawthorn, G. A. Botton, and J. Y. Wei, Synthesis of high-oxidation Y-Ba-Cu-O phases in superoxygenated thin films, Physical Review Materials 2, 33803 (2018).
  * Paturi  _et al._ [2004] P. Paturi, M. Peurla, K. Nilsson, and J. Raittila, Crystalline orientation and twin formation in YBCO thin films laser ablated from a nanocrystalline target, Superconductor Science and Technology 17, 564 (2004).
  * McCoy  _et al._ [2019] S. J. McCoy, E. Y. Cho, H. Li, and S. A. Cybart, Ho-Ba-Cu-O Thin Films for Superconductive Electronics, IEEE International Superconductive Electronics Conference , 1 (2019).
  * Vovk  _et al._ [2013] R. V. Vovk, Z. F. Nazyrov, I. L. Goulatis, and A. Chroneos, Metal-to-insulator transition in Y1-xPrxBa 2Cu3O7-$\delta$ single crystals with various praseodymium contents, Physica C: Superconductivity and its Applications 485, 89 (2013).
  * Antognazza  _et al._ [1990] L. Antognazza, J. M. Triscone, O. Brunner, L. Mieville, M. G. Karkut, and O. Fischer, Layer by layer preparation of YPBCO alloy thin films, Journal of the less common metals 165, 344 (1990).
  * Peng  _et al._ [1989] J. L. Peng, P. Klavins, R. N. Shelton, H. B. Radousky, P. A. Hahn, and L. Bernardez, Upper critical field and normal-state properties of single-phase YPBCO compounds, Physical Review B 40, 4517 (1989).
  * Prokscha  _et al._ [2008] T. Prokscha, E. Morenzoni, K. Deiters, F. Foroughi, D. George, R. Kobler, A. Suter, and V. Vrankovic, The new $\mu$e4 beam at psi: A hybrid-type large acceptance channel for the generation of a high intensity surface-muon beam, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 595, 317 (2008).
  * Suter and Wojek [2012] A. Suter and B. Wojek, Musrfit: A free platform-independent framework for $\mu$sr data analysis, Physics Procedia 30, 69 (2012), 12th International Conference on Muon Spin Rotation, Relaxation and Resonance ($\mu$SR2011).
  * Yaouanc and Dalmas de Réotier [2011] A. Yaouanc and P. Dalmas de Réotier,  _Muon Spin Rotation, Relaxation, and Resonance_ (Oxford University Press, 2011).


# Data-driven battery operation for energy arbitrage using rainbow deep
reinforcement learning

Daniel J. B. Harrold Email<EMAIL_ADDRESS>School of Computing and
Mathematics, Keele University, United Kingdom  Jun Cao  School of Geography,
Geology, and the Environment, Keele University, United Kingdom  Zhong Fan
School of Computing and Mathematics, Keele University, United Kingdom

###### Abstract

As the world seeks to become more sustainable, intelligent solutions are
needed to increase the penetration of renewable energy. In this paper, the
model-free deep reinforcement learning algorithm Rainbow Deep Q-Networks is
used to control a battery in a small microgrid to perform energy arbitrage and
more efficiently utilise solar and wind energy sources. The grid operates with
its own demand and renewable generation based on a dataset collected at Keele
University, as well as using dynamic energy pricing from a real wholesale
energy market. Four scenarios are tested including using demand and price
forecasting produced with local weather data. The algorithm and its
subcomponents are evaluated against two continuous control benchmarks with
Rainbow able to outperform all other method. This research shows the
importance of using the distributional approach for reinforcement learning
when working with complex environments and reward functions, as well as how it
can be used to visualise and contextualise the agent’s behaviour for real-
world applications.

Keywords— Actor-Critic Methods, Deep Q-Networks, Demand Response, Microgrids,
Renewable Energy

##  1 Introduction

As the world seeks to reduce its carbon emissions, increasing the penetration
of renewable energy sources (RES) is a necessity to move towards a more
sustainable future [1]. However, despite RES being some of the most cost-
effective forms of energy generation [2], the intermittent nature of solar
photovoltaic (PV) and wind turbine (WT) generation means that non-renewable
sources such as natural gas and coal make up the majoirty of our energy
networks [3]. Therefore, intelligent solutions are required to manage energy
supply and demand more resourcefully.

Smart energy networks combine power engineering with information technologies
to manage utilities for improving the efficiency, portability, stability, and
security of the grid through the use of demand response (DR) to manipulate
energy usage [4]. The primary objectives of DR are to increase the utilisation
of renewable energy sources (RES) and facilitate energy arbitrage using energy
storage systems (ESS), in addition to performing load levelling, peak shaving,
and frequency regulation [5]. Although difficult to enable DR at a utility
grid scale due to scalability, it is much more feasible to implement in a
localised microgrid with its own ESSs and RESs [6]. However, this will require
an intelligent energy management system to control the ESS able to observe the
fluctuating energy demand, intermittent RESs, and volatile wholesale energy
pricing markets.

For this management system, we propose reinforcement learning (RL), a branch
of machine learning in which an agent learns to interact with its environment
to find an optimal control policy [7]. The agent iteratively chooses an action
based on observations of the environment to receive a reward with the goal to
maximise its total future reward. In addition, artificial neural networks
(ANN) can be used to combine RL with deep learning for solutions to more
complex problems. The model-free nature of RL algorithms makes them incredibly
versatile and easy to implement in situations where there is little
information or data available as the agent learns effectively through trial-
and-error.

In this paper, we present the use of a state-of-the-art RL algorithm to
control an ESS for energy arbitrage in a microgrid. The aim of the agent is to
reduce energy costs and consequently increase the utilisation of RES.

###  1.1 Related Literature

Research into applying RL to DR and smart energy networks management has
increased significantly over the past decade [8, 9].

Before Mnih et al. [10] first introduced modern deep RL, basic tabular methods
were the most widely used. Kuznetsova et al. [11] used Q-learning [12] to
control an ESS connected to WT generation to increase its utilisation under a
number of scenarios. However, an issue with tabular RL is that the observation
values must be from a discrete state-space. Kofinas et al. [13] use fuzzy
logic to adapt Q-learning to a continuous state and action space, but more
modern deep RL algorithms can fulfil this too.

Authors generally have moved away from tabular methods to deep RL methods
using ANNs as functional approximators, such as Deep Q-Networks (DQN) [10].
Examples of using DQN for energy arbitrage and energy cost reduction include
François-Lavet et al. [14] with a hybrid ESS setup in a solar microgrid, Bui
et al. [15] with a variant of DQN considering uncertainities in the
environment, and Cao et al. [16] with another more advanced variant based on a
realistic lithium-ion battery degradation model. Although value-function
methods such as Q-learning and DQN are simple algorithms with robust learning
properties, a limitation is that the agent must select from a discrete set of
actions which could be considered unsuitable for continuous control problems;
such as ESS operation.

As a result, many authors have instead investigated actor-critic methods.
These are a branch of RL where one ANN is used to choose a continuous value
for the action and a second ANN is used to evaluate performance. Examples
include Mocanu et al. [17] using deep deterministic policy gradient (DDPG)
[18] to schedule a building’s energy management system, Zhang et al. [19] used
proximal policy optimisation (PPO) [20] for energy cost reduction using ESS
control in a AC/DC hybrid microgrid, and Pinto et al. [21] used soft actor-
critic [22] for cost reduction and load levelling across multiple smart
buildings. However, actor-critic methods are significantly more sensitive to
hyperparameter tuning than value-function methods and can suffer from unstable
and unpredictable learning [7]. Therefore, the correct balance must be found
between control resolution and robust learning.

###  1.2 Contributions

The contributions of this work are as follows:

  * •

The paper presents the state-of-the-art deep RL algorithm Rainbow DQN [23] for
the control agent. To the best of our knowledge, this is the first time the
algorithm has been used for a DR application and is tested against two
continuous control benchmarks: an actor-critic method and a linear programming
model.

  * •

The paper assesses the distributional approach of deep RL where the agent
evaluates performance using a value distribution, presenting the benefits of
using this approach for DR applications.

  * •

The methods are evaluated under four scenarios with different sizes of state
and action-spaces to evaluate how different algorithms utilise forecasting as
well as having a greater resolution of control.

  * •

The demand and weather data used in the simulations were collected from the
same site at Keele University Campus as part of the Smart Energy Network
Demonstrator project [24], as well as using dynamic energy prices from a
wholesale energy market [25].

###  1.3 Structure

The rest of this paper is organised into a description of the microgrid
environment in Section 2, the background behind the RL methodology in Section
3, outlining the case study in Section 4, results and discussions in Section 5
with final conclusions in Section 6.

##  2 Microgrid Description

Figure 1: Basic schematic of the microgrid.

Microgrids are localised energy networks with their own demand and RES that
are able to trade energy with the main utility grid or operate as an
independent island [6]. This section will describe the microgrid used in the
future case studies and how it is modelled as a Markov decision process (MDP)
for RL.

###  2.1 Environment Description

The microgrid is fitted with a large lithium-ion battery ESS as well as both
PV and WT generation. It is connected to the main utility grid that sets
dynamic energy prices from which the microgrid can import energy from or
export to. AC and DC power lines are present connected via inverters, as well
as transformers between the microgrid and both the utility grid and the WT.
The basic schematic of the environment is shown in Figure 1.

####  2.1.1 Demand and Pricing

The demand data collected at Keele University Campus ranges from 00:00 January
1st 2014 up to and including 23:00 December 31st 2017. The raw data is
separated into different residential, industrial, and commercial sites as well
as readings for key individual buildings. The demand used for this simulation
is from the three main incomer substations into the campus with the half-
hourly readings summed to match the frequency of the hourly weather data. All
loads are AC.

The microgrid operates under a dynamic energy pricing scheme. The energy
prices from the utility grid are from a real day-ahead energy trading market
covering the UK [25] but also with a set maximum price of $250\text{GBP/MWh}$.
The price for exporting to the grid is the same as the import price.

####  2.1.2 Renewable Energy Sources

The output for the renewable generation is modelled using weather data
collected at the Keele University weather station. The WT’s output
$\text{X}^{\text{WT}}$ uses power curve modelling with wind speed $v$:

$\text{X}^{\text{WT}}=\begin{cases}0&v<v_{\text{ci}}\text{ or
}v>v_{\text{co}}\\\ \frac{1}{2}\rho\pi r^{2}c_{p}v^{3}&v_{\text{ci}}\leqslant
v<v_{r}\\\ \text{X}_{r}^{\text{WT}}&v_{r}\leqslant v\leqslant
v_{\text{co}}\end{cases}$ (1)

where $v_{ci}=3\text{ms}^{-1}$, $v_{r}=12\text{ms}^{-1}$, and
$v_{co}=25\text{ms}^{-1}$ are the cut-in, rated, and cut-out wind-speeds
respectively. These values are comparable to a 1MW turbine with a blade radius
of $r=30\text{m}$ and a power coefficient $c_{p}=0.4$ [26]. This microgrid
considers two of these turbines for a combined maximum wind capacity of
$2\text{MW}$. The same power curve approach is used by both Kuznetsova et al.
[11] and Zhang et al. [19] for their RES output. The solar output is
calculated using hourly solar radiation data and scaled to model a solar farm
with a maximum capacity of 5MW. The output of the WT is AC while the output of
the PV is DC.

####  2.1.3 ESS

(a) Battery efficiency when charging.

(b) Battery efficiency when discharging.

(c) Transformer and inverter efficiencies.

Figure 2: Efficiency profiles of different components of the microgrid.

The charge of the ESS is bound between 0 and its maximum capacity
$\text{C}_{\text{max}}=5\text{MWh}$, with the amount of energy the battery can
charge or discharge each step bound between 0 and its maximum power
$\text{X}_{\text{max}}=2\text{MW}$. The charge of the ESS $c_{t}$ is
calculated using:

$c_{t}=(x_{t}\eta^{\text{ESS}}_{t}+c_{t-1})\eta^{\text{SDC}}_{t}$ (2)

The ESS’s charging and discharging efficiencies are nonlinear and modelled as
a function of the current charge and operating power using the same
approximation as Morstyn et al. [27] with the charging and discharging
profiles shown in Figures 2(a) and 2(b) respectively. For this, the 5MWh ESS
is modelled as 50 100kWh lithium-ion batteries which charge and discharge in
unison. The ESS also has a linear self-discharge efficiency of
$\eta_{\text{SDC}}=0.999$ such that it loses $0.1\%$ of its charge every hour.
The output of the ESS is DC.

####  2.1.4 Transformer and Inverter

The case study microgrid is a hybrid AC/DC network with inverters between the
two power lines and transformers to step-down the voltage to the microgrid.
The efficiency profiles are nonlinear and given as a function of the load
factor; the current load divided by the rated load. The rated loads of each of
the transformers and inverters can be found in Figure 1 with the efficiency
profiles modelled using a polynomial estimation, shown in Figure 2(c).
Although the real efficiency of these devices at low load factor would
approach 0%, the minimum efficiency is capped at 10% to prevent divisions by 0
during environment calculations.

Note that the ESS and PV generation are connected to both 2MW and 5MW rated
inverters as well as to each other. When energy is to pass between the AC and
DC lines, the control agent will automatically use whichever inverter that
would result in lower power losses. Mbuwir et al. [28] used a similar profile
for inverters [29] for RL in a microgrid with solar energy.

###  2.2 Markov Decision Process

For RL to operate the microgrid environment, the control process must be
modelled as an MDP. This means the control process must satisfy the Markov
property and can be generalised to a a state-space $\mathcal{S}$, action-space
$\mathcal{A}$, reward function $\mathcal{R}$, and a state transition function.
The rest of this section will explain the states, actions, and rewards for
this case study.

####  2.2.1 States

The simulations in Section 5 use agents with two different state-spaces. The
basic agents makes a total of 8 observations: the ESS charge $c_{t}$, grid
demand $\text{X}^{\text{D}}$, energy price P, PV generation
$\text{X}^{\text{PV}}$, WT generation $\text{X}^{\text{WT}}$, as well as the
hour in the day, week, and if it was a workday or not. The forecasting agents
also observe forecasts calculated using ANNs based on previous values and
weather data. This means the advanced agents observe the same 8 observations
as well as predictions for the next hour of demand, solar and wind generation,
and price for 12 total observations. As the agent does not observe every
aspect of the underlying state, the process is treated as a partially
observable MDP.

Each of the values are scaled and bound between 0 and 1 to improve the
training stability of the ANNs. The charge, demand, PV and WT generation are
all divided by the maximum capacity of the battery 5MWh. The dynamic energy
price is divided by the set price limit of 250GBP/MWh. The hourly values begin
at 0 and end at 1 for the day or week. The workday value is set to 0 for
weekends or UK bank holidays; otherwise 1.

####  2.2.2 Actions

The action the agent selects controls the ESS’s power $x$ at that step. The
agent can select from several discrete actions where the ESS can charge
($x>0$), discharge ($x<0$), or can remain idle ($x=0$). Two different actions-
spaces of different size are considered: the first with 5 actions and the
second with 9 actions. The actions are evenly spaced between the maximum
charging and discharging power such that the action-space of size 5 is:

$\mathcal{A}=[2\text{MW},1\text{MW},0\text{MW},-1\text{MW},-2\text{MW}]$ (3)

Although a real ESS would be able to charge or discharge by any continuous
value from its minimum to its maximum power, restricting the size of the
action-space improves learning and performance in value-function RL methods as
the agent has to explore fewer state-action pairs [7].

Due to the configuration of the grid shown in Figure 1 the ESS will always
charge from RES first before importing from the utility grid to minimises
energy losses from the transformer and inverter.

####  2.2.3 Reward

The reward at each step is equal to the energy cost savings made by the agent.
As the ESS will always look to charge from the RES generation first, the
excess PV generation $\text{X}^{\text{PV+}}$ and excess WT generation
$\text{X}^{\text{WT+}}$ represent the amount of RES left after the charging
has been considered. From this, the demand of the DC portion of the grid
$\text{X}^{\text{dc}}$ is calculated as:

$\text{X}^{\text{dc}}=(x\eta^{\text{ESS}}-\text{X}^{\text{PV+}})\eta^{\text{inv}}_{\text{[2,5]MW}}$
(4)

This is then used to calculate the total amount of energy being imported from
the utility grid $\text{X}^{\text{in}}$:

$\text{X}^{\text{in}}=(\text{X}^{\text{dc}}+\text{X}^{\text{D}}-\text{X}^{\text{WT}}\eta^{\text{tra}}_{\text{2MW}})\eta^{\text{tra}}_{\text{2MW}}$
(5)

During preliminary studies it was found that the agent would fail to learn an
optimal policy properly from this value alone because the agent would be
rewarded if the demand was low with high renewable generation and punished in
the reverse case. This is an issue because these are environmental factors
which the agent has no control over so leads to problems with credit
assignment. Therefore, a base demand $\text{X}^{\text{base}}$ is used to
normalise the reward:

$\text{X}^{\text{base}}=(\text{X}^{\text{D}}-\text{X}^{\text{PV}}\eta^{\text{inv}}_{\text{[2,5]MW}}-\text{X}^{\text{WT}}\eta^{\text{tra}}_{\text{2MW}})\eta^{\text{tra}}_{\text{2MW}}$
(6)

The agent should also be discouraged from taking an action that would exceed
the ESS’s charge boundaries. In this event, the agent receives a negative
reward $R_{\text{punish}}$ equal to the squared amount it would have exceeded
the boundaries divided by the maximum power $\text{X}_{\text{max}}$. This
gives a combined reward function:

$\mathcal{R}=\frac{\text{P}_{t}(\text{X}^{\text{base}}_{t}-\text{X}^{\text{in}}_{t})}{\text{P}_{\text{max}}\text{C}_{\text{max}}}-R_{\text{punish}}$
(7)

The energy savings value is divided by the product of the maximum price
$\text{P}_{\text{max}}$ and the ESS capacity $\text{C}_{\text{max}}$ so that
the reward largely remains between -1 and 1 for stable learning. This reward
function is different to using pure energy arbitrage in that the nonlinear
efficiencies mean the demand and RES are taken into account as well. If the
efficiencies were constant then the reward function would cancel them out and
the agent would be performing straight arbitrage. However, in this case the
agent considers both arbitrage as well as how to optimally utilise RES taking
into account the battery, inverter, and transformer efficiencies.

##  3 Reinforcement Learning Methodology

This section will explore the theory behind RL and introduce the specific
algorithms used in the simulations later in the paper.

The key benefit of RL is that the methods are tpically model-free meaning the
state transition function and reward function do not need to be known. This
makes them incredibly versatile and easy to implement as the agents learn
entirely from its own experience; effectively through trial-and error [7].
Also, model-based solutions are only as accurate as the model they are based
on.

###  3.1 Fundamentals

Each time-step $t$ in RL, the agent observes the current state of the
environment $s_{t}$ from its state-space $\mathcal{S}$ and selects an action
$a_{t}$ from an action-space $\mathcal{A}$ following a learnt policy. The
agent then receives a reward $r_{t}$ from a reward function $\mathcal{R}$ and
transitions to a new state $s_{t+1}$ following the state transition function.
This is shown in Figure 3. The goal of the agent is to learn an optimal policy
that maximises its total discounted future reward, scaled with a discount
factor $\gamma$.

Figure 3: Reinforcement learning control process.

$R_{t}=r_{t}+\gamma
r_{t+1}+\gamma^{2}r_{t+2}+\ldots=\sum\limits_{k=0}^{\infty}\gamma^{k}r_{t+k}$
(8)

Value-function methods evaluate performance by assigning value to each state
or state-action pair. The state-value $V(s)$ represents how valuable it is to
be in a given state and the action-value $Q(s,a)$ how valuable it is to take a
particular action in that state. The control problem is solved recursively
using the Bellman equation [30]:

$Q(s_{t},a_{t})=\mathbb{E}\mathcal{R}(s_{t},a_{t})+\gamma\mathbb{E}Q(s_{t+1},a_{t+1})$
(9)

In these methods, the agent will follow a deterministic policy where the agent
will always choose the action it predicts has the highest action-value. For
exploration, value-function methods typically follow an $\epsilon$-greedy
policy where each time-step the agent has a probability of $\epsilon$ to
choose a random action and $(1-\epsilon)$ to select the action with the
highest action-value. The $\epsilon$ value will typically decrease over time
where maximising the return becomes more important than exploration.

###  3.2 Algorithms

This section will outline the methods used in the simulation. Each one is a
valuation function method using an ANN as a functional approximator so are
able to process continuous observation values, rather than tabular RL methods
which require discrete state and action-spaces.

####  3.2.1 Deep Q-Networks (DQN)

Q-learning [12] is a tabular temporal difference (TD) learning method and had
been the most widely used RL algorithm prior to modern deep RL methods. DQN
[10] is an adaptation of Q-learning using a functional approximator to predict
the action-values and then updating the network weights $\theta$ using
stochastic gradient descent. TD methods such as these learn using the TD-error
$\delta$ calculated as the difference between the current state-value and a
target value:

$\delta_{i}\doteq r_{i}+\gamma Q(s_{i+1},a_{i+1})-Q(s_{i},a_{i})$ (10)

It is common practice to use a stable target network to calculate the target
action-value $\hat{Q}(s,a)$ with network weights $\hat{\theta}$ frozen from
the evaluation network at a previous step and then update
$\hat{\theta}\leftarrow\theta$ every $C$ steps. A target value $y_{i}$ for the
Bellman update is calculated by:

$y_{i}=r_{i}+\gamma\max_{a\in\mathcal{A}}\hat{Q}_{\hat{\theta}}(s_{i+1},a)$
(11)

Experience replay [31] is also used to improve learning through the use of a
replay buffer. Each step, the transition $\langle
s_{t},a_{t},r_{t},s_{t+1}\rangle$ is stored in the buffer and then the
evaluation network is trained by sampling a random batch of those tuples. The
evaluation network is updated using the mean squared error loss which
corresponds to the TD-error $\delta$:

$\mathcal{L}_{i}(\theta)=\mathbb{E}\left[(y_{i}-Q_{\theta}(s_{i},a_{i}))^{2}\right]$
(12)

Many environments will reach a terminal state and then reset for the next
episode. In this microgrid environment, there are no terminal states so the
algorithm continues until finished. The pseudocode for DQN is shown in
Algorithm 1.

Algorithm 1 Deep Q-Networks

Initialise network $Q_{\theta}(s,a)$ with random weights $\theta$

Initialise target network $\hat{Q}_{\hat{\theta}}(s,a)$ and set
$\hat{\theta}\leftarrow\theta$

Initialise replay buffer

for $t=0,T$ do

Observe state $s_{t}$

if probability $p<\epsilon$, choose random action $a_{t}$

else $a_{t}=\operatorname*{arg\,max}_{a}Q_{\theta}(s_{t},a_{t})$

Execute $a_{t}$, observe reward $r_{t}$ and next state $s_{t+1}$

Store transition $\langle s_{t},a_{t},r_{t},s_{t+1}\rangle$

Sample batch $\langle s_{i},a_{i},r_{i},s_{i+1}\rangle$ randomly

Calculate targets
$y_{i}=r_{i}+\gamma\max_{a}\hat{Q}_{\hat{\theta}}(s_{i+1},a_{i})$

Loss
$\mathcal{L}_{i}(\theta)=\mathbb{E}\left[(y_{i}-Q_{\theta}(s_{i},a_{i}))^{2}\right]$

Perform gradient descent on $\theta$ using loss

if $\epsilon>\epsilon_{min}$, $\epsilon\leftarrow(1-\epsilon_{decay})\epsilon$

Every $C$ steps, set $\theta^{-}\leftarrow\theta$

end for

####  3.2.2 Double Deep Q-Networks (DDQN)

An issue with standard DQN is that the target network selects the next-best
action with the maximisation step which is then itself used for evalualtion.
This can lead to overestimations in the value calculations causing the agent
to prefer those particular actions. DDQN [32] decouples this by removing the
maximisation step and instead using the evaluation network to predict the
next-best action:

$a^{*}=\operatorname*{arg\,max}_{a\in\mathcal{A}}Q_{\theta}(s_{i+1},a)$ (13)

The benefit of this is that the agent is not automatically assumed to always
take the most valuable action from the next state. The updated target is then
calculated by:

$y_{i}=r_{i}+\gamma\hat{Q}_{\hat{\theta}}(s_{i+1},a^{*})$ (14)

Bui et al. [15] use DDQN to manage a community battery ESS in a microgrid,
showing how the algorithm can perform well despite the uncertainties in
microgrid environments.

####  3.2.3 Double Dueling Deep Q-Networks (D3QN)

Using only $Q(s,a)$ for the Bellman update does not take into account the
inherit value of the state $V(s)$. Therefore, the advantage function $A(s,a)$
separates the two to give greater merit for taking more valuable actions when
in low value states. Dueling DQN [33] builds on DQN by using dueling ANN
architecture using a shared encoder with two output streams: one to estimate
$A(s,a)$ for each of the actions with weights $\psi$ and the other to estimate
$V(s)$ with weights $\xi$. From this, $Q(s,a)$ can be calculated with:

$Q_{\theta}(s_{i},a_{i})=V_{\xi}(\phi_{i})+A_{\psi}(\phi_{i},a_{i})-\bar{A}_{\psi}(\phi_{i})$
(15)

using the mean advantage value across all actions $\bar{A}(\phi_{i})$ for
identifiability. Here, $\phi$ is the output of the shared encoder which are
the hidden layers of the functional approximator that both output streams
share:

$\phi_{i}=f_{\zeta}(s_{i})$ (16)

Therefore, the weights of $Q(s,a)$ are $\theta=\\{\zeta,\xi,\psi\\}$ in their
concatenation.

In previous work [34], we used D3QN to control a battery in a simple microgrid
setup with RES under a dynamic wholesale energy price where it outperformed
both DQN and DDQN.

####  3.2.4 Prioritised Experience Replay DQN (PER-DQN)

Random sampling of the replay buffer can be inefficient as newer transitions
will not be sampled as frequently as older transitions. A prioritised
experience replay (PER) [35] buffer samples transitions by assigning a
priority $p$ to each transition such that more important transitions are
sampled quicker and more frequently. Each priority is equal to the absolute
TD-error $\delta$ of its most recent sample plus a small offset with a
probability $P$ of that transition being sampled:

$P(i)=\frac{p_{i}^{\alpha}}{\sum_{k}p_{k}^{\alpha}}$ (17)

with a priority scale $\alpha$. An increasing importance scale $\beta$ is then
used to calculate the importance weight $w$ which scales the loss in training
to give more attention to transitions with larger errors:

$w_{i}=\left(\frac{1}{N}\times\frac{1}{P(i)}\right)^{\beta}/\max_{j}w_{j}$
(18)

Yang et al. [36] used D3QN with PER to control an ESS in a wind farm to
account for uncertainties in the wind power generation and wholesale energy
price.

####  3.2.5 Multistep DQN (MS-DQN)

In some environments, taking the return over multiple steps for the target
calculation can increase learning speed [37] as the agent follows the current
learnt policy for $n$-steps and then performs the backup over that time-
period. The value of $n$ should be tuned for different environments. The new
target $y_{i}$ is then calculated using:

$y_{i}=R_{t:t+n-1}+\gamma^{n}\max_{a\in\mathcal{A}}\hat{Q}\left(s_{t+n},a;\theta^{-}\right)$
(19)

Kuznetsova et al. [11] used a 3-step return for tabular Q-learning for 2-step-
ahead battery scheduling with scenario-based wind power generation. However,
to our knowledge, there is no case of MS-DQN being used for DR.

####  3.2.6 NoisyNet DQN (NN-DQN)

This uses a pure-greedy policy but replaces the linear layers of the DQN ANN
with noisy layers that inject noise into the $Q(s,a)$ calculation. The benefit
of this is to an $\epsilon$-greedy policy is that the functional approximator
can learn to ignore the noise over time in different parts of the state-space
rather than having a policy-wide exploration rule. Using Gaussian noise
parameters $\epsilon^{w}$ and $\epsilon^{b}$ and an initial sigma value
$\sigma_{0}$, the output of the noisy layer is calculated as:

$y=\left(\mu^{w}+\sigma^{w}\odot\epsilon^{w}\right)x+\mu^{b}+\sigma^{b}\odot\epsilon^{b}$
(20)

Cao et al. [16] used NN-DDQN for energy arbitrage using a realistic lithium-
ion battery model and found faster learning properties than both DQN and D3QN.

####  3.2.7 Categorical DQN (C51)

An issue with using the expectation value for $Q(s,a)$ is that it can over-
generalise in environments with a non-deterministic reward function [38].
Categorical DQN removes the expectation step in the Bellman equation and
instead estimates the value distribution $Z(s,a)$:

$Z(s_{t},a_{t})\stackrel{{\scriptstyle D}}{{=}}\mathcal{R}(s_{t},a_{t})+\gamma
Z(s_{t+1},a_{t+1})$ (21)

A number of atoms are assigned to each action for $N_{\text{actions}}\times
N_{\text{atoms}}$ output neurons with the probability distribution of the
return $d(s,a)$ calculated using a softmax applied separately across the atoms
for each action. The algorithm C51 [38] uses 51 atoms across each action. The
probability mass on each atom is equal to:

$d_{i}(s,a)=\frac{\exp(\theta_{i}(s,a))}{\sum_{j}\exp(\theta_{j}(s,a))}$ (22)

Each of these probability masses can then be projected onto a support $z$ with
a value equally spaced between a minimum return $v_{\text{min}}$ and a maximum
return $v_{\text{max}}$:

$z_{i}=v_{\text{min}}+(i-1)\frac{v_{\text{max}}-v_{\text{min}}}{N_{\text{atoms}}-1}$
(23)

The sum of these probability masses on each support is equal to $Q(s,a)$. The
target used for training is the projected target distribution calculated using
the categorical algorithm [38]:

$(\Phi\hat{\mathcal{T}}Z_{\hat{\theta}}(s,a))_{i}=\sum^{N_{\text{atoms}}}_{j=0}\left[1-\frac{|[\hat{\mathcal{T}}z_{j}]^{v_{\text{max}}}_{v_{\text{min}}}-z_{i}|}{\Delta
z}\right]^{1}_{0}\hat{d}_{j}(s_{i+1},a^{*})$ (24)

where $\Phi$ is the projection of the distribution onto $z$ and $\mathcal{T}$
is the distributional Bellman operator. As this generates a discrete
probability distribution, the Kullback-Leibler divergence between the
predicted value and the projected target is used for the loss function:

$\mathcal{L}_{i}(\theta)=D_{\text{KL}}\left(\Phi\hat{\mathcal{T}}Z_{\hat{\theta}}(s_{i},a_{i})\|Z_{\theta}(s_{i},a_{i})\right)$
(25)

To the best of our knowledge, no authors have used C51 or a similar
distributional RL algorithm for DR applications.

####  3.2.8 Rainbow DQN

Rainbow [23] combines all of the components just described into a single
algorithm. The noisy linear layer has $N_{\text{actions}}\times
N_{\text{atoms}}$ output neurons for the advantage stream and
$N_{\text{atoms}}$ for the state-value stream to implement the dueling
architecture into the value distribution calculation. The new distribution
probability on an atom is then calculated by:

$d_{i}(s,a)=\frac{\exp(\xi_{i}(\phi)+\psi_{i}(\phi,a)-\bar{\psi}_{i}(\phi))}{\sum_{j}\exp(\xi_{j}(\phi)+\psi_{j}(\phi,a)-\bar{\psi}_{j}(\phi))}$
(26)

The multistep return is used with the next-best action to calculate the target
distribution used for the Kullback-Leibler divergence loss function. The
transitions are stored in a PER buffer but instead uses the divergence error
to update the priority with a smaller offset.

As with C51, we cannot find an example of authors using the complete form of
Rainbow DQN for DR. The pseudocode for Rainbow is shown in Algorithm 2.

Algorithm 2 Rainbow Deep Q-Networks

Initialise network $Z_{\theta}(s,a)$ with random weights $\theta$

Initialise target network $Z_{\hat{\theta}}(s,a)$ and set
$\hat{\theta}\leftarrow\theta$

Initialise PER buffer with initial priority $p_{0}=1$

for $t=0,T$ do

Observe state $s_{t}$

Choose action $a_{t}=\operatorname*{arg\,max}_{a}\sum_{i}z_{i}d_{i}(s_{t},a)$

Execute $a_{t}$, observe reward $r_{t}$ and next state $s_{t+1}$

Store transition $\langle s_{t},a_{t},r_{t},s_{t+1}\rangle$ with
$p_{t}=\max_{i<t}p_{i}$

Sample batch $\langle s_{i},a_{i},r_{i},s_{i+1}\rangle$ with
$P(i)=p_{i}^{\alpha}/\sum_{j}p_{j}^{\alpha}$

Calculate weights $w_{i}=(N\cdot P(i))^{-\beta}/\max_{j}w_{j}$

Calculate projected targets
$\Phi\hat{\mathcal{T}}Z_{\hat{\theta}}(s_{i},a_{i})$

Loss $\mathcal{L}_{i}(\theta)=w_{i}\cdot
D_{\text{KL}}(\Phi\hat{\mathcal{T}}Z_{\hat{\theta}}(s_{i},a_{i})\|Z_{\theta}(s_{i},a_{i}))$

Perform gradient descent on $\theta$ using loss

Update $p_{i}\leftarrow
D_{\text{KL}}(\Phi\hat{\mathcal{T}}Z_{\hat{\theta}}(s_{i},a_{i})\|Z_{\theta}(s_{i},a_{i}))$

Every $C$ steps, set $\hat{\theta}\leftarrow\theta$

end for

##  4 Case Study

This section outlines the full methodology and setup used in the case study
simulations. Figure 4 shows the process at each time-step. The prediction ANNs
provide forecasted values for the state observations which are then passed to
the agent to output the action following the learnt policy. The environment
then executes that action and the transition is stored in the replay buffer,
which is sampled for learning.

Figure 4: The control process at each step of the environment.

###  4.1 Microgrid Simulation

The simulations are run across 4 years of data from 2014 up to 2018 with
hourly readings for energy consumption, price, and weather parameters. Each
time-step represents one hour with each episode representing one week or 168
time-steps. The simulations are run over 200 episodes with the first 100 are
used for training and the second 100 for evaluation. As the environment has no
terminal states, the episodes are only relevant for reviewing performance.

Both the forecasting ANNs and algorithms are implemented using Tensorflow. For
fairness and repeatability, each ANN is set to the same Tensorflow and Numpy
random seed. For the evaluation episodes, the $\epsilon$ value of
$\epsilon$-greedy policy methods was set to 0 and the noise was removed from
agents with noisy layers. However, the agents will continue to learn during
the evaluation steps so are always learning online.

Each method is tested under four different scenarios with different sizes of
state and action-space: with and without the forecast predictions as well as
with a smaller and larger action-space. Therefore, the agents are tested on
how they effectively utilise the forecasted values as well as the resolution
of control using a set of discrete actions.

A notable challenge in this case study is the limitation of data. RL agents
are often tested using the Arcade Learning Environment [39] or similar
simulated game environments where agents are able to train indefinitely. Due
to the limited amount of data available in this case study, it is important
that the algorithms are able to train quickly and efficiently.

###  4.2 Hyperparamter Testing

Table 1: Agent Hyperparameters Parameter | Value
---|---
ANN Learning Rate | 0.001
Discount Factor $\gamma$ | 0.99
Batch Size | 64
Epsilon $\epsilon$ |  1 $\rightarrow$ 0.01
PER Priority Scale $\alpha$ | 0.7
PER Importance Scale $\beta$ | 0.5
MS-DQN Backup Steps | 2
NN-DQN Initial Noise $\sigma_{0}$ | 0.5
C51 Atoms | 51
C51 Boundaries $v_{\text{min}}$, $v_{\text{max}}$ | -10, 10
DDPG Noise | 0, 0.1

Different hyperparameters were tested in preliminary simualtions with the
final key values shown in Table 1. The PER uses a relitavely high priority
scale and importance value with a low offset, all of which promote aggressive
sampling of the more important transitions [35]. Various MS-DQN back-up sizes
were tested but performance decreased the larger the step size so Rainbow only
uses the standard 1-step backup. Gradient clipping has been suggested to
improve the performance of D3QN [33] and reward clipping for C51 [38].
However, these either had no effect or even lowed performance in those
preliminary studies as the rewards are typically small between -1 and 1
anyway.

###  4.3 Forecasting

(a) Demand real and predicted values across one week.

(b) Price real and predicted values across one week.

Figure 5: The prediction ANN forecasts for Week 200.

The scenarios with the large state-space use forecasted values calculated
using a regression ANN each for demand, price, PV generation, and WT
generation. The RES forecasts use basic ANNs provided with the previous 24
hours worth of the values as input. The demand and price prediction ANNs are
more advanced and use several weather parameters: dry bulb temperature, air
humidity, and sky clearness index. Using the amount of rainfall was also
tested but was found to increase the errors due to the high amount of zero-
data. The first 100 weeks of the data are used for training and the second 100
for evaluation so that it coincides with the training and evaluating periods
of the simulations themselves.

Several different architectures of ANN were tested for the prediction networks
including convolutional and recurrent long short-term memory layers for each
of the input parameters. However, it was found that a basic ANN had comparable
mean absolute percentage errors (MAPE) to the other architectures but often
with a significantly faster training time. The MAPE for the used demand
predictions is 5.01% and price predictions is 11.62%. An example of the values
for the demand and price forecasts for the final week given in Figure 5(a) and
5(b) respectively.

###  4.4 Benchmarks

Two benchmarks are used to assess the performance of the DQN agents: an
advanced actor-critic method and a linear programming model. Both methods have
the advantage of continuous control rather than the discrete set of actions
for DQN.

####  4.4.1 Deep Deterministic Policy Gradients (DDPG)

The first benchmark used is DDPG [18], an actor-critic algorithm with many
similarities to DQN. The key difference is as well as having a the critic
network that estimates the action-value, an actor network is used to calculate
a deterministic action with a noise process for exploration.

The main benefit of DDPG over DQN is the continuous action-space meaning it
can charge or discharge by any valid amount, whereas DQN must select from one
of the discrete actions. However, actor-critic methods typically have much
more unstable learning properties than value-function methods. An example of
its use for DR is Zhang et al. [19] as a benchmark for control in a hybrid
energy microgrid to compare against a different actor-critic method.

####  4.4.2 Linear Programming Model

The second benchmark is a linear programming model, similar to that used by
Cao et al. [16]. As there are no effective methods for optimally solving
general nonlienar programming problems [40], the model uses a simplified
environment with constant efficiencies for the ESS, transformer, and inverter.
The battery charging and discharging efficiency is generalised to
$\eta^{\text{ESS}}=0.95$ with the transformer and inverter efficiency combined
and generalised to $\eta^{\text{Grid}}=0.92$.

As the efficiencies are applied differently if the battery is charging or
discharging, the power $x_{t}$ is separated into the charging power
$x^{\text{ch}}_{t}$ and discharging power $x^{\text{dis}}_{t}$. The model is
bound by the following constraints:

$\text{Constraints}:\begin{cases}0\leqslant
x^{\text{ch}}_{t}\leqslant\text{X}^{\text{max}}\\\ 0\leqslant
x^{\text{dis}}_{t}\leqslant\text{X}^{\text{max}}\\\ 0\leqslant
c_{t}\leqslant\text{C}_{\text{max}}\\\ c_{0}=0\\\
c_{t}=c_{t-1}\eta^{\text{SDC}}+x^{\text{ch}}_{t}\eta^{\text{ESS}}-x^{\text{dis}}_{t}\end{cases}$
(27)

The model performs energy arbitrage using the predicted value of the wholesale
energy price at the next step $\text{P}^{\text{pred}}_{t+1}$. The objective of
the solver is to maximise the revenue over the 100 evaluation weeks:

$\text{Obj}=\max_{x}\sum\limits_{i=0}\text{P}^{\text{pred}}_{i}\left(x^{\text{dis}}_{i}\eta^{\text{Grid}}\eta^{\text{ESS}}-\frac{x^{\text{ch}}_{i}}{\eta^{\text{Grid}}}\right)$
(28)

The sequence of $x$ values the produced by the model are saved and input into
the complete environment for comparison against the agents.

##  5 Results and Discussion

Table 2: ESS Savings Results (1000 GBP) and percentage difference against DQN Algorithm |  |  Basic ($\mathcal{S}$=8) |  |  Forecasting ($\mathcal{S}$=12)
---|---|---|---|---
|  $\mathcal{A}$=5 | vs DQN |  $\mathcal{A}$=9 | vs DQN |  |  $\mathcal{A}$=5 | vs DQN |  $\mathcal{A}$=9 | vs DQN
DQN |  | 75.00 | - | 68.58 | - |  | 78.91 | - | 76.1 | -
DDQN |  | 71.46 | -4.72% | 70.88 | 3.35% |  | 81.02 | 2.67% | 77.97 | 2.46%
D3QN |  | 75.11 | 0.15% | 73.36 | 6.97% |  | 80.69 | 2.26% | 79.83 | 4.90%
PER-DQN |  | 76.69 | 2.25% | 72.29 | 5.41% |  | 81.55 | 3.35% | 77.76 | 2.18%
MS-DQN |  | 70.55 | -5.93% | 64.95 | -5.29% |  | 72.56 | -8.05% | 76.79 | 0.91%
NN-DQN |  | 76.47 | 1.96% | 73.26 | 6.82% |  | 79.65 | 0.94% | 75.98 | -0.16%
C51 |  | 79.25 | 5.67% | 81.56 | 18.93% |  | 85.28 | 8.07% | 83.30 | 9.46%
Rainbow |  | 79.46 | 5.95% | 80.00 | 16.65% |  | 89.33 | 13.20% | 91.49 | 20.22%
DDPG |  | 77.14 | 2.85% | 77.14 | 12.48% |  | 82.93 | 5.09% | 82.93 | 8.98%
Model |  | 79.54 | 6.05% | 79.54 | 15.98% |  | 79.54 | 0.80% | 79.54 | 4.52%
Figure 6: Cumulative energy savings as a percentage relative to DQN across the
evaluation episodes. Figure 7: Smoothed episodic rewards (over 5 episodes) of
the DRL algorithms across all episodes.

In these simulation, the DQN algorithms introduced in Section 3 are used to
reduce the energy costs of the microgrid environment outlined in Section 2.
The results for four different scenarios using different state and action-
spaces can be found in 2. The energy cost savings across the 100 evaluation
weeks relative to DQN can be found in Figure 6. The episodic rewards of the
different algorithms across all 200 weeks can be found in Figure 7.

This section will analyse the results in the context of the different
algorithms used, how they utilise the large state and action-spaces, and
finally looking into the behaviour of the distributional agents.

###  5.1 Algorithms

In the most basic simulation, DQN with the smaller action-space performs
largely the same as most of the other DQN variants but performs noticeably
worse when either the state or action-space size is increased. It also is
never able to achieve a greater saving than the linear programming model. This
suggests it is sufficient for simple environments but alternatives should be
used for more complex problems.

D3QN and PER-DQN perform better than DQN in all scenarios meaning there is
sufficient merit in using the advantage function and a PER respectively. NN-
DQN performs better than DQN in the basic scenarios but performs largely the
same in the forecasting scenarios. Therefore, an $\epsilon$-greedy policy is
sufficient for this application but NN-DQN may complement a PER as the
exploration method is tied directly to the learning. However, these methods
are only able to outperform the linear programming model in one scenario and
cannot outperform DDPG in any. MS-DQN performs worse than all the other
methods almost every time suggesting it is inappropriate for this application
because the agent cannot see the immediate outcome of its action.

The greatest improvements come from C51 and Rainbow. Both are able to
outperform DDPG in every scenario and are only bested by the model solution in
the basic environment with the smaller action-space. Rainbow is able to
achieve a saving of £91.49k in the forecasting scenario with larger action-
space, a 20.22% increase over DQN’s score. This suggests the various
subcomponents of the Rainbow algorithm all complement each other well as the
result is much greater than the model, DDPG, and C51.

The distributional approach appears to be the most significant advancement to
DQN, as was also found by Hessel et al. [23] for Atari games. For this
application, we suggest the distributional approach is better able to
interpret the punishments for exceeding the ESS’s capacity limits which would
affect the expectation for a direct action-value calculation, as well as
learning how to navigate the nonlinear efficiency profiles of the
environment’s components.

###  5.2 Action-Space and Forecasting

Figure 8: Control process at each step of Episode 200 using Rainbow DQN.

Most of the DQN variants all had reduced performance when the action-space
$\mathcal{A}$ was increased from 5 to 9, meaning the agent could not explore
enough state-action pairs as effectively when there are more actions to select
from. The two notable exceptions to this are C51 and Rainbow where their
performance either stays largely the same or even slightly increased. This is
particularly interesting as it suggests the distributional approach to value
calculation is more effective in larger action-spaces and can give finer
control of the environment than regular DQN. Despite the model and DDPG have a
continuous action-space, the superior learning properties of Rainbow and C51
proved more valuable.

All of the agents performed better when provided with future forecasted values
but the most notable improvement was in the Rainbow agent. The various
subcomponents were able to combine well and proved effective together in
interpreting and learning from the forecasted values with Rainbow performing
significantly better than all other DQN variants and the two benchmarks. An
example of the control process of the Rainbow agent for the final week of the
forecasted state-space and larger action-space can be found in Figure 8.

###  5.3 Value Distribution

(a) Episode 1, Day 2.

(b) Episode 50, Day 6.

(c) Episode 100, Day 4.

(d) Episode 170, Day 7.

(e) Episode 190, Day 2.

(f) Episode 200, Day 4.

Figure 9: Value distributions across different episodes, with the current step
marked by the black dotted line.

Figure 9 shows how the Rainbow agent learns the value distribution over time
at different steps. In each of the subfigures, the plot on the left shows the
probability of getting a particular return while the plot on the right shows
the state of the environment at that step. For context, the value distribution
of less valuable actions appear to the left with more valuable actions to the
right. Despite the figures being taken in the larger action-space scenario,
the distributions of only five actions have been shown in the interest of
legibility.

At episode 1 (9(a)) the distribution is random as the agent has not yet begun
learning. By episode 50 (9(b)), the agent has begun to understand the
environment but there is still a probability of taking a negative action; such
as charging or discharging when the ESS is full or empty respectively. By
episode 100 (9(c)), the agent has learnt that it can almost always take a
positive action in any state.

Figure 9(d) shows an example of a step where the ESS has a lot of charge
before a period of higher energy prices with the higher return of the
discharging actions showing the agent has learnt to discharge in that event.
Figure 9(e) show a case where the ESS is empty but prices are relatively high
so remaining idle is the most valuable action. Figure 9(f) shows a case at the
end of the day as the prices begin to drop where the agent has learnt its most
valuable to either remain idle or begin charging more slowly.

These findings are important as it shows the benefits of using the
distributional RL approach. Whereas most methods learn the expectation of the
action-value, Rainbow and C51 can better understand non-deterministic reward
functions and under what conditions those irregularities occur; such as the
punishment received for exceeding the ESS’s boundaries in this case study. The
distributional approach also provides a means to visualise how the agent is
learning, which is valuable for contextualising and interpreting the behaviour
of RL agents.

###  5.4 Improvements and Further Work

One improvement that could be made to this case study is more rigorous testing
of the forecasting methods with the breadth of historical weather data
available. However, the accuracy of the prediction ANNs were more than
sufficient in this case study as a supplement to the RL agents.

In the future, this case study will be built upon using a multi-agent system.
Here, different agents will be able to control multiple different types of ESS
collaboratively alongside forms of distributed generation. This may also
include competitive agents with selfish interests for other stakeholders; such
as home owners with their own smaller ESSs or for the charging of electric
vehicles.

##  6 Conclusion

In this paper, the deep RL algorithm Rainbow DQN was used to operate an ESS in
a microgrid with its own demand, renewable energy, and dynamic energy prices.
The various subcomponents of Rainbow were tested and found that the value
distribution approaches - C51 and Rainbow - outperformed the other algorithms
by a notable margin. Basic DQN was sufficient in the simpler scenarios but
struggled in the more complex cases. The components of Rainbow were able to
effectively compliment each other and allowed the agent to effectively use and
learn from the future forecasting values, outperforming the other DQN variants
and the benchmarks showing the superior learning capabilities of Rainbow. By
plotting the value distribution, it was also found to be an effective means to
contextualise and understand the agent’s behaviour.

## Acknowledgements

This work was supported by the Smart Energy Network Demonstrator project
(grant ref. 32R16P00706) funded by ERDF and BEIS. This work is also supported
by the EPSRC EnergyREV project (EP/S031863/1) and the Royal Society Research
Grant (RGS/R1/191395).

## References

  * [1] D. Bogdanov, M. Ram, A. Aghahosseini, A. Gulagi, A. S. Oyewo, M. Child, U. Caldera, K. Sadovskaia, J. Farfan, L. De Souza Noel Simas Barbosa, M. Fasihi, S. Khalili, T. Traber, C. Breyer, Low-cost renewable electricity as the key driver of the global energy transition towards sustainability, Energy 227 (2021) 120467.  doi:10.1016/j.energy.2021.120467.
  * [2] M. Ram, M. Child, A. Aghahosseini, D. Bogdanov, A. Lohrmann, C. Breyer, A comparative analysis of electricity generation costs from renewable, fossil fuel and nuclear sources in G20 countries for the period 2015-2030, Journal of Cleaner Production 199 (2018) 687–704.  doi:10.1016/j.jclepro.2018.07.159.
  * [3] Department for Business, Energy, & Industrial Strategy, Digest of UK Energy Statistics 2020, Tech. rep., The Government of the United Kingdom (Jul. 2020).
  * [4] P. Palensky, D. Dietrich, Demand Side Management: Demand Response, Intelligent Energy Systems, and Smart Loads, IEEE Transactions on Industrial Informatics 7 (3) (2011) 381–388.  doi:10.1109/TII.2011.2158841.
  * [5] S. Vazquez, S. M. Lukic, E. Galvan, L. G. Franquelo, J. M. Carrasco, Energy Storage Systems for Transport and Grid Applications, IEEE Transactions on Industrial Electronics (2010) 3881–3895.
  * [6] F. Katiraei, R. Iravani, N. Hatziargyriou, A. Dimeas, Microgrids management, IEEE Power and Energy Magazine 6 (3) (2008) 54–65.  doi:10.1109/MPE.2008.918702.
  * [7] R. S. Sutton, A. G. Barto, Reinforcement learning: an introduction, 2nd Edition, Adaptive computation and machine learning series, The MIT Press, Cambridge, MA, 2018.
  * [8] J. R. Vázquez-Canteli, Z. Nagy, Reinforcement learning for demand response: A review of algorithms and modeling techniques, Applied Energy 235 (2019) 1072–1089.
  * [9] A. Perera, P. Kamalaruban, Applications of reinforcement learning in energy systems, Renewable and Sustainable Energy Reviews 137 (2021) 110618.  doi:10.1016/j.rser.2020.110618.
  * [10] V. Mnih, K. Kavukcuoglu, D. Silver, A. A. Rusu, J. Veness, M. G. Bellemare, A. Graves, M. Riedmiller, A. K. Fidjeland, G. Ostrovski, S. Petersen, C. Beattie, A. Sadik, I. Antonoglou, H. King, D. Kumaran, D. Wierstra, S. Legg, D. Hassabis, Human-level control through deep reinforcement learning, Nature 518 (7540) (2015) 529–533.
  * [11] E. Kuznetsova, Y.-F. Li, C. Ruiz, E. Zio, G. Ault, K. Bell, Reinforcement learning for microgrid energy management, Energy 59 (2013) 133–146.  doi:10.1016/j.energy.2013.05.060.
  * [12] C. J. C. H. W. Watkins, Learning from Delayed Rewards, Ph.D. thesis, University of Cambridge, UK (1989).
  * [13] P. Kofinas, A. Dounis, G. Vouros, Fuzzy Q-Learning for multi-agent decentralized energy management in microgrids, Applied Energy 219 (2018) 53–67.  doi:10.1016/j.apenergy.2018.03.017.
  * [14] V. François-Lavet, D. Taralla, D. Ernst, R. Fonteneau, Deep Reinforcement Learning Solutions for Energy Microgrids Management, European Workshop on Reinforcement Learning (EWRL 2016).
  * [15] Y.-H. Bui, A. Hussain, H.-M. Kim, Double Deep Q-Learning-Based Distributed Operation of Battery Energy Storage System Considering Uncertainties, IEEE Transactions on Smart Grid 11 (1) (2020) 457–469.  doi:10.1109/TSG.2019.2924025.
  * [16] J. Cao, D. Harrold, Z. Fan, T. Morstyn, D. Healey, K. Li, Deep Reinforcement Learning-Based Energy Storage Arbitrage With Accurate Lithium-Ion Battery Degradation Model, IEEE Transactions on Smart Grid 11 (5) (2020) 4513–4521.  doi:10.1109/TSG.2020.2986333.
  * [17] E. Mocanu, D. C. Mocanu, P. H. Nguyen, A. Liotta, M. E. Webber, M. Gibescu, J. G. Slootweg, On-line Building Energy Optimization using Deep Reinforcement Learning, IEEE Transactions on Smart Grid (2018) 1–1doi:10.1109/TSG.2018.2834219.
  * [18] T. P. Lillicrap, J. J. Hunt, A. Pritzel, N. Heess, T. Erez, Y. Tassa, D. Silver, D. Wierstra, Continuous control with deep reinforcement learning, arXiv:1509.02971 [cs, stat].
  * [19] G. Zhang, W. Hu, D. Cao, W. Liu, R. Huang, Q. Huang, Z. Chen, F. Blaabjerg, Data-driven optimal energy management for a wind-solar-diesel-battery-reverse osmosis hybrid energy system using a deep reinforcement learning approach, Energy Conversion and Management 227 (2021) 113608.
  * [20] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, O. Klimov, Proximal Policy Optimization Algorithms, arXiv:1707.06347 [cs]ArXiv: 1707.06347.
  * [21] G. Pinto, M. S. Piscitelli, J. R. Vázquez-Canteli, Z. Nagy, A. Capozzoli, Coordinated energy management for a cluster of buildings through deep reinforcement learning, Energy 229 (2021) 120725.  doi:10.1016/j.energy.2021.120725.
  * [22] T. Haarnoja, A. Zhou, P. Abbeel, S. Levine, Off-Policy Maximum Entropy Deep Reinforcement Learning with a Stochastic Actor (2018) 10.
  * [23] M. Hessel, J. Modayil, H. van Hasselt, T. Schaul, G. Ostrovski, W. Dabney, D. Horgan, B. Piot, M. Azar, D. Silver, Rainbow: Combining Improvements in Deep Reinforcement Learning, AAAI Conference on Artificial Intelligence (2018) 3215–3222.
  * [24] Keele University, Smart Energy (2020).
URL tinyurl.com/keeleSEND

  * [25] Nord Pool, Historical Market Data (2020).
URL nordpoolgroup.com/historical-market-data/

  * [26] C. Carrillo, A. Obando Montaño, J. Cidrás, E. Díaz-Dorado, Review of power curve modelling for wind turbines, Renewable and Sustainable Energy Reviews 21 (2013) 572–581.
  * [27] T. Morstyn, B. Hredzak, R. P. Aguilera, V. G. Agelidis, Model Predictive Control for Distributed Microgrid Battery Energy Storage Systems, IEEE Transactions on Control Systems Technology 26 (3) (2018) 1107–1114.  doi:10.1109/TCST.2017.2699159.
  * [28] B. V. Mbuwir, D. Geysen, F. Spiessens, G. Deconinck, Reinforcement learning for control of flexibility providers in a residential microgrid, IET Smart Grid 3 (1) (2019) 98–107.  doi:10.1049/iet-stg.2019.0196.
  * [29] A. Driesse, P. Jain, S. Harrison, Beyond the curves: Modeling the electrical efficiency of photovoltaic inverters, in: 2008 33rd IEEE Photovoltaic Specialists Conference, 2008, pp. 1–6.  doi:10.1109/PVSC.2008.4922827.
  * [30] R. Bellman, The theory of dynamic programming, Bulletin of the American Mathematical Society 60 (6) (1954) 503–516.
  * [31] L.-J. Lin, Self-improving reactive agents based on reinforcement learning, planning and teaching, Machine Learning 8 (1992) 293–321.
  * [32] H. van Hasselt, A. Guez, D. Silver, Deep Reinforcement Learning with Double Q-Learning, AAAI Conference on Artificial Intelligence (2016) 2094–2100.
  * [33] Z. Wang, T. Schaul, M. Hessel, H. van Hasselt, M. Lanctot, N. de Freitas, Dueling Network Architectures for Deep Reinforcement Learning, AAAI Conference on Artificial Intelligence.
  * [34] D. J. B. Harrold, J. Cao, Z. Fan, Battery Control in a Smart Energy Network using Double Dueling Deep Q-Networks, in: 2020 IEEE PES Innovative Smart Grid Technologies Europe (ISGT-Europe), 2020, pp. 106–110.  doi:10.1109/ISGT-Europe47291.2020.9248785.
  * [35] T. Schaul, J. Quan, I. Antonoglou, D. Silver, Prioritized Experience Replay, arXiv:1511.05952 [cs].
  * [36] J. Yang, M. Yang, M. Wang, P. Du, Y. Yu, A deep reinforcement learning method for managing wind farm uncertainties through energy storage system control and external reserve purchasing, International Journal of Electrical Power & Energy Systems 119 (2020) 105928.  doi:10.1016/j.ijepes.2020.105928.
  * [37] R. S. Sutton, Learning to predict by the methods of temporal differences, Machine Learning 3 (1) (1988) 9–44.  doi:10.1007/BF00115009.
  * [38] M. G. Bellemare, W. Dabney, R. Munos, A Distributional Perspective on Reinforcement Learning, arXiv:1707.06887 [cs, stat].
  * [39] M. G. Bellemare, Y. Naddaf, J. Veness, M. Bowling, The Arcade Learning Environment: An Evaluation Platform for General Agents, Journal of Artificial Intelligence Research 47 (2013) 253–279.  doi:10.1613/jair.3912.
  * [40] S. P. Boyd, L. Vandenberghe, Convex optimization, Cambridge University Press, Cambridge, UK ; New York, 2004.