# On two non-ergodic reversible cellular automata, one classical, the other
quantum

Tomaž Prosen  Faculty of Mathematics and Physics, University of Ljubljana,
Jadranska 19, SI-1000 Ljubljana, Slovenia<EMAIL_ADDRESS>
###### Abstract

We propose and discuss two variants of kinetic particle models — cellular
automata in 1+1 dimensions, which have some appeal due to their simplicity and
intriguing properties which could warrant further research and applications.
The first model is a deterministic and reversible automaton describing two
species of quasiparticles: stable massless _matter_ particles moving with
velocity $\pm 1$ and unstable, standing (zero velocity) _field_ particles. We
discuss two distinct continuity equations for three conserved charges of the
model. While the first two charges and the corresponding currents have support
three (3) lattice sites and represent a lattice analogue of conserved energy-
momentum tensor, we find an additional conserved charge and current with
support of nine (9) sites, implying non-ergodic behaviour and potentially
signalling integrability of the model with a highly nested R-matrix structure.
The second model represents a quantum (or stochastic) deformation of a
recently introduced and studied _charged hardpoint lattice gas_ , where
particles of different binary charge ($\pm 1$) and binary velocity ($\pm 1$)
can nontrivially mix upon elastic collisional scattering. We show that while
the unitary evolution rule of this model does not satisfy the full Yang-Baxter
equation, it still satisfies an intriguing related identity which gives birth
to an infinite set of local conserved operators, the-so-called glider
operators.

To Prof. Giulio Casati, on the occasion of his 80th birthday

##  1 Introduction

This article is dedicated to my friend and mentor, Giulio Casati. For that
reason I will also take the liberty to write the paper from a very personal
perspective. One of the most important things that Giulio taught me, was to
passionately appreciate extremely simple models of dynamics. Constructing and
solving simple models of (lately most often, many-body) dynamics with
carefully chosen physical properties has thus become my personal obsession
throughout my career. The other thing which I owe Giulio was a pragmatic but
deep appreciation of ergodic theory. Together we have been responsible for a
few endeavours in experimental mathematics, which sometimes rose attention of
experts. For instance, we pointed out intriguing ergodic properties of
classical and quantum polygonal billiards and related non-hyperbolic and non-
integrable dynamical systems [CP99, CP00, WCP14, LCP22].

Crudely speaking, ergodic theory [CFS82] divides dynamical systems into
ergodic and non-ergodic ones. The former are characterized by the property
that time averages of (relevant) observables can be computed in terms of
ensemble average over the entire space of states (phase space), while the
latter exhibit memory of the initial condition in typical (say, physically
relevant) observables. Additionally, ergodic dynamics can have various degree
of dynamical complexity or chaos, while the non-ergodic ones, can conform to
some algebraic tools of exact solvability, or integrability. However, it is
not clear to what extent breaking of ergodicity is connected to any form of
integrability as it can precisely be defined in terms of Lax pairs (in
classical, deterministic systems) or Yang-Baxter equation (in quantum or
stochastic setting). In physics literature, these questions have been
extensively discussed recently in the context of many-body dynamics, where
mathematics of ergodic theory is much less developed, say on the lattice and
having local interactions, both in the classical and quantum realm (see e.g.
some recent reviews [AKPR16, AABS19, MBR22, S23]). Besides integrability, many
different forms of ergodicity breaking have been suggested, such as phase
space or Hilbert space fragmentation due to various forms of kinetic
constraints, the so-called many-body localization due to static disorder,
many-body scarring due to hidden weakly broken non-abelian symmetries, etc.

In this paper we propose two very simple many-body locally interacting
dynamical systems defined on a discrete 1+1 space-time lattice. The first
model, which we call a matter-field automation for reasons that will hopefully
appear clear to the reader, is an example of a deterministic particle
kinematics with two distinct types of degrees of freedom (matter&field), which
is non-ergodic for a nontrivial reason: existence of a highly nontrivial but
local conserved charge. Inspired by brute-force empirical computer-algebra
calculations, we suggest that the model may be integrable, but if it is, the
corresponding R-matrix should be highly nontrivial. The second model, which is
a quantum deformation of the previously studied charged hardpoint lattice gas
[MKP17, MPPRV19, KSPIP22], obeys a remarkable reduced Yang-Baxter-like
identity for a non-abelian spectral parameter. This again ensures manifest
ergodicity breaking in the model through an existence of an infinite set of
local translationally invariant conserved operators, the so-called gliders.

##  2 Matter-Field Automaton

###  2.1 Definition of the automaton

Let us define a deterministic reversible cellular automaton on
$\\{0,1\\}^{\mathbb{Z}}\ni{\underline{s}}$ as follows. Even lattice sites
correspond to matter variable $s_{2x}\in\\{0,1\\}$ while odd lattice sites
correspond to field variable $s_{2x+1}\in\\{0,1\\}$, $x\in\mathbb{Z}$. The
instantaneous system configuration is thus specified by an infinite binary
sequence ${\underline{s}}$, on which we specify deterministic dynamics.

The dynamics is defined in a staggered fashion, in even and odd time layers.
Specifically, for so-called _even_ time steps we update the matter-field-
matter triples $(s_{4x},s_{4x+1},s_{4x+2})$, for all $x\in\mathbb{Z}$, with a
deterministic, reversible rule
$(ss^{\prime}s^{\prime\prime})\longrightarrow(rr^{\prime}r^{\prime\prime})$,
while for _odd_ time-steps, the triples $(s_{4x-2},s_{4x-1},s_{4x})$, are
updated by the same rule. The rule is specified as

$\displaystyle(000)\longleftrightarrow(000),$
$\displaystyle(001)\longleftrightarrow(100),$
$\displaystyle(011)\longleftrightarrow(110),$ (1)
$\displaystyle(101)\longleftrightarrow(010),$
$\displaystyle(111)\longleftrightarrow(111),$

or, graphically and perhaps more intuitively, by the following diagrams

0 0 0 1 0 0 1 0 1 1 1 0 0 0 1 0 1 0 0 1 1 1 1 1

0 0 0 0 0 1 0 1 0 0 1 1 1 0 0 1 0 1 1 1 0 1 1 1

which indicates mappings of all eight lower matter-field-matter triples to
corresponding upper matter-field-matter triples, or vice versa due to time
reversal symmetry. The black and red lines thus represent the worldlines of
matter and field particles, respectively. In Fig. 1 we show an example of a
piece of a many-body trajectory $\\{s_{x}^{t};x,t\in\mathbb{Z}\\}$, where time
$t$ runs in vertical direction (upwards). The construction of the trajectory
can be understood as a checkerboard tiling of the rule’s yellow squares which
is uniquely prescribed by the lowest two lines of tiles (the initial
condition).

Note that this dynamical system has a rich set of three independent
$\mathbb{Z}^{2}$ symmetries, specifically C,P, and T:

  1. 1.

“Charge” conjugation C, $s\to 1-s$, i.e., if $s^{t}_{x}$ is a valid
trajectory, then $1-s^{t}_{x}$ is a valid trajectory,

  2. 2.

Parity P, $x\to-x$, i.e., if $s^{t}_{x}$ is a valid trajectory, then
$s^{t}_{-x}$ is a valid trajectory,

  3. 3.

Time-reversal T, $t\to-t$, i.e., if $s^{t}_{x}$ is a valid trajectory, then
$s^{-t}_{x}$ is a valid trajectory.

In Figs. 2-6 we depict some interesting characteristic trajectories of the
matter-field automaton. In Fig. 2 we prepare a maximum entropy initial state
in a finite box of size $L$, where both field and matter are sampled randomly
and uniformly in the set $\\{0,1\\}$ and plot an emerging deterministic
trajectory. Since the initial condition vanishes outside the box,
$s^{t=0}_{x}=0$, for $x<0$ or $x>L$, we see a slow ‘evaporation’ of matter
particles from the edges of the box. In Fig. 3 we plot a similar trajectory,
where initially we have only matter (at maximum entropy) and no field, while
in Fig. 4 we plot the reverse situation where initially there is only the
field at maximum entropy and no matter. The fact that spatio-temporal patterns
of these trajectory seem qualitatively distinct depending on the initial
condition, and in particular the relative ratio of field versus matter
drastically differs, suggests that dynamics may be non-ergodic. We shall
establish this more precisely below by constructing non-trivial algebraic
conservation laws. While previous plots represent cartoons of equilibrium
dynamics, we plot a typical trajectory representing far-from-equilibrium
dynamics in Fig. 5: a head on collision of two streams of matter (with no
field in the initial data), which experience a non-trivial scattering process
via formation of a complex field pattern, finally resulting in two oppositely
moving scattered streams of matter.

Finally, we demonstrate that matter-field system can be understood as a very
general automaton which can also mimic the dynamics of popular Rule 54
reversible cellular automation [BKP21] for a dynamically closed (invariant)
subset of initial configurations. More generally, we can show it can simulate
dynamics for classes of negative-length hard-rod systems with commensurate
initial data. For instance, if we take initial condition ${\underline{s}}$
where only every sixth site can be occupied, i.e. $s_{x}\equiv 0$ for $x\neq
0\pmod{6}$, then one can show that no field worldline can be crossed by
another matter particle, and hence each field particle decays at immediate
next time step (the field is maximally short-lived). This means that
scattering of matter particles always experiences a fixed phase-shift, which
makes dynamics qualitatively similar (yet not precisely equivalent) to Rule 54
on a certain reduced lattice (see Fig. 6 for a demonstration). In Rule 54
dynamics, for instance, the minimal distance between movers/kinks (=2) is
twice the time-tag in kink-scattering (=1), while in this case, the minimal
distance between (left/right) movers (=6) is three times a scattering time-lag
(=2), i.e. free field life-time. Similar commensurate negative-length hard-rod
dynamics is obtained for the subset of initial data $\\{{\underline{s}}\\}$ of
the form $s_{x^{\prime}}=0$, $x^{\prime}\neq 0\pmod{2(2k+1)}$, for any
$k=1,2\ldots$, where the ratio between minimal distance between parallel kinks
to scattering-lag is $2k+1$.

Figure 1: Example of an automaton’s trajectory (time runs in vertical
direction, say upwards). Yellow squares denote application of the local rule,
or equivalently, of the three-bit permutation gate $Y$. Black/red lines denote
the worldlines of matter/field quasiparticles.

###  2.2 Time evolution over abelian algebra of local observables

Let us now reinterpret dynamics of matter-field automaton in terms of an
abelian $C^{*}$-algebraic dynamical system, and then study its conservation
laws. We begin by defining a quasilocal algebra of observables, i.e. real
functions over space of configurations $\\{0,1\\}^{\mathbb{Z}}$. For more
details on a fully analogous but more precise construction in a related model
of Rule 54 reversible cellular automaton see Ref. [KMPV19]. Let $[0]_{2x}$,
$[1]_{2x}$ define a local basis of observables, checking if $s_{2x}=0$ (or
$1)$ for the matter degrees of freedom, and correspondingly
$[0]_{2x+1},[1]_{2x+1}$ for the field:

