# Stripe Domains and First Order Phase Transition in the Vortex Matter of
Anisotropic High-temperature Superconductors

V. K. Vlasko-Vlasov1, J. R. Clem2, A. E. Koshelev1, U. Welp1, W. K. Kwok1 1
Materials Science Division, Argonne National Laboratory, Argonne, Illinois
60439  2 Department of Physics and Astronomy, Iowa State University, Ames, IW
50011-3160,USA

###### Abstract

We report the direct imaging of a novel modulated flux striped domain phase in
a nearly twin-free YBCO crystal. These domains arise from instabilities in the
vortex structure within a narrow region of tilted magnetic fields at small
angles from the in-plane direction. By comparing the experimental and
theoretically derived vortex phase diagrams we infer that the stripe domains
emerge from a first order phase transition of the vortex structure. The size
of domains containing vortices of certain orientations is controlled by the
balance between the vortex stray field energy and the positive energy of the
domain boundaries. Our results confirm the existence of the kinked vortex
chain phase in an anisotropic high temperature superconductor and reveal a
sharp transition in the state of this phase resulting in regular vortex
domains.

In type II superconductors (SC) the magnetic field penetrates in the form of
vortices - magnetic flux tubes each carrying a single flux quantum 1 .
Interactions between vortices and their coupling to crystal defects result in
a rich variety of vortex structures. They include a triangular or square
vortex lattice that can subsequently melt at high temperatures to a vortex
liquid, a pinned vortex glassy state that can sustain high current carrying
capacity, and in highly anisotropic superconductors, 2D pancake vortices that
can interact with in-plane Josephson vortices to create 1D vortex chain
states. These states of vortex matter define all transport and magnetic
properties of applied superconductors.

Here we report on an unusual vortex domain structure that arises following a
first order phase transition in the vortex state in nearly untwined YBCO
crystal under tilted magnetic fields. Typically, vortices repulse each other
when the distance between them is much smaller than the penetration depth
$\lambda$, the scale over which the magnetic field of the vortices decays. In
this case, the vortices can arrange into a uniformly spaced lattice or form
smooth density gradients defined by the critical current, the maximum current
that the superconductor can sustain before reverting to the normal state. In
contrast, at larger distances the interactions between vortices can be
attractive, resulting in the formation of vortex chains and bundles (see
review 2 ). For example, periodic stripe domains and circular bundles of
vortices interspersed with flux free Meissner regions occur in thin niobium
disks with a low Ginzburg-Landau parameter $\kappa=\lambda/\xi\sim 1$, where
$\xi$ is the coherence length setting the size of the vortex core. These
vortex domains are similar to the alternating normal (N) and superconducting
Meissner (M) domains found in the intermediate state of type I SCs with a
nonzero demagnetizing factor. In the latter, the lamellae M/N domains are
known to minimize the magnetostatic energy of the normal regions at the
expense of the positive energy of the boundaries between the normal and SC
phases.

Vortex attraction can also appear in large $\kappa$, anisotropic type II SCs
when an applied magnetic field is tilted from the anisotropy axis 3 ; 4 ; 5 .
Here, the circulating supercurrents that form the vortex tend to flow
perpendicular to the anisotropy axis rather than to the vortex line. This
creates an inversion of the magnetic field at some distance from the vortex
core resulting in an attraction of the neighboring vortex within the tilt
plane 3 ; 5  and formation of dense vortex chains coexisting with the dilute
lattice of Abrikosov vortices, as observed in decoration, Hall probe, and
electron microscopy experiments in layered high-Tc superconductors (see review
6 ).

Formation of patterns including stripe structures is general for a wide
variety of condensed matter and colloidal systems with competing attractive
and repulsive interactions, which are usually extended at different length
scales 7 ; 8 ; 9 . In fact, long range attraction and short range repulsion of
vortices in superconductors with two order parameters [10] as also in
multilayers of different superconductors 11 ; 12  can result in formation of
vortex stripes. Experimentally, stripes of vortices intermittent with Meissner
phase were observed in $MgB_{2}$ single crystals 13 ; 14  and were treated as
a possible consequence of the two different SC condensates in the 2-band
magnesium diborate.

In this work we discovered a new vortex structure of regular stipe domains
with alternating flux density, $B_{z}$, generated by magnetic fields tilted
from the ab-plane in a nearly twin-free YBCO platelet crystal. Unlike
vortex/Meissner domains and vortex chain structures, our vortex domains are
created via an instability in the vortex state 15 ; 16 ; 17 ; 18 ; 19 ; 20 ;
21 ; 22  due to a first order phase transition characterized by an abrupt jump
in the vortex tilt angle 23 . A possibility of such first order transition in
the vortex orientation was first noticed by Buzdin and Simonov in 24 . During
the phase transition, an intermediate state with domains of different vortex
phases (orientations) emerges. We propose that similar to the case of type I
SCs and ferromagnets, the regular structure of domains with alternating
density of $B_{z}$ reduces the magnetostatic energy of the vortex stray fields
perpendicular to the surface of the sample. Between the domains with different
flux orientations vortices gradually rotate forming a domain boundary.
Responsible for this rotation supercurrents yield a positive domain wall
energy which limits the refinement of domains caused by magnetostatics. We
suggest that the attractive coupling of tilted vortices into chains along the
vortex tilt plane 3 ; 4 ; 5  and attraction between the neighboring chains
with mutually shifted vortices stabilizes domains with larger $B_{z}$. Our
calculated vortex structure phase diagram in tilted fields defines the
boundaries of the angular vortex instability in YBCO in qualitative agreement
with the experiment.

The $YBa_{2}Cu_{3}O_{7-\delta}$ rectangular platelet crystal (1130x340x20 $\mu
m^{3}$) used in our experiments was grown using a flux method. The crystal had
a few narrow twin lamellae near one end but most of the area was twin-free.
The onset temperature of the superconducting transition and the width were
$T_{c}$ = 92.4K and $\Delta T$ = 0.3K respectively, as determined from
magnetization measurements. We imaged the perpendicular magnetic fields
($B_{z}$) on the sample surface using a magneto-optical (MO) indicator
technique 25 . Data were obtained after cooling the sample in applied magnetic
fields tilted from the ab-plane of the crystal as well as during
remagnetization of the sample in tilted fields. The resulting intensity of the
MO images depicts the strength and polarity of the normal field $B_{z}$,
Fig.1a delineates the spatial field pattern at T = 60K after cooling the
sample in a field of $H_{//}$ = 1388 Oe parallel to the long edge of the
sample and then switching off the field. The dark (negative $B_{z}$) and
bright (positive $B_{z}$) image contrast on the right and left edges of the
sample demonstrates that trapped vortices induced by the applied field lie
parallel to the sample plane and their ends produce stray fields that diverge
up at the left edge and converge down at the right edge as shown in Fig.1d.
Similar images were observed when $H_{//}$ was tilted by less than $\sim
0.5^{\circ}$ from the sample plane.

At tilt angles $>0.5^{\circ}$ a qualitatively new stripe pattern emerges at
the center of the sample. Fig.1b shows such a pattern formed after cooling the
sample in a field $H_{//}$ = 1388 Oe with a tilt angle of $\sim 1.60^{\circ}$
and switching off the field. Similar to Fig.1a there is bright and dark
contrast at the short sample ends and a new dark contrast appears at the top
and bottom long edges of the crystal. The latter is associated with negative
stray field around the sample due to the positive $B_{z}$ trapped in the
central region of the sample. It is in this region the periodic stripes with
alternating $B_{z}$ values reside. Fig.1e shows the $B_{z}$ profile measured
across the stripes. It clearly reveals periodic oscillations of Bz
superimposed on a monotonic spatial induction gradient. The latter is typical
of perpendicular trapped flux in a superconducting plate and corresponds to
the average critical current circulating in the ab-plane of the sample 26 ; 27
; 28 . We associate the observed stripe pattern with a domain structure of
vortices with different tilt angles or, equivalently, with different density
of in-plane and $\it{c}$-axis vortex segments within the vortex staircase
structure inherent for cuprates 29 ; 30 ; 19 ; 23  (see Fig.1f). This novel
domain structures exist in a narrow range of tilt angles between $\theta\sim
0.5^{\circ}$ and $2.40^{\circ}$. Within this range, the width and length of
the trapped $B_{z}$ region and thus the length of the stripe domains were
slightly shorter for smaller angles. At a fixed angle the trapped $B_{z}$
region was shorter for smaller fields. Fig. 2 depicts the phase diagram of the
observed stripe domain structures obtained by cooling the sample to 60K under
various angles of applied fields $H_{//}$ ranging from $\sim$280 to 1388 Oe.

Figure 1: (a) Trapped flux (H=0) after field cooling from $T>T_{c}$ to 60K in
$H_{//}$ = 1388 Oe. (b) Same after cooling in H = 1388 Oe tilted by
$1.6^{\circ}$. (c) Same after cooling in H = 1388 Oe tilted by $3.5^{\circ}$.
(d) Profile of $B_{z}$ between two arrows shown in Fig.1b. (d) Scheme of the
stray fields corresponding to (a). (e) $B_{z}$ profile across the sample
between arrows shown in (b). (f) Scheme of the staircase vortices in cuprates.
(g) Scheme of the vortex domains. Currents (short arrows) in the domain wall
flow along vortices in force free configuration.

To track the initial stage of the domain nucleation we imaged the flux
structure during gradual reduction of the field after initially field cooling
the sample to T = 60K with $H_{//}$ = 1388 Oe tilted by $1.6^{\circ}$ (see MO
pictures in Supplemental Info). Before reducing $H_{//}$, the $B_{z}$
distribution is homogeneous over the entire imaged area. Decreasing field to
$\sim$ 900 Oe results in the appearance of a dark edge contrast due to the
stray fields coming from the trapped positive $B_{z}$. Simultaneously, stripes
of alternating bright and dark contrast (modulated $B_{z}$) emerge along
$H_{//}$ near short ends of the sample. Further decreasing the magnetic field,
extends the stripes over the entire length of the sample. At even smaller
$H_{//}\sim$150 Oe, vortices turn into the plane near the perimeter of the
sample and a region of trapped $B_{z}$ with stripe domains forms inside the
crystal. This configuration remains after reducing the field to zero. For
angles larger than $2.4^{\circ}$ the field cooling resulted in a larger region
of stronger trapped $B_{z}$ but without any stripe domain features (Fig.1c).

Figure 2: Experimental $H_{z}-H_{x}$ domain existence map.

The nucleation of domains could be also initiated if after cooling in
$H_{//}$=1388 Oe the field was first increased by a few hundred gauss and then
decreased back to the initial value. After field cooling the MO picture was
homogeneous. But after increasing and decreasing field back to 1388 Oe,
stripes appeared at the short ends of the sample. They extended over the
central area during further reduction of $H_{//}$ as described above.

Similar observations were also performed in tilted fields
$H_{\perp}(\theta_{\perp})$ perpendicular to the long sample side. In this
case the stripe domains formed along the width of the sample and were wider
and somewhat less regular than in $H_{//}$ (Fig.4a). They also exist in a
narrow range of tilt angles ($\theta_{\perp}\sim$ 0.4 to $3^{\circ}$) and
become shorter with decreasing $\theta_{\perp}$ (Fig.4b). At smaller angles
the total trapped $B_{z}$ is weaker and black and white contrast at the long
edges of the crystal appears due to the stronger in-plane flux $B_{y}$. The
domains emerge in the shape of narrow wedges stretching along $H_{\perp}$ from
the long edges during application of large enough field. At decreasing field
they detach from the edges and remain in the central region after switching
$H_{\perp}$ off.

So far there were no experimental reports of vortex domains in high-$\kappa$
SCs. The coexistence of vortex lattices with different tilt angles could be
expected from theoretical predictions of vortex angle instabilities. These
instabilities were calculated for some range of angles for (i) a single tilted
vortex line in anisotropic three-dimensional superconductors 15 ; 16 ; 17 ; 20
; 21 ; 22 , (ii) a single tilted vortex line in layered superconductor with
purely magnetic coupling19 . and (iii) a chain of tilted vortex lines 23 .
Below, we derive the vortex phase diagram for YBCO under tilted magnetic
fields by taking into account the energy of individual vortices and energies
of their interactions and their coupling to external fields. We then analyze
the stray fields in thin SC plates and show that they define the width of the
stripe domains and discuss vortex attraction as a possible mechanism for
stabilizing domains with larger $B_{z}$.

Figure 3: Vortex domains at 60K after application and switching off
$H_{\perp}$=2100 Oe along the short side tilted from the ab-plane by $\sim
1.20^{\circ}$ (a) and $\sim 0.40^{\circ}$(b)

The vortex phase diagram in tilted fields is constructed using analytic
results for vortex-chain energies derived in Ref. 23  (see 31  for details).
Our calculations for a bulk sample show abrupt changes in the vortex angle
which are mostly driven by the energy of an isolated vortex chain formed in
the tilt plane while interactions between the chains may be neglected. The
total energy of a vortex chain consists of isolated vortex line and inter-
vortex coupling contributions, $E_{TV}=E_{TV}^{s}+E_{TV}^{i}$. Both terms
depend on the in-plane and $\it{c}$-axis magnetic field components and on the
SC materials parameters 31 . The structure of vortices tilted by a small angle
from the in-plane direction, $\theta<1/\gamma$, can be described by the
addition of $\it{c}$-axis pancake vortices (kinks) to the Josephson vortex
line lattice in purely in-plane field. The equilibrium state is obtained by
minimizing the sum of $E_{TV}$ and the energy of interaction of the
$\it{c}$-axis flux with thermodynamic magnetic field $H_{z}$. We found that in
some range of parameters the total energy has minima at two values of $B_{z}$
and the system undergoes a first-order phase transition at the critical value
of $H_{z}$ resulting in a jump in $B_{z}$. Roughly, the transition is located
in the in-plane field range $B_{x}<\gamma\Phi_{0}^{2}/\lambda_{ab}^{2}$ and
bypassed range of $B_{z}$ corresponds to tilt angles around $\theta\sim
1/\gamma$. Fig. 4 shows numerically computed phase diagram obtained for
typical YBCO parameters, the penetration depth $\lambda_{ab}$ =200 nm,
distance between $CuO_{2}$ layers $\it{s}$= 1.2 nm, and anisotropy $\gamma$ =
7. In the region between two lines in Fig.4 the vortex structure experiences
an angular instability and the c-axis component of the flux, $B_{z}$, jumps
from $B_{z}^{min}$ to $B_{z}^{max}$ via a first order phase transition.
Although the sample shape and vortex pinning are not accounted for in this
theory, the calculated phase diagram is in reasonable agreement with our
experimental $H_{z}-H_{x}$ data of the emergence of the stripe domains
(Fig.2). The experimental $H_{z}$ fields limiting the domain existence are
somewhat smaller compared to $B_{z}$ of the theoretical instability lines. We
attribute this difference to the demagnetization factor of the sample and
infer the angular vortex instability accompanying the vortex angle transition
as the main reason for the appearance of domains.

Figure 4: Calculated phase diagram of the vortex state in YBCO under tilted
fields. With increasing tilt, the $B_{z}$ component jumps from the bottom to
top line. The intermediate state domains should form in the area between the
lines.

In large SC samples the sum of energies of individual vortices and bulk
vortex-vortex interactions is the main factor defining the flux state.
However, in thin SC plates the external stray fields induced by vortices come
into play. First, they define the long-range vortex interactions near the
surface, which can be treated as a repulsion of magnetic monopoles 32 .
Second, they can introduce modulated flux structures, which would reduce the
energy of the vortex stray fields around the superconductor ($E_{s}$). In type
I SCs and low-$\kappa$ materials this energy is responsible for the emergence
and size of the SC/N or Meissner/Shubnikov lamella domain structures 2 . We
conjecture that the regularity of our vortex stripe domain structure can also
be explained by the effect of the stray field energy, $E_{s}$. Similar to the
case of ferromagnetic domains, $E_{s}$ can be minimized by periodic
oscillations of $B_{z}$. In relatively small fields, we can consider our
system as comprised of up and down magnetized domains superimposed onto an
averaged trapped field background. The latter is not important for the
variational problem. Formulas for the stripe domains give $E_{s}\sim
1.7\mu_{0}^{-1}(\Delta B_{z}/2)^{2}D$ 33 . Here $\Delta B_{z}$ is the
difference of $B_{z}$ in neighboring domains and $D$ is the domain width. An
important factor limiting the refinement of domains is the positive energy at
their boundaries $\sigma_{B}$. The balance between the reduction of the stray
field energy and the increase of the boundary energy yields the equilibrium
size of the domains. Positive $\sigma_{B}$ is produced by currents responsible
for the rotation of vortices in the boundary layer between domains with
different flux direction. For example, let us consider two infinite blocks of
vortices of the same density $B_{0}$ turned by an angle $\alpha_{0}$
($\alpha_{1}-\alpha_{2}$ in Fig.1f) with a $x_{0}$-thick boundary between them
where $B_{0}$ rotates in a force-free way ($J(x)//B(x)$), so that
$B(x)=\hat{y}B_{0}sin(\alpha(x))+\hat{z}B_{0}cos(\alpha(x))$. The current
$J(x)$ flows only within the boundary layer and vanishes in the homogeneous
domains. In the boundary region the current density is
$J(x)=\mu_{0}^{-1}B_{0}(d\alpha/dx)$. For a linear variation of the vortex
angle $\alpha(x)=x(\alpha_{0}/x_{0})$ , the force-free current $J(x)$ has some
maximum density $J_{c//}$ with the kinetic energy density
$\mu_{0}\lambda^{2}J_{c//}^{2}/2$.The resulting energy per unit area of the
boundary will be $\sigma_{B}=J_{c//}B_{0}\alpha_{0}\lambda^{2}/2>0$.
Obviously, the real structure of the boundary is more complex, but the energy
due to the boundary currents should still be positive. Minimizing the energy
per unit area of the plate $E_{S}+E_{B}=1.7\mu_{0}^{-1}(\Delta
B_{z}/2)^{2}D+\sigma_{B}(d/D)$, where $d$ is the thickness of the sample,
yields the equilibrium size of the stripe domains,
$D=(4\mu_{0}d\sigma_{B}/1.7\Delta B_{z}^{2})^{1/2}$. Remarkably, this simple
formula gives the domain width $D=30\mu m$ (for $d=20\mu
m,J_{c//}=10^{7}A/cm^{2},\alpha_{0}=2^{\circ},\lambda=200nm$ and $\Delta
B_{z}$=10 G) very close to the observed values ($D\sim 20\mu m$).

One could expect that in the increased Bz domains there is an enhanced
repulsion of vortices. However, we suggest that the vortex attraction in the
chains along the tilt plane and attraction between the chains can stabilize
these domains. For a linear chain of tilted vortices the intrachain attraction
emerges through the pancake components of vortices and has a dipolar character
34 . The same dipolar attraction can occur between the chains if pancake
stacks in the neighboring chains are shifted by half-a-period along the chain.
A possible confirmation of this scenario could be zigzag vortex chains
observed in $NbSe_{2}$ in fields close to the basal plane 35 . This mechanism
can also account for the formation of chain bundles in BSCCO in fields
strongly tilted from the c-axis 36 .

In summary, we report the first direct imaging of a novel vortex stripe domain
phase in a nearly twin-free YBCO crystal. The domains with alternating $B_{z}$
arise from instabilities in the vortex state within a narrow region of tilted
magnetic fields at small angles from the ab-plane. By comparing the
experimental and theoretically derived vortex phase diagrams we infer that the
stripe domains emerge due to a first order phase transition of the vortex
structure. We show that the size of domains with different vortex angles is
defined by the balance of the vortex stray field energy and positive energy of
domain boundaries. Our results confirm the emergence of the kinked vortex
chain phase in an anisotropic high temperature superconductors and reveal
sharp transformations in their state.

Acknowledgements This work was supported by the US Department of Energy DOE
BES under Contract No. DE-AC02-06CH11357. One of us (JR) acknowledges support
from the Center for Emergent Superconductivity, an Energy Frontier Research
Center funded by the US Department of Energy, Office of Science, Office of
Basic Energy Sciences.

## References

  * (1) A. A. Abrikosov, J. Phys. Chem. Sol. 2, 199 (1957).
  * (2) E. H. Brandt, M. P. Das, J Supercond. Nov. Magn. 24, 57 (2011).
  * (3) A. M, Grishin, A.Y. Martynovich, S.V. Yampolskii , Zh. Eksp. Teor. Fiz. 97 1930 (1990); JETP 70, 1089 (1990); Physica B165/166, 1103 (1990).
  * (4) A. I. Buzdin, A. Yu. Simonov Pisma Zh. Eksp. Teor. Fiz. 51, 18 (1990); JETP Lett. 51, 191 (1990 ).
  * (5) V. G. Kogan, N. Nakagawa, and S. L. Thiemann, Phys. Rev. B 42, 2631 (1990).
  * (6) S. J. Bending and M. J. W. Dodgson, J. Phys.: Cond. Matt. 17, R955 (2005).
  * (7) G. Malescio and G. Pellicane, Nature Mater. 2, 97 (2003); Phys. Rev. E 70, 021202 (2004).
  * (8) C. J. Olson Reichhardt, C. Reichhardt, and A. R. Bishop, Phys. Rev. E 82, 041502 (2010).
  * (9) K. Nelissen, B. Partoens, and F. M. Peeters, Phys. Rev. E 71, 066204 (2005).
  * (10) E. Babaev and M. Speight, Phys. Rev. B 72, 180502(R) (2005).
  * (11) L. Komendova, M.V. Milosevic, and F.M. Peeters, Phys. Rev. B 88, 094515 (2013).
  * (12) C. N. Varney, K. A. H. Sellin, Q.-Z. Wang, H. Fangohr and E. Babaev, J. Phys. Cond. Mat., 25, 415702 (2013).
  * (13) V. V. Moshchalkov, M. Menghini, T. Nishio, Q. H. Chen, A. V. Silhanek, V. H. Dao, L. F. Chibotaru, N. D. Zhigadlo, and J. Karpinski, Phys. Rev. Lett. 102, 117001 (2009).
  * (14) J. Gutierrez, B. Raes, A. V. Silhanek, L. J. Li, N. D. Zhigadlo, J. Karpinski, J. Tempere, and V. V. Moshchalkov, Phys. Rev. B 85, 094511 (2012).
  * (15) A. Sudbo, E. H. Brandt, and D.A. Huse, Phys. Rev. Lett. 71, 1451 (1993).
  * (16) L. L. Daemen , L. J. Campbell , A. Yu. Simonov, and V. G. Kogan, Phys. Rev. Lett. 70, 2943 (1993).
  * (17) E. Sardella and M. A. Moore, Phys. Rev. B 48, 9664 (1993).
  * (18) G. Preosti and P. Muzikar, Phys. Rev. B 48, 9921 (1993).
  * (19) M. Benkraouda and J. R. Clem, Phys. Rev B 53, 438 (1996).
  * (20) A. K. Nguyen and A. Sudbo, Phys. Rev. B 53, 843 (1996).
  * (21) A.M. Thompson and M. A. Moore, Phys. Rev. B 55, 3856 (1997).
  * (22) P. Muzikar, J. Low Temp. Phys. 115, 135 (1999).
  * (23) A. E. Koshelev, Phys. Rev. B 71, 174507 (2005).
  * (24) A. I. Buzdin and Yu. A. Simonov, Supercond. Sci. Technol. 5, s284 (1992).
  * (25) V.K. Vlasko-Vlasov, U. Welp, G.W. Crabtree, and V.I. Nikitenko, in Physics and Materials Science of Vortex States, Flux Pinning and Dynamics, Vol. 356 of NATO Advanced Study Institute, Series E: Applied Sciences, edited by R. Kossowsky et al. Kluwer, Dordrecht, 1999, pp. 205-237.
  * (26) P. N. Mikheenko and Yu. E. Kuzovlev, Physica C 204, 229 (1993).
  * (27) E. H. Brandt and M. Indenbom, Phys. Rev. B 48, 12 893 (1993).
  * (28) E. Zeldov, J. R. Clem, M. McElfresh, and M. Darwin, Phys. Rev.B 49, 9802 (1994).
  * (29) B. I. Ivlev, Y. N. Ovchinnikov, and V. L. Pokrovskii, Mod. Phys. Lett. 5, 73 (1991).
  * (30) L. N. Bulaevskii, M. Ledvij, and V. G. Kogan, PRB 46, 366 (1992).
  * (31) Details of calculation of the phase diagram in tilted fields can be found in supplementary materials.
  * (32) G. Carneiro and E. H. Brandt, Phys. Rev. B 61, 6370 (2000).
  * (33) L. D. Landau and E. M. Lifshitz, Electrodynamics of Continuous Media (Pergamon, Oxford, 1984).
  * (34) A. Buzdin and I. Baladi Phys. Rev. Lett. 88, 147002 (2002).
  * (35) H. F. Hess, C. A. Murray, and J. V. Waszczak, Phys. Rev. Lett. 69, 2138 (1992).
  * (36) T. Tamegai, H. Aoki, M. Matsui, and M. Tokunaga, Pramana J. Phys. 66 , 271 (2006).

##  Appendix A Phase diagram for tilted chains

Supplemental Info for ”Stripe Domains and First Order Phase Transition in the
Vortex Matter of Anisotropic High-temperature Superconductors” by V. K.
Vlasko-Vlasov1, J. R. Clem2, A. E. Koshelev1, U. Welp1, W. K. Kwok1

1 Materials Science Division, Argonne National Laboratory, Argonne, Illinois
60439

2 Department of Physics and Astronomy, Iowa State University, Ames, IW
50011-3160,USA

Several theoretical studied demonstrated that in anisotropic superconductors
some orientations of dilute tilted vortex lattice are unstable [1s, 2s]. As a
consequence, there is a first order phase transition where the lattice
orientation angle jumps. We analyze behavior of the equilibrium c-axis
induction $B_{z}$ for uniformly tilted vortex lattice and find range of
parameters where $B_{z}$ has a jump. We consider tilted lattice composed of
tilted chains separated by the distance $c_{y}$. The in-plane distance between
tilted vortices in a chain is $a$ and the $c$-axis distance is $c_{z}$, see
Fig. 5, so that the tilting angle with respect to the $c$-axis is given by
$\tan\theta=a/c_{z}$. The components of magnetic induction are given by
$B_{x}=\Phi_{0}/c_{z}c_{y}$ and $B_{z}=\Phi_{0}/ac_{y}=B_{x}c_{z}/a$.

Figure 5: Scheme of tilted vortices.

To facilitate calculations, we make several simplifying assumptions. We
consider the range of small c-axis fields $B_{z}\sim B_{x}/\gamma$. This
allows us to neglect weak $B_{z}$ dependence of separation between the chains
and take the lattice parameters $c_{z}$ and $c_{y}$ to be the same as for the
Josephson vortex lattice for the field aligned with the layers. Two aligned
configurations of the Josephson vortex lattice are possible which are
degenerate within the London model. The $c$-axis lattice parameter is
connected with the in-plane magnetic field by relation
$c_{z}=\sqrt{\beta\Phi_{0}/\left(\gamma B_{x}\right)}$ with $\beta=2\sqrt{3}$
and $2/\sqrt{3}$ for these two configurations. Analysis of commensurability
oscillations at high fields [3s] suggests that the first configuration is
preferable. Our analysis shows that discontinuous behavior exists in the range
$c_{z}\gtrsim\lambda_{ab}/\gamma$. In this range the in-plane lattice
parameter $c_{y}$ exceeds the London penetration depth and interaction between
the tilted chains gives only small contribution to the total energy. The
ground-state c-axis vortex density is mostly determined by the energy of an
isolated tilted chain and we focus on the analysis of this energy [2s].

Figure 6: Phase diagram of the vortex state in YBCO in tilted fields. Between
the curves vortices are unstable.

The energy of an isolated tilted chain is composed of two parts, the energy of
isolated tilted vortices and the interaction term,
$E_{TV}=E_{TV}^{s}+E_{TV}^{i}$. To facilitate numerical evaluation, we present
the dimensionless energy $\tilde{E}_{TV}=\lambda_{ab}E_{TV}/\varepsilon_{0}$
with $\varepsilon_{0}=\Phi_{0}^{2}/(4\pi\lambda_{ab})^{2}$ as function of the
reduced planar density of pancake vortices $\tilde{n}=\lambda_{ab}/a$ and the
lattice parameter $\tilde{c}=c_{z}\gamma/\lambda_{ab}$. The field components
are related to these reduced parameters as

$B_{x}=\frac{\beta\gamma\Phi_{0}}{\lambda_{ab}^{2}\tilde{c}^{2}};\
B_{z}=\frac{\beta\Phi_{0}\tilde{n}}{\lambda_{ab}^{2}\tilde{c}}.$ (1)

Analytical expression for the energy of an isolated tilted vortex is only
available for two asymptotic regimes [2s]:

$\tilde{E}_{TV}^{s}\approx\tilde{E}_{PS}^{s}+\genfrac{\\{}{.}{0.0pt}{0}{\tilde{E}_{JV}^{s}+\tilde{n}\left(\ln\frac{\gamma}{l_{s}}-0.81+\frac{\tilde{n}\tilde{c}}{2}\left[\ln\left(l_{s}\tilde{n}\tilde{c}\right)-\frac{3}{2}\right]\right)\text{
for
}\tan\theta\gg\gamma}{\tilde{n}\ln\frac{\sqrt{1+\left(\gamma/\tilde{n}\tilde{c}\right)^{2}}+1}{2}+\frac{\ln\left(8l_{s}\tilde{n}\tilde{c}\right)-\gamma_{E}}{2\tilde{n}\tilde{c}^{2}}\text{
for }\tan\theta\ll\gamma},$ (2)

where $l_{s}=\lambda_{ab}/s$,
$\tilde{E}_{PS}^{s}=\tilde{n}\left[\ln(\lambda_{ab}/\xi)+0.497\right]$ is the
energy of aligned pancake stacks and
$\tilde{E}_{JV}^{s}=\frac{1}{\tilde{c}}\left(\ln l_{s}+1.54\right)$ is the
energy of an isolated Josephson vortices. The first asymptotics corresponds to
the regime of kinked lines. The interaction energy is given by

$\tilde{E}_{TV}^{i}\\!\approx\\!\frac{\pi\gamma}{\tilde{c}^{2}}\\!+\\!\pi\tilde{n}^{2}\\!-\\!\sqrt{\tilde{n}^{2}+\tilde{c}^{-2}}\left[\\!\ln\left(4\pi\gamma\sqrt{\tilde{n}^{2}+\tilde{c}^{-2}}\right)\\!-\\!\gamma_{E}\right]+\\!\tilde{n}\ln\frac{\gamma\left(\tilde{n}+\sqrt{\tilde{n}^{2}+\tilde{c}^{-2}}\right)}{\tilde{n}+\sqrt{\tilde{n}^{2}+\left(\gamma/\tilde{c}\right)^{2}}}$
(3)

with $\zeta(3)\approx 1.202$. This result is valid for
$\tan\theta<c\gamma/(2\pi s)$ covering the relevant angular range.

To find equilibrium $\tilde{n}$ for fixed reduced magnetic field strength
$\tilde{h}$, one has to minimize
$F_{TV}(\tilde{n})=E_{TV}(\tilde{n})-\tilde{n}\tilde{h}$. In the range of
large $\tilde{c}$ above the certain critical value the function
$F_{TV}(\tilde{n})$ has two minima and the equilibrium pancake density
experiences a jump between the two values $\tilde{n}_{\mathrm{dn}}$ and
$\tilde{n}_{\mathrm{up}}$ at the critical value of field . These values are
converted to the magnetic inductions using Eq. (1). Figure 6 shows example of
the phase diagram obtained assuming $\lambda_{ab}=200$ nm, $s=1.2$ nm and
$\gamma=7$. We also assumed $\beta=2\sqrt{3}$ in Eq. (1).

##  Appendix B Imaging formation of the vortex domains

To track the initial stage of the domain nucleation we imaged the flux
structure during gradual reduction of the field after initially field cooling
the sample to T = 60K with $H_{//}$ = 1388 Oe tilted by $1.6^{o}$ (see Fig.
7). Before reducing $H_{//}$, the $B_{z}$ distribution (not shown) is
homogeneous over the entire imaged area with very faint darker contrast near
the edges of the sample due to weak Meissner flux expulsion. Decreasing field
to $H_{//}$ = 832 Oe results in the appearance of a dark edge contrast due to
the stray fields coming from the trapped positive $B_{z}$. Simultaneously,
stripes of alternating bright and dark contrast (modulated $B_{z}$) emerge
along $H_{//}$ near the left and right ends of the sample (Fig. 7a). Further
decreasing the magnetic field, extends the stripes over the entire length of
the sample (Fig. 7b). At even smaller $H_{//}$=129 Oe, vortices turn into the
plane near the perimeter of the sample (darker periphery) and a region of
trapped $B_{z}$ with stripe domains forms inside the crystal as depicted by
the bright contrast in Fig. 7c. This configuration remains after reducing the
field to zero (Fig. 7d).

References

[1s] L. L. Daemen, L. J. Campbell, A. Yu. Simonov, and V. G. Kogan Phys. Rev.
Lett. 70, 2948 (1993); A. Sudbø, E. H. Brandt, and D. A. Huse Phys. Rev. Lett.
71, 1451 (1993); E. Sardella, Physica C 257, 231 (1997).

[2s] A. E. Koshelev, Phys. Rev. B 71, 174507 (2005).

[3s] S. N. Gordeev, A. A. Zhukov, P. A. J. de Groot, A. G. M. Jansen, R.
Gagnon, and L. Taillefer, Phys. Rev. Lett. 85, 4594 (2000).

Figure 7: Nucleation of domains during reduction of $H_{//}$ at $1.60^{o}$.
The sample has been cooled in 1388 Oe to 60K. (a)-(d) $H_{//}$=832, 555, 129,
and 0 Oe. The zig-zag features in (d) are magnetic domains in the MO
indicator.


# Thermostatic control for demand response using distributed averaging and
deep neural networks

Kshitij Singh<EMAIL_ADDRESS>Pratik Bajaria
<EMAIL_ADDRESS>Student/Research Scholar, Veermata Jijabai
Technological Institute, Mumbai - 400 019, India

###### Abstract

Smart buildings are the need of the day with increasing demand-supply ratios
and deficiency to generate considerably. In any modern non-industrial
infrastructure, these demands mainly comprise of thermostatically controlled
loads (TCLs), which can be manoeuvred. TCL loads like air-conditioner, heater,
refrigerator, are ubiquitous, and their operating times can be controlled to
achieve desired aggregate power. This power aggregation, in turn, helps
achieve load management targets and thereby serve as ancillary service (AS) to
the power grid. In this work, a distributed averaging protocol is used to
achieve the desired power aggregate set by the utility using steady-state
desynchronization. The results are verified using a computer program for a
homogeneous and heterogeneous population of TCLs. Further, load following
scenario has been implemented using the utility as a reference. Apart from
providing a significant AS to the power grid, the proposed idea also helps
reduce the amplitude of power system oscillations imparted by the TCLs.
Hardware-based results are obtained to verify its implementation feasibility
in real-time. Additionally, we extend this idea to data-driven paradigm and
provide comparisons therein.

††preprint: AIP/123-QED

##  I Introduction

Sustainability is the process of maintaining change in a balanced environment,
in which the exploitation of resources, the direction of investments, the
orientation of technological development, and institutional changes are all in
harmony and enhance both current and future potential to meet human needs and
aspirations. With the advent of civilization, to serve electrical and
transportation needs, fossil fuel reserves have been over-utilized. Excessive
use of fossil fuels and increasing human settlements, has brought imbalance to
the environment. Thus, researchers have moved towards building settlements
that can serve individual requirements, thereby mean no harm to the
environment utilizing least natural resources; hence are sustainable. To
support the central power grid, ancillary services (ASs) - a bidirectional
power flow architecture is employed.

ASs support the transmission of electric power from seller to purchaser, given
the obligation of control areas and transmitting utilities within the control
area, maintaining reliable operation of the interconnected transmission system
Hirst and Kirby (1997). Load following, reactive power voltage regulation,
system protection services, loss compensation, system control, load dispatch
services, and energy imbalance services are some ASs. Load following is a
service program that monitors real-time consumption about scheduled loads to
consume what is estimated. In contrast, load dispatch services take real-time
action in case of any discrepancies. System protective services are energy
reserve that compensates for any discrepancies in the scheduled loads and
dispatches where the loss compensation scheme guides its use. The energy
imbalance programs take all these actions using system control service.
Reactive power voltage regulation refers to services that mitigate
perturbation in the voltage online. As observed, these services either fall
into demand-side or supply-side management, but with increasing demand to
supply ratio, it is advisable to optimize loads to match available resources.

ASs have gained attention with increasing penetration of renewable energy
sources (RESs) in the course of replacing fossil fuels. Although renewable
energy sources (RESs) are abundantly available in nature, they suffer from
intermittencies. These intermittencies are significant enough to render a
power grid unstable and hence require additional regulating circuits to avoid
direct injection Malik  _et al._ (2017); Yan  _et al._ (2017); Kara, Bergés,
and Hug (2015); Meyn  _et al._ (2015). One of the ways of injection is to
store it in a battery and then provide it to the grid using a battery
management system (BMS) with a stable discharge rate. Although reliable, this
strategy bears a high capital expenditure for installation and hence is
commercially less feasible Siano (2014); Pavić, Capuder, and Kuzle (2016);
Teng  _et al._ (2017); Malik and Ravishankar (2018). Therefore, research has
moved in the direction of devising methods to regulate the loads with techno-
commercial feasibility. Demand Response (DR) programs have been developed to
ally with several handshaking criteria and thereby maintain power grid
stability.

DR programs allow consumers to adjust their electricity consumption in
response to energy prices or incentive payments. ASs form a major part of
these DR programs, and optimal use of these services can come to an advantage
in terms of adjusting loads as per utility demands. As per a survey in U.S.
(2012) (US), cooling and heating loads account for the largest annual
residential electricity consumption. And the scenario remains the same
globally, thereby requiring optimization of these loads to attain required
power aggregation. Fundamentally, all heating/cooling loads (or TCLs) comprise
a temperature transducer, a thermostat, and other assemblies like compressor,
and condenser. It can be noted, though, that TCLs follow a scheme whereby they
switch between ON and OFF states to attain required setpoint temperatures. ON
states consume full capacity power of the TCL, and OFF states imply dormancy.
Thus, adjusting these operation states can help achieve the required power
aggregate and thus provide load-following ancillary service.

In literature Lakshmanan  _et al._ (2016a, b), field experiments are conducted
using domestic household refrigerators to quantify the flexibility of TCLs by
advancing or delaying operation times, thereby achieve required aggregate
power. To avoid synchronization, parameter randomization of TCLs by
compartmentalization of compressor cycles is also an effective scheme Malik
and Ravishankar (2018). An automatic generation control signal (AGC - a
demand-supply balancing technique) tracking Ma  _et al._ (2017) and a two-
level scheduling protocol Zhou  _et al._ (2017) achieves flexible TCL
operating cycles in an intra-day electricity market making it a valuable asset
to DR programs. However, with increasing heterogeneity in system parameters,
the performance worsens in all the schemes mentioned above. On the other hand,
a multi-objective model predictive control (MPC) can also be employed for
residential heating (or using heat as an energy carrier) with heat pumps
Baeten, Rogiers, and Helsen (2017); Hu  _et al._ (2017). Again, the number of
parameters used makes it difficult to be implemented for a large population of
TCLs.

Aggregated models for TCLs can carry information about the overall power
consumption of TCLs and can be framed using discrete state-space dynamics. For
instance, a centralized control scheme on an aggregated model can provide
effective frequency regulation Hu  _et al._ (2017); Wu  _et al._ (2018). A
Markovian probabilistic framework Zhao  _et al._ (2018) or a load shedding
strategy Luo  _et al._ (2017) can also be used to achieve similar objectives.
Concepts from optimization theory related to the tracking device, convex
polytopes, and binary multi-swarm, particle optimization, are also used to
control the desynchronization of TCLs Iacovella  _et al._ (2017); Zhao  _et
al._ (2017); Ma  _et al._ (2018). Batch reinforcement learning based on Monte
Carlo methods is used to allow TCLs to participate in day-ahead markets,
thereby achieve DR objectives Ruelens  _et al._ (2017). A nano-grid solution
to store thermal energy in TCLs for later use has also been experimented
Burmester, Rayudu, and Seah (2017). This approach avoids wastage by regulating
maximum power point tracking (MPPT) of a photovoltaic panel. In literature
Ghanavati and Chakravarthy (2018), a partial differential equation (PDE) based
Fokker-Plank model for demand-side management is also proposed, and as known,
requires highly sophisticated tools for solving. Additionally, an event-
triggered control of TCLs on a hierarchical architecture is implemented over a
network for optimization of communication cost Meng  _et al._ (2017). All of
the methods mentioned above are probabilistic, discretized, and involve
considerable computation cost.

Complexity of population of TCLs increases with the number of TCLs, and with
the reference aggregate power set by the utility. The undesired
synchronization of TCL duty cycles by the control policies limits the time
scale and capacity of ASs that a population can provide. By sequentially
changing the set point temperature, safe protocols have been developed to
minimize the unwanted power oscillations Sinitsyn, Kundu, and Backhaus (2013);
Mehta  _et al._ (2014). Other ways of providing power regulation include
variable deadband strategy Sinitsyn, Kundu, and Backhaus (2013), injecting a
delay of $M$-minutes Sinitsyn, Kundu, and Backhaus (2013), and deploying state
stacking technique based on priority measures Hao  _et al._ (2015); Vrettos,
Koch, and Andersson (2012). These measures are defined as the distance from
the switching boundaries to reduction in aggregate power. TCLs in the ON stack
closer to switching boundary turn ON first, for an increase in aggregate
power; TCLs in the OFF stack closer to the switching boundary are turned OFF
first and so on. On similar lines, a phase response curve (PRC) based
technique effectively helps to manoeuvre an ensemble of TCLs by using a
control input to modulate duty cycles as well as induce delay/advance and
thereby changing the set point temperature Bomela, Zlotnik, and Li Jr (2018).
The switching dynamics of TCLs are similar to oscillators, and the phase-
dependent model has been realized to achieve the desired objective. A
continuum model for studying the synchronization dynamics of TCLs to
understand and mitigate potential risks of a decentralized response provider
offers deeper insights to TCL dynamics Webborn (2018). All these methods are
effective but involve high computation costs and challenges for practical
implementation.

This study extends the Kuramoto phase oscillator model Bajaria and Singh
(2019) used to control the phase delay/advance of the operating cycles of the
TCLs which helps to achieve permissible load requirements set by the utility.
We provide an extensive study on the possibilities of practical implementation
of the model discussed in the previous work Bajaria and Singh (2019). Based on
our findings, we propose a distributed averaging protocol to achieve similar
(and in a better way) objectives obtained using Kuramoto model Bajaria and
Singh (2019). Hardware-based outcomes, in conjunction with computer
simulations, provide validity to the proposed idea.

Additionally, we extend the idea of phase desynchronization using machine
learning (ML) fundamentals to control overall aggregate power consumption of
TCLs. The major concern being reduction in computational cost of deterministic
approaches to data-driven ones. Some of the articles from literature justify
using ML for control of TCLs Kara (2014); Kara  _et al._ (2014); Costanzo  _et
al._ (2016); Ruelens  _et al._ (2019); Almalaq and Edwards (2017); Kazmi  _et
al._ (2019). In particular, our views are concurrent with inclusions of deep
neural network (DNN) to learn the right phase thereby achieve required power
aggregate Ruelens  _et al._ (2019). We exploit the characteristics of delay
calculator (explained later in the manuscript) and extend it to DNNs to reduce
computations every time the network architecture/parameters gets perturbed.

The paper is organized as follows. In section II, we lay down fundamentals to
develop a platform to propose the mentioned idea. The motivation behind the
proposed paper has been discussed briefly in section III. A novel model for
the control of TCLs namely a distributed averaging protocol, is discussed in
section IV. In section V, we provide software and hardware outcomes using LEDs
as acting TCLs, which simulate similar characteristics, thereby validating the
proposed idea. Section VI presents a competitive benchmark, where we compare
the proposed model with those available in the literature and provide its
competencies using peer-to-peer reviews. In the last section, we provide a
deeper insight to the characteristics of delay calculator and resolves its
“curse of dimensionality" problem using a DNN.

##  II Prerequisites

###  II.1 Behaviour of TCLs

A TCL switches ON until it reaches the upper limit of the set threshold
temperature, post which it turns OFF until reaching the lower limit of the
threshold. Usually, TCLs employ a thermostat that measures the surrounding
temperature and thereby generates switching signals to regulate the TCL.
Typically, a hybrid model is employed to express the dynamics of TCL
mathematically. Incorporating the thermal resistance, capacitance, energy
transfer rate, and switching signal, the hybrid model describes the
temperature dynamics. It can be stated as follows,

---
(a)
(b)
Figure 1: Hybrid model. $\left(a\right)$ Temperature variation of a single TCL
with respect to time. $\left(b\right)$ Phase portrait of a TCL system
dynamics.

$\displaystyle\dot{T}(t)$ $\displaystyle=-\frac{1}{RC}[T(t)-T_{a}+s(t)PR],$
(1) $\displaystyle s(t)$ $\displaystyle=\begin{cases}0&\text{if}\
T(t)<T_{min},\\\ 1&\text{if}\ T(t)>T_{max},\\\
s(t)&\text{otherwise},\end{cases}$

Table 1: List of parameters Parameter | Meaning | Value | Unit
---|---|---|---
$T_{a}$ | Ambient temperature | $32$ |  ∘C
$\delta$ | Thermostat deadband | $1$ |  ∘C
R | Thermal resistance | $2$ |  ∘C/kW
C | Thermal capacitance | $10$ | kWh/∘C
P | Energy transfer rate | $14$ | kW
K | Coupling strength | $0.267$ | Mhz

where $T_{a}$ and $P$ are ambient temperature and average energy transfer rate
of a TCL respectively. The minimum and maximum temperature of TCLs are
$T_{min}=\left(T_{s}-\delta/2\right)$ and
$T_{max}=\left(T_{s}+\delta/2\right)$ respectively, $T_{s}$ and $\delta$ being
the thermostat set point temperature and deadband respectively. $s(t)$ is the
switching signal of the TCLs. FIG. 1$\left(a\right)$ depicts the temperature
dynamics where the TCLs oscillates between $19.5^{\circ}C$ and $20.5^{\circ}C$
while FIG. 1$\left(b\right)$ is phase portrait of the TCL obtained using
hybrid model (1) and data from TABLE LABEL:tab1.

Thus, for $N$ TCLs, the aggregate power consumption is given by,

$\displaystyle P_{agg}=\sum^{N}_{j=1}P\eta_{j}s_{j}(t),$ (2)

where $\eta_{j}$ is the coefficient of performance of $j$th TCL. For
simplicity of calculation, we consider $\eta_{j}=\eta=1$ for all TCLs.

######  Remark 1.

All the parameters used henceforth are taken from TABLE LABEL:tab1 unless
otherwise specified.

---
Figure 2: (From top to bottom) Pendulums coupled by a spring, oscillating in
two equilibrium modes; Oscillations in in-phase mode (Top), Oscillations in
anti-phase mode (Bottom).

##  III Motivation

###  III.1 Coupled Pendulums

To further understand the effect of TCL switching oscillations on aggregate
power consumption, we correlate the periodic behavior of TCL switching with a
setup of spring coupled pendulums. It is well known that the Kuramoto model
can be well understood using a system of coupled oscillators, and various
modes of oscillations can be analyzed. Consider two pendulums coupled with a
spring, as shown in FIG. 2, where $m_{i}$ are masses, $\phi_{1}$ and
$\phi_{2}$ the phases and $K$ is the spring constant. For the sake of
simplicity, assume both the pendulums to be identical in terms of system
parameters and dissipative losses to be zero. As known, the phase of coupled
pendulums exhibit ‘in-phase’ (phase difference $=0$ radians), ‘anti-phase’
(phase difference $=\pi$ radians) and ‘quasi in-phase’ ($0<$ phase difference
$<\pi$ radians) type oscillatory behavior Strogatz (2008). Let all the phases
measured in the anti-clockwise direction be positive and assume a constant
external force that drives these oscillations. These external forces can be
impulsive or persistent in time. It can be seen analogous to the addition of
an external delay to maintain required steady-state behavior (see FIG. 2).

Now, all external forces (delay $=0$ radians) that keep aggregate switching
amplitudes (a mathematical sum of phases) to its maximum can be considered to
be oscillating in ‘in-phase’ modes. All forces (delay $=\pi$ radians) that
keep aggregate switching amplitudes to its minimum can be considered to be
oscillating in ‘anti-phase’ modes. All the other states oscillate between
these two extrema can be termed as ‘quasi in-phase’ modes. We portray ‘in-
phase’ and ‘anti-phase’ oscillation modes pictorially, as shown in FIG. 2.
Further, using the concept of clustering, we show how these modes come to
benefit in achieving the required aggregate power.

###  III.2 Boolean phase oscillator model for TCLs

In the previous work Bajaria and Singh (2019), the author makes use of the
Boolean phase oscillator model to incorporate the system dynamics and thereby
control the aggregate power consumed ($P_{agg}$). Visualizing the TCLs as
oscillators Bajaria and Singh (2019), the author employs the Kuramoto
framework to desynchronize TCLs. Beginning with the conventional Kuramoto
equation Bajaria and Singh (2019), to accommodate the digital switching
signals, uses the Heaviside function, thereby imparting the Kuramoto framework
the capability to couple digital signals. The standard Kuramoto equation can
be expressed as follows,

$\displaystyle\dot{\phi}_{i}$
$\displaystyle=\omega_{i}+\frac{\epsilon}{N}\sum^{N}_{j=1,j\neq
i}\left(sin(\phi_{j}-\phi_{i})\right),$ (3)

where, $\omega_{i}$, $\phi_{i}$, $\epsilon$ and $N$ are natural frequency,
phase of $i$th oscillator, gain and number of oscillators respectively,
$\forall i\in\mathbb{Z^{+}}$.

---
(a)
(b)
(c)
Figure 3: Boolean phase oscillator model. (a) Angular separation for $N=4$
TCLs. (b) TCL temperature and switching signal dynamics. (c) $P_{agg}$ for
$N=4$ TCLs at a given time.

The inherent coupling dynamics of (3) enables the oscillators to establish
communication between all TCLs. As mentioned earlier, the author using
Heaviside function modifies (3) Bajaria and Singh (2019) and the resultant
equation can be given as,

$\displaystyle\dot{\phi}_{i}=\omega_{i}+K\sum^{N}_{j=1,j\neq
i}\left|\Theta\left[sin(\phi_{j})\right]-\Theta\left[sin(\phi_{i}+\alpha_{ij})\right]\right|,$
(4)

where $\Theta$ is Heaviside function given by,

$\displaystyle\Theta[\sin{(.)}]=\begin{cases}\text{$0$}&\quad\text{$if$
$\sin{(.)}<0$},\\\ \text{$1$}&\quad\text{$if$ $\sin{(.)}\geq 0$},\end{cases}$
(5)

and $K$ is the coupling between oscillators. The phase separation
$\alpha_{ij}$ in (4) induces a phase separation between the $i_{th}$ and
$j_{th}$ oscillator, or in this case TCLs as shown in FIG. 3$\left(a\right)$.
Following equations characterize temperature dynamics and stabilization of
aggregate power consumption as depicted in FIG. 3$\left(b\right)$ and
3$\left(c\right)$ respectively,

$\begin{split}\dot{T}_{i}&=-\frac{1}{[RC]_{i}}\left[T_{i}(t)-T_{a}+s_{i}(t)[PR]_{i}\right],\\\
\dot{\phi}_{i}&=\omega_{i}+K\sum^{N}_{j=1,j\neq
i}\left|\Theta\left[sin(\phi_{j})\right]-\Theta\left[sin(\phi_{i}+\alpha_{ij})\right]\right|,\\\
s_{i}(t)&=\Theta\left[sin(\phi_{i})-s_{i,0}\right],\\\
s_{i,0}&=sin\left[\frac{\pi-T_{ON,i}}{2}\right],\\\ \end{split}$ (6)

Where $T_{ON,i}$ is ON time of the $i$th switching signal (mapped to its
equivalent phase in radians), $\omega_{i}$ is the frequency of switching, $K$
is the coupling constant. A bias, ($s_{i,0}$), helps the generated signal map
the duty cycle in the Kuramoto framework.

---
Figure 4: Kuramoto based model. Aggregate power consumption of $N=100$
heterogeneous TCLs.

In the previous work Bajaria and Singh (2019), the idea is further supported
by simulations results. It was thus conclusively stated that the power
oscillations and aggregation are minimum when the phase delay ($\alpha_{ij}$)
is $2\pi/N$.

To delineate this theory, the author Bajaria and Singh (2019) takes a case of
$N=4$ homogeneous TCLs, with the natural frequencies as $\omega_{i}=0.55$Hz.
Furthermore, the TCLs desynchronize at their respective phase differences with
their common frequency being $\omega_{FFT}=0.298$Hz, where $\omega_{FFT}$ is
calculated using fast Fourier transform $\left(FFT\right)$ of the switching
signal. Finally, the Kuramoto model desynchronizes TCLs at a common frequency
lesser than their mean frequency. FIG. 3$\left(c\right)$ shows the power vs.
time plot for $N=4$ TCLs.

The author then takes a case of $N=100$ heterogeneous TCLs with their natural
frequencies in the range of $\pm 5\%$ of its nominal value. In FIG. 4, it is
observed that the load fluctuations are minimized to within $\pm 2\%$ after
applying Kuramoto based model.

######  Remark 2.

Frequency calculated in the previous work Bajaria and Singh (2019) are scaled
down by hours, i.e., $\omega_{i}$ is calculated as $1/T$ where, $T$ is in
hours. However, in this paper we consider time in seconds.

###  III.3 Shortcomings of Boolean phase oscillator model

Although the Boolean phase oscillator model Bajaria and Singh (2019) provides
motivating results, some shortcomings in the model proposed can be listed as
follows:

  1. 1.

It must be noted that the natural frequency $\omega_{i}$ used in (6) can be
easily calculated using mathematical manipulations as follows,

$\begin{split}\omega_{i}=&\ \omega_{FFT}-K\sum^{N}_{j=1,j\neq
i}\left|\Theta\left[sin(\phi_{j})\right]-\Theta\left[sin(\phi_{i}+\alpha_{ij})\right]\right|,\end{split}$
(7)

However, it is observed that as the number of TCLs increases, the value of
$\omega_{i}$ needs correction. For instance, $N=16$ requires $\omega_{i}=8.81$
Hz and $N=100$ needs $\omega_{i}=62.18$ Hz, to achieve $\omega_{FFT}=0.298$
Hz. Moreover, as the number of TCLs increases, $\omega_{i}$ computed from (7)
requires an additional correction factor $e$ to give desired results. Failing
to find the correct $e$ leads to erroneous deadband values. Also, since the
system under consideration consists of mechanical devices such as compressors
and switches, such high frequencies are undesirable. Additionally, the TCLs
are frequency sensitive devices, and drastic changes in the frequency must be
avoided to prolong the product life Mehta  _et al._ (2014). In a practical
scenario, where system parameters (e.g., $N,\alpha_{ij}$) change dynamically,
computation of $\omega_{i}$ in real-time would seem difficult. Thus, rendering
the Boolean phase oscillator model difficult to be practically realized.

  2. 2.

$K$ is a parameter that is computed using dynamical equations governing the
system parameters Dörfler, Chertkov, and Bullo (2013). In case of TCLs, the
value of $K$ can be calculated from the induced time delays in a programmable
frame Bajaria and Singh (2019). However, the inter-dependence of $K$ and
$\omega_{i}$ in (6) adds complexity to its accurate computation.

---
(a)
(b)
Figure 5: Frequency dynamics. (a) FFT of switching signal at $t$ = $t_{0}$.
(b) FFT of switching signal at $t$ = $t_{f}$.

##  IV Distributed Averaging Model for the TCLs

Consider an undirected graph $G=\\{V,E\\}$ with set of nodes $N$ in $V$ and
edge set $E$. Each edge of the graph $E\in G$ is an unordered pair of distinct
nodes. A real scalar quantity $x_{i}$(t) is associated with node $i$ at any
given time $t$. The distributed averaging consensus algorithm computes the
average $\frac{1}{N}\sum^{N}_{i=1}x_{i}(0)$ iteratively thus, facilitating a
local communication link between the nodes. The nodes are updated continuously
based on their states as well as the states of all other neighboring nodes.

A widely accepted and used linear iterative algorithm can be stated as
follows,

$\displaystyle\dot{x}_{i}(t)$ $\displaystyle=\sum^{N}_{j=1,j\neq
i}W_{ij}\Big{(}x_{j}(t)-x_{i}(t)\Big{)},$ (8)

for $i\in\\{1,\ldots,n\\}$ and $t\in\mathbb{R^{+}}$ and $W_{ij}$ being the
weight associated with each edge $E_{ij}$. Since the nodes, in the presented
case, are undirected, the weights associated with them are symmetric i.e.,
$W_{ij}=W_{ji}$.

$\displaystyle W=w_{ij}(\mathbbm{1}_{N\times N}-diag(\mathbbm{1}_{N})),$ (9)

where, $w_{ij}$ is the weight of edge $E_{ij}$ and $diag(\mathbbm{1}_{N})$ is
diagonal matrix with $\mathbbm{1}_{N},\mathbbm{1}_{N\times N}$ denoting a
vector and matrix of ones, respectively. Thus, the matrix $W$ is the
communication matrix that describes the strength of the link between two
nodes. Analogous to frequency synchronization seen in the Kuramoto framework
(6) we extend the distributed theory to the system of TCLs. The nodes are
similarized to TCLs and the scalar quantity to the frequency of switching
signal as follows,

$\begin{split}\dot{f_{i}}&=W\sum^{N}_{j=1,j\neq i}(f_{j}-f_{i}),\\\
\end{split}$ (10)

where $W$ is the tuning parameter of synchronization, left to the control
designer. (10) synchronizes the frequencies randomly chosen from distribution
to their mean value, as shown in FIG. 5. Thus, (6) can be re-framed as
follows,

$\begin{split}\dot{T}_{i}&=-\frac{1}{[RC]_{i}}\left[T_{i}(t)-T_{a}+s_{i}(t)[PR]_{i}\right],\\\
\dot{f_{i}}&=W\sum^{N}_{j=1,j\neq i}(f_{j}-f_{i}),\\\
s_{i}(t)&=\Theta\left[sin(2\pi f_{i}t+\alpha_{ij})-s_{i,0}\right],\\\
s_{i,0}&=sin\left[\frac{\pi-T_{ON,i}}{2}\right].\\\ \end{split}$ (11)

It should be noted that (11) bifurcates the onus of controlling the frequency
and angular separation.

---
Figure 6: $P_{agg}$ vs time showing the aggregate power consumed at a given
time for $N=100$ heterogeneous TCLs.

For the matter of comparison, a similar ensemble of heterogenous TCLs, as in
FIG. 4 (natural frequencies in the same range as of $\pm 5\%$ of
$\omega_{i}=0.271$Hz), are taken. Additionally, the value of weight($W$)
chosen as $0.06$. By using (11), the system desynchronizes, and FIG. 6 depicts
the results obtained. Similar results can be noted by comparing the two plots,
however, as opposed to Kuramoto based model in FIG. 6, high frequency
switching and computational hindrance have been avoided.

##  V Software and Hardware Implementation

This section presents the application of the distributed averaging protocol in
the form of software simulations and later describes its practical
implementation.

Table 2: List of Parameters (at set point temperature = $24^{\circ}C$) Parameter | Meaning | Value | Unit
---|---|---|---
$f$ | frequency of TCL | 0.0027 | Hz
$\delta$ | thermostat deadband | 2 | ${}^{\circ}C$
$d$ | duty cycle of TCL | 0.5083 | -
Table 3: List of Parameters (at set point temperature = $27^{\circ}C$) Parameter | Meaning | Value | Unit
---|---|---|---
$f$ | frequency of TCL | 0.0031 | Hz
$\delta$ | thermostat deadband | 3 | ${}^{\circ}C$
$d$ | duty cycle of TCL | 0.452 | -

###  V.1 Case Study

To make the two cases put forth as realistic as possible, real time data of a
TCL has been used to simulate the cases. For this, a TCL of a certain make was
observed and the data listed in the TABLE 2 and TABLE 3 were collected. This
TCL had a power rating of $1.66$kW and capacity of $1.5$tons. The values
$R=9.8^{\circ}$C/kW and $C=0.0746$kW/∘C were computed using (1), TABLE 2 and
TABLE 3.

  1. 1.

Case 1: $N=1000$

For $N=1000$, heterogeneous TCLs are defined having a duty cycle in the range
of $0.422$ to $0.482$, $\delta$ as $3.0^{\circ}C$ and $f$ between $0.0029$Hz
and $0.0033$Hz respectively. The TCLs have been set to work at a set point
temperature of $27^{\circ}C$.

Following the distributed averaging model, TCLs settle at their respective
phase difference, i.e., $\alpha_{ij}=\alpha=2\pi/N$, which in this case is
$\alpha=\pi/500$. As depicted in FIG. 7, the power fluctuation for $N=1000$
heterogeneous TCLs dampen over time and move towards steady state. The
aggregate power at which the TCLs settle is $750.4472$kW and fluctuations in
the power are $\pm 2.67\%$.

Figure 7: $P_{agg}$ vs time showing the aggregate power consumed at a given
time for $N=1000$ heterogeneous TCLs.
Figure 8: $P_{agg}$ vs time showing the aggregate power consumed at a given
time for $N=10000$ heterogeneous TCLs.

  2. 2.

Case 2: $N=10000$

For $N=10000$, heterogeneous TCLs are defined having a duty cycle in the range
of $0.4812$ to $0.5354$, $\delta$ as $2.0^{\circ}C$ and $f$ between $0.0026$Hz
and $0.0036$Hz respectively. The TCLs have been set to work at a set point
temperature of $24^{\circ}C$.

Following the distributed averaging model, TCLs settle at their respective
phase difference, i.e., $\alpha_{i}=\alpha=2\pi/N$, which in this case is
$\alpha=\pi/5000$. As depicted in FIG. 8, the power fluctuation for $N=10000$
heterogeneous TCLs dampen over time and move towards steady state. The
aggregate power at which the TCLs settle is $8.4453$MW, and fluctuations in
power are $\pm 2.5\%$. In FIG. 9, phase portraits of $100$ TCLs chosen
randomly from the population are shown. It is clear that the portrait
dissimilarity is minimal, although setpoint temperatures remain unaffected.

Figure 9: Phase portrait for $N=10000$ TCLs.

(a)

(b)

Figure 10: (a) $P_{norm}(\%)$ vs time showing variation of $P_{norm}(\%)$ with
$\alpha$. (b) $P_{agg}$ vs time plot showing variation of $P_{agg}$ with
$\alpha$.

###  V.2 Load Following

Next, the authors show how a load-following case can be emulated, providing a
valuable AS to the utility. In the previous work Bajaria and Singh (2019), a
function mapping RMS (root mean square) aggregate power of a population of
TCLs with the effective delay induced was proposed. Using the same idea, we
choose different utility demands, and the corresponding delay of $\alpha_{ij}$
is calculated. The effective power so obtained is normalized using total RMS
power capacities; i.e.,

$\displaystyle
P_{norm}(\%)=\frac{(P_{rms,agg}-P_{rms,\alpha})}{P_{rms,agg}}\times 100,$ (12)

where $P_{rms,agg}$ is the RMS of aggregate power, $P_{rms,\alpha}$ is the RMS
of aggregate power at steady-state at $\alpha_{ij}=\alpha$ and $P_{norm}$ is
the normalized power of all the TCLs along a given period respectively. These
results can be seen from FIG. 10$\left(a\right)$ where utility demands for
random reduction/increase in the total allowable aggregate power and the
corresponding change in the delay induced causes controlled desynchronization
and hence the required load following scenario is achieved. As seen from FIG.
10$\left(b\right)$($N=100,P=14kW$) it can be inferred that the resultant
aggregate power is a significant AS to the centralized power grid.

(a)

(b)

(c)

Figure 11: (a) $P_{agg}$ vs time plot for random and desynchronized state for
$N=100$. (b) $P_{red}(\%)$ vs time plot showing the reduction in fluctuations
for $N=100$. (c)$P_{red}(\%)$ vs time plot showing the reduction in
fluctuations for $N=10000$.

###  V.3 Effects in a Power Grid

It should be noted that apart from the control of power aggregation, the
proposed algorithm benefits the grid by reducing the power system
fluctuations. To compare/visualize the same, we take a small population of
$N=100$ TCLs. To exhibit random behavior i.e., depicting as-is characteristics
without applying control), the frequency and phases of TCLs are chosen
randomly from a distribution. Next, the same set of heterogeneous TCLs is
acted upon by the distributed averaging protocol achieving effective
desynchronization. As shown in FIG. 11$\left(a\right)$, the system without
application of distributed averaging strategy induce higher oscillations in
the grid as compared to after application of the proposed scheme. In FIG.
11$\left(b\right)$ the quantum of reduction is observed to be $\sim 40\%$
($560kW$) of random case, showing a significant reduction($P_{red}(\%)$)
where,

$\displaystyle
P_{red}(\%)=\frac{(P_{random}-P_{desynchronized})}{P_{random}}\times 100.$
(13)

Similarly, a case of $N=10000$ TCLs is also shown in FIG. 11$\left(c\right)$.
It can be seen a reduction of $\sim 2\%$ can be noticed which translates to
$2.8MW$. Hence, apart from providing load following objectives, the proposed
scheme can also help reduce the amplitude of power system oscillations.

###  V.4 Hardware Implementation

To validate the presented theory, the authors have implemented the model on a
hardware system. The authors use the circuit diagram as shown in FIG. 12,
which practically looks as seen in FIG. 13.

Light-emitting diodes (LEDs) provide close resemblance to a switching device,
analogous to TCLs. Hence, for the demonstration purpose, LEDs have been used
in place of an actual TCLs. Since the number of LEDs was limited to four, a
micro-controller with necessary computational capabilities and decent clock
frequency was sufficient.

Figure 12: Circuit diagram used for practical implementation.

####  V.4.1 Hardware

The components used for the implementation include 4 LEDs, micro-controller, a
breadboard, 5V DC supply, wires for connection, and probes to display the
output on the oscilloscope. The digital storage oscilloscope (DSO) used has
the following specifications,

  * •

Product Range: TBS1000B-EDU

  * •

Scope Channel: $2$

  * •

Bandwidth : $70$MHz

  * •

Sampling Rate :$1$ GSPS

  * •

Display Memory Depth : $2.5$kpts

  * •

Scope Display Type: WVGA LCD Colour

  * •

Plug Type: EU, SWISS, UK

The micro-controller used is Arduino Uno with following specifications,

  * •

Micro-controller: ATmega$328$

  * •

Operating Voltage: $5$V

  * •

Digital I/O Pins: $14$ (of which $6$ are PWM pins)

  * •

Analog Input Pins: $6$

  * •

DC Current per I/O Pin: $40$ mA

  * •

DC Current for $3.3$V Pin: $50$ mA

  * •

Clock Speed: $16$ MHz

Figure 13: Practical setup used for implementation of distributed averaging
model.
Figure 14: Oscilloscope output for $N=4$ depicting random, synchronized and
desynchronized period of LEDs.

####  V.4.2 Implementation

The setup for the practical implementation of the distributed averaging model
is shown in the FIG. 13. Based on the distributed averaging model, as
described by (11), the logic is coded and burned on Arduino. The Arduino
generates switching signals for the LEDs. The output of the LEDs is observed
on the oscilloscope.

As it can be seen, the output shown in FIG. 14 is similar to that of FIG. 6.

In the beginning, i.e., between $t=0s$ to $t=6.25s$, FIG. 14, depicts the
condition where the LEDs emulate random behaviour of the TCLs where no control
is enforced.

Figure 15: Transition of LEDs from the random period to synchronized state.
Figure 16: Transition of LEDs from the synchronized period to desynchronized
state.

After $t=6.25s$, the phases of the LEDs get synchronized. Following (10), LEDs
oscillate with their mean frequency. During this period, the LEDs have no
phase difference between them. FIG. 15 shows the synchronized state of the
LEDs. The synchronized state of the TCLs is a highly undesirable condition in
a power system.

Further, succeeding $t=12s$, the LEDs induce phase difference of
$\alpha=\pi/2$ between them. As seen from FIG. 16, the output of the LEDs are
minimized. Additionally, the fluctuations are also minimized.

The results of the hardware implementation were in conjunction with the
results of computer simulations. Thus, the results of the proposed theory are
validated.

##  VI Competitive Benchmark

In this section, the authors compare the results obtained from the proposed
model with the models mentioned in literature. From FIG. 7, the $P_{agg}$ is
$750$kW ($45.18\%$ of maximum power consumption). The fluctuations are
quantified using root mean square error ($RMSE\%$) given as follows,

$\displaystyle RMSE\%$
$\displaystyle=\sqrt{\frac{1}{T}\frac{\int_{0}^{T}\left(P_{ref}(t)-P_{agg}(t)\right)^{2}dt}{P_{base}^{2}}}\times
100.$ (14)

For a heterogeneous ensemble of TCLs, the relative error computed is within
$2.49\%$, and the $RMSE\%$ is $4.8685\%$. Juxtaposing these results alongside
with the minimum variance law (MVC) Callaway (2009), the latter achieves
relative error lesser than $5\%$.

TABLE 4 shows the efficacy of a distributed averaging model against several
state-of-the-art methods. The relative error obtained in reference to PRC
based technique Bomela, Zlotnik, and Li Jr (2018) were $4\%$ to $8\%$. Whereas
compared to the statistical model presented in the literatureKoch, Mathieu,
and Callaway (2011), the $RMSE\%$ of the presented model is marginally higher.
As compared to the Boolean phase oscillator model Bajaria and Singh (2019),
the relative errors are similar. Thus, the distributed averaging model
achieves competent results while avoiding computational complexities.

Table 4: Comparison with state-of-the-art methods Parameter | Relative error$\left(\%\right)$ | RMSE $\left(\%\right)$
---|---|---
Callaway (2009) | $5\%$ | $1\%$
Vrettos, Koch, and Andersson (2012) | $30-50\%$ | $1.18-8\%$
Bomela, Zlotnik, and Li Jr (2018) | $4-8\%$ | $2.31-5.89\%$
Koch, Mathieu, and Callaway (2011) | - | $2.27\%$
Bajaria and Singh (2019) | $2\%$ | $6.3\%$
Distributed Averaging | $2.49\%$ | $4.8685\%$

##  VII Data-driven Approach to Phase Desynchronization of TCLs

As discussed in our previous work Bajaria and Singh (2019), delay calculator
(DC) is the underlying theory behind choosing the right phase value so as to
achieve controlled desynchronization. In this section, we further exploit the
benefits of DC by analysing it for higher dimensions. It can be noted though,
DC provides accurate results to deploy load following (refer section V) with a
sectional phase distribution. But, what if DC needs to be evaluated for all
the possible permutations and combinations? It is quite obvious that, as the
number of dimensions increase, the volume of data blow up exponentially making
it difficult for analysis, commonly referred to as “curse of dimensionality
(COD)" Bellman (1966). In order to avoid COD and provide a different
perspective to DC for controlled desynchronization of TCLs, we propose a newer
machine learning (ML) based model for the same. Out of the available sets of
ML techniques, we try implementing DNN for the discussed case. An obvious
choice is to use a multi-variate regressor, but to our analysis DNN proves to
be more accurate. Also, DNN must be a right choice taking into consideration
COD and making use of the available model/data for a simple application like
TCL. There is an immense amount of literature available on DNNs, in order to
understand the concepts and terminologies of DNNs (out of scope of this
context, but the reader can refer the citation) Goodfellow, Bengio, and
Courville (2016).

---
Figure 17: Variation of $P_{norm}$ vs $\alpha$ across $N$.

###  VII.1 Generating Dataset

In order to train the supervised learning model, for predicting the value of
$\alpha$, we generate a dataset using parameters mentioned in TABLE 2. For
emulating a practical environment, we choose the value of $P$ and dutycycle
from a random distribution. The number of TCLs vary from $N=10$ to $N=500$.
Additionally as the value of $P_{agg}$ varies from a maximum at $\alpha=0$
radians, to a minimum at $\alpha$ = $2\pi/N$ radians, we compute $P_{norm}$
using equation 14 for $100$ equally spaced values of $\alpha$ between $0$ and
$2\pi/N$ radians. As evident from FIG. 17, the value of $P_{norm}$ varies
between maximum and minimum.

This dataset encompasses the parameters of most commercial and household TCLs
and gives us a sufficiently large data to train our DNN for predicting values
of $\alpha$ for any $N$. The dataset contains $49490$ data points with two
independent variables $N$ and $P$, and one target variable $\alpha$.

###  VII.2 Data Preprocessing

The dataset generated using feature engineering generally do not have missing
value or $NaN$ and so goes for the data generated using DC theory. However, we
need to normalize the dataset prior to model building and training for
improved accuracy.

For normalization, we use minmax normalization to linearly transform the
dataset. Minmax transformation transforms each data point based on the minimum
and maximum value along a particular feature. It can be expressed as,

$\displaystyle X_{norm}=\frac{(X-X_{min})}{(X_{max}-X_{min})}$ (15)

Where $X_{norm}$, $X_{min}$, and $X_{max}$ are the normalized value of the
data point $X$, minimum and maximum value across the feature $X$,
respectively. To validate our model we segregate the dataset into training and
testing data in the ratio $7:3$.

###  VII.3 Model Training and Results

Next, we train the DNN using dataset generated in previous sections. In order
to optimize computational cost and standardization of code, we make use of
scikit-learn Pedregosa  _et al._ (2011) and keras Chollet  _et al._ (2015)
libraries in Python $3.7.3$. A simple $4$-layered DNN with $2$ dense hidden
layers comprising of $2$ and $4$ nodes, respectively has been implemented.
With an appropriate choice of activation functions, epochs and batch sizes a
RMSE of $3\%$ is attained. The model evaluation metrics are mentioned in TABLE
LABEL:tb:tab4 showing competitive results and accuracies.

Table 5: Model Evaluation Metric | Value
---|---
Root Mean Squared Error | $3.06\%$
Mean Squared Error | $0.09\%$
Mean Absolute Error |  $0.0214^{\circ}$ ($3.7{e^{-4}}$ radians)
---
Figure 18: Root Mean Squared Error(RMSE) variation with iterations.
---
Figure 19: Mean Squared Error(MSE) variation with iterations.
---
Figure 20: $\alpha_{true}$ vs $\alpha_{pred}$.
---
Figure 21: Plot of $\alpha_{pred}$ compared to $\alpha_{true}$.

The variation of RMSE, MSE with epochs are shown in FIG. 18, FIG. 19;
depicting DNN getting trained along the iterations. It can be observed nearly
after $~{}15$ iterations, the accuracy of the DNN starts to saturate.

We compare the values of $\alpha_{pred}$, $\alpha_{true}$ using a correlation
map as shown in FIG. 20. Additionally, as seen from FIG. 21 and TABLE
LABEL:tb:tab4 the value of $\alpha_{pred}$ is within $0.1^{\circ}$($1.7e^{-3}$
radians) of the true value showing the efficiency of predictions.

##  VIII Conclusion

The major highlight of this work is the application of a distributed averaging
protocol for the desynchronization of a population of TCLs. The system of TCLs
are visualized as oscillators, and their power fluctuations have been
minimized. A significant cornerstone of the proposed work is that user comfort
is unimpaired. The parameters governing the dynamics of the system are then
articulated as a set of differential equations.

Apart from power fluctuations, aggregation, and user comfort, another
breakthrough is the simplicity of the presented model and reduced computation.
Load following helps the utility by improving the usability of the generated
power, thereby supporting in load-frequency control and maintaining voltage
profile. A reduction in power fluctuation and aggregation would ease off the
demand pressure from the industries directly or indirectly associated with
power generation, distribution, protection, and maintenance. The results of
both simulations, as well as hardware, support the findings of the presented
works. The DNN model provides a convenient way to compute the value of the
delay further saving us computational cost and optimizing the system. The
results obtained are motivating and show great potential to be implemented as
a plug-n-play device and be made available in the market.

## References

  * Hirst and Kirby (1997) E. Hirst and B. Kirby,  _Creating competitive markets for ancillary services_ (Oak Ridge National Laboratory, 1997).
  * Malik  _et al._ (2017) M. Malik, M. Tabone, D. P. Chassin, E. C. Kara, R. V. Guha, and S. Kiliccote, “A common data architecture for energy data analytics,” in  _Smart Grid Communications (SmartGridComm), 2017 IEEE International Conference on_ (IEEE, 2017) pp. 417–422.
  * Yan  _et al._ (2017) J. Yan, Y. Zhai, P. Wijayatunga, A. M. Mohamed, and P. E. Campana, “Renewable energy integration with mini/micro-grids,” Applied Energy 201, 241–244 (2017).
  * Kara, Bergés, and Hug (2015) E. C. Kara, M. Bergés, and G. Hug, “Impact of disturbances on modeling of thermostatically controlled loads for demand response,” IEEE Transactions on Smart Grid 6, 2560–2568 (2015).
  * Meyn  _et al._ (2015) S. P. Meyn, P. Barooah, A. Busic, Y. Chen, J. Ehren, _et al._ , “Ancillary service to the grid using intelligent deferrable loads.” IEEE Trans. Automat. Contr. 60, 2847–2862 (2015).
  * Siano (2014) P. Siano, “Demand response and smart grids—a survey,” Renewable and sustainable energy reviews 30, 461–478 (2014).
  * Pavić, Capuder, and Kuzle (2016) I. Pavić, T. Capuder, and I. Kuzle, “Low carbon technologies as providers of operational flexibility in future power systems,” Applied Energy 168, 724–738 (2016).
  * Teng  _et al._ (2017) F. Teng, Y. Mu, H. Jia, J. Wu, P. Zeng, and G. Strbac, “Challenges on primary frequency control and potential solution from evs in the future gb electricity system,” Applied energy 194, 353–362 (2017).
  * Malik and Ravishankar (2018) A. Malik and J. Ravishankar, “A hybrid control approach for regulating frequency through demand response,” Applied Energy 210, 1347–1362 (2018).
  * (10) E. I. A. (US),  _Annual Energy Outlook 2012: With Projections to 2035_ (Government Printing Office, 2012).
  * Lakshmanan  _et al._ (2016a) V. Lakshmanan, M. Marinelli, J. Hu, and H. W. Bindner, “Provision of secondary frequency control via demand response activation on thermostatically controlled loads: Solutions and experiences from denmark,” Applied energy 173, 470–480 (2016a).
  * Lakshmanan  _et al._ (2016b) V. Lakshmanan, M. Marinelli, A. M. Kosek, P. B. Nørgård, and H. W. Bindner, “Impact of thermostatically controlled loads’ demand response activation on aggregated power: A field experiment,” Energy 94, 705–714 (2016b).
  * Ma  _et al._ (2017) K. Ma, C. Yuan, J. Yang, Z. Liu, and X. Guan, “Switched control strategies of aggregated commercial hvac systems for demand response in smart grids,” Energies 10, 953 (2017).
  * Zhou  _et al._ (2017) Y. Zhou, C. Wang, J. Wu, J. Wang, M. Cheng, and G. Li, “Optimal scheduling of aggregated thermostatically controlled loads with renewable generation in the intraday electricity market,” Applied energy 188, 456–465 (2017).
  * Baeten, Rogiers, and Helsen (2017) B. Baeten, F. Rogiers, and L. Helsen, “Reduction of heat pump induced peak electricity use and required generation capacity through thermal energy storage and demand response,” Applied Energy 195, 184–195 (2017).
  * Hu  _et al._ (2017) J. Hu, J. Cao, M. Z. Chen, J. Yu, J. Yao, S. Yang, and T. Yong, “Load following of multiple heterogeneous tcl aggregators by centralized control,” IEEE Transactions on Power Systems 32, 3157–3167 (2017).
  * Wu  _et al._ (2018) X. Wu, J. He, Y. Xu, J. Lu, N. Lu, and X. Wang, “Hierarchical control of residential hvac units for primary frequency regulation,” IEEE Transactions on Smart Grid 9, 3844–3856 (2018).
  * Zhao  _et al._ (2018) H. Zhao, Q. Wu, S. Huang, H. Zhang, Y. Liu, and Y. Xue, “Hierarchical control of thermostatically controlled loads for primary frequency support,” IEEE Transactions on Smart Grid 9, 2986–2998 (2018).
  * Luo  _et al._ (2017) F. Luo, Z. Y. Dong, K. Meng, J. Wen, H. Wang, and J. Zhao, “An operational planning framework for large-scale thermostatically controlled load dispatch,” IEEE Transactions on Industrial Informatics 13, 217–227 (2017).
  * Iacovella  _et al._ (2017) S. Iacovella, F. Ruelens, P. Vingerhoets, B. Claessens, and G. Deconinck, “Cluster control of heterogeneous thermostatically controlled loads using tracer devices,” IEEE Transactions on Smart Grid 8, 528–536 (2017).
  * Zhao  _et al._ (2017) L. Zhao, W. Zhang, H. Hao, and K. Kalsi, “A geometric approach to aggregate flexibility modeling of thermostatically controlled loads,” IEEE Transactions on Power Systems 32, 4721–4731 (2017).
  * Ma  _et al._ (2018) K. Ma, C. Yuan, X. Xu, J. Yang, and Z. Liu, “Optimising regulation of aggregated thermostatically controlled loads based on multi-swarm pso,” IET Generation, Transmission & Distribution 12, 2340–2346 (2018).
  * Ruelens  _et al._ (2017) F. Ruelens, B. J. Claessens, S. Vandael, B. De Schutter, R. Babuška, and R. Belmans, “Residential demand response of thermostatically controlled loads using batch reinforcement learning,” IEEE Transactions on Smart Grid 8, 2149–2159 (2017).
  * Burmester, Rayudu, and Seah (2017) D. Burmester, R. Rayudu, and W. Seah, “Use of maximum power point tracking signal for instantaneous management of thermostatically controlled loads in a dc nanogrid,” IEEE Transactions on Smart Grid (2017).
  * Ghanavati and Chakravarthy (2018) M. Ghanavati and A. Chakravarthy, “Demand-side energy management by use of a design-then-approximate controller for aggregated thermostatic loads,” IEEE Transactions on Control Systems Technology 26, 1439–1448 (2018).
  * Meng  _et al._ (2017) K. Meng, D. Wang, Z. Y. Dong, X. Gao, Y. Zheng, and K. P. Wong, “Distributed control of thermostatically controlled loads in distribution network with high penetration of solar pv,” CSEE Journal of Power and Energy Systems 3, 53–62 (2017).
  * Sinitsyn, Kundu, and Backhaus (2013) N. A. Sinitsyn, S. Kundu, and S. Backhaus, “Safe protocols for generating power pulses with heterogeneous populations of thermostatically controlled loads,” Energy Conversion and Management 67, 297–308 (2013).
  * Mehta  _et al._ (2014) N. Mehta, N. A. Sinitsyn, S. Backhaus, and B. C. Lesieutre, “Safe control of thermostatically controlled loads with installed timers for demand side management,” Energy Conversion and Management 86, 784–791 (2014).
  * Hao  _et al._ (2015) H. Hao, B. M. Sanandaji, K. Poolla, and T. L. Vincent, “Aggregate flexibility of thermostatically controlled loads,” IEEE Transactions on Power Systems 30, 189–198 (2015).
  * Vrettos, Koch, and Andersson (2012) E. Vrettos, S. Koch, and G. Andersson, “Load frequency control by aggregations of thermally stratified electric water heaters,” in  _Innovative Smart Grid Technologies (ISGT Europe), 2012 3rd IEEE PES International Conference and Exhibition on_ (IEEE, 2012) pp. 1–8.
  * Bomela, Zlotnik, and Li Jr (2018) W. Bomela, A. Zlotnik, and S. Li Jr, “A phase model approach for thermostatically controlled load demand response,” (2018).
  * Webborn (2018) E. Webborn, “The dynamics of thermostatically controlled loads for power system frequency control,” (2018).
  * Bajaria and Singh (2019) P. Bajaria and N. Singh, “Dynamic dispatch of thermostatically controlled loads using phase oscillator model,” Sustainable Energy, Grids and Networks , 100220 (2019).
  * Kara (2014) E. C. Kara,  _Data-Driven Approaches to Demand Response: Studies on Thermostatically Controlled Loads and Electric Vehicles_ , Ph.D. thesis, figshare (2014).
  * Kara  _et al._ (2014) E. C. Kara, M. D. Tabone, J. S. MacDonald, D. S. Callaway, and S. Kiliccote, “Quantifying flexibility of residential thermostatically controlled loads for demand response: a data-driven approach,” in  _Proceedings of the 1st ACM Conference on Embedded Systems for Energy-Efficient Buildings_ (ACM, 2014) pp. 140–147.
  * Costanzo  _et al._ (2016) G. T. Costanzo, S. Iacovella, F. Ruelens, T. Leurs, and B. J. Claessens, “Experimental analysis of data-driven control for a building heating system,” Sustainable Energy, Grids and Networks 6, 81–90 (2016).
  * Ruelens  _et al._ (2019) F. Ruelens, B. J. Claessens, P. Vrancx, F. Spiessens, and G. Deconinck, “Direct load control of thermostatically controlled loads based on sparse observations using deep reinforcement learning,” CSEE Journal of Power and Energy Systems 5, 423–432 (2019).
  * Almalaq and Edwards (2017) A. Almalaq and G. Edwards, “A review of deep learning methods applied on load forecasting,” in  _2017 16th IEEE International Conference on Machine Learning and Applications (ICMLA)_ (IEEE, 2017) pp. 511–516.
  * Kazmi  _et al._ (2019) H. Kazmi, J. Suykens, A. Balint, and J. Driesen, “Multi-agent reinforcement learning for modeling and control of thermostatically controlled loads,” Applied energy 238, 1022–1035 (2019).
  * Strogatz (2008) S. Strogatz,  _Sync: How Order Emerges from Chaos in the Universe, Nature, and Daily Life_ (Paw Prints, 2008).
  * Dörfler, Chertkov, and Bullo (2013) F. Dörfler, M. Chertkov, and F. Bullo, “Synchronization in complex oscillator networks and smart grids,” Proceedings of the National Academy of Sciences 110, 2005–2010 (2013).
  * Callaway (2009) D. S. Callaway, “Tapping the energy storage potential in electric loads to deliver load following and regulation, with application to wind energy,” Energy Conversion and Management 50, 1389–1400 (2009).
  * Koch, Mathieu, and Callaway (2011) S. Koch, J. L. Mathieu, and D. S. Callaway, “Modeling and control of aggregated heterogeneous thermostatically controlled loads for ancillary services,” (Citeseer, 2011).
  * Bellman (1966) R. Bellman, “Dynamic programming,” Science 153, 34–37 (1966).
  * Goodfellow, Bengio, and Courville (2016) I. Goodfellow, Y. Bengio, and A. Courville,  _Deep learning_ (MIT press, 2016).
  * Pedregosa  _et al._ (2011) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, _et al._ , “Scikit-learn: Machine learning in python,” Journal of machine learning research 12, 2825–2830 (2011).
  * Chollet  _et al._ (2015) F. Chollet _et al._ , “Keras,” (2015).
  * Nolan and O’Malley (2015) S. Nolan and M. O’Malley, “Challenges and barriers to demand response deployment and evaluation,” Applied Energy 152, 1–10 (2015).
  * Ihara and Schweppe (1981) S. Ihara and F. C. Schweppe, “Physically based modeling of cold load pickup,” IEEE Transactions on Power Apparatus and Systems , 4142–4150 (1981).
  * Rosin (2014) D. P. Rosin,  _Dynamics of complex autonomous Boolean networks_ (Springer, 2014).
  * Abrams and Strogatz (2004) D. M. Abrams and S. H. Strogatz, “Chimera states for coupled oscillators,” Phys. Rev. Lett. 93, 174102 (2004).
  * Acebrón  _et al._ (2005) J. A. Acebrón, L. L. Bonilla, C. J. P. Vicente, F. Ritort, and R. Spigler, “The kuramoto model: A simple paradigm for synchronization phenomena,” Reviews of modern physics 77, 137 (2005).