# Dipolar Ising Model: Phases, Growth Laws and Universality

Shikha Kumari  Sanjay Puri  School of Physical Sciences, Jawaharlal Nehru
University, New Delhi - 110067, India  Varsha Banerjee  Department of Physics,
Indian Institute of Technology Delhi, New Delhi - 110016, India

###### Abstract

The behavior of many magnetic and dielectric solids, and the more contemporary
magnetic super-lattices, is governed by dipolar interactions. They are
anisotropic and long-ranged, having varied consequences ranging from ground
states with complicated magnetic order to the presence of glassy dynamics
characterized by a plethora of relaxation times. These systems are well-
captured by the dipolar Ising model (DIM) with nearest-neighbor exchange
interactions ($J$) and long-range dipolar interactions ($D$). Depending on the
relative interaction strength $\Gamma=J/D$, there are four phases of distinct
magnetic order and symmetry. Using Monte Carlo simulations, we perform deep
quenches to study domain growth or coarsening in the $d=3$ DIM. This important
non-equilibrium phenomenon has not been addressed as dipolar interactions are
notoriously difficult to handle theoretically. Our study reveals that, in
spite of the anisotropy in interactions and diversity in ground state
configurations, we observe universality in the ordering dynamics of all
phases.

##  I Introduction

Dipolar interactions are prevalent in magnetic and dielectric solids composed
of rare-earths and transition metals Kretschmer and Binder (1979); Merz
(1954); Molnar et al. (1997); Chakraborty et al. (2004); Chamberlain and
Corruccini (2005); Biltmo and Henelius (2009); Wu et al. (2017); Cai et al.
(2019). They are anisotropic and long-ranged, and arise from nuclear magnetic
moments in alkali hydrides and solid ${}^{3}He$, electron magnetic moments in
rare earth fluorides, chlorides and hydroxides, electric dipole moments in
ferroelectric and anti-ferroelectric structures, etc. A large class of these
compounds is well-described by the nearest-neighbor (nn) Ising model with
dipolar interactions or the dipolar Ising model (DIM). For $N$ Ising spins on
a $d$-dimensional cubic lattice with sites labeled by $i$, the Hamiltonian is
given by

$\mathcal{H}=-J\sum_{\langle
ij\rangle}\sigma_{i}\sigma_{j}-D\sum_{\begin{subarray}{c}{i,j}\\\ i\neq
j\end{subarray}}\left(\frac{3\cos^{2}{\theta}_{ij}-1}{r_{ij}^{3}}\right)\sigma_{i}\sigma_{j},\quad\sigma_{i}=\pm
1.$ (1)

The first term on the right-hand-side is the contribution from the nn exchange
energy. The parameter $J$ represents the interaction strength, and favors
ferromagnetic (antiferromagnetic) alignment of spins for $J>0$ ($J<0$). The
second term is the contribution from the dipole-dipole interactions, whose
strength is given by $D=\mu^{2}/a^{3}$. Here, $\mu$ is the dipole moment of
the spin, and $a$ is the spacing between nn sites of the underlying lattice.
Further, $\vec{r}_{ij}$ is a vector joining sites $i$ and $j$ in units of
lattice spacing $a$, and $\theta_{ij}$ is the angle made by $\vec{r}_{ij}$
with the Ising axis ($z$-axis). The presence of $r_{ij}^{3}$ in the
denominator makes dipolar interactions long-ranged due to which spin-spin
interactions are significant up to multiple lattice spacings. Consequently,
the sum in the second term extends over all spin pairs. In addition to being
long-ranged, the dipolar interactions are anisotropic, fluctuating in sign and
strongly influenced by the underlying lattice structure. The
$\theta_{ij}$-dependence implies that dipolar interactions can be zero,
positive or negative – depending on the positions of the spins $i$ and $j$.
For a reference spin $i$, the interaction with spin $j$ is zero for
$3\cos^{2}{\theta}_{ij}-1=0$, antiferromagnetic for
$$55\text{\,}\mathrm{\SIUnitSymbolDegree}$<\theta_{ij}<$125\text{\,}\mathrm{\SIUnitSymbolDegree}$$,
and ferromagnetic for other values. Therefore, a ferromagnetic alignment of
spins is favored along the $z$-direction, but domain walls are preferred along
the $xy$-plane.

A phase diagram of the DIM was obtained by Kretschmer and Binder for a simple
cubic lattice $(L^{3},\ L=6,8)$ using Monte Carlo (MC) simulations Kretschmer
and Binder (1979). The interplay of the nn exchange interactions and the
complicated dipolar interactions reveals rich phase behavior. The phase
diagram in Fig. 1 has been obtained by varying $J,T$ with fixed $D=1$.
Depending on the relative interaction strength $\Gamma=J/D$, the system
exhibits four phases with distinct magnetic order and symmetry:
(I) For $\Gamma<-1.338$, the nn exchange interaction with $J<0$ dominates over
the dipolar term. As expected, the ground state (GS) is an antiferromagnet
(AFM).
(II) For $-1.338<\Gamma<0.127$, the dipolar term dominates over the nn
exchange interaction. The ground state is anisotropic and consists of
ferromagnetic columns along the $z$-axis arranged antiferromagnetically in the
$xy$-plane. We refer to this GS as a columnar antiferromagnet (CAFM).
(III) For $0.127<\Gamma<0.164$, the nn exchange favors ferromagnetic alignment
and the dipolar interaction continues to dominate. The GS structure changes to
a layered antiferromagnet (LAFM) with a two-component order parameter.
(IV) For $\Gamma>0.164$, the exchange interactions are dominant and the GS is
a ferromagnet (FM).

It is interesting to note the systematic transition from the AFM phase to the
FM phase via the onset of ferromagnetic alignment – first along a single-site
one-dimensional column (CAFM), then along a two-dimensional layer (LAFM) and
finally . The lower panel depicts the corresponding ground state (GS)
configurations for each phase. The symmetry of the Hamiltonian under
$\sigma_{i}\rightarrow-\sigma_{i}$ spins results in degenerate GS. The
degeneracy is two for Phases I, II, IV, and four for Phase III. The
equilibrium states for a given value of $\Gamma$ at non-zero $T$ will be
combinations of the corresponding degenerate GS, separated by frustrated spins
due to conflicting interactions. The free energy landscape is complex with
several local minima, especially in the dipolar-dominated regimes.

A large number of experimental realizations of the DIM indeed exhibit a non-
trivial organization of dipoles in the (so-called) equilibrium phase. For
example, magnetic measurements of rare earth hydroxides reveal FM ordering for
$\text{Tb(OH)}_{3}$, $\text{Dy(OH)}_{3}$ and $\text{Ho(OH)}_{3}$, and a
complex AFM ordering for $\text{Nd(OH)}_{3}$ and $\text{Gd(OH)}_{3}$ Wolf et
al. (1968). Amongst rare-earth fluorides, LiHoF4 and LiTbF4 exhibit FM
ordering, while LiErF4 shows AFM ordering Battison et al. (1975); Biltmo and
Henelius (2009); Holmes et al. (1973); Beauvillain et al. (1977). Similar
observations have been made in the context of rare earth ethyl sulphates Cooke
et al. (1959), rare earth perovskites Wu et al. (2017, 2019), and rare earth
garnets Cai et al. (2019). The above examples have captured much attention in
modern solid state physics and materials science because they exhibit
interesting magnetic, multiferroic, and optical effects. Some of these
systems, due to their low critical temperatures, are also realizations of
quantum Hamiltonians and exhibit phase transitions driven by quantum
fluctuations Bitko et al. (1996); Schechter and Stamp (2005). The presence of
a glassy phase due to frustrated moments arising from conflicting FM and AFM
interactions has also been seriously contemplated Piatek et al. (2013);
Andresen et al. (2014); Alonso and Allés (2017). Recent inclusions in the
family of dipolar solids are self-assembled super-lattices of mono-disperse
magnetic nanoparticles such as ferric oxide (Fe2O3) or magnetite (Fe3O4) Bae
et al. (2007); Singh et al. (2014). They exhibit a rich phase diagram
revealing a variety of stable structures such as hexagonal close-packed, face
centered cubic, spherical, cylindrical, etc Talapin et al. (2007); Boles et
al. (2016). Complex and almost perfect geometric ordering of the nano-dipoles
has also been observed in microscopy snapshots Taheri et al. (2015). Both
these observations have been interpreted in the context of dipolar
interactions between individual particles.

Laboratory experiments generally require application of external fields that
drive the system out-of-equilibrium. The system re-equilibrates, and the
approach to equilibrium critically depends on the free energy landscape. An
important non-equilibrium study in this context is the kinetics of domain
growth or phase ordering, initiated by a sudden quench of the system from the
disordered phase to the ordered phase Puri and Wadhawan (2009). The subsequent
domain growth or coarsening, characterized by a growing length scale
$\ell(t)$, is monitored with time. The growth law reveals details of the free
energy landscape and relaxation time-scales Lai et al. (1988). For example,
pure and isotropic systems with non-conserved kinetics such as the Ising model
($D=0$) obey the Lifshitz-Allen-Cahn (LAC) law: $\ell(t)\sim t^{1/2}$ Lifshitz
(1962); Allen and Cahn (1979). It is characteristic of systems with no energy
barriers to coarsening and a unique relaxation time-scale. Systems with
complicated free energy landscape due to disorder and competing interactions
have a plethora of relaxation time-scales Binder and Young (1986); Fischer and
Hertz (1991); Janke (2008); Ros et al. (2019). The interfaces are usually
rough fractals, and the barriers to coarsening grow as a power law of the
domain size Villain (1984); Shrivastav et al. (2011, 2014). Domain growth in
these systems exhibits logarithmic behavior in the asymptotic limit Corberi et
al. (2012).

In spite of its wide-ranging presence, there are only a few studies of the DIM
in $d=3$, and non-equilibrium phenomena are even less addressed Bupathy et al.
(2017, 2018). This is primarily because handling long-ranged interactions is
notoriously difficult, both analytically and computationally. Thus challenged,
we were motivated to develop theoretical techniques to study the DIM. In a
benchmarking study, we investigated coarsening via large-scale Monte Carlo
(MC) simulations on cubic lattices $(L=128)$ using Ewald summation procedures
to accurately handle the long-range dipolar interactions Bupathy et al.
(2017). We investigated their effect on growth laws in the FM phase for a
limited range of $\Gamma$ values $(\gtrsim 0.16)$. Encouraged by the unusual
observation of anisotropic growth laws in this simplest phase, we now
undertake the task of exploring the entire phase diagram. Our investigations
are guided by the following questions: Are the growth laws distinct in the
four phases or are they universal? What information can we obtain about the
energy landscapes of the four phases? Is the system characterized by a
universal scaling behavior or is it phase specific? Does the complexity and
frustration introduced by dipolar interactions yield rough fractal interfaces?

There are two important results in our paper. Our first finding is that,
although the equilibrium states in the four phases have distinct symmetry and
the interactions are anisotropic, the system obeys the LAC domain growth law
$\ell(t)\sim t^{1/2}$ across phases and directions. The growth exponent $1/2$
is universal for the DIM. The second important result is that the two-point
equal-time correlation function that quantifies domain growth exhibits
generalized scaling in all phases:
$C\left(x,y,z,t\right)=g(\rho/\ell_{\rho},z/\ell_{z})$, where $\rho$ is the
radial coordinate in the $xy$-plane. The coarsening system is thus
characterized by unique but distinct length-scales along $xy$ and $z$
directions. The scaling function in both directions is universal, and can be
approximated by the well-known Ohta-Jasnow-Kawasaki (OJK) function Ohta et al.
(1982); Puri and Wadhawan (2009). The rest of the paper is organized as
follows. In Section 2, we present detailed numerical results on coarsening in
phases 1-4 of the DIM. Section 3 provides a summary and discussion.

##  II Coarsening Studies

We now proceed to present detailed results from our studies of coarsening.
These studies use local single-spin flip moves which are computationally
expensive in contrast to the cluster algorithms often used for equilibrium
studies. Further, in systems with long-range dipolar interactions, spins
separated by multiple lattice spacings also contribute to the energy. An
explicit evaluation of the latter is only possible for small systems as the
number of computations scales as $O(N^{2})$. A successful procedure for
dealing with long-range interactions is the Ewald summation technique Frenkel
and Smit (1996). The basic idea here is to write the potential in two parts
using the identity:

$\dfrac{1}{r^{3}}=\dfrac{f(r)}{r^{3}}+\dfrac{1-f(r)}{r^{3}}.$ (2)

Here, $f(r)$ is an appropriate splitting function with the following
properties: (i) The first part is negligible beyond a certain cutoff $r_{\rm
cut}$ so that the summation up to this cutoff is a good approximation of this
contribution to the total energy. (ii) The second part is a slowly varying
function for all $r$ so that its Fourier transform can be represented by only
a few $k$-vectors with $|k|\leq k_{\rm cut}=2\pi/r_{\rm cut}$. Therefore, this
part can be efficiently evaluated in reciprocal space. There are many choices
of $f(r)$ that satisfy the above two conditions, but the usual choice is a
complementary error function:

$\text{erfc}(r)=\dfrac{2}{\sqrt{\pi}}\int_{r}^{\infty}e^{-t^{2}}dt.$ (3)

In particular, we use $f(r)=\text{erfc}(\sqrt{\eta}~{}r)$, where the Ewald
splitting parameter $\eta$ decides the relative weights of the real and
Fourier terms. By a suitable choice of $\eta$, we can optimize the computation
time for a specified error bound. An excellent discussion of this procedure is
found in Ref. Holm (2004).

Using cubic systems of up to $128^{3}$ spins, we performed deep quenches to
$T=0.5\ T_{c}(\Gamma)$ in the ordered phase. The quench locations are marked
in the phase diagram in Fig. 1. The initial state for all quenches was chosen
to be a random configuration of $\sigma_{i}=\pm 1$ corresponding to the
disordered (paramagnetic) phase. The Ewald summation technique with metallic
boundary conditions was used to compute the dipolar term of Eq. (1) Frenkel
and Smit (1996) . We chose the Ewald splitting parameter $\eta=0.032$, which
yields an error of $\delta=10^{-3}$ in the evaluation of the dipolar term. We
also performed simulations designed to yield an error of $\delta=10^{-4}$ on
smaller lattices. This did not alter the growth laws obtained with
$\delta=10^{-3}$ for larger lattices.

The system evolution was studied using spin-flip Glauber dynamics with the
standard Metropolis procedure up to 1024 MC steps (MCS). In each phase, the
acceptance rate (fraction of spins flipped in 1 MCS) decreases exponentially
with time. In Table 1, we show typical values of the acceptance rate in the 4
quenches studied in this paper.

Table 1: Typical values of the acceptance rate at different times in the four quenches studied here. $\Gamma$ | $-10.0$ | 0.0 | 0.14 | $\infty$
---|---|---|---|---
$t=10^{0}$ | 0.345 | 0.358 | 0.335 | 0.537
$t=10^{1}$ | 0.052 | 0.080 | 0.066 | 0.147
$t=10^{2}$ | 0.014 | 0.024 | 0.027 | 0.055
$t=10^{3}$ | 0.001 | 0.003 | 0.003 | 0.012

All statistical quantities have been averaged over $10$ different initial
conditions. This results in error bars for the correlation function and length
scale data which are smaller than the symbol sizes we use subsequently for
these quantities.

In Fig. 2, we present typical snapshots of coarsening morphologies in the four
phases. The top row shows the snapshots at $t=8$ MCS, and the second row at
$t=64$ MCS. During the coarsening process, the degenerate equilibrium states
compete with one another, and the corresponding domains are separated by
interfacial defects. As time evolves, the defects annihilate and the system
selects one of the ground states. The different colors in each snapshot
represent domains of one of the degenerate ground states of that phase, as
shown in the keys below. To identify the domains in phases I-III, a standard
prescription is used to define the staggered magnetization:
(i) In phase I, the morphology has AFM spin arrangement along $x$, $y$ and $z$
directions. The staggered spin variables are obtained using:

$\psi_{xyz}=(-1)^{x+y+z}\sigma_{xyz}.$ (4)

The green (blue) domains in Fig. 2(a) correspond to correlated regions of up
(down) spins in the staggered representation. The key below also shows the
corresponding GS configurations in terms of the original spin variables
$\sigma_{xyz}$.
(ii) In phase II, the morphology is AFM along the $x$ and $y$ directions, and
FM along the $z$ direction. So the spins are staggered using:

$\psi_{xyz}=(-1)^{x+y}\sigma_{xyz}.$ (5)

The green (blue) domains in Fig. 2(b) correspond to correlated regions of up
(down) spins in the staggered representation of the GS.
(iii) In phase III, the morphology comprises of FM $xz$ (or $yz$) planes with
an AFM arrangement in the $y$ ($x$) direction. In such morphologies, the
system is characterized by a two-component order parameter
$\vec{\psi}_{xyz}\equiv\left(\psi^{1}_{xyz},\psi^{2}_{xyz}\right)$ Krinsky and
Mukamel (1977); Sadiq and Binder (1984) with

$\displaystyle\psi^{1}_{xyz}$ $\displaystyle=$
$\displaystyle(-1)^{x}\frac{1}{4}\left[\sigma_{xyz}+\sigma_{x(y+1)z}-\left(\sigma_{(x+1)yz}+\sigma_{(x+1)(y+1)z}\right)\right],$
(6) $\displaystyle\psi^{2}_{xyz}$ $\displaystyle=$
$\displaystyle(-1)^{y}\frac{1}{4}\left[\sigma_{xyz}+\sigma_{(x+1)yz}-\left(\sigma_{x(y+1)z}+\sigma_{(x+1)(y+1)z}\right)\right].$
(7)

In the morphologies in Fig. 2(c), a site is assigned: green for
$\vec{\psi}_{xyz}=(1,0)$; blue for $\vec{\psi}_{xyz}=(-1,0)$; red for
$\vec{\psi}_{xyz}=(0,1)$; orange for $\vec{\psi}_{xyz}=(0,-1)$. Each of the
colored domains represents correlated regions corresponding to one of the four
GS, as shown in the key below.
(iv) No staggering is required for the morphologies in Fig. 2(d) in the FM
phase. For uniformity of notation, we assign $\psi_{xyz}=\sigma_{xyz}$.
Clearly, the domains in each phase grow with time, although the growth is
faster in the $z$ direction as compared in the $xy$ plane. However, the growth
in the $xy$ plane for Fig. 2(c) is arrested as the system gets stuck in
metastable states. The problem of metastability is widespread in phase III due
to the presence of four competing GS.

The length scales associated with the evolving domains can be evaluated from
the two-point equal-time correlation function:

$C(\vec{r},t)=\langle\psi(\vec{r}_{1},t)\psi(\vec{r}_{2},t)\rangle-\langle\psi(\vec{r}_{1},t)\rangle\
\langle\psi(\vec{r}_{2},t)\rangle.$ (8)

Here, $\psi$ is the appropriate order parameter,
$\vec{r}=\vec{r}_{1}-\vec{r}_{2}$, and $\langle\cdot\cdot\cdot\rangle$
represents an average over independent runs Puri and Wadhawan (2009). In the
case of isotropic domain growth, the dynamical scaling ansatz assumes the
presence of a single length scale $\ell(t)$. This is demonstrated post facto
by the scaling collapse of the correlation function: $C(\vec{r},t)=g(r/\ell)$.
The validity of this has been shown in many experiments and simulations. In
the anisotropic case considered here, we propose (following Ref. Puri et al.
(1992, 1994)) the simplest anisotropic generalization of this ansatz. This
also has to be verified post facto by the dynamical scaling of the correlation
function. Thus, we introduce $C(\vec{r},t)\equiv C(\vec{\rho},z;t)$ where
$\vec{\rho}=(x,y)$. In the case of unique length scales $\ell_{\rho}$ and
$\ell_{z}$ characterizing domain growth in the $xy$ and $z$-directions, the
correlation function should exhibit generalized dynamical scaling:
$C(\vec{\rho},z,t)=g(\rho/\ell_{\rho},z/\ell_{z})$.

In the isotropic case, an approximate analytical form of the correlation
function for a system described by a scalar nonconserved order parameter has
been obtained by Ohta et al. (OJK) by studying the interfacial defect dynamics
Ohta et al. (1982); Puri and Wadhawan (2009). The OJK function is given by

$C(r,t)=\frac{2}{\pi}\sin^{-1}\gamma,\quad\gamma=\exp\left(-r^{2}/\ell^{2}\right).$
(9)

In Fig. 3, we show the scaled correlation functions for phase II with
$\Gamma=0$: $C(\vec{\rho},0,t)$ vs. $\rho/\ell_{\rho}$ in Fig. 3(a), and
$C(0,z,t)$ vs. $z/\ell_{z}$ in Fig. 3(b). The length scales are self-
consistently defined from the data as the distance over which the relevant
correlation function decays to 0.2 times its maximum value. We emphasize that
no fitting parameters have been used to observe the data collapse. Both data
sets exhibit dynamical scaling, indicating that the system is characterized by
unique but distinct length-scales along $xy$ and $z$-directions. More striking
is their agreement with the OJK function (solid line), which shows that the
defect dynamics is robust across directions in spite of the inherent
anisotropy of the dipolar interactions. Fig. 3(c) shows $C(\vec{\rho},0,t)$
vs. $\rho/\ell_{\rho}$ for $t=64$ in phases I-IV for the specified values of
$\Gamma$. Fig. 3(d) shows the corresponding data for $C(0,z,t)$ vs.
$z/\ell_{z}$ for $t=64$. The data collapse is excellent in both Figs.
3(c)-(d), and is well-approximated by the OJK form. These observations reveal
that the coarsening morphologies are robust across diverse phases of the DIM.

Finally, we discuss the growth laws in the four phases of the DIM. In this
context, it is useful to reiterate observations from earlier studies Bray and
Rutenberg (1994) of domain growth in the Ising model with isotropic long-range
interactions:

$J(r_{ij})\sim r_{ij}^{-(d+\mu)}.$ (10)

Using energy scaling arguments, Bray-Rutenberg Bray and Rutenberg (1994)
predicted the following domain growth laws:

$\ell(t)\sim\left\\{\begin{array}[]{ll}t^{1/1+\mu},&\mu<1,\\\ (t\ln
t)^{1/2},&\mu=1,\\\ t^{1/2},&\mu>1.\\\ \end{array}\right.$ (11)

A recent work by Christiansen et al. Christiansen et al. (2019) has studied
coarsening in the $d=2$ Ising model with additional long-range interactions as
in Eq. (10). Using efficient numerical schemes, the authors confirmed the
predictions in Eq. (11) for many values of $\mu$.

The angular dependence of the dipolar interactions yields anisotropic
morphologies due to ferromagnetic interactions along the $z$-direction and
antiferromagnetic interactions in the $xy$-plane. Roughly speaking, the system
behaves as a long-range Ising model along the $z$-axis with $d=1,\mu=2$. Then,
the expected growth law from Eq. (11) is $\ell_{z}(t)\sim t^{1/2}$. On the
other hand, in the $xy$-plane, we have $d=2,\mu=1$ in Eq. (10). Again, we
expect the growth to be predominantly diffusive with logarithmic corrections.

In Fig. 4, we present the growth laws $\ell_{\rho}(t)$ vs. $t$ and
$\ell_{z}(t)$ vs. $t$ on a log-log scale for phases I-IV. We have not shown
$\ell_{\rho}$ vs. $t$ for the LAFM phase as the growth in the $xy$ plane is
almost frozen in that phase Kretschmer and Binder (1979). Further, as the
growth in the $z$-direction is very rapid, we see finite-size effects in the
late-time data for the FM phase. The error bars for each data point are
smaller than the symbol sizes used except in the regime where finite-size
effects are seen. The dashed lines with the expected slope 1/2 have been
plotted for reference. It is striking that the growth in both directions, in
spite of the diversity of phases and the inherent anisotropy of the dipolar
interactions, obeys the universal LAC law $\ell(t)\sim t^{1/2}$! Their
presence is thus confined to the prefactors alone. There are no length-scale-
dependent barriers to coarsening in the DIM, and the non-equilibrium evolution
is characterized by a unique relaxation time-scale. (In our earlier study, the
domain growth in the FM phase is studied for smaller values of $\Gamma\gtrsim
0.16$ that are close to the phase boundary separating the LAFM and FM phase
Bupathy et al. (2017). Here, longer timescales are required to observe the
$t^{1/2}$ law and the anisotropy is again contained in the pre-factors.)

Before concluding, it is relevant to ask how the finite size of the system
would affect the growth laws. In contrast to the equilibrium case, various
physical quantities do not change systematically with the system size. Rather,
the domains grow in a power-law manner until the domain scale becomes some
significant fraction of the lateral system size. After that, finite-size
effects are seen via a crossover to a regime with flattening and saturation of
the domain growth law. Typically, the data sets for different system sizes are
numerically coincident until they encounter finite-size effects. This is shown
in Fig. 4(b) for the case with $\Gamma=0$.

##  III Summary and Discussion

We end this paper with a summary of our results, their implications, and
future directions. We have performed large-scale Monte Carlo (MC) simulations
to study coarsening dynamics in the dipolar Ising model (DIM) which
encompasses short-range exchange interactions as well as long-range dipolar
interactions. This model is characterized by four distinct phases: (i)
antiferromagnet (AFM), (ii) columnar antiferromagnet (CAFM), (iii) layered
antiferromagnet (LAFM), (iv) ferromagnet (FM). The dipolar interactions lead
to diverse ground state (GS) configurations with strong anisotropy. Yet the
non-equilibrium dynamics is characterized by universality. Our main
observations are (i) the spatial correlation function exhibits universal
scaling; (ii) The domain growth law obeys the universal Lifshitz-Allen-Cahn
law $\ell(t)\sim t^{1/2}$.

Thermal quenches are a starting point for many non-equilibrium studies in the
laboratory. In these experiments, the system often accesses long-lived
metastable states that are encountered in our simulations. So our observations
could provide a fresh outlook to interpret relaxation phenomena in dipolar
solids Chamberlain and Corruccini (2005); Biltmo and Henelius (2009); Wu et
al. (2019). Our work also suggests novel experimental investigations in the
more contemporary self-assembled super-lattices Ku et al. (2011). In the
latter system, the constituent magnetic nanoparticles are usually
functionalized with an insulating surfactant layer to prevent aggregation. The
surfactant thickness can be adjusted to manipulate dipole-dipole interactions
to tailor spin morphologies dictated by applications. For example, a large
class of self-assembled lattices finds applications for spintronic devices
that require an AFM arrangement of the nano-dipoles for efficient operation
Black et al. (2000); Kechrakos and Trohidou (2005); Dugay et al. (2011);
Usmani et al. (2018); Bupathy et al. (2019); Déjardin (2011); Anand et al.
(2019). An improved understanding of the interplay of short-range and long-
range interactions would yield better strategies to achieve such challenges.

To the best of our knowledge, this is the first study of coarsening in ($d=3$)
dipolar solids. The developed methodologies provide a basis for relaxation
studies in systems with anisotropic and long-range interactions in general. We
have identified unexpected dynamical universalities in the DIM, which is
representative of a large class of microscopic and mesoscopic systems.
However, there are many puzzles that still need to be understood. For example
the basis for universality, the role of lattice geometry and other forms of
long-range interactions, and the possibility of the glassy state in the LAFM
phase are open questions. Another non-equilibrium phenomenon of great
relevance in experimental systems is that of aging Binder and Young (1986);
Fischer and Hertz (1991). This property probes the history-dependent evolution
of correlation and response functions when the system is driven out-of-
equilibrium by, e.g., a thermal quench or the application of a magnetic field.
There are a few such studies for the Ising model with long-range interactions
in $d=2$, and these have suggested novel violations, phases and exponents
Stariolo and Cannas (1999); Cannas et al. (2008); Christiansen et al. (2020).
It will be interesting to study aging phenomena in the anisotropic DIM. The
consequence of special directions on the relationship between waiting times
and relaxation times is an open question. We hope that our study will motivate
investigations to seek these answers, which are important for fundamental
understanding as well as experimental interpretations.

### Acknowledgements

SK acknowledges the HPC facility at the Jawaharlal Nehru University, New
Delhi, and useful discussions with Arunkumar Bupathy. VB acknowledges DST,
India for a MATRICS research grant.

## References

  * Kretschmer and Binder (1979) R. Kretschmer and K. Binder, Z. Phys. B 34, 375 (1979).
  * Merz (1954) W. J. Merz, Phys. Rev. 95, 690 (1954).
  * Molnar et al. (1997) A. A. Molnar, Y. M. Vysochanskii, A. A. Horvat, and Y. S. Nakonechnii, Ferroelectrics 192, 137 (1997).
  * Chakraborty et al. (2004) P. B. Chakraborty, P. Henelius, H. Kjønsberg, A. W. Sandvik, and S. M. Girvin, Phys. Rev. B 70, 144411 (2004).
  * Chamberlain and Corruccini (2005) S. L. Chamberlain and L. R. Corruccini, Phys. Rev. B 71, 024434 (2005).
  * Biltmo and Henelius (2009) A. Biltmo and P. Henelius, Europhys. Lett. 87, 27007 (2009).
  * Wu et al. (2017) L. S. Wu, S. E. Nikitin, M. Frontzek, A. I. Kolesnikov, G. Ehlers, M. D. Lumsden, K. A. Shaykhutdinov, E.-J. Guo, A. T. Savici, Z. Gai, et al., Phys. Rev. B 96, 144407 (2017).
  * Cai et al. (2019) Y. Cai, M. N. Wilson, J. Beare, C. Lygouras, G. Thomas, D. R. Yahne, K. Ross, K. M. Taddei, G. Sala, H. A. Dabkowska, et al., Phys. Rev. B 100, 184415 (2019).
  * Wolf et al. (1968) W. P. Wolf, H. Meissner, and C. A. Catanese, J. Appl. Phys. 39, 1134 (1968).
  * Battison et al. (1975) J. E. Battison, A. Kasten, M. J. M. Leask, J. B. Lowry, and B. M. Wanklyn, J. Phys. C 8, 4089 (1975).
  * Holmes et al. (1973) L. M. Holmes, T. Johansson, and H. J. Guggenheim, Solid State Commun. 12, 993 (1973).
  * Beauvillain et al. (1977) P. Beauvillain, J. P. Renard, and P. E. Hansen, J. Phys. C 10, L709 (1977).
  * Cooke et al. (1959) A. H. Cooke, D. T. Edmonds, C. B. P. Finn, and W. P. Wolf, Proc. Phys. Soc. 74, 791 (1959).
  * Wu et al. (2019) L. S. Wu, S. E. Nikitin, M. Brando, L. Vasylechko, G. Ehlers, M. Frontzek, A. T. Savici, G. Sala, A. D. Christianson, M. D. Lumsden, et al., Phys. Rev. B 99, 195117 (2019).
  * Bitko et al. (1996) D. Bitko, T. F. Rosenbaum, and G. Aeppli, Phys. Rev. Lett. 77, 940 (1996).
  * Schechter and Stamp (2005) M. Schechter and P. C. E. Stamp, Phys. Rev. Lett. 95, 267208 (2005).
  * Piatek et al. (2013) J. O. Piatek, B. Dalla Piazza, N. Nikseresht, N. Tsyrulin, I. Živković, K. W. Krämer, M. Laver, K. Prokes, S. Mataš, N. B. Christensen, et al., Phys. Rev. B 88, 014408 (2013).
  * Andresen et al. (2014) J. C. Andresen, H. G. Katzgraber, V. Oganesyan, and M. Schechter, Phys. Rev. X 4, 041016 (2014).
  * Alonso and Allés (2017) J. J. Alonso and B. Allés, J Phys. Condens. Matter 29, 355802 (2017).
  * Bae et al. (2007) C. J. Bae, S. Angappane, J.-G. Park, Y. Lee, J. Lee, K. An, and T. Hyeon, Appl. Phys. Lett. 91, 102502 (2007).
  * Singh et al. (2014) G. Singh, H. Chan, A. Baskin, E. Gelman, N. Repnin, P. Král, and R. Klajn, Science 345, 1149 (2014).
  * Talapin et al. (2007) D. V. Talapin, E. V. Shevchenko, C. B. Murray, A. V. Titov, and P. Král, Nano Lett. 7, 1213 (2007).
  * Boles et al. (2016) M. A. Boles, M. Engel, and D. V. Talapin, Chem. Rev. 116, 11220 (2016).
  * Taheri et al. (2015) S. M. Taheri, M. Michaelis, T. Friedrich, B. Förster, M. Drechsler, F. M. Römer, P. Bösecke, T. Narayanan, B. Weber, I. Rehberg, et al., PNAS 112, 14484 (2015).
  * Puri and Wadhawan (2009) S. Puri and V. Wadhawan, _Kinetics of Phase Transitions_ (CRC Press, Boca Raton, 2009).
  * Lai et al. (1988) Z. W. Lai, G. F. Mazenko, and O. T. Valls, Phys. Rev. B 37, 9481 (1988), URL https://link.aps.org/doi/10.1103/PhysRevB.37.9481.
  * Lifshitz (1962) I. M. Lifshitz, Sov. Phys. JETP 15, 939 (1962).
  * Allen and Cahn (1979) S. M. Allen and J. W. Cahn, Acta Metall. 27, 1085 (1979).
  * Binder and Young (1986) K. Binder and A. P. Young, Rev. Mod. Phys. 58, 801 (1986), URL https://link.aps.org/doi/10.1103/RevModPhys.58.801.
  * Fischer and Hertz (1991) K. H. Fischer and J. A. Hertz, _Spin Glasses_ , Cambridge Studies in Magnetism (Cambridge University Press, 1991).
  * Janke (2008) W. Janke, _Rugged Free-Energy Landscapes - An Introduction: In Lecture Notes in Physics_ (Springer Heidelberg, 2008).
  * Ros et al. (2019) V. Ros, G. Ben Arous, G. Biroli, and C. Cammarota, Phys. Rev. X 9, 011003 (2019).
  * Villain (1984) J. Villain, Phys. Rev. Lett. 52, 1543 (1984).
  * Shrivastav et al. (2011) G. P. Shrivastav, S. Krishnamoorthy, V. Banerjee, and S. Puri, Europhys. Lett. 96, 36003 (2011).
  * Shrivastav et al. (2014) G. P. Shrivastav, M. Kumar, V. Banerjee, and S. Puri, Phys. Rev. E 90, 032140 (2014).
  * Corberi et al. (2012) F. Corberi, E. Lippiello, A. Mukherjee, S. Puri, and M. Zannetti, Phys. Rev. E 85, 021141 (2012).
  * Bupathy et al. (2017) A. Bupathy, V. Banerjee, and S. Puri, Phys. Rev. E 95, 060103 (2017).
  * Bupathy et al. (2018) A. Bupathy, V. Banerjee, and S. Puri, Europhys. Lett. 122, 36002 (2018).
  * Ohta et al. (1982) T. Ohta, D. Jasnow, and K. Kawasaki, Phys. Rev. Lett. 49, 1223 (1982).
  * Frenkel and Smit (1996) D. Frenkel and B. Smit, _Understanding Molecular Simulations: From Algorithms to Applications_ (Academic, San Diego, 1996).
  * Holm (2004) C. Holm, Computational Soft Matter: From Synthetic Polymers to Proteins 23, 195 (2004).
  * Krinsky and Mukamel (1977) S. Krinsky and D. Mukamel, Phys. Rev. B 16, 2313 (1977).
  * Sadiq and Binder (1984) A. Sadiq and K. Binder, J. Stat. Phys. 35, 517 (1984).
  * Puri et al. (1992) S. Puri, K. Binder, and S. Dattagupta, Phys. Rev. B 46, 98 (1992).
  * Puri et al. (1994) S. Puri, N. Parekh, and S. Dattagupta, J. Stat. Phys. 77, 839 (1994).
  * Bray and Rutenberg (1994) A. J. Bray and A. D. Rutenberg, Phys. Rev. E 49, R27 (1994).
  * Christiansen et al. (2019) H. Christiansen, S. Majumder, and W. Janke, Phys. Rev. E 99, 011301 (2019), URL https://link.aps.org/doi/10.1103/PhysRevE.99.011301.
  * Ku et al. (2011) J. Ku, D. M. Aruguete, A. P. Alivisatos, and P. L. Geissler, J. Am. Chem. Soc. 133, 838 (2011).
  * Black et al. (2000) C. T. Black, C. B. Murray, R. L. Sandstrom, and S. Sun, Science 290, 1131 (2000).
  * Kechrakos and Trohidou (2005) D. Kechrakos and K. N. Trohidou, Phys. Rev. B 71, 054416 (2005).
  * Dugay et al. (2011) J. Dugay, R. P. Tan, A. Meffre, T. Blon, L.-M. Lacroix, J. Carrey, P. F. Fazzini, S. Lachaize, B. Chaudret, and M. Respaud, Nano lett. 11, 5128 (2011).
  * Usmani et al. (2018) S. Usmani, M. Lepesant, A. Bupathy, T. Blon, L.-M. Lacroix, V. Banerjee, B. Chaudret, and J. Carrey, Phys. Rev. B 98, 104433 (2018).
  * Bupathy et al. (2019) A. Bupathy, V. Banerjee, and J. Carrey, Phys. Rev. B 100, 064420 (2019).
  * Déjardin (2011) P. M. Déjardin, J. Appl. Phys. 110, 113921 (2011).
  * Anand et al. (2019) M. Anand, V. Banerjee, and J. Carrey, Phys. Rev. B 99, 024402 (2019).
  * Stariolo and Cannas (1999) D. A. Stariolo and S. A. Cannas, Phys. Rev. B 60, 3013 (1999), URL https://link.aps.org/doi/10.1103/PhysRevB.60.3013.
  * Cannas et al. (2008) S. A. Cannas, M. F. Michelon, D. A. Stariolo, and F. A. Tamarit, Phys. Rev. E 78, 051602 (2008), URL https://link.aps.org/doi/10.1103/PhysRevE.78.051602.
  * Christiansen et al. (2020) H. Christiansen, S. Majumder, M. Henkel, and W. Janke, Phys. Rev. Lett. 125, 180601 (2020), URL https://link.aps.org/doi/10.1103/PhysRevLett.125.180601.

Figure 1: Schematic phase diagram of the DIM model in the
$(\Gamma,k_{B}T/D)$-plane, where $\Gamma=J/D$. We set $D=1$ and vary $J,T$.
There are four distinct phases: (I) Antiferromagnet (AFM); (II) Columnar
antiferromagnet (CAFM); (III) Layered antiferromagnet (LAFM); (IV) Ferromagnet
(FM). The corresponding ground state configurations are shown below. The
quench locations are marked by crosses, and satisfy $T=0.5T_{c}(\Gamma)$. The
corresponding parameter values are $(\Gamma,k_{B}T/D)=(-10,1.75)$ for AFM;
$(0,1.19)$ for CAFM; $(0.14,0.825)$ for LAFM; and $(\infty,\infty)$ for FM.
Figure 2: Evolution snapshots of the staggered magnetization at $t=8$ MCS (top
row) and $t=64$ MCS (second row) for the four phases: (a) AFM, (b) CAFM, (c)
LAFM, and (d) FM. The key shows the color code for the degenerate GS in each
phase, as described in the text. Figure 3: Scaled correlation functions at
$t=16,64,128$ for $\Gamma=0$ in phase II: (a) $C(\rho,0,t)$ vs.
$\rho/\ell_{\rho}$, and (b) $C(0,z,t)$ vs. $z/\ell_{z}$. Scaled correlation
functions at $t=64$ for specified $\Gamma$-values: (c) $C(\rho,0,t)$ vs.
$\rho/\ell_{\rho}$, and (d) $C(0,z,t)$ vs. $z/\ell_{z}$. The solid line in
(a)-(d) represents the OJK function. Figure 4: Characteristic length scales,
$\ell_{\rho}(t)$ vs. $t$ and $\ell_{z}(t)$ vs. $t$ on a log-log scale. The
panels correspond to (a) $\Gamma=-10$; (b) $\Gamma=0$; (c) $\Gamma=0.14$; (d)
$\Gamma=\infty$. In (b), we also include data for systems of size $32^{3}$ and
$64^{3}$ to show the finite-size effects. In (c), we have not shown
$\ell_{\rho}$ vs. $t$ as there is almost no growth in that direction. In (d),
the data sets for $\ell_{\rho}$ and $\ell_{z}$ are numerically
indistinguishable as the growth is isotropic. The dashed line in each panel
denotes the LAC law: $\ell(t)\sim t^{1/2}$.


# Co-designing Intelligent Control of Building HVACs and Microgrids

Rumia Masburah∗,1, Sayan Sinha∗,1, Rajib Lochan Jana∗,1, Soumyajit Dey1, and
Qi Zhu2
1Indian Institute of Technology Kharagpur, 2Northwestern University
{rumiamasburah, sayan.sinha<EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS>∗ denotes equal contribution. Ack: The work was funded
by MHRD and Department of Power, Govt. of India, under the IMPRINT project no.
6158.

###### Abstract

Building loads consume roughly 40% of the energy produced in developed
countries, a significant part of which is invested towards building
temperature-control infrastructure. Therein, renewable resource-based
microgrids offer a greener and cheaper alternative. This communication
explores the possible co-design of microgrid power dispatch and building HVAC
(heating, ventilation and air conditioning system) actuations with the
objective of effective temperature control under minimised operating cost. For
this, we attempt control designs with various levels of abstractions based on
information available about microgrid and HVAC system models using the Deep
Reinforcement Learning (DRL) technique. We provide control architectures that
consider model information ranging from completely determined system models to
systems with fully unknown parameter settings and illustrate the advantages of
DRL for the design prescriptions.

##  I Introduction

Microgrids with renewable power sources have emerged as a viable alternative
to the growing energy demand (about 40%) [1] of building loads in developed
countries. In large smart building loads, the heating, ventilation, and air
conditioning (HVAC) systems account for the highest percentage ($70\%$) of the
overall building load [24]. Given this huge energy requirement, an important
aim of Building Energy Management Systems (BEMS) is to optimise the HVAC
energy consumption, while maintaining the comfort needs of building occupants.
This has an immediate impact on energy cost expense reduction for the
building.

Recent research in the area of electrical cost minimisation for a network of
smart buildings can be broadly divided into two major classes. One direction
of research focuses on optimal control of building microgrids to decrease
dependency on the utility grid and balance energy supply and demand by
utilising storage systems efficiently [15, 7]. The other direction focuses on
the optimal control for building HVAC systems and load management by the BEMS
to optimise the building energy consumption for reducing electricity costs [5,
13]. Thus one may observe that while intelligent control of microgrid systems
reduces the energy dependency for building networks (i.e. a collection of
community building loads) on the utility grid, the BMS tries to optimise in-
building actuators for reducing building power-demands. These two avenues of
power dispatch, control and system operations can be made to operate more
efficiently if cross-layer, coordinated optimisations are enabled among them.
This idea forms the premise of our work. The operational cost optimisation of
a building network depends on the supply from the microgrid, demand of BEMS
and real-time price of the utility grid, and different operational constraints
imposed by the microgrid and BEMS. So the control actions for the two systems,
i.e. microgrid and BEMS, may be intertwined for reducing the total electrical
cost. It leads to co-designing the control of the two systems (microgrid and
the BEMS).

Model predictive control (MPC) shows its efficiency for optimal control design
with underlying complex constrained multivariable optimisation problems in
many domains [5]. The MPC-based strategies require accurate prediction on the
power consumption of building loads like HVACs, which depend on the accuracy
of building thermal models [12]. Developing an accurate building thermal model
from scratch is a time-consuming, and costly process [5]. It is often
intractable to get accurate dynamics of the temperature inside the building,
which is affected by many factors. Moreover, thermal models are building-
specific due to heterogeneity present in building design, the material used
and operations. For this reason, an alternative control method is to use
learning-based algorithms that directly learn the control strategies from
sensor data. However, learning-based techniques require a significantly large,
effective and diverse dataset depending on building size and internal
structure. In the case of building models, such datasets are not readily
available [11].

In order to address these issues, Reinforcement Learning (RL) [22] has been
used as a possible real-time control technique in building scenarios.
Leveraging Deep Reinforcement learning (DRL) for continuous large datasets
[10] generates near-optimal actions without having a supervisor, via trial-
and-error paradigm, for handling a complex control problem. DRL based HVAC
control works [26] mainly focus on calculating energy schedules in BEMS using
predicted day-ahead data rather than real-time experience and information. In
retrospect, the present work is the first to apply DRL technique for co-
optimising the control actions for microgrid and HVAC systems with the
objective of minimising the electrical expenses of building networks while
maintaining the comfort needs of the building occupants at different thermal
zones.

However, such co-optimisation of two structurally independent but functionally
related systems involves a very large state-space, making parameter learning
and real-time control challenging. Based on this observation, the present work
analyses three different possible state information and control scenarios and
performs the co-optimisation. Scenario-1: Parameters of grid components and
the building thermal models are all available. In this case, we employ an
integrated MPC based coordinated control of the grid and building loads.
Scenario-2: Parameters of the grid components are not known, but the building
thermal model is available with an MPC based temperature control. In this
case, we design a DRL based microgrid controller, which additionally
supervises the parameters of the HVAC MPC.
Scenario-3: All the parameters of the grid, as well as the building thermal
models, are unavailable. So, we design a pure DRL based approach for both grid
and HVAC control surfaces. We formulate the underlying sequential decision
making problem as an MDP (Markov Decision process) with suitable environment
states, actions and reward function. The main contributions of our work are
summarised as follows:

  * •

We propose a framework for co-designing the control operations of a microgrid
system with building loads. The control operations include the energy storage
system charging and discharging, power exchange between the microgrid and the
utility grid, real-time price of utility grid, comfortable temperature range
for the building, HVAC power adjustment and weather uncertainties.

  * •

We solve the control problem assuming three different abstraction levels as
discussed earlier. In that way, we provide a systematic approach for the
application of DRL techniques in this problem, assuming different levels of
system parameter information. We compare the performance of the proposed DRL-
based algorithms ( Scenario-2 and 3) with the pure model-based control
approach (Scenario-1).

##  II System Description and Problem Statement

The system model for this work is depicted in Fig. 1. It comprises a set of
building electrical loads and a microgrid (coupled with the utility grid) that
supplies the electricity to the building loads. Brief description of the
different components of our studied system are illustrated in this section.

Figure 1: Control Framework

###  II-A Microgrid System

The microgrid consists of PV panels as uncontrollable generation units, DC
diesel generators as controllable generation units and batteries as the
electrical storage system (ESS). The units are connected to the microgrid
through controllable voltage source converters (VSC). We adopt the standard
method for VSC control with inner loop current controller and outer loop
voltage controller, which regulates the VSC output voltage to a desired
reference [20].

Let us assume, the sampling interval of the microgrid is $t_{s}$. At some
instance $k$, the PV power generation is governed by the following dynamical
equation [25],

$P_{solar}(k)=A*y*r(k)*R$ (1)

where, we denote $A$ as the total PV panel area, $y$ as the solar panel yield,
$r$ as the solar radiation and $R$ as the performance ratio. For the storage
system i.e. the battery model, the charging/discharging dynamics of the
battery is given as,

$\begin{split}SOC(k+1)=(1-\delta_{ESS}t_{s})SOC(k)-u_{ESS}(k)\frac{P_{ESS}(k)}{\eta_{dis}{P^{max}_{ESS}}}\\\
+(1-u_{ESS}(k))\frac{P_{ESS}(k)}{\eta_{ch}{P^{max}_{ESS}}}\end{split}$ (2)

following [25], where $SOC(k),P_{ESS}(k)$ are the current state of charge
(SOC) of the battery and the power supply from/to the battery during
discharge/charge respectively. Also $\delta_{ESS}$ is the self-discharge rate
of storage. The symbol $\eta_{dis}$/$\eta_{ch}$ represents
discharging/charging efficiency and $P^{max}_{ESS}$ is maximum storage
capacity of the battery. The $u_{ESS}$ is the charging/discharging command,
$u_{ESS}(k)=1$ indicates the discharging mode and $u_{ESS}(k)=0$ indicates the
charging mode. Power output dynamics ($P_{DG}$) of the diesel generator can be
given as [23],

${P}_{DG}(k+1)=-\frac{1}{\tau_{DG}}P_{DG}(k)+\frac{u_{DG}(k)}{\tau_{DG}}$ (3)

where, $\tau_{DG}$ is the average delay between the delivery and power command
and $u_{DG}$ is the power commanded to the generator by the controller.

###  II-B HVAC system & Control Framework

In this work, we consider each building of the building network to be equipped
with an HVAC system. We use the hierarchical economic Model Predictive Control
(eMPC) framework from [12] for HVAC systems to satisfy the comfort need of the
in-building environment while minimising the energy consumption. The higher-
level eMPC of HVAC calculates the optimal airflow rate into each zone. The
input to the building for temperature control is, in fact, the air mass flow
rate that should enter each currently used thermal zones through the ducts.
These optimal input values are given to the lower-level controllers (PID
controller) as the setpoints for air mass flow. The output of the lower-level
controllers is basically control signals that adjust the angle of dampers
which regulate the amount of air blown into the thermal zones. In this work,
we are interested in the electrical load of the building that affects the
control of microgrid power dispatch. We assume that other electrical loads
except HVAC are constant. The power requirement of the HVAC system is mainly
governed by the power requirement of the variable speed driven fan and
chiller. The power demand of the fan can be modelled as follows [12]:

$P_{f}(t)=k_{f}(\dot{m_{i}})^{2}$ (4)

where $k_{f}$ is a parameter that captures both the fan efficiency and the
duct pressure losses while $\dot{m_{i}}$ is the supply air mass flow rate of
the fan. The chiller power demand can be modelled as

$P_{c}(t)=\dfrac{c_{p}}{\eta}\dot{m_{i}}(T_{s}-T_{c})$ (5)

following [12] where, $\eta$ is coefficient of performance for chiller,
$c_{p}$ is specific heat capacity of air, $\dot{m_{i}}$ is the supply air mass
flow rate, $T_{s}$ is the temperature of air entering the chiller and $T_{c}$
is the temperature of air going out from the chiller.

Let us assume, prediction horizon and sampling period are $H^{p}$ and $t_{s}$,
respectively. In order to calculate the optimal airflow rate, the higher-level
eMPC in [12] solves an optimisation problem at each sampling time $k$. The
objective of the problem is to trade off the energy consumption and comfort of
the building that is captured in the following objective function.

$w_{E}E_{H}(k)+w_{c}\sum_{i=0}^{n}CR_{i}*\dfrac{1}{CF_{i}(k)}$ (6)

where, $E_{H}(k)=\int_{kt_{s}}^{(k+H_{p})t_{s}}[P_{f}(t)+P_{c}(t)]dt$ is the
energy consumption of the HVAC used in [12] within the time interval
$[kt_{s},(k+H_{p})t_{s}]$. The term
$CF_{i}(k)=\sum_{j=k}^{k+H_{p}}\dfrac{1}{|T^{desire}_{i}-T_{i}(j)|}$ is the
comfort factor of $i^{th}$ zone ($1\leq i\leq n$) within the time interval
$[kt_{s},(k+H_{p})t_{s}]$, where $T^{desire}_{i}$ is the desired temperature
(set point) and $T_{i}(j)$ is the temperature of the $i^{th}$ zone at $j$-th
sampling instant.The comfort factor implies the deviation of current
temperature from the set point and the higher value of it indicates that the
zone temperature is closer to the desired value. The term $CR_{i}\in(0,1)$
indicates the air conditioning criticality of the $i^{th}$ zone. The
criticality of zones sets the relative priority for maintaining high comfort
among the zones. The $w_{E},w_{c}\in(0,1)$ are the weight parameters for
energy consumption and comfort need respectively. The energy consumption of
the HVAC system depends on the choice of system parameters $w_{E}$ and
$w_{c}$. These two parameters are needed to be suitably tuned for deciding the
trade-off between the relative priorities of comfort and energy consumption.

##  III Proposed HVAC-aware MEMS Framework

In this section, we describe our proposed control framework, depicted in Fig.
1 for the microgrid energy management system (MEMS). We use a two-level
hierarchical control design, similar to [18], for optimally managing the
microgrid operation to minimise the electrical expenses of a building network.
At each sampling time $k$, our high-level supervisory controller senses the
state $x_{m}(k)$ of the microgrid and generates the optimal schedule for the
ESS and diesel generators (DGs) of microgrid along with a suitable choice of
HVAC system input $\mathcal{H}(k)$. The schedule for the storage,
$\mathcal{S}_{ESS}(k)$, indicates the amount of power to be fed to or drawn
from the microgrid for a specific time duration $[k,k+1]$ during discharging
or charging respectively. Similarly, the schedule of diesel generators
$\mathcal{S}_{DG}(k)$ indicates the amount of power to be fed to the microgrid
for the duration $[k,k+1]$. The input $\mathcal{H}(k)$ of HVAC system is
either parameters $w_{E}$, $w_{c}$ or control inputs $\dot{m}$ of HVAC. As
mentioned earlier, the $w_{E}$ and $w_{c}$ are the design parameters for HVAC
control, and the energy consumption of the HVAC system depends on the choice
of them [12]. As a result, the suitable values for $w_{E}$ and $w_{c}$ are
chosen by the supervisory controller to minimise the electrical expenses and
the model-based HVAC control system generates the control inputs based on
these choices of values. Therefore, the supervisory controller generates
$w_{E}$ and $w_{c}$, when the model of HVAC is known and generates control
inputs $\dot{m}$ of HVAC when the model of HVAC is unknown. The operation of
the supervisory controller can be formally described by

$\\{\mathcal{S}_{ESS}(k),\mathcal{S}_{DG}(k),\mathcal{H}(k)\\}=f(x_{m}(k),\mathcal{B}(k),\mathcal{W}(k),v(k))$
(7)

where, $\mathcal{B}(k),\mathcal{W}(k),v(k)$ represent the power demand of the
building network, weather data forecast and real-time price of grid
respectively for the interval $[k,k+1]$.

For designing this supervisory controller, we explore a DDPG-based DRL
technique and reduce the dependability of the historical data. We describe a
systematic approach for applying such a learning technique to design the MEMS
supervisory controller for Scenario-2 and 3, for achieving the optimisation
goal, i.e. the co-optimisation of two subsystems (microgrid and the building
HVAC system).

##  IV DRL Based Controller Design for Partially Known System model

In microgrids, the operations of ESS and HVAC systems are intertwined, which
means that the current decision in any of the systems would affect the future
decision of the other system. So, we formulate the sequential decision making
problem in the microgrid as an MDP (Markov Decision Process). In this case,
instead of using the air mass flow rate for regulating the temperature, we
predict the values of the MPC parameters($w_{E}$ and $w_{c}$) within a range
for minimising the electrical consumption of the building load further and
also maximising the comfort factor.

MDP Formulation: A MDP is formally defined as a four-tuple $M=(s,a,P,R)$,
where $s$ is the set of environment states, $a$ is the set of actions, $P$ is
the state transition probability function and $R$ is a reward function.

In this work, the agent is defined as the learner and decision-maker (i.e.,
microgrid with HVAC controller), while the environment is comprised of objects
like renewable generators, ESS, HVAC system, utility grid, indoor/outdoor
temperature etc. The agent observes environment state $s_{t}$ and takes action
$a_{t}$. Then, environment state becomes $s_{t+1}$ and the reward $R_{t+1}$ is
returned. In the following parts, we will design the key components of the
MDP, including environment state, action and reward function.
The optimal solution of the given problem can be solved by observing the
current system states. If $SOC$ denotes the State of Charge (SOC) of ESS and
$P^{max}_{ESS}$ is the maximum capacity of ESS. For preventing overcharging
and over-discharging of ESS, we introduce a factor $\gamma$ as a safety
factor. Then we can write

$\gamma.P^{max}_{ESS}\leq SOC\leq(1-\gamma).P^{max}_{ESS}$ (8)

We define the net power demand $P_{net}^{t}$ of the building load (assuming
HVAC load $P_{H}^{t}$ as the controllable load and other loads $P_{const}^{t}$
as constant) at a time step $t$ as

$P_{net}^{t}=P_{H}^{t}+P_{const}^{t}-P_{solar}^{t}$ (9)

We also assume that there is no loss in ESS while charging, discharging or
staying idle. $P_{solar}^{t}$ is the predicted PV power for time step $t$. The
time is repeated with a period of 24 hours, the state of the time
stamp($\widehat{t}$) is considered to be repeated over 24 hours.

Environmental State: We denote the real-time or dynamic electricity price of
the utility grid as $v_{t}$ and the parameters related to two contradictory
factors (energy consumption and comfort) of the MEMS as $w_{E}$ and $w_{c}$.
The state of the Microgrid with the building load can be written as
$s_{t}=[SOC^{t}\quad P_{net}^{t}\quad v_{t}\quad T_{in}\quad T_{out}\quad
w_{c}\quad\widehat{t}]\in$ $S$ where $S$ is the state space of the energy
management system of the building.

Action Space: Let us denote the power exchange between the microgrid and
utility grid as $P_{ut}$ and ESS charging/discharging power as $P_{ESS}^{t}$.
The range of $P_{ESS}^{t}$ is [$-d_{max},c_{max}$] where $d_{max}$ is the
maximum discharging power and $c_{max}$ is the maximum charging power. For
maintaining the power balance, the power supply should be equal to the power
demand. Then we can have

$P_{ut}+P_{solar}^{t}-P_{ESS}^{t}=P_{const}^{t}+P_{H}^{t}$ (10)

Again we know $P_{H}^{t}$ is a function of $w_{E}$ and $w_{c}$. For reducing
the overall electrical consumption, we have to tune these two parameters. Let
the action related to the tuning parameters be $T_{p}$. In this case, the
action space $A$ is: $A=\\{P_{ESS}\quad T_{p}\\}$.

Reward: Let $R(s_{t},a_{t})$ be the reward function that gives the cost of the
energy that the microgrid sell to the utility company. The reward function can
be defined as :

$\displaystyle R(s_{t},a_{t})$ $\displaystyle=r^{b}+r^{net}\quad\text{where},$
(11) $\displaystyle r^{b}$
$\displaystyle=-[SOC^{t}.u_{t}+\epsilon.(P^{max}_{ESS})],\quad
SOC^{t}\leq\gamma.P^{max}_{ESS}$ (12) $\displaystyle r^{b}$
$\displaystyle=-[SOC^{t}.u_{t}+(1-\epsilon).P^{max}_{ESS}],SOC^{t}\geq\gamma^{\prime}.P^{max}_{ESS}$
(13) $\displaystyle r^{net}$
$\displaystyle=-[P_{net}^{t}.u_{t}+\kappa^{up}(w_{c}-w_{c}^{max})],\quad
w_{c}>w_{c}^{max}$ (14)
$\displaystyle=-[P_{net}^{t}.u_{t}+\kappa^{low}(w_{c}^{min}-w_{c})],\quad
w_{c}<w_{c}^{min}$ (15) $\displaystyle=-[P_{net}^{t}.u_{t}],\quad otherwise.$
(16)

Here $\epsilon$ is the penalty factor for overcharging or overdischarging the
ESS and $\kappa^{up}$, $\kappa^{low}$ are penalty factors for the MPC based
HVAC parameter values if they are out of range and
$\gamma^{\prime}=(1-\gamma)$. The values of the parameters are in the range
[0,1]. Selling electricity price $u_{t}$ can be defined as $u_{t}=\sigma
v_{t}$ where $\sigma$ is a constant.
Action-Value Function: The estimation of the function $f$ illustrates the
relationship between the parameters, the total power consumption $P_{total}$
and comfort factor:

$\widehat{f}(w_{c},w_{E},P_{t-1}^{total})=P_{t}^{total}$ (17)

where $\widehat{f}$ represents the approximated function that explains the
relationship between the input data.

##  V DRL Based Controller Design for unknown System model

In this section, we try to optimise the MEMS with renewable power sources,
energy storage system, HVAC systems without having any knowledge about the
system. To be explicit, our goal is to minimise the electricity cost by
maintaining a comfortable indoor temperature. However, due to the
unavailability of accurate dynamics of indoor temperature, which is affected
by many parameters [24], it is tough to track the statistical distributions of
all the random system parameters like renewable generation output, outdoor
temperature and electricity price.

Unlike the partially known system method, here we consider the whole building
network along with the microgrid to be unknown along with HVAC controller
dynamics. This motivates the case of a pure DRL based approach where a DRL
agent is responsible for controlling both the HVAC control surfaces and the
microgrid settings.

MDP Formulation: Many control tasks like building energy management systems,
having the stochastic idea of the environment (e.g., occupants, weather)
frequently include continuous actions (e.g., temperature setpoint, supply air
mass flow rate). A mere discretisation of the continuous action space rapidly
gets intractable because of the dimensionality problem [6]. As discussed
earlier, we present a Deep Deterministic Policy Gradient (DDPG) based DRL
model to represent the controller of the HVAC system. The Markov Decision
Process (MDP) model for the DDPG is given as explained below.:

Figure 2: The agent-environment interaction in the MDP.

Environmental State: The state space $S$ of the MDP can be represented as a
five dimensional vector $s_{t}=[SOC^{t}\quad v_{t}\quad T_{all}(i)\quad
T_{\delta}\quad T_{out}]$ at a given time step $t$. Here, $e_{t}^{b}$ is the
amount of charge on the battery, $T_{all}(i)$ is the temperature of all the
zones for ith room and $T_{\delta}$ is the difference in temperature from the
comfortable range and $T_{out}$ is the outside air temperature. The vector
$T_{all}$ is cycled in an orderly fashion depending on the room it is
representing such that the temperature of the current room is represented by
the first dimension. For instance, if for room $1$, $T_{all}(1)$ is given by
$[T_{1}\quad T_{2}\quad T_{3}\quad...\quad T_{n}$], then for room $2$, the
same would be given by $[T_{n}\quad T_{1}\quad T_{2}\quad...T_{n-1}]$. Here
$n$ is the total number of rooms.

Action Space: The action space $A$ for the model is given by $a_{t}$=
$[P_{ESS}^{t},\quad$ṁ$]$, where $P_{ESS}^{t}$ at a given time step $t$, is the
amount of charging/discharging power of the battery and ṁ is the output air
mass flow rate.

Reward: The reward function is defined as

$r=-(power\textunderscore cost)-\lambda T_{\delta}$ (18)

where $power\textunderscore cost$ is the cost of the power given the current
pricing scheme. The cost is taken to be directly proportional to the air mass
flow rate. $T_{\delta}$ is the difference in temperature of the room and the
comfortable range of temperature and $\lambda$ here is a hyper-parameter.

Action-Value Function: While co-optimizing the ESS and the HVAC system at time
step t, the MEMS agent intends to maximize the total reward collected over the
future [22]. Let $Q(s,a)$ be the action-value function under the DDPG policy
$\pi$, which represents the expected return if action $a_{t}=a$ is taken in
state $s_{t}=s$ under the policy $\pi$ . Then, the optimal action-value
function $Q^{\ast}(s,a)$ is max $Q(s_{t},a_{t})$ which can be evaluated using
Bellman equation, i.e.,

$\displaystyle Q^{\ast}(s,a)$ $\displaystyle=\mathbb{E}[r_{t+1}+\gamma
max_{a\textquotesingle}Q^{\ast}(s_{t+1},a\textquotesingle)|s_{t}=s,a_{t}=a]$
(19)
$\displaystyle=\sideset{}{{}_{s\textquotesingle,r}}{\sum}P(s\textquotesingle,r|s,a)[r+\gamma
max_{a\textquotesingle}Q^{\ast}(s\textquotesingle,a\textquotesingle)$ (20)

where
$s\textquotesingle\in\mathcal{S},r\in\mathcal{R},a\textquotesingle\in\mathcal{A}$
and $P\in\mathcal{P}$.
The state transition probabilities $P(s\textquotesingle,r|s,a)$ are required
for evaluating $Q^{\ast}(s,a)$. Due to available noises and disturbances, it
is difficult to obtain these transition probabilities accurately as mentioned
earlier. To support the case with continuous system states, we use DDPG based
algorithm.

Figure 3: DDPG Based HVAC-Aware Microgrid Control Framework.

##  VI DDPG Based HVAC-Aware Microgrid Control Algorithm

As discussed earlier, for solving this MDP problem, we propose a DDPG-based
algorithm instead of DQN [17]. In the proposed Deep neural network model, the
agent learns how the changes in MPC-based HVAC controller parameters affect
the total power consumption of the building.

1

2Initialize online critic network $C:Q(s_{t},a_{t}|\theta^{Q})$ and actor
network $A:\mu(s_{t}|\theta^{\mu})$ with random weights $\theta^{Q}$ and
$\theta^{\mu}$ respectively.

3Initialize the target critic network
$C^{\prime}:Q^{\prime}(s_{t},a_{t}|\theta^{Q^{\prime}})$ and target actor
network $A^{\prime}:\mu^{\prime}(s_{t}|\theta^{\mu^{\prime}})$ with
$\theta^{Q^{\prime}}\leftarrow\theta^{Q}$ and
$\theta^{\mu^{\prime}}\leftarrow\theta^{\mu}$ respectively.

4Initialise replay buffer $B$ of size $N$.

5for  _episode m=1,2,…….,M_ do

6 Obtain the initial system states $s_{0}$ for a random day in the training
set.

7 Initialize a random exploration noise $n_{t}$.

8 for  _time step(in hours) t=1,2,……,T_ do

9 Obtain control action $a_{t}$ according to equation(21)

10 Execute this control action $a_{t}$ in the environment, observe the reward
$r_{t}$ after transition to the new system state $s_{t+1}$ at the end of time
slot $t$.

11 Store the experience $(s_{t},a_{t},r_{t},s_{t+1})$ in the replay buffer
$B$.

12 Randomly sample a mini-batch of $K$ transitions
$(s_{i},a_{i},r_{i},s_{i+1})$ from replay buffer $B$, $1\leq i\leq K$.

13 Estimate reward for each sampled transition using equation (22)

14 Update $C$ using equation (23) over the sampled mini-batch.

15 Update $A$ using sampled policy gradient as shown in equation (24)

16 Update target networks $C^{\prime}$ and $A^{\prime}$ using equations (25)
and (26) respectively.

17 end

18

19 end

20

Algorithm 1 Training Procedure with DDPG

The detailed DDPG based training procedure is described in Algorithm 1. The
DDPG based policy adopts an actor-critic framework based on Deterministic
Policy Gradient (DPG) [21]. The actor accepts the current state and predicts
the actions as per our action space. It is trained using a mean squared error
(MSE) loss. The output of the critic network is scalar, and it predicts the
reward. It tries to map the reward to the reward function that we have
provided for the environment. This network too is trained using MSE loss. The
actor-network is defined as $a_{t}=\mu(s_{t}|\theta^{\mu})$, where $a_{t}$
represents the control action, $s_{t}$ is the system state and $\theta^{\mu}$
represents the weights of the actor-network. The critic network is defined as
$Q(s_{t},a_{t}|\theta^{Q})$, where $a_{t}$ is the specified control action by
the actor-network and $\theta^{Q}$ represents the weights of the critic
network. The action-value function, $Q(s_{t},a_{t}|\theta^{Q})$ defines the
expected reward by taking action $a_{t}$ at state $s_{t}$ by following the
policy.

DDPG uses a replay buffer to perform a random sampling of experiences. It
updates the network parameters accordingly. All the features of the
experiences are saved $(s_{i},a_{i},r_{i},s_{i+1})$ and stored in a finite-
sized (say $N$) cache referred as a “replay buffer.” Then, random mini-batches
of experiences are sampled while the value and policy networks are updated.
The training for the DDPG network while interacting with the HVAC-aware MEMS
environment. At each time slot t, we first obtain the system states, and send
these states as input to the policy network, which finally give the output as
the control action. We need to explore the state space so that the policy does
not converge to local optimal solutions during the training. Therefore, for
exploration we add a random noise to the obtained control action,

$a_{t}=\mu(s_{t}|\theta^{\mu})+n_{t}$ (21)

where $n_{t}$ is the random noise. In our work, we use an Ornstein-Uhlenbeck
process [9] for generating the noise, $n_{t}$, for exploration. After that the
obtained $a_{t}$ is applied to our system. At the end of time slot $t$, the
new thermal state $s_{t+1}$ is obtained and then the total reward $r_{t}$
during that time slot. Each transition $(s_{t},a_{t},r_{t},s_{t+1})$ will be
stored in the replay buffer B for training the actor network and the policy
network. After that, we will randomly sample N transitions from the replay
buffer for training the network. For each transition
$(s_{i},a_{i},r_{i},s_{i+1})\in N$, the estimated reward is calculated as
follows where $\Gamma$ is the discount factor for model update,

$r_{i}=r_{i}+\Gamma
Q^{\prime}(s_{i+1},\mu^{\prime}(s_{i+1}|\theta^{\mu^{\prime}})|\theta^{Q^{\prime}})$
(22)

Then the critic network is updated by minimising the MSE error(i.e. loss)
between the estimated reward ($r_{0}$ calculated from Eq.(22)) and the reward
predicted by the critic network $Q(s_{i},a_{i})$ over the sampled mini-batch
as shown below,

$Loss(L_{\theta})=(1/N)\sum_{i}(y_{i}-Q(s_{i},a_{i}|\theta^{Q}))^{2}$ (23)

After that, the actor network is updated using the policy gradient [16] as
follows,

$\Delta_{\theta^{\mu}}=\Delta_{a}Q(s,a)\Delta_{\theta^{\mu}}\mu(s|\theta^{\mu})$
(24)

Finally, the target critic and actor networks are updated using the following
equations respectively,

$\theta^{Q^{\prime}}=\tau\theta^{Q}+(1-\tau)\theta^{Q^{\prime}}$ (25)
$\theta^{\mu^{\prime}}=\tau\theta^{\mu}+(1-\tau)\theta^{\mu^{\prime}}$ (26)

After the completion of the training, only the actor network can make control
action. The parameters i.e. the weights of the actor and critic networks
become fixed after the training as illustrated in Algorithm 1. Now the trained
DDPG model can be deployed to provide the real-time energy management strategy
for the HVAC-aware microgrid system. In Algorithm 2, the decision-making
process is outlined. For a certain day, the MEMS agent observes the current
system state $s_{t}$, and sets the control actions according to the learned
policy by DDPG at each time step t. For the partially known system model, the
actions are battery charging/discharging power and tuning the parameters of
the system. The actions are battery charging/discharging power and air mass
flow rate for the completely unknown system model. Then the obtained control
actions are mapped to energy schedules of the physical different units in the
MEMS (Section III) along with HVAC of the building.

1

2Load the parameter $\theta^{\mu}$ of the online actor-network after the
training procedure as illustrated in Algorithm 1.

3for  _day =1,2,……., $N^{test}$_ do

4 Obtain the initial system states $s_{0}$ for the day.

5 for  _time step(in hours) t=1,2,……,T_ do

6 Obtain control action $a_{t}$ according to equation(21).

7 Execute this control action $a_{t}$ in the MEMS environment, observe the
reward $r_{t}$ after transition to the new system state $s_{t+1}$ at the end
of time slot $t$. Use $s_{t+1}$ for finding $a_{t+1}$.

8 end

9

10 end

11

Algorithm 2 DDPG Based HVAC-aware MEMS Control Strategy

Figure 4: RMSE (%) in Temperature: subfigures a-g for each zone, subfigure h
for all zones.

##  VII Simulation Setup

We implement the MEMS along with the building loads with the HVAC system for
capturing the dynamics of the overall system. The main components of the
system are illustrated in Figure-1. We use OpenModelica [2] platform to
simulate the HVAC system and the thermal dynamics in the building. In the
simulations, we have assumed that there is no change in building occupancy or
its thermal load due to radiative forcing. The experiments are conducted on a
Linux system, with 8-core $3.4$ GHz each and $8$ GiB of RAM. These experiments
are implemented using Pytorch $1.5.1$ [19]. We start with an MPC based
controller for HVAC (designed using Matlab MPC toolbox). In this simulation,
we worked with seven zones. Zones represent different segments of a building
that have negligible thermal coupling. This controller takes two variables as
input - a state vector for the temperature of each zone, and the scalar
parameter $w_{c}$ (as defined earlier), which is set manually (to an initial
value of $0.5$). The output produced by the controller is the air mass flow
rate for each zone. We have used the weather data [4] of Kolkata, India, for
the experiment. The electricity buy price for our system follows the real-time
price structure provided in [3].

For the DDPG implementation, we use two DNNs (i.e. online and target) for the
actor and the critic, respectively. The actor network and the critic network
have two hidden layers with 128 neurons in each layer. We use the ReLU
activation function and batch normalisation in each layer. We adopt Adam [14]
for gradient-based optimization and the learning rate is 0.001. We use the
target network’s soft update parameter, $\tau$ as 0.001 and the batch size
128. The duration of each time slot is 30 minutes, and each episode consists
of 48 time slots. The initial exploration noise scale $\epsilon$ is taken as
0.35. Model parameter $\gamma$ is set to 0.05, $\sigma$ to 0.3, and dmax is
set to 0.9, along with cmax. We take the value of $\kappa_{up}$ as well as
$\kappa_{low}$ as 0.5. Here, $\lambda$ is assigned a value of 0.15. The
partially known model(Combo) is trained for 50 epochs, whereas the completely
unknown model(DRL) is trained for 100 epochs.

##  VIII Results & Discussions

Figure 5: Subfigures (a-g) denote the Temperature for each zone, subfigure (h)
denotes the Average Temperature for all zones.

We present comparative plots among the three modes of abstraction. Figure-5
gives the average temperature for each zone, along with an overall comparison.
Figure-4 shows the Root Mean Square Error of temperature with respect to the
desired temperatures for each zone. The average power usage of all zones has a
similar curve as Figure-5. Here the DRL model (scenario-3) consistently
performs at par (and even better at times) with that of the MPC model
(scenario-1). However, the MPC model optimised with a DRL model (the combo
scheme, i.e. scenario-2) does not give better results in comparison with the
other two. As mentioned earlier, the day has been divided into 48 parts at
half an hour intervals. All three plots validate similar inference. RMSE is
computed at intervals of 3 hours, leading to eight such points in a 24 hour
period.

Table-I provides the quantitative results of our work at various levels of
abstraction of the system. The word “Combo” is used to designate the partially
known model. As evident from the values, the mean values of the completely
unknown (DRL) model are comparable (and at times better) w.r.t. those of the
completely known (MPC) model. Hence, in this specific case, even in the
presence of abstractions, AI-based learning algorithms can successfully mimic
the scenario. The Table-I also highlights that the standard deviation of
average temperature is less in the case of DRL compared to MPC. This is a
positive result owing to the fact that drastic changes in temperature, even
within the comfortable range of temperature, is undesirable. This is because
the DRL learns from experience and avoids unnecessary fluctuations. However,
the standard deviation is less for the partially known model though the
model’s power consumption was significantly higher.

| Power Consumption |  Temperature (∘C)
---|---|---
(kW) | Average | RMSE
| Mean
MPC | 1.496 | 25.971 | 1.602
DRL | 1.384 | 25.907 | 1.506
Combo | 2.035 | 26.885 | 2.180
| Standard Deviation
MPC | 4.170 | 6.025 | 3.840
DRL | 4.208 | 5.909 | 3.853
Combo | 5.199 | 5.870 | 4.591
TABLE I: Mean and Standard Deviation of Power Consumption and Temperature for
various levels of abstractions of the system model

Adapting an MPC model for a new building environment requires additional
parameters while incorporating precise environment features into the model,
leading to an increase in the chances of errors. However, using DRL, an
existing model can be retrained and installed, which leads to saving efforts
and cost towards deployment.

##  IX Conclusion

In this work, we attempt to develop two levels of abstraction for the
completely known model of MEMS. The first level of abstraction deals with an
MPC-based method amalgamated with DRL. The DRL tunes certain parameters of the
MPC model instead of it being directly derived from the environment. The next
layer of abstraction replaces the completely known model with a black-box
model, which is completely simulated using a DRL architecture. We find that
this model performs at par with our specific test cases and exhibits less
standard deviation w.r.t. the MPC-based MEMS model. Using a DRL model has
certain design advantages leading to efficient deployment. One of the future
work would be using Twin Delayed Deep Deterministic Policy Gradients (TD3)[8]
instead of DDPG in our problem.

## References

  * [1] Energy.gov http://www.energy.gov/news2009/7550.html, 2009.
  * [2] Openmodelica https://www.openmodelica.org/, 2017.
  * [3] Real-time market of electricity in india https://www.iexindia.com/marketdata/rtm_areaprice.aspx, 2018.
  * [4] Weather data (kolkata) https://energyplus.net/weather-location/asia_wmo_region_2/IND//IND_Kolkata.428090_ISHRAE, 2018\.
  * [5] Abdul Afram and Farrokh Janabi-Sharifi.  Theory and applications of hvac control systems-a review of model predictive control.  Building and Environment, pages 343–355, 2014.
  * [6] D. P Bertsekas.  Dynamic programming and optimal control.  Belmont, MA: Athena scientific, 1(2), 1995.
  * [7] Jongwoo Choi, Youngmee Shin, and Choi.  Robust control of a microgrid energy storage system using various approaches.  IEEE Transactions on Smart Grid, pages 2702–2712, 2018.
  * [8] Scott Fujimoto, Herke Hoof, and David Meger.  Addressing function approximation error in actor-critic methods.  In International Conference on Machine Learning, pages 1587–1596. PMLR, 2018.
  * [9] L. S. Ornstein G. E. Uhlenbeck.  On the theory of the brownian motion, 1930.
  * [10] Shixiang Gu, Timothy Lillicrap, Ilya Sutskever, and Sergey Levine.  Continuous deep q-learning with model-based acceleration.  JMLR.org, 2016.
  * [11] Achin Jain and Derek Nong.  Digital twins for efficient modeling and control of buildings: An integrated solution with scada systems.  In 2018 Building Performance Analysis Conference and SimBuild, 2018\.
  * [12] Rajib Lochan Jana, Soumyajit Dey, and Pallab Dasgupta.  A hierarchical hvac control scheme for energy-aware smart building automation.  ACM Transactions on Design Automation of Electronic Systems (TODAES), 2020.
  * [13] Anthony Kelman and Francesco Borrelli.  Bilinear model predictive control of a hvac system using sequential quadratic programming.  IFAC Proceedings Volumes, 44(1):9869–9874, 2011.
  * [14] Diederik Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  International Conference on Learning Representations, pages 1–15, 2015.
  * [15] Elizaveta Kuznetsova, Yan-Fu Li, Carlos Ruiz, Enrico Zio, Graham Ault, and Keith Bell.  Reinforcement learning for microgrid energy management.  Energy, 59:133–146, 2013.
  * [16] T. Lillicrap, J. Hunt, A. Pritzel, N. Heess, T. Erez, Y. Tassa, D. Silver, and Daan Wierstra.  Continuous control with deep reinforcement learning.  CoRR, abs/1509.02971, 2016.
  * [17] V. Mnih.  Human-level control through deep reinforcement learning.  Nature, 518:529–541, 2015.
  * [18] Thomas Morstyn, Branislav Hredzak, Ricardo P Aguilera, and Vassilios G Agelidis.  Model predictive control for distributed microgrid battery energy storage systems.  IEEE Transactions on Control Systems Technology, 26(3):1107–1114, 2017.
  * [19] Adam Paszke, Sam Gross, Francisco Massa, and Lerer.  Pytorch: An imperative style, high-performance deep learning library.  2019\.
  * [20] Nagaraju Pogaku, Milan Prodanovic, and Timothy C Green.  Modeling, analysis and testing of autonomous operation of an inverter-based microgrid.  IEEE Transactions on power electronics, 22(2):613–625, 2007.
  * [21] David Silver, Guy Lever, Nicolas Heess, Thomas Degris, Daan Wierstra, and Martin Riedmiller.  Deterministic policy gradient algorithms.  JMLR.org, 2014.
  * [22] R.S. Sutton and A.G. Barto.  Reinforcement learning: an introduction.  The MIT Press, London, England, 2018.
  * [23] Kasemsak Uthaichana and DeCarlo.  Hybrid optimal theory and predictive control for power management in hybrid electric vehicle.  2018\.
  * [24] Tianshu Wei, Yanzhi Wang, and Qi Zhu.  Deep reinforcement learning for building hvac control.  In Proceedings of the 54th Annual Design Automation Conference 2017, pages 1–6, 2017.
  * [25] Hongbin Wu, Xingyue Liu, and Ming Ding.  Dynamic economic dispatch of a microgrid: Mathematical models and solution algorithm.  International Journal of Electrical Power & Energy Systems, 2014\.
  * [26] Shichao Xu, Yixuan Wang, Yanzhi Wang, Zheng O’ Neill, and Qi Zhu.  One for many: Transfer learning for building hvac control.  2020\.