$[r]_{x}({\underline{s}})=\delta_{r,s_{x}}.$ (2)

We have $[0]_{x}+[1]_{x}=\mathbbm{1}$ (unit element in the algebra of
observables), and denote $[\bullet]_{x}=[1]_{x}-[0]_{x}$, for $x$ either even
or odd. we note that $[0]_{x},[1]_{x}$ (or $\mathbbm{1},[\bullet]_{x}$)
generate the entire (quasi)local algebra under multiplication, addition (and
closure). We denote local observables with larger support as
$[s_{0},s_{1},\ldots,s_{r}]_{x}=[s_{0}]_{x}[s_{1}]_{x+1}\cdots[s_{r}]_{x+r}$.

The fundamental three site operator – local propagator – in the algebra of
observables implementing the rules

$\displaystyle[000]_{2x}\longleftrightarrow[000]_{2x},$
$\displaystyle[001]_{2x}\longleftrightarrow[100]_{2x},$
$\displaystyle[011]_{2x}\longleftrightarrow[110]_{2x},$ (3)
$\displaystyle[101]_{2x}\longleftrightarrow[010]_{2x},$
$\displaystyle[111]_{2x}\longleftrightarrow[111]_{2x},$

is encoded in terms of $8\times 8$ permutation matrix

$Y=\left(\begin{array}[]{cccccccc}1&0&0&0&0&0&0&0\\\ 0&0&0&0&1&0&0&0\\\
0&0&0&0&0&1&0&0\\\ 0&0&0&0&0&0&1&0\\\ 0&1&0&0&0&0&0&0\\\ 0&0&1&0&0&0&0&0\\\
0&0&0&1&0&0&0&0\\\ 0&0&0&0&0&0&0&1\\\ \end{array}\right)\,.$ (4)

The full time step of dynamical automorphism on the algebra of observables
$a^{t}({\underline{s}})=a({\underline{s}}^{-t})$,

$a^{2(t+1)}=Ua^{2t},\qquad a^{2t+1}=U_{\rm e}a^{2t},\quad a^{2t+2}=U_{\rm
o}a^{2t+1},$ (5)

is constructed as follows

$U=U_{\rm o}U_{\rm e}\,,$ (6)

where even and odd time steps are generated as

$\displaystyle U_{\rm e}=\prod_{x}Y_{4x,4x+1,4x+2},\quad U_{\rm
o}=\prod_{x}Y_{4x-2,4x-1,4x}\,,$ (7)

where three indices denote the positions in the string, or in tensor product
of local algebras ($\bigotimes_{x}{\rm span}\\{[0]_{x},[1]_{x}\\}$), where the
three site operator (4) acts nontrivially.

Figure 2: Example of a typical trajectory in a maximum entropy initial state,
where $s_{x}^{t=0}$ are sampled randomly and uniformly in $\\{0,1\\}$, for a
section of $0\leq x\leq 600$ and duration $t_{\rm max}=150$, $0\leq t\leq
t_{\rm max}$. Note the evaporation of matter particles from the sides, and
reduction of the field intensity, which is due to initial vacuum in the
complement region $s^{t=0}_{x^{\prime}}=0$, $x^{\prime}<0$ or
$x^{\prime}>600$. Figure 3: Example of a typical trajectory from initial state
where the matter is in maximal entropy state, i.e. $s_{2x}$ are sampled
randomly and uniformly for $0\leq 2x\leq 600$, and there is _no_ field, i.e.
$s_{2x+1}\equiv 0$, of duration $t_{\rm max}=150$, $0\leq t\leq t_{\rm max}$.
Figure 4: Example of a typical trajectory from initial state where the field
is in maximal entropy state, i.e. $s_{2x+1}$ are sampled randomly and
uniformly for $0\leq 2x+1\leq 600$, and there is _no_ matter, i.e.
$s_{2x}\equiv 0$, of duration $t_{\rm max}=150$, $0\leq t\leq t_{\rm max}$.
Figure 5: Head-on collision of a random (maximum entropy) sample of left
movers with random sample of right movers. The total width of each initial
stream of left/right movers is $800$ and duration is $t_{\rm max}=600$. Figure
6: Trajectory starting from a random initial condition where only $s_{6x}$
($0\leq x\leq 180$) are chosen random, uniform in $\\{0,1\\}$, while all other
$s_{x^{\prime}}=0$. Such trajectory features only maximally short-lived field
particles which provide ‘unit’ phase shifts for left/right movers upon
scattering and corresponds to dynamics of negative-length hard-rods, like e.g.
in the Rule 54 system [BKP21].

###  2.3 Conservation laws

Let $a_{[x,y]}$ be any local observable supported in the interval of sites
$[x,y]\subset\mathbb{Z}$. We define a shift automorphism $\Pi$

$\Pi a_{[x,y]}=a_{[x+2,y+2]}\,,$ (8)

which shifts an observable by ‘one lattice unit’, so that the identity of
matter/field sites is preserved. Note that

$\Pi U^{\rm e}=U^{\rm o}\Pi,\quad\Pi U^{\rm o}=U^{\rm
e}\Pi,\quad[\Pi^{2},U^{\rm e}]=0,\quad[\Pi^{2},U^{\rm o}]=0.$ (9)

A local observable $q$ is called a conserved charge, if together with an
appropriate conserved current $p$, satisfies a conservation law (space-time
discrete continuity equation):

$U_{\rm e}q-U_{\rm o}q+\Pi p-\Pi^{-1}p=0.$ (10)

Then, extensive local observable is globally conserved

$Q=\sum_{x\in\mathbb{Z}}q_{x},\quad q_{x}:=\Pi^{2x}q,\qquad UQ=Q.$ (11)

Limiting to local observables on clusters of _five_ sites (beginning end
ending with a matter site), covering one ’unit cell’ (periodicity block of 4
sites) of the automaton, we find (using exact computer algebra) exactly two
conservation laws

$\displaystyle q^{(1)}_{x}$ $\displaystyle=$
$\displaystyle[\bullet]_{4x}-[\bullet]_{4x+2},$ (12) $\displaystyle
p^{(1)}_{x}$ $\displaystyle=$
$\displaystyle\frac{1}{2}[\bullet]_{4x}+\frac{1}{2}[\bullet]_{4x+1}+\frac{1}{2}[\bullet]_{4x+2}+\frac{1}{2}[\bullet]_{4x}[\bullet]_{4x+1}[\bullet]_{4x+2},$
(13) $\displaystyle q^{(2)}_{x}$ $\displaystyle=$
$\displaystyle[\bullet]_{4x+1}+[\bullet]_{4x+2}+[\bullet]_{4x+3},$ (14)
$\displaystyle p_{x}^{(2)}$ $\displaystyle=$ $\displaystyle-[\bullet]_{4x+2},$
(15)

with local densities and current supported on up to three (3) consecutive
sites. $q^{(1)}_{x}$ can be interpreted as the momentum density while
$q^{(2)}_{x}$ is the analog of the energy density. Indeed, writing
$P=\sum_{x}q^{(1)}_{x}$, $E=\sum_{x}q^{(2)}_{x}$ we have an intuitive
expressions for total momentum and total energy of a configuration:

$P({\underline{s}})=\sum_{x}(-1)^{x}s_{2x}\,,\quad
E({\underline{s}})=\sum_{x}(s_{2x}+2s_{2x+1})\,,$ (16)

meaning that the field particle caries no momentum and has an energy content
(rest energy) equal to the energy of two matter particles.

If we then increase the support size to clusers of _nine_ ($2\times 4+1$)
consecutive sites, we obtain precisely one new additional conservation law:
(using compact notation: $[\circ]\equiv\mathbbm{1}$,
$[r_{0}r_{1}r_{2}\ldots]_{4x}\equiv[r_{0}r_{1}r_{2}\ldots]$)

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!q^{(3)}=2[\circ\circ\circ\circ\circ\circ\bullet\circ\bullet]-2[\circ\circ\circ\circ\circ\circ\bullet\bullet\circ]-4[\circ\circ\circ\circ\circ\bullet\circ\circ\bullet]-2[\circ\circ\circ\circ\circ\bullet\bullet\circ\circ]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+2[\circ\circ\circ\circ\circ\bullet\bullet\bullet\bullet]-8[\circ\circ\circ\circ\bullet\circ\circ\circ\bullet]-4[\circ\circ\circ\circ\bullet\circ\circ\bullet\circ]+2[\circ\circ\circ\circ\bullet\circ\bullet\circ\circ]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!-2[\circ\circ\circ\circ\bullet\circ\bullet\bullet\bullet]+4[\circ\circ\circ\circ\bullet\bullet\circ\bullet\bullet]-2[\circ\circ\circ\circ\bullet\bullet\bullet\circ\bullet]+2[\circ\circ\circ\circ\bullet\bullet\bullet\bullet\circ]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!-2[\circ\circ\circ\bullet\circ\circ\circ\circ\bullet]-2[\circ\circ\circ\bullet\circ\circ\circ\bullet\circ]+2[\circ\circ\circ\bullet\circ\circ\bullet\circ\circ]-2[\circ\circ\circ\bullet\circ\circ\bullet\bullet\bullet]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!-4[\circ\circ\circ\bullet\circ\bullet\circ\circ\circ]+2[\circ\circ\circ\bullet\bullet\bullet\circ\circ\bullet]+2[\circ\circ\circ\bullet\bullet\bullet\circ\bullet\circ]-2[\circ\circ\circ\bullet\bullet\bullet\bullet\circ\circ]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+2[\circ\circ\circ\bullet\bullet\bullet\bullet\bullet\bullet]+[\circ\circ\bullet\circ\circ\circ\circ\circ\bullet]+[\circ\circ\bullet\circ\circ\circ\circ\bullet\circ]-[\circ\circ\bullet\circ\circ\circ\bullet\circ\circ]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+[\circ\circ\bullet\circ\circ\circ\bullet\bullet\bullet]+2[\circ\circ\bullet\circ\circ\bullet\circ\circ\circ]-[\circ\circ\bullet\circ\bullet\bullet\circ\circ\bullet]-[\circ\circ\bullet\circ\bullet\bullet\circ\bullet\circ]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+[\circ\circ\bullet\circ\bullet\bullet\bullet\circ\circ]-[\circ\circ\bullet\circ\bullet\bullet\bullet\bullet\bullet]+[\circ\circ\bullet\bullet\circ\bullet\circ\circ\bullet]+[\circ\circ\bullet\bullet\circ\bullet\circ\bullet\circ]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!-[\circ\circ\bullet\bullet\circ\bullet\bullet\circ\circ]+[\circ\circ\bullet\bullet\circ\bullet\bullet\bullet\bullet]-[\circ\circ\bullet\bullet\bullet\circ\circ\circ\bullet]-[\circ\circ\bullet\bullet\bullet\circ\circ\bullet\circ]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+[\circ\circ\bullet\bullet\bullet\circ\bullet\circ\circ]-[\circ\circ\bullet\bullet\bullet\circ\bullet\bullet\bullet]-2[\circ\circ\bullet\bullet\bullet\bullet\circ\circ\circ]-[\circ\bullet\circ\circ\circ\circ\circ\circ\bullet]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!-[\circ\bullet\circ\circ\circ\circ\circ\bullet\circ]+[\circ\bullet\circ\circ\circ\circ\bullet\circ\circ]-[\circ\bullet\circ\circ\circ\circ\bullet\bullet\bullet]-2[\circ\bullet\circ\circ\circ\bullet\circ\circ\circ]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+[\circ\bullet\circ\circ\bullet\bullet\circ\circ\bullet]+[\circ\bullet\circ\circ\bullet\bullet\circ\bullet\circ]-[\circ\bullet\circ\circ\bullet\bullet\bullet\circ\circ]+[\circ\bullet\circ\circ\bullet\bullet\bullet\bullet\bullet]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!-[\circ\bullet\circ\bullet\circ\bullet\circ\circ\bullet]-[\circ\bullet\circ\bullet\circ\bullet\circ\bullet\circ]+[\circ\bullet\circ\bullet\circ\bullet\bullet\circ\circ]-[\circ\bullet\circ\bullet\circ\bullet\bullet\bullet\bullet]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+[\circ\bullet\circ\bullet\bullet\circ\circ\circ\bullet]+[\circ\bullet\circ\bullet\bullet\circ\circ\bullet\circ]-[\circ\bullet\circ\bullet\bullet\circ\bullet\circ\circ]+[\circ\bullet\circ\bullet\bullet\circ\bullet\bullet\bullet]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+2[\circ\bullet\circ\bullet\bullet\bullet\circ\circ\circ]-[\bullet\circ\circ\circ\circ\circ\circ\circ\bullet]-[\bullet\circ\circ\circ\circ\circ\circ\bullet\circ]+[\bullet\circ\circ\circ\circ\circ\bullet\circ\circ]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!-[\bullet\circ\circ\circ\circ\circ\bullet\bullet\bullet]-2[\bullet\circ\circ\circ\circ\bullet\circ\circ\circ]+[\bullet\circ\circ\circ\bullet\bullet\circ\circ\bullet]+[\bullet\circ\circ\circ\bullet\bullet\circ\bullet\circ]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!-[\bullet\circ\circ\circ\bullet\bullet\bullet\circ\circ]+[\bullet\circ\circ\circ\bullet\bullet\bullet\bullet\bullet]-[\bullet\circ\circ\bullet\circ\bullet\circ\circ\bullet]-[\bullet\circ\circ\bullet\circ\bullet\circ\bullet\circ]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+[\bullet\circ\circ\bullet\circ\bullet\bullet\circ\circ]-[\bullet\circ\circ\bullet\circ\bullet\bullet\bullet\bullet]+[\bullet\circ\circ\bullet\bullet\circ\circ\circ\bullet]+[\bullet\circ\circ\bullet\bullet\circ\circ\bullet\circ]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!-[\bullet\circ\circ\bullet\bullet\circ\bullet\circ\circ]+[\bullet\circ\circ\bullet\bullet\circ\bullet\bullet\bullet]+2[\bullet\circ\circ\bullet\bullet\bullet\circ\circ\circ]-[\bullet\bullet\bullet\circ\circ\circ\circ\circ\bullet]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!-[\bullet\bullet\bullet\circ\circ\circ\circ\bullet\circ]+[\bullet\bullet\bullet\circ\circ\circ\bullet\circ\circ]-[\bullet\bullet\bullet\circ\circ\circ\bullet\bullet\bullet]-2[\bullet\bullet\bullet\circ\circ\bullet\circ\circ\circ]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+[\bullet\bullet\bullet\circ\bullet\bullet\circ\circ\bullet]+[\bullet\bullet\bullet\circ\bullet\bullet\circ\bullet\circ]-[\bullet\bullet\bullet\circ\bullet\bullet\bullet\circ\circ]+[\bullet\bullet\bullet\circ\bullet\bullet\bullet\bullet\bullet]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!-[\bullet\bullet\bullet\bullet\circ\bullet\circ\circ\bullet]-[\bullet\bullet\bullet\bullet\circ\bullet\circ\bullet\circ]+[\bullet\bullet\bullet\bullet\circ\bullet\bullet\circ\circ]-[\bullet\bullet\bullet\bullet\circ\bullet\bullet\bullet\bullet]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+[\bullet\bullet\bullet\bullet\bullet\circ\circ\circ\bullet]+[\bullet\bullet\bullet\bullet\bullet\circ\circ\bullet\circ]-[\bullet\bullet\bullet\bullet\bullet\circ\bullet\circ\circ]+[\bullet\bullet\bullet\bullet\bullet\circ\bullet\bullet\bullet]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+2[\bullet\bullet\bullet\bullet\bullet\bullet\circ\circ\circ]$
(17)

with the current

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!p^{(3)}=2[\circ\circ\circ\circ\circ\bullet\bullet\circ\circ]-4[\circ\circ\circ\circ\bullet\circ\circ\circ\bullet]-4[\circ\circ\circ\circ\bullet\circ\circ\bullet\circ]+2[\circ\circ\circ\circ\bullet\circ\bullet\circ\circ]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!-4[\circ\circ\circ\circ\bullet\circ\bullet\bullet\bullet]-2[\circ\circ\circ\bullet\circ\circ\circ\circ\bullet]-2[\circ\circ\circ\bullet\circ\circ\circ\bullet\circ]-2[\circ\circ\circ\bullet\circ\circ\bullet\bullet\bullet]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+2[\circ\circ\circ\bullet\bullet\bullet\circ\circ\bullet]+2[\circ\circ\circ\bullet\bullet\bullet\circ\bullet\circ]+2[\circ\circ\circ\bullet\bullet\bullet\bullet\bullet\bullet]+[\circ\circ\bullet\circ\circ\circ\circ\circ\bullet]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+[\circ\circ\bullet\circ\circ\circ\circ\bullet\circ]+[\circ\circ\bullet\circ\circ\circ\bullet\bullet\bullet]-[\circ\circ\bullet\circ\bullet\bullet\circ\circ\bullet]-[\circ\circ\bullet\circ\bullet\bullet\circ\bullet\circ]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!-[\circ\circ\bullet\circ\bullet\bullet\bullet\bullet\bullet]+[\circ\circ\bullet\bullet\circ\bullet\circ\circ\bullet]+[\circ\circ\bullet\bullet\circ\bullet\circ\bullet\circ]+[\circ\circ\bullet\bullet\circ\bullet\bullet\bullet\bullet]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!-[\circ\circ\bullet\bullet\bullet\circ\circ\circ\bullet]-[\circ\circ\bullet\bullet\bullet\circ\circ\bullet\circ]-[\circ\circ\bullet\bullet\bullet\circ\bullet\bullet\bullet]-[\circ\bullet\circ\circ\circ\circ\circ\circ\bullet]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!-[\circ\bullet\circ\circ\circ\circ\circ\bullet\circ]-[\circ\bullet\circ\circ\circ\circ\bullet\bullet\bullet]+[\circ\bullet\circ\circ\bullet\bullet\circ\circ\bullet]+[\circ\bullet\circ\circ\bullet\bullet\circ\bullet\circ]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+[\circ\bullet\circ\circ\bullet\bullet\bullet\bullet\bullet]-[\circ\bullet\circ\bullet\circ\bullet\circ\circ\bullet]-[\circ\bullet\circ\bullet\circ\bullet\circ\bullet\circ]-[\circ\bullet\circ\bullet\circ\bullet\bullet\bullet\bullet]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+[\circ\bullet\circ\bullet\bullet\circ\circ\circ\bullet]+[\circ\bullet\circ\bullet\bullet\circ\circ\bullet\circ]+[\circ\bullet\circ\bullet\bullet\circ\bullet\bullet\bullet]-[\bullet\circ\circ\circ\circ\circ\circ\circ\bullet]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!-[\bullet\circ\circ\circ\circ\circ\circ\bullet\circ]-[\bullet\circ\circ\circ\circ\circ\bullet\bullet\bullet]+[\bullet\circ\circ\circ\bullet\bullet\circ\circ\bullet]+[\bullet\circ\circ\circ\bullet\bullet\circ\bullet\circ]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+[\bullet\circ\circ\circ\bullet\bullet\bullet\bullet\bullet]-[\bullet\circ\circ\bullet\circ\bullet\circ\circ\bullet]-[\bullet\circ\circ\bullet\circ\bullet\circ\bullet\circ]-[\bullet\circ\circ\bullet\circ\bullet\bullet\bullet\bullet]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+[\bullet\circ\circ\bullet\bullet\circ\circ\circ\bullet]+[\bullet\circ\circ\bullet\bullet\circ\circ\bullet\circ]+[\bullet\circ\circ\bullet\bullet\circ\bullet\bullet\bullet]-[\bullet\bullet\bullet\circ\circ\circ\circ\circ\bullet]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!-[\bullet\bullet\bullet\circ\circ\circ\circ\bullet\circ]-[\bullet\bullet\bullet\circ\circ\circ\bullet\bullet\bullet]+[\bullet\bullet\bullet\circ\bullet\bullet\circ\circ\bullet]+[\bullet\bullet\bullet\circ\bullet\bullet\circ\bullet\circ]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+[\bullet\bullet\bullet\circ\bullet\bullet\bullet\bullet\bullet]-[\bullet\bullet\bullet\bullet\circ\bullet\circ\circ\bullet]-[\bullet\bullet\bullet\bullet\circ\bullet\circ\bullet\circ]-[\bullet\bullet\bullet\bullet\circ\bullet\bullet\bullet\bullet]$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+[\bullet\bullet\bullet\bullet\bullet\circ\circ\circ\bullet]+[\bullet\bullet\bullet\bullet\bullet\circ\circ\bullet\circ]+[\bullet\bullet\bullet\bullet\bullet\circ\bullet\bullet\bullet].$
(18)

If we furthermore increase to 3-cell cluster of $13=3\times 4+1$ consecutive
sites, we obtain only one additional conservation law which, however, is
trivially related to the previous one:

$q^{(4)}_{x}=\Pi q^{(3)}_{x}=q^{(3)}_{x+2},\qquad p^{(4)}_{x}=-\Pi
p^{(3)}_{x}=-p^{(3)}_{x+2}.$ (19)

Of course, it holds in general, that if $(q,p)$ is a local conservation law,
then $(\Pi q,-\Pi p)$ is also a conservation law. Scanning for existence of
higher conservation laws with densities supported over four or more unit cells
seems impossible with brute-force computer algebra.

Also note that the conservation law property (10) is invariant under local
space-time discrete gauge transformation, i.e. taking any local observable $a$

$\displaystyle q$ $\displaystyle\longrightarrow q+\Pi a-\Pi^{-1}a,$ (20)
$\displaystyle p$ $\displaystyle\longrightarrow p-U_{\rm e}a+U_{\rm o}a.$ (21)

In the explicit expressions listed above, the gauge has been fixed by right-
alignment of the densities.

One might expect that matter-field automaton hosts an infinite tower of local
conserved quantities of increasing support size, however these may be very
difficult to find empirically. This hypothesis would suggest also that the
matter-field automaton may be a completely integrable system in a similar
spirit as Rule 54 [GP22b], as they also display similar negative-length hard
rod dynamics for a class of initial data, but preliminary attempts to find a
Yang-Baxter or Lax structure failed. The other option is that matter-field
system has only a few, perhaps finitely many (or an incomplete set) of
conserved local charges and would then represent a paradigmatic theory between
integrable and ergodic dynamics.

Deterministic cellular automata on discrete phase space (space of
configurations) often admit continuous, single (or few) parameter
modifications, which render the models quantum or stochastic (depending on the
type of deformation which turns a deterministic evolution to a unitary or
markov map). For instance, in the definition of matter-field automaton, one
may in each step allow for the tripple $(010)$ to remain $(010)$ with
probability $\alpha$ and to map to $(101)$ with probability $1-\alpha$.
Similarly, $(101)$ remains in $(101)$ with probability $\beta$ and maps to
$(010)$ with probability $1-\beta$. All the definitions of the previous
section remain valid, except that the permutation matrix $Y$ is now replaced
by a stochastic (markov) matrix:

$Y_{\rm stoch}=\left(\begin{array}[]{cccccccc}1&0&0&0&0&0&0&0\\\
0&0&0&0&1&0&0&0\\\ 0&0&\alpha&0&0&1-\beta&0&0\\\ 0&0&0&0&0&0&1&0\\\
0&1&0&0&0&0&0&0\\\ 0&0&1-\alpha&0&0&\beta&0&0\\\ 0&0&0&1&0&0&0&0\\\
0&0&0&0&0&0&0&1\\\ \end{array}\right)\,.$ (22)

Interestingly, the first two conservation laws, the momentum and energy (16),
see also (12,13), survive the stochastic deformation, while, as indicated by
computer algebra, the third conservation law $(q^{(3)},p^{(3)})$ is _broken_.
Hence, the deformation of matter-field automaton likely no longer posses
nontrivial ergodicity breaking.

##  3 Deformed (quantized) hardpoint lattice gas

Nevertheless, we show briefly in this section how continuous deformations of
some other reversible deterministic cellular automata can lead to interesting
non-ergodic dynamics. For this purpose we use the so-called charged hardpoint
lattice gas [MKP17] (also referred to as XXC model in [GP22a] since it is a
deterministic limit of models introduced and studied by Maassarani [M98]).

In comparison to matter-field automaton, the model has a richer local
configuration space, i.e. three states per site, and a simpler local
interaction map, applying only to a pair of neighbouring sites instead of
three. In physical terms, the model describes dynamics of free point particles
which move with speeds $\pm 1$ with initial positions placed at commensurate
(multiples of integer) coordinates, but which may carry an internal (say
binary) degree of freedom (charge). Denoting the three local states as
$s\in\\{0,+,-\\}$, 0 referred to as a vacancy, the local reversible and
deterministic rule, which we apply sequentially to pairs $(2x,2x+1)$ and
$(2x+1,2x)$ reads

$(00)\longleftrightarrow(00)\,,\quad(0\alpha)\longleftrightarrow(\alpha
0)\,,\quad(\alpha\beta)\longleftrightarrow(\alpha\beta)\,,\quad\alpha,\beta\in\\{+,-\\}.$
(23)

The model allows for a remarkable set of exact results on equilibrium and non-
equilibrium dynamics, ranging from dynamical correlations in equilibrium and
the transport coefficients [MKP17], to analytic _matrix-product-state_
description of time dependent quenches and non-equilibrium steady states in
boundary driven setup [MPPRV19], to recent exact solution of full counting
statistics [KSPIP22]. One may wander if the model still allows for some degree
of solvability, if hard-core collision is relaxed to a general stochastic, or
unitary scattering

$(+-)\longrightarrow u_{11}(+-)+u_{12}(-+),\qquad(-+)\longrightarrow
u_{21}(+-)+u_{22}(-+),$ (24)

where $u$ may be a $2\times 2$ stochastic, or unitary matrix. In the latter
case, we think of quantum lattice gas model and quantum superposition states
as elaborated precisely below. We note that this system is closely related to
a semiclassically quantized sine-Gordon model recently studied in [KVZ22].

Figure 7: Diagrammatic proof of glider property (31) in terms of equivalence
of quantum circuits, where gates $U[u],U[\mathbbm{1}],U[u^{-1}]$ are indicated
by blue boxes. The proof of (32) is analogous (via left-right reflection of
the above diagrams).

Consider a local Hilbert space $\mathcal{H}_{1}=\mathbb{C}^{3}$ with states
${|s\rangle}_{s=0,+,-}$. We refer to states ${|\pm\rangle}$ as charged
particles and we denote them by greek index ${|\alpha\rangle}_{\alpha=\pm}$,
and to the state ${|0\rangle}$ as vacancy. Let us define a local propagator
over $\mathcal{H}_{1}\otimes\mathcal{H}_{1}$ as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!U[u]$
$\displaystyle=$ $\displaystyle{|00\rangle}{\langle
00|}+\sum_{\alpha\in\\{\pm\\}}({|0\alpha\rangle}{\langle\alpha 0|}+{|\alpha
0\rangle}{\langle 0\alpha|}+{|\alpha\alpha\rangle}{\langle\alpha\alpha|})+$
$\displaystyle+$ $\displaystyle
u_{11}{|+-\rangle}{\langle+-|}+u_{12}{|+-\rangle}{\langle-+|}+u_{21}{|-+\rangle}{\langle+-|}+u_{22}{|-+\rangle}{\langle-+|},$

where $u$ can be in principle any invertible matrix, $u\in GL(2)$, for the
properties that will be discussed below, while for quantum physics
applications we will think of $u$ as a unitary mixing matrix. Such $U[u]$
satisfies two remarkable identities, valid for any $u$:

$\displaystyle
U_{12}[u]U_{23}[\mathbbm{1}]U_{12}[u^{-1}]=U_{23}[u^{-1}]U_{12}[\mathbbm{1}]U_{23}[u],$
(26) $\displaystyle U[u]U[u^{-1}]=\mathbbm{1},$ (27)

which can be checked by straightforward computation. Note that identity (26)
bears some similarity to the braid group form of Yang-Baxter equation if $u$
is considered as a non-abelian spectral parameter. Specifically,
$U[\mathbbm{1}]$ obeys the braid relation – it constitutes the so-called Yang-
Baxter map – and provides a deterministic update rule for the charged
hardpoint lattice gas. However $U[u]$ does not have the full Yang-Baxter
property, but _just enough_ to render an existence of an infinite set of
conserved local operators, as shown below.

For simplicity we can now assume a finite system of even number of sites $L$
and define its many-body Hilbert space as a tensor product
$\mathcal{H}=\mathcal{H}_{1}^{\otimes L}$. Then we define a locality
preserving discrete time dynamical system (a quantum cellular automaton) over
it in terms of a brickwork circuit, i.e.

$\mathcal{U}=\mathcal{U}_{\rm o}\mathcal{U}_{\rm e},\quad\mathcal{U}_{\rm
e}=\prod_{x=1}^{L/2}U_{2x-1,2x}[u],\quad\mathcal{U}_{\rm
o}=\prod_{x=1}^{L/2}U_{2x,2x+1}[u],\quad L+1\equiv 1,$ (28)

and write evolution of observables as

$a^{2t}=\mathcal{U}_{\rm e}a^{2t-1}\mathcal{U}_{\rm e}^{-1},\quad
a^{2t+1}=\mathcal{U}_{\rm o}a^{2t}\mathcal{U}_{\rm o}^{-1}.$ (29)

Let us define the following pair of local operators (note periodic boundary
conditions if needed):

$\displaystyle g^{(+)}_{x}$ $\displaystyle=$ $\displaystyle
U_{x,x+1}[u]U_{x-1,x}U[\mathbbm{1}]U_{x,x+1}[u^{-1}],$ (30) $\displaystyle
g^{(-)}_{x}$ $\displaystyle=$ $\displaystyle
U_{x-1,x}[u]U_{x,x+1}U[\mathbbm{1}]U_{x-1,x}[u^{-1}].$

Observing identities (26,27), we can straightforwardly derive a few remarkable
properties of these operators:

$\displaystyle\mathcal{U}_{\rm e}g^{(+)}_{2x}\mathcal{U}^{-1}_{\rm e}$
$\displaystyle=$ $\displaystyle g^{(+)}_{2x+1},\qquad\mathcal{U}_{\rm
o}g^{(+)}_{2x-1}\mathcal{U}^{-1}_{\rm o}=g^{(+)}_{2x},$ (31)
$\displaystyle\mathcal{U}_{\rm e}g^{(-)}_{2x}\mathcal{U}^{-1}_{\rm e}$
$\displaystyle=$ $\displaystyle g^{(-)}_{2x-1},\qquad\mathcal{U}_{\rm
o}g^{(-)}_{2x+1}\mathcal{U}^{-1}_{\rm o}=g^{(-)}_{2x},$ (32)

which justify the term glider operators [GUWZ10] (see Fig. 7 for a simple
graphical proof of the properties (31,32) in terms of quantum circuit
representation). We note that our gliders would be trivial
$g^{(\pm)}_{x}=\mathbbm{1}$ for integrable trotterization [VZP18] where
$U(\mathbbm{1})=\mathbbm{1}$.

The gliders are particular kind of quantum conservation laws, where the formal
charge and current densities are identical. Consequently, we can build
products of such operators, separated by at least two pairs of sites, which
again behave in the same way

$g^{(\pm)}_{x_{1},x_{2}\ldots
x_{m}}=\prod_{n=1}^{m}g^{(\pm)}_{2x_{n}},\quad\mathcal{U}g^{(\pm)}_{x_{1},x_{2}\ldots
x_{m}}\mathcal{U}^{-1}=g^{(\pm)}_{x_{1}\pm 1,x_{2}\pm 1\ldots x_{m}\pm 1}$
(33)

under constraints $x_{j+1}-x_{j}\geq 2$. We thus have a large set of extensive
charges

$\displaystyle Q^{(\pm)}_{r_{1},r_{2}\ldots r_{m}}$ $\displaystyle=$
$\displaystyle\sum_{x=1}^{L/2}g^{(\pm)}_{x,x+r_{1}\ldots x+r_{m}}\,,\quad
r_{1}\geq 2,\;r_{m+1}-r_{m}\geq 2,$ $\displaystyle
Q^{(\pm)}_{{\underline{r}}}$ $\displaystyle\equiv$
$\displaystyle\mathcal{U}Q^{({\pm})}_{{\underline{r}}}\mathcal{U}^{-1}.$ (34)

These conserved operators concisely encapsulate the non-ergodic properties of
the deformed hardpoint lattice gas. The non-ergodicity should not be
considered surprising as the pattern of left and right movers is evolving
regularly according to ballistic (free) dynamics. However, internal (charge)
degrees of freedom undergo interacting dynamics which is nontrivially
constrained by conserving $Q^{\pm}_{{\underline{r}}}$.

The fact that the number of distinct glider charges
$Q^{(\pm)}_{{\underline{r}}}$ increases exponentially with the size of
support, say $\ell$ such that all $r_{j}\leq\ell$, is reminiscent of
superintegrable quantum cellular automata [BKP21, GP22b, GP22a], yet the
absence of exact Yang-Baxter property suggests that the dynamical system is
not completely integrable (i.e. there should be no Bethe ansatz!?). We are
thus suggesting a new class of non-ergodic many-body quantum dynamics. It
would be interesting to understand if there is any connection to a class of
non-ergodic dual unitary circuits (for definition and basic properties see
[BKP19]) possessing glider operators without any (even partial) Yang-Baxter
property (introduced and studied in Ref. [BP22]).

##  4 Conclusion

The dynamical systems introduced and discussed in this short contribution have
been constructed for the sole amusement of the author. Yet, we hope these
constructions may at some point find some non-trivial applications in the
context of non-equilibrium statistical mechanics. Specifically, one of the
long-standing key goals in the field is to derive irreversible macroscopic
transport laws, such as Fick’s law of diffusion or Fourier’s law of heat
transport, from microscopic reversible equations of motion. This question has
also been among the main research focuses of Giulio Casati [CFVV84] as well as
stimulating the beginning of my research career [PR92], and some of our joint
work investigated the Fourier’s law in systems very similar to the models
studied here: hardpoint colliding masses in one dimension [CP03].

The coexistence of diffusive and ballistic transport has been rigorously
demonstrated in Rule 54 automaton (reviewed in Ref. [BKP21], see also Ref.
[GHKV18]), and one may hope that such results can be extended to the matter-
field automaton, which seems richer and perhaps more generic (for example,
Rule 54 dynamics is super-integrable with exponentially growing number of
conserved local charges, while the matter-field automaton is perhaps at most
integrable as it has much rarer nontrivial local charges).

I acknowledge enjoying numerous discussions with Balazs Pozsgay on related
problems. This work has been supported by the Slovenian Research Agency (ARRS)
under the Program P1-0402 and grants N1-0233 and N1-0219.

## References

## References

  * [CP99] Casati, G.; Prosen, T. Mixing property of triangular billiards, Phys. Rev. Lett. 1999, 83, 4729–4732.
  * [CP00] Casati, G.; Prosen, T. Triangle map: A model of quantum chaos, Phys. Rev. Lett. 2000, 85, 4261-4264.
  * [WCP14] Casati, G.; Prosen, T. Nonergodicity and localization of invariant measure for two colliding masses, Phys. Rev. E 2014, 89, 042918.
  * [LCP22] Lozej, Č.; Casati, G.; Prosen, T. Quantum chaos in triangular billiards, Phys. Rev. Res. 2022, 4, 013138.
  * [CFS82] Cornfeld, I. P.; Fomin, S. V.; Sinai, Ya. G. Ergodic theory , Springer-Verlag New York, USA, 1982.
  * [AKPR16] D’Alessio, L.; Kafri, Y.; Polkovnikov, A.; Rigol, M. From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics, Adv. Phys. 2016, 65, 239-362.
  * [AABS19] Abanin, D. A.; Altman E.; Bloch, I.; Serbyn, M. Many-body localization, thermalization, and entanglement, Rev. Mod. Phys. 2019, 91, 021001.
  * [MBR22] Moudgalya, S.; Bernevig, B. A.; Regnault, N. Quantum many-body scars and Hilbert space fragmentation: a review of exact results, Rep. Prog. Phys. 2022, 85, 086501.
  * [S23] Spohn, H. Hydrodynamic scales of integrable many-particle systems, preprint 2023, arXiv:2301.08504.
  * [MKP17] Medenjak, M.; Klobas, K.; Prosen, T. Diffusion in Deterministic Interacting Lattice Systems, Phys. Rev. Lett. 2017, 119, 110603.
  * [MPPRV19] Medenjak, M.; Popkov, V.; Prosen, T.; Ragoucy, R.; Vanicat, M. Two-species hardcore reversible cellular automaton: matrix ansatz for dynamics and nonequilibrium stationary state, SciPost Phys. 2019, 6, 074.
  * [KSPIP22] Krajnik, Ž.; Schmidt, J.; Pasquier, V.; Ilievski, I.; Prosen, T. Exact Anomalous Current Fluctuations in a Deterministic Interacting Model, Phys. Rev. Lett. 2022, 128, 160601.
  * [GP22a] Gombor, T.; Pozsgay, B. Superintegrable cellular automata and dual unitary gates from Yang-Baxter maps, SciPost Phys. 2022, 12, 102.
  * [M98] Maassarani, Z. The XXC models, Phys. Lett. A 1998, 244, 160–164.
  * [KVZ22] Kormos, M.; Vörös, D.; Zaránd, G. Finite temperature dynamics in gapped 1D models in the sine-Gordon family, Phys. Rev. B 2022, 106, 205151.
  * [KMPV19] Klobas, K.; Medenjak, M.; Prosen, T.; Vanicat, M. Time-Dependent Matrix Product Ansatz for Interacting Reversible Dynamics, Commun. Math. Phys. 2019, 371, 651–688.
  * [BKP21] Buča, B.; Klobas, K.; Prosen, T. Rule 54: Exactly solvable model of nonequilibrium statistical mechanics, J. Stat. Mech. 2021, 074001.
  * [GP22b] Gombor, T.; Pozsgay, B. Integrable deformations of superintegrable quantum circuits, preprint 2022, arXiv:2205.02038 .
  * [GUWZ10] Gütschow, J.; Uphoff, S.; Werner, R. F.; Zimborás, Z. Time Asymptotics and Entanglement Generation of Clifford Quantum Cellular Automata, J. Math. Phys. 2010, 51, 015203.
  * [VZP18] Vanicat, M.; Zadnik, L.; Prosen, T. Integrable Trotterization: Local Conservation Laws and Boundary Driving, Phys. Rev. Lett. 2018, 121, 030606.
  * [BKP19] Bertini, B.; Kos, P.; Prosen, T. Exact Correlation Functions for Dual-Unitary Lattice Models in 1+1 Dimensions, Phys. Rev. Lett. 2019, 123, 210601.
  * [BP22] Borsi, M.; Pozsgay, B. Construction and the ergodicity properties of dual unitary quantum circuits, Phys. Rev. B 2022, 106, 014302.
  * [CFVV84] Casati, G.; Ford, J.; Vivaldi, F.; Visscher, W. M. One-Dimensional Classical Many-Body System Having a Normal Thermal Conductivity, Phys. Rev. Lett. 1984, 52, 1861.
  * [PR92] Prosen, T.; Robnik, M. Energy transport and detailed verification of Fourier heat law in a chain of colliding harmonic oscillators, J. Phys. A: Math. Gen. 1992, 25, 3449.
  * [CP03] Casati, G.; Prosen, T. Anomalous heat conduction in a one-dimensional ideal gas, Phys. Rev. E 2003, 67, 015203(R).
  * [GHKV18] Gopalakrishnan, S.; Huse, D. A.; Khemani, V.; Vasseur, R. Hydrodynamics of operator spreading and quasiparticle diffusion in interacting integrable systems, Phys. Rev. B 2018, 98, 220303(R).


© 20xx IEEE. Personal use of this material is permitted.Permission from IEEE
must be obtained for all other uses, in any current or future media, including
reprinting/republishing this material for advertising or promotional purposes,
creating new collective works, for resale or redistribution to servers or
lists, or reuse of any copyrighted component of this work in other works.

# Modified Auto Regressive Technique for Univariate Time Series Prediction of
Solar Irradiance

Umar Marikkar, A. S. Jameel Hassan, Mihitha S. Maithripala
Roshan I. Godaliyadda, Parakrama B. Ekanayake, Janaka B. Ekanayake
_Department of Electrical and Electronic Engineering, Faculty of Engineering_
_University of Peradeniya, Peradeniya 20400, Sri Lanka_
{e14219, jameel.hassan.2014<EMAIL_ADDRESS>
{roshangodd, mpb.ekanayake<EMAIL_ADDRESS>
###### Abstract

The integration of renewable resources has increased in power generation as a
means to reduce the fossil fuel usage and mitigate its adverse effects on the
environment. However, renewables like solar energy are stochastic in nature
due to its high dependency on weather patterns. This uncertainty vastly
diminishes the benefit of solar panel integration and increases the operating
costs due to larger energy reserve requirement. To address this issue, a
Modified Auto Regressive model, a Convolutional Neural Network and a Long
Short Term Memory neural network that can accurately predict the solar
irradiance are proposed. The proposed techniques are compared against each
other by means of multiple error metrics of validation. The Modified Auto
Regressive model has a mean absolute percentage error of 14.2%, 19.9% and
22.4% for 10 minute, 30 minute and 1 hour prediction horizons. Therefore, the
Modified Auto Regressive model is proposed as the most robust method,
assimilating the state of the art neural networks for the solar forecasting
problem.

###### Index Terms:

Univariate irradiance forecast, Auto Regression, Neural networks

##  I Introduction

Over the years, penetration of renewable energy has vastly increased in the
electricity grid network. The increase in energy demand, adverse effects of
fossil fuel generation and awareness towards climate change has advanced the
use of renewable resources [1]. The growing concern towards environmental
pollution has rendered the sustainable development goal #7 aiming to develop
cleaner energy sources by the United Nations environment program [2].

The increased intervention by renewable resources such as solar and wind pose
a highly volatile problem due to the intermittent nature of power generation.
In comparison, solar power poses a bigger problem than wind power due to
higher fluctuations due to cloud covers and effects of instant weather
changes. There is a critical need in improving the real time solar irradiance
forecast since the core operations of the utility and Independent System
Operators (ISO) depends on the power generation capacity. For instance, the
day ahead prediction of solar irradiance will make the Unit Commitment (UC)
more efficient whereas the improvement in short term forecast will reduce
errors due to fluctuations in solar power and minimise the strain on the grid.
This indicates that a higher prediction accuracy is useful in many timescales.

The integration of renewable energy is inevitable in the ever growing energy
demand and is a huge contribution to the smart grid evolution. For a viable
progress of the smart grid, this integration needs to be addressed in terms of
its stakeholders; the ISO and Independent Power Producers (IPP). In the energy
market, the ISO can dispatch and share resources amongst IPPs for the energy
bidding process. In order to optimally utilise the power production, the IPPs
need to minimise their error of power production forecast. A lower forecast
error, gives the IPP a larger window of bidding during the intra-day energy
bidding process in light of the solar power prediction. This intra-day bidding
ranges in window time from 15 minutes to a couple of hours. However, failure
to produce the bid power will be settled in terms of cumbersome penalties for
the energy bidder [3]. Furthermore, these predictions are seldom deployed in-
house which provide the specific forecast horizon depending on the necessity
of the power producer. The critical feature for short-term prediction schemes
is the time constraint in data acquisition and processing [4]. Having a very
narrow window needs to ensure that the error is minimised by a very quick
operating forecast technique. Thus, the robustness of the forecast scheme will
depend not only on its accuracy of prediction but also on the ease of
deployability having an instant processing capability.

The solar prediction scheme can be categorised into long term, short-term and
very short-term predictions. In the literature a day ahead solar irradiance
prediction has been performed using neural networks considering the ground
sensor data and weather information as input in [5]. A combined neural network
approach using Convoluted Neural Networks (CNN) and Long Short Term Memory
(LSTM) is used for day ahead prediction in [6]. A day ahead probabilistic
Photo Voltaic (PV) power forecasting method based on auantile CNN based on
feature extraction is discussed in [7]. A novel CNN framework using genetic
algorithm and particle swarm optimisation for the hyper parameter tuning is
presented in [8] using multiple meteorological data. Meanwhile, very short-
term prediction ranging from few minutes to 6 hours is performed using
multiple techniques including deep learning techniques, sky image based
methods and statistical methods for time series data such as Auto Regression
(AR), Auto Regression and Moving Average (ARMA) and Auto Regression Integrated
Moving Average (ARIMA) [9, 10, 11]. Statistical methods incorporating AR,
ARIMA, ARMA have been extensively exercised in time series data in various
other applications in literature [12, 13, 14]. The statistical models AR and
ARIMA are compared to LSTM in [15]. The LSTM model is reported to outperform
the statistical methods here. A hybrid model incorporating a discrete wavelet
transform, ARMA and a Recurrent Neural Network (RNN) is implemented for 1
minute ahead prediction in [16] showing considerable improvement in precision
of prediction.

In the literature, very short-term prediction of solar irradiance using deep
learning techniques is consolidated to outperform others; specifically the
LSTM neural network. However, a simplistic statistical forecast technique will
be more robust in deployment due to the time constraint in very short-term
prediction requirements. Therefore, in comparison to the deep learning
techniques the statistical methods are convenient in computational expense.
Furthermore, considering the above contributions, the usage of multiple data
sources extensively is a major drawback in the ease of in-house implementation
of the forecast schemes due to the scarcity of resources. In this paper, we
propose a Modified Auto Regressive (MAR) approach, a CNN model and LSTM model
for univariate solar prediction. The three models are compared across multiple
error metrics for validation. Finally we propose the MAR as the best approach
as it is able to assimilate its performance to the LSTM model for multiple
prediction horizons as verified in this paper for 10 minute, 30 minute and 1
hour horizons using only the past irradiance measurements as inputs. This
ensures a highly robust model that is easily deployable in-house, for real
time very short-term predictions.

##  II Data Preparation

###  II-A Study Area

The solar irradiance data was obtained from the PV plant stationed at the
Faculty of Engineering in University of Peradeniya in Sri Lanka. The city of
Peradeniya is located surrounded by the hills of Hantana with a tropical
climate. This results in fluctuations of the solar irradiance curve rather
than yielding the typical _"bell"_ shaped curve. This setting gives a more
challenging data set which highly reflects the volatile nature of solar
irradiance in contrast to data sets often encountered in the literature. The
data is collected for a period of one year with data points at every 10 minute
interval.

###  II-B Training/Testing Split and Data Standardisation

For all forecasting models, the training/testing data split is divided as
70/30$\%$ considering conventional deep learning practice. As the collected
data spans a whole year, this gives a sufficiently large data-set
($\approx$110 days) for testing.

For an efficient training and forecast performance, the input data is
standardised as in equation (1) as a pre processing step, and de-standardised
in the post processing stage.

$z=\frac{x-\mu}{\sigma}$ (1)

where,
$z$ = Normalised signal value $x$ = Irradiance level at each timestamp $\mu$ =
Mean of the dataset $\sigma$ = Variance of the dataset

##  III Methodology

The short-term prediction of solar irradiance is implemented for the time
horizons of 10 minute, 30 minute and 1 hour intervals. The forecasting schemes
are developed using machine learning techniques in terms of Convolutional
Neural Networks (CNN), Long-Short Term Memory (LSTM) networks, and in addition
a Modified Auto Regressive (MAR) forecast model is implemented. Out of the
three techniques, the MAR approach is highlighted as the best model for solar
prediction.

###  III-A Convolutional Neural Network (CNN)

CNNs are a type of neural networks most prominently deployed in image
processing and classification problems. The key notion of CNN is its ability
to learn abstract features and recognise them during training by means of
kernels in the network [17]. Therefore, in this paper, the CNN has been
employed in a time series model to identify the temporal abstract level
features in order to predict the next time step.

In order to encapsulate the complex modelling of features, the CNN utilises
three separate layers, namely: the convolutional layer, the pooling layer and
the fully connected layer. The convolutional layer is responsible for the
identification of relationships between the inputs in the locality of the
convolution operation that takes place between inputs and the kernels. The
pooling layer performs a down-sampling of the output from the convolution
operation. This is then fed to a fully connected layer which is responsible
for predicting the output depending on the features. A series of convolution-
pooling layers can be used if necessary.

In this paper, one convolution layer and an average pooling layer is used.
These layers are designed to extract the feature on the input variables, which
is the past 4 samples (selected as in section III-D2) of the time series
sequence, as in equation (2).

$h^{k}_{ij}=g((W^{k}*x)_{ij}+b_{k})\\\ $ (2)

where, $W^{k}$ is the weight of the kernel connected to $k^{th}$ feature map,
$g$ is the activation function and $b_{k}$ is the bias unit.

The Rectified Linear Unit (ReLU) function is used as the activation function
after evaluating the performance against other activation functions. The ReLU
function is defined by equation (3).

$g(x)=max(0,x)$ (3)

The Adam optimisation algorithm is used as the training function which is an
efficient implementation of the gradient descent algorithm [18]. Finally, two
dense - fully connected layers are implemented following the pooling layer
with one final dense layer with a single neuron that outputs the prediction.
An abstraction of the CNN architecture implemented is shown in Fig 1. The
hyper parameters of the model are chosen by optimisation of a grid search
algorithm as highlighted in Table I.

####  III-A1 Pre processing and Post processing stages

The solar irradiance curve has a trend of the $"bell"$ shape to it. In order
to remove this trend in the input data, pre processing is performed at the
input stage. In addition to the data standardisation described in II-B, a
difference transform of lag 1 is performed to the input signal after
standardisation. The transformed input is fed to the CNN and the predicted
signal is obtained. The predicted signal is passed in a post processing stage
to reconstruct the solar irradiance curve as predicted. The pre processing
difference transform and post processing reconstruction equations are given in
equation (4).

$\begin{split}X&=[x_{0},x_{1},\dots,x_{n}]\\\
\tilde{X}&=[(x_{0}-0),(x_{1}-x_{0}),\dots,(x_{n}-x_{n-1})]\\\
\tilde{Y}&=[\tilde{y_{0}},\tilde{y_{1}},\dots,\tilde{y_{n}}]\\\
Y&=[(\tilde{y_{0}}+0),(\tilde{y_{1}}+x_{0}),\dots,(\tilde{y_{n}}+x_{n-1})\\\
\end{split}$ (4)

here,
$X$ = Normalised signal value $\tilde{X}$ = Difference transformed input
$\tilde{Y}$ = Predicted signal $Y$ = Reconstructed predicted signal value

Figure 1: CNN architecture

###  III-B Long-Short Term Memory Neural Network (LSTM)

The LSTM network is a type of Recurrent Neural Networks (RNN), used for time
series prediction. A major drawback of RNN is the inability to capture long-
term dependencies in a signal, due to memory constraints. The LSTM cell has a
selective storage of trends in its memory, hence it ignores repetitive
information. The cell state is defined by which information is stored in or
discarded. This is controlled by means of three gates; the input gate $i_{t}$,
output gate $O_{t}$ and forget gate $f_{t}$. The output of the LSTM networks
depends on the current input and the cell state [19]. The working mechanism
and cell architecture of the LSTM network is shown in Fig. 2 and Fig. 3
respectively.

Figure 2: Working mechanism of LSTM cells Figure 3: LSTM cell architecture

At time $t$, the inputs to the network are the sequence vector $X_{t}$, the
hidden state output $h_{t-1}$ and the cell state $C_{t-1}$. The outputs of the
network are, the LSTM hidden state $h_{t}$ and the cell state $C_{t}$. The
forget gate, input gate and output gate are calculated using equations (5),
(6) and (7). Here, $i_{t}$ is the input gate and $O_{t}$ is the output gate.
The forget gate $f_{t}$ is used to update, maintain or delete the cell state
information.

$f_{t}=\sigma(W_{f}\times[h_{t-1},x_{t}]+b_{f})\\\ $ (5)
$i_{t}=\sigma(W_{i}\times[h_{t-1},x_{t}]+b_{i})\\\ $ (6)
$O_{t}=\sigma(W_{O}\times[h_{t-1},x_{t}]+b_{o})\\\ $ (7)

The current candidate cell state $\bar{C}$ is calculated by equation (8), and
is updated to produce the output cell state $C_{t}$ as in equation (9). Using
the output cell state, the current hidden state $h_{t}$ is calculated by
equation (10).

$\bar{C_{t}}=tanh(W_{C}\times[h_{t-1},x_{t}]+b_{c})\\\ $ (8)
$C_{t}=f_{t}\times C_{t-1}+i_{t}\times\bar{C_{t}}\\\ $ (9) $h_{t}=O_{t}\times
tanh(C_{t})\\\ $ (10)

$W_{f}$, $W_{i}$, $W_{O}$, and $b_{f}$, $b_{i}$ , $b_{o}$ and $b_{c}$ are
weights and bias parameters of each gate. $\sigma$ is the sigmoid activation
function.

TABLE I: Hyper-parameter optimization for implemented networks Network Model |  Model Hyper-Parameter Names |  Search Space for Optimal Hyper-Parameters
---|---|---
CNN |  Optimizer |  Adam
|  Learning rate ($\alpha$) |  [0.1 0.01 0.005 0.001]
|  Convolution Layers |  [1 2]
|  Fully connected Layers |  [1 2 3]
|  Hidden Layer Neurons |  [8 16 50 150 300]
|  Number of Kernels |  [3 16 80 150]
|  Kernel Size |  [3 2]
|  Batch Size |  [16 32 64 128 256 512]
|  Pooling Size |  [1 2]
|  Epochs |  [10 30 100 500]
LSTM |  Optimizer |  Adam
|  Initial learning rate ($\alpha$) |  [0.1 0.05 0.01 0.005]
|  Learning rate drop period |  [10 30 100 300]
|  LSTM Layers |  [1 2]
|  LSTM Neurons |  [8 16 32 64 128 256]
|  Fully Connected Layers |  [1 2 3]
|  Hidden Layer Neurons |  [8 16 32 64 128]
|  Epochs |  [50 100 300 500]

###  III-C Network Design for Deep-Learning Models

All simulations are run on an Intel core-i7 @4.5GHz computer. Implemented
deep-learning networks are designed using MATLAB deep learning toolbox.

Neural networks, if poorly trained, leads to over-fitting or under-fitting of
the training data, resulting in disparity between training data prediction and
actual prediction performance. Similarly, bad design of an neural network
architecture could lead to error propagation, high computational cost, or
simply overkill.

Hyper-parameter optimization plays an important role in choosing the optimal
neural network architecture and training parameters. Brute force methods such
as grid search, probabilistic models such as bayesian optimization and random
searches are widely used. As high computational power is available for
training, grid search algorithm was implemented. Initially, a coarse search
was carried out on a large search space as shown in Table I. Then, a fine
search was implemented on a smaller search space. As all hyper-parameters were
well optimized throughout the smaller search space, coarse search hyper-
parameters were chosen, as highlighted in Table I.

###  III-D Modified Auto-Regressive Model (MAR)

In the AR model, the predicted signal value at the next time step is linearly
dependent on observed values at a set number of previous time steps. However,
our proposed model does not work on the standardized irradiance measurements,
but ensemble deducted values, as described in section III-D1. The AR model
equation relating predicted value to the previously observed values is given
by equation (11).

$x_{n,pred}=\sum_{k=1}^{m}w_{k}\times x_{n-k}$ (11)

where,
$m$ = order of the AR model $x_{n,pred}$ = predicted signal value for next
timestamp $w_{k}$ = model weights $x_{n-k}$ = past signal values

####  III-D1 Feature Engineering by Ensemble Deduction

Prior to the prediction, the expected value of $m$ number of past signal
values at each timestamp is deducted from its corresponding irradiance
measurement, as shown in equation (12). This ensures that the periodic nature
of the days, governed by the bell shape curve, is unaffected at the time of
prediction. The ensemble deduction in a given day to predict the 20th
timestamp of the day is illustrated in Fig. 4.

$x_{n-i,ens}=x_{n-i}-E[x_{n-i}]\\\ $ (12)

where,
$i$ = [1, …,$m$] $n$ = prediction timestamp $x_{n-i,ens}$ = ensemble deducted
signal value at $n-i$ $x_{n}$ = actual standardized signal value at $n-i$
$E[.]$ = statistical expectation operator

Figure 4: Ensemble Deduction to predict 20th timestamp

####  III-D2 Parameter Optimisation

The order of the model ($m$) depends on the Partial Auto-Correlation Function
(PACF) of the given data. The PACF provides the correlation between a fixed
time series value $x_{n}$ and its lagged values $x_{n-\tau}$ relative to the
fixed value. The equation to compute the PACF is described in equation (13).
Fig. 5 shows a graphical representation of equation (13). As observed, $m=4$
was chosen as the optimal order.

$R_{\tau}=E[x_{n-\tau}\cdot x_{n}]$ (13)

where,
$E[.]$ = statistical expectation operator $R_{\tau}$ = correlation function
$\tau$ = lag of the previous timestamp

Figure 5: Partial Auto-Correlation Function (PACF)

The prediction error $x_{n,pred}-x_{n,real}$ is chosen to calculate model
parameters. They are calculated using optimisation; where a positive,
monotonically increasing error function is minimized. A squared error function
as given by equation (14) exhibits these characteristics. Therefore, the Yule-
Walker equation given by equation (15) is used to calculate model parameters.

$f(e_{n})=(x_{n,pred}-x_{n,real})^{2}$ (14)

where,
$f(e_{n})$ = error function $e_{n}$ = error at a given time step $n$
$x_{n,pred}$ = predicted value at $n$ $x_{n,real}$ = observed value at $n$

$W=(X^{T}X)^{-1}X^{T}Y$ (15)

where,
$W$ = weights matrix $X$ = design matrix (dependent on order $m$) $Y$ = output
matrix ($X_{real}$)

The design matrix $X$ contains the training examples as its rows, and features
for each example as its columns. The number of columns depends on the order
$m$. After optimizing the model parameters, a finite loop is run for each time
step of the day, predicting the signal value $x_{n,pred}$ at the next time
step. To calculate predicted solar irradiance, $x_{n,pred}$ is de-
standardized.

##  IV Results and Discussion

Irradiance prediction for two randomly chosen successive days for 10 minute,
30 minute and 1 hour prediction horizons are shown in Fig. 6, Fig. 7, Fig. 8
respectively. The deep learning models and the MAR model are designed with one
specific model to forecast across all time horizons discussed in the paper. As
observed, when the prediction horizon increases, the tendency for predicted
curves to follow sudden changes in irradiance is less. Therefore, the
dependency on the ensemble (bell-shaped) features of data increases.

TABLE II: Error comparison of Prediction Models Error Model | Horizon | CNN | AR | LSTM | MAR
---|---|---|---|---|---
| 10 min | 113.62 | 115.62 | 114.79 | 110.38
RMSE /$Wm^{-2}$ | 30 min | 164.63 | 170.05 | 146.50 | 148.25
| 1 h | 181.98 | 182.17 | 161.40 | 158.56
| 10 min | 64.66 | 74.52 | 69.70 | 68.21
MAE /$Wm^{-2}$ | 30 min | 102.82 | 123.16 | 98.34 | 99.06
| 1 h | 124.52 | 138.07 | 111.98 | 112.09
| 10 min | 14.56 | 16.02 | 14.71 | 14.20
MAPE /% | 30 min | 21.81 | 23.67 | 19.36 | 19.94
| 1 h | 24.59 | 27.45 | 22.18 | 22.42

The performance of the four models are evaluated with respect to the metrics
of measure Root Mean Square Error (RMSE), Mean Absolute Error (MAE) and Mean
Absolute Percentage Error (MAPE) in Table II. It can be observed that the
error increases for all models when the prediction horizon increases. However,
the performance of the CNN and conventional AR model deteriorates faster than
the other two. It is noteworthy that, the MAR model, being a simplistic
implementation with pre processing, consolidates a robust performance with the
time horizon change while matching the performance of a deep learning LSTM
model in all aspects; both errors and increased time horizons.

##  V Conclusion

In this paper we propose three models of solar prediction; a Modified Auto
Regressive (MAR) model, two deep learning models each based on CNN and LSTM
neural networks. The performance of the models are quantified by the error
metrics RMSE, MAE and MAPE, and it affirms that the MAR model fits best for
the case of very-short term prediction of solar irradiance.

In a system such as a tropical environment, variability of irradiance at a
given timestamp is high, reducing the correlation between consecutive samples.
Hence, deep neural networks tend to mostly capture the bell-shaped nature of
solar irradiance, as intra-day variations are highly uncorrelated. By means of
the ensemble mean curve deduction the MAR, having the least computational
cost, is capable of predicting solar irradiance with a performance similar to
LSTM- the state of the art prediction scheme- across all tested prediction
horizons.

Existing prediction models use multi-sensory data; such as temperature,
humidity, cloud cover and irradiance. The proposed MAR uses a single sensor
measurement as input for the prediction sufficing in performance for most use
cases, with an MAPE of less than 15% for 10 minute prediction, and less than
20% for 30 minute prediction. This enables an easy acquisition of data, which
facilitates an easily deployable

Figure 6: Observed vs. Predicted irradiance for 10 minute prediction horizon
extracted for 1 day from the predicted dataset Figure 7: Observed vs.
Predicted irradiance for 30 minute prediction horizon extracted for 2
consecutive days from the predicted dataset Figure 8: Observed vs. Predicted
irradiance for 1 hour prediction horizon extracted for 2 consecutive days from
the predicted dataset

forecast system. Thus, taking into account the aforementioned conditions, MAR
is chosen as the optimal solar irradiance prediction model.

## References

  * [1] D. Abbot, “Keeping the energy debate clean: How do we supply the world’s energy needs?” _Proceedings of the IEEE_ , vol. 98, no. 1, pp. 42–66, Jan 2010.
  * [2] UNDP, “UNDP strategic plan, 2018–2021. Special session (executive board of the United Nations Development Programme, the United Nations Population Fund and the United Nations office for project services),” New York: United Nations; 2017.
  * [3] A. Kaur, L. Nonnenmacher, H. Pedro, and C. Coimbra, “Benefits of solar forecasting for energy imbalance markets,” _Renewable Energy_ , vol. 86, pp. 819–830, Feb 2016.
  * [4] D. U. F. Malte Zieher, Dr. Matthias Lange, “Varaiable Renewable Energy Forecasting- Integration into Electricity Grid and Markets- A Best Practice Guide,” Federal Ministry for Economic Corporation and Development, Germany, Tech. Rep., 2015.
  * [5] Y. Wang, Y. Shen, S. Mao, X. Chen, and H. Zou, “LASSO and LSTM integrated temporal model for short-term solar intensity forecasting,” _IEEE Internet of Things Journal_ , vol. 6, no. 2, pp. 2933–2944, April 2019\.
  * [6] W. Lee, K. Kim, J. Park, J. Kim, and Y. Kim, “Forecasting solar power using Long-Short Term Memory and Convolutional Neural Networks,” _IEEE Access_ , vol. 6, pp. 73 068–73 080, 2018.
  * [7] Q. Huang and S. Wei, “Improved Quantile Convolutional Neural Network with two-stage training for daily-ahead probabilistic forecasting of photovoltaic power,” _Energy Conversion and Management_ , vol. 220, p. 113085, Sep 2020.
  * [8] Dong _et al._ , “A novel Convolutional Neural Network framework based solar irradiance prediction method,” _International Journal of Electrical Power & Energy Systems_, vol. 114, p. 105411, Jan 2020.
  * [9] D. M. L. H. Dissawa _et al._ , “Cross-correlation based cloud motion estimation for short-term solar irradiation predictions,” in _2017 IEEE International Conference on Industrial and Information Systems (ICIIS)_ , 2017, pp. 1–6.
  * [10] M. Cervantes, H. Krishnaswami, W. Richardson, and R. Vega, “Utilization of low cost, sky-imaging technology for irradiance forecasting of distributed solar generation,” in _2016 IEEE Green Technologies Conference (GreenTech)_ , April 2016, pp. 142–146.
  * [11] A. Moreno-Munoz, J. J. G. de la Rosa, R. Posadillo, and F. Bellido, “Very short term forecasting of solar radiation,” in _IEEE Photovoltaic Specialists Conference_ , San Diego, CA, USA, May 2009.
  * [12] M. Rafi, M. T. Wahab, M. Bilal Khan, and H. Raza, “Atm cash prediction using time series approach,” in _2020 3rd International Conference on Computing, Mathematics and Engineering Technologies (iCoMET)_ , 2020, pp. 1–6.
  * [13] H. Matsila and P. Bokoro, “Load forecasting using statistical time series model in a medium voltage distribution network,” in _IECON 2018 - 44th Annual Conference of the IEEE Industrial Electronics Society_ , 2018, pp. 4974–4979.
  * [14] N. Chen-xu and W. Jie-sheng, “Auto Regressive Moving Average (ARMA) prediction method of bank cash flow time series,” in _2015 34th Chinese Control Conference (CCC)_ , 2015, pp. 4928–4933.
  * [15] W. Ji, C. Chan, J. Loh, F. Choo, and L. Chen, “Solar radiation prediction using statistical approaches,” in _2009 7th International Conference on Information, Communications and Signal Processing (ICICS)_ , 2009, pp. 1–5.
  * [16] H. Nazaripouya, B. Wang, Y. Wang, P. Chu, H. R. Pota, and R. Gadh, “Univariate time series prediction of solar power using a hybrid wavelet-ARMA-NARX prediction method,” in _2016 IEEE/PES Transmission and Distribution Conference and Exposition (T D)_ , 2016, pp. 1–5.
  * [17] S. Oehmcke, O. Zielinski, and O. Kramer, “Input quality aware Convolutional LSTM networks for virtual marine sensors,” _Neurocomputing_ , Nov 2017\.
  * [18] L. Wen, K. Zhou, S. Yang, and X. Lu, “Optimal load dispatch of community microgrid with deep learning based solar power and load forecasting,” _Energy_ , vol. 171, 03 2019.
  * [19] H. Zhou, Y. Zhang, L. Yang, Q. Liu, K. Yan, and Y. Du, “Short-term photovoltaic power forecasting based on Long Short Term Memory neural network and attention mechanism,” _IEEE Access_ , vol. 7, pp. 78 063–78 074, 2019.


# A Generic Bundle Forwarding Interface

Felix Walter *  | D3TN GmbH
Dresden, Germany
<EMAIL_ADDRESS>
###### Abstract

A generic interface for determining the next hop(s) for a DTN bundle is a
valuable contribution to DTN research and development as it decouples the
topology-independent elements of bundle processing from the topology-dependent
forwarding decision. We introduce a concept that greatly increases flexibility
regarding the evaluation and deployment of DTN forwarding and routing
techniques and facilitates the development of software stacks applicable to
heterogeneous topologies.

##  I Motivation

A universal deployment of _Delay- and Disruption-tolerant Networking (DTN)_
will consist of a huge number of diverse, heterogeneous (sub-)networks, such
as networks with:

  * •

persistent low-latency end-to-end connectivity,

  * •

short or long-latency, asymmetric or unidirectional, but precisely-scheduled
_contacts_ between the nodes,

  * •

probabilistic encounters (e.g., if nodes move somewhat randomly but stay in
local vicinity most of the time),

  * •

or totally random connectivity.

Further variations of these characteristics can be imagined. To prevent the
need to develop specialized implementations of the DTN protocols for every
(sub-)network, we advocate for a flexible interface attached to common core
infrastructure.

##  II Generic Interface Concept

As sketched in Figure 1, our concept introduces a dedicated component
integrating with the _Bundle Protocol Agent (BPA)_ , which we call the _Bundle
Dispatcher Module (BDM)_. Communication occurs via a network or IPC socket
using a low-overhead data exchange format with wide compatibility such as
Protobuf or Cap’n Proto.

Figure 1: High-level overview of the Bundle Dispatcher Module concept as
extension to Figure 2 from the DTN Bundle Protocol specification, RFC 9171
[1]. The concept will be implemented in µD3TN [2].

###  II-A Event-based Bundle Processing

To achieve a loose coupling of the components and increased flexibility with
respect to the BDM implementation, the interface leverages an event-based
approach: The BPA posts information about changes in connectivity, bundles
that require a forwarding decision, bundles expiring, etc. via a _publish-
subscribe_ interface. A module subscribing to these _events_ can then decide
a) whether or not to act on them and b) invoke behavior to realize the
forwarding decision. This approach further provides the flexibility to attach
multiple _subscribers_ beside the BDM, e.g., for monitoring purposes.

To pass information about bundles to the BDM, a _Bundle Metadata_ data
structure is defined, containing the bundle headers, extension blocks, and
further metadata like timestamps, but excluding the payload for performance
reasons.

###  II-B Bundle Processing Actions

As counterpart to the event-based information flow from the BPA to the BDM, a
_Remote Procedure Call (RPC)_ mechanism is provided by the BPA, which offers
functions for updating the stored bundles and associated forwarding decisions.
In this context, we propose to attach a list of _Bundle Processing Actions_ to
each bundle, that defines in-order what should be done next with the given
bundle and can be updated by the BDM. At least two actions must be supported
by every BPA:

  * •

SendTo(node): Forward the bundle to the specified next-hop node. Note that
this action may occur multiple times or contain a multicast identifier.

  * •

Drop: Remove the “Forward pending” retention constraint if the previous action
was successful.

This list can get implementation-dependent extensions, e.g., actions to
fragment bundles, control their storage, and so on. Like in some Software-
Defined Networking implementations, there might be a function to query the
supported actions for announcing such extensions in a straight-forward manner.

The BPA executes the action lists beginning with the bundle that was updated
first, which allows the BDM to control the order of bundle forwarding. It
should be noted that an empty action list is a valid state: In this case, the
BPA will keep the bundle in memory until its lifetime expires or the BDM
triggers an update on the action list. For maximum flexibility, the action
list can have a configurable default assigned to bundles upon reception, thus,
even allowing BPA operation without a dedicated BDM in simple cases (e.g., if
there is a persistent “default gateway” contact).

###  II-C BPA-BDM Interaction

The two complementary mechanisms provide for a flexible interaction between
BPA and BDM and the implementation of basically arbitrary forwarding
techniques. Figures 2 and 3 depict how two possible example cases would be
handled.

Figure 2: Example sequence diagram for injecting a new Application Data Unit
(ADU) addressed to node Z for which the next hop is node Y, assuming
opportunistic single-copy forwarding in the BDM (simplified). Figure 3:
Example sequence diagram for receiving and forwarding a bundle addressed to
node Z for which the next hop is node Y, assuming contact-based routing in the
BDM (simplified).

###  II-D Inter-BDM Communication

Some forwarding techniques require the exchange of specific data, e.g.,
PRoPHET [8] needs to distribute delivery predictability values between nodes.
A BDM implementation has three options to exchange such data: a) attach it to
bundles as extension blocks (the used metadata structure contains them), b)
inject bundles on its own by registering as a BP application, or c) use an
independent mechanism or channel.

##  III Related Work

The presented concept enables a DTN software stack that can adapt to
heterogeneous topologies, for which some alternative approaches exist. _DTN2_
[3] implements bundle forwarding in dedicated C++ classes that can be
triggered by over 60 different events and can access the bundle data structure
in memory. The _Interplanetary Overlay Network (ION)_ [4] consists of multiple
daemons that access a shared data structure. The forwarding decision is made
in a separate daemon, which is chosen depending on the destination endpoint
scheme. _ProgDTN_ [5] is a novel approach executing a JavaScript program from
within the BPA that represents a function returning the forwarding decision
and can be flexibly exchanged. Overall, however, none of these implementations
provides as much flexibility to implement and exchange the forwarding
component as the approach introduced here.

##  IV Summary and Outlook

This paper outlined a clean, low-overhead, event-driven socket interface to
support a dedicated bundle forwarding component that can be flexibly exchanged
and, thus, provides adaptability to heterogeneous topologies. At the time of
writing, a combined Rust (of the BPA) and Python 3 (of the BDM) implementation
in µD3TN [2] is almost complete and the integration of techniques such as
Schedule-Aware Bundle Routing (SABR) [6] and IP-based Neighbor Discovery
(IPND) [7] is ongoing. The pending evaluation of this interface will be based
on an extended field test using drones and a satellite link to emulate a
heterogeneous DTN internetwork.

## Acknowledgments

The presented concept stems from a long series of discussions in which most of
the team at D3TN was involved. Specifically, without the extensive
contributions of the following people, the concept would not be as advanced:
Marius Feldmann, Juan Andres Fraire, Tobias Nöthlich, and Georg Alexander
Murzik. The concept is developed as part of the REDMARS2 project that is
funded by Germany’s Federal Ministry of Education and Research (FKZ16KIS1356).

## References

  * [1] S. Burleigh et al., “Bundle Protocol Version 7,” _RFC 9171_ , IETF, 2022.
  * [2] “µD3TN,” [online] Available: https://d3tn.com/ud3tn.html.
  * [3] “DTN Reference Implementation (DTN2),” [online] Available: https://github.com/delay-tolerant-networking/DTN2.
  * [4] “Interplanetary Overlay Network (ION),” [online] Available: http://sourceforge.net/projects/ion-dtn/.
  * [5] M. Sommer et al., “ProgDTN: Programmable Disruption-tolerant Networking,” in _International Conference on NETworked SYStems_ , Online, 2022. Available: https://dtn7.github.io/assets/sommer2022progdtn.pdf
  * [6] “Schedule-Aware Bundle Routing,” _CCSDS 734.3-B-1_ , Washington, DC, USA, 2019.
  * [7] D. Ellard et al., “DTN IP Neighbor Discovery,” _Internet-Draft_ , IRTF, 2015. https://datatracker.ietf.org/doc/html/draft-irtf-dtnrg-ipnd-03
  * [8] A. Lindgren et al., “Probabilistic routing protocol for intermittently connected networks,” _RFC 6693_ , IRTF, 2012.