# Quantum gravity phenomenology and the blackbody radiation

R. Turcati<EMAIL_ADDRESS>Centro Brasileiro de Pesquisas Físicas, Rua Dr.
Xavier Sigaud, 150, URCA, Rio de Janeiro CEP 22290-180, RJ, Brazil  I. Soares
<EMAIL_ADDRESS>Centro Brasileiro de Pesquisas Físicas, Rua Dr. Xavier Sigaud,
150, URCA, Rio de Janeiro CEP 22290-180, RJ, Brazil  S. B. Duarte<EMAIL_ADDRESS>Centro Brasileiro de Pesquisas Físicas, Rua Dr. Xavier Sigaud, 150, URCA, Rio
de Janeiro CEP 22290-180, RJ, Brazil

###### Abstract

We analyze the problem of blackbody radiation in the presence of quantum
gravity effects encoded in modified dispersion relations. In this context, the
spectral radiance and the generalized Stefan-Boltzmann law are studied.
Furthermore, the regime of low temperatures is contemplated as well, where
features related to the blackbody thermal laws and the thermodynamics
quantities such as energy, pressure, entropy, and specific heat are obtained.
Possible implications in compact objects such as neutron stars are also
discussed.

##  I Introduction

The nature of spacetime is perhaps the main challenge in modern physics. It is
commonly accepted that at the Planck scale, i.e., for energies near
$10^{19}GeV$, our physical theories cease to be predictive, and quantum
gravity effects should become relevant Gross and Mende (1988); Maggiore
(1993); Addazi  _et al._ (2022). Nevertheless, such an energy scale is far
beyond our technological limitations, which has hampered any direct
observation of these relic signatures.

On the other hand, in the past years, it has been suggested that quantum
gravity effects may induce small deviations in the standard predictions, which
could be observable, at least in principle, through astronomical observations
and ground experiments Amelino-Camelia (2000); Addazi  _et al._ (2022). The
possibility of observing these phenomena includes neutrino oscillations, gamma
rays, laser interferometry, and CPT violation, among others Kostelecky (2004).

The framework where these phenomena can be investigated is the so-called
Quantum Gravity Phenomenology. This formalism includes several tests such as
prints on the initial cosmological perturbations, black holes related to extra
dimensions, and Planck scale spacetime fuzziness. Features related to
violations of fundamental symmetries like the Lorentz and CPT symmetries are
also contemplated in these scenariosKostelecky (2004).

Among the several ways to probe these signatures arising from the quantum
gravity, a major role is played through modifications of the dispersion
relationAmelino-Camelia  _et al._ (1998); Ellis  _et al._ (1999). Modified
dispersion relations together with higher-order derivatives field theories are
common aspects of several phenomenological approaches to quantum gravity.
Indeed, in the framework of quantum gravity phenomenology, modifications of
the relativistic dispersion relation associated with free particles are
expected at the Planck scale, which can shed some light on the underlying
quantum theory of gravity.

On the other hand, it is expected that the modifications of the relativistic
dispersion relations be suppressed by the Planck scale. In this vein, one
might wonder if these effects are really important and why consider them after
all. From the cosmology perspective, trans-Planckian effects should have
affected the early-stage evolution of the universe, which could have led to
fingerprints in the cosmic microwave background. Furthermore, the thermal
emission of black holes is associated with modes of arbitrarily high energies
near the horizon, which have been based on modified dispersion relations, and
have provided remarkable insights with regards to the Hawking
radiationJacobson (1991). In addition, Lorentz violation effects in
perturbative quantum field theories could be enhanced in the mechanism of
regularization and renormalization.

Here, we will adopt a pragmatic point of view and try to understand the
eventual consequences of modified dispersion relations in the blackbody
radiation laws and the related thermodynamic quantities.

The structure of this paper is organized as follows. In Sec. (2) modified
dispersion relations are introduced and the grand-canonical formalism is
employed for these dispersion relations. In Sec. (III), the spectral radiance
and the Stefan-Boltzmann law in the present context are derived. The blackbody
thermal laws and the thermodynamic quantities in the regime of low
temperatures are contemplated in Secs (IV) and (IV.1). Our final remarks and
further perspectives can be found in Sec. (V).

##  II Modified dispersion relation and the statistical mechanics

Motivated by the idea that quantum gravity effects may induce modifications in
the dispersion relation Amelino-Camelia  _et al._ (2003); Alfaro and Palma
(2003); Amelino-Camelia (2005); Alfaro  _et al._ (2004); Magueijo and Smolin
(2005); Amelino-Camelia  _et al._ (2005), we consider the general dispersion
relation of the following form:

$\displaystyle
c^{2}\mathbf{p}^{2}=E^{2}-m^{2}c^{4}+\sum_{n=1}^{\infty}\lambda_{n}E^{(n+2)},$
(1)

where $\lambda_{n}$ are factors with the dimension of the inverse of energy,
where its value depends on the specific quantum gravity model, and we have
assumed a Taylor expansion for $E\ll 1/\lambda_{n}$. In addition, whenever
$\lambda_{n}\to 0$, one recovers the standard relativistic dispersion
relation.

Since we are considering the gauge sector of the electromagnetism, we will
then assume a zero rest mass. Furthermore, we will focus our attention on the
first order correction, which reduces (1) to

$\displaystyle c^{2}\mathbf{p}^{2}=E^{2}+\lambda_{1}E^{3}.$ (2)

At this point, some comments are in order.

  1. 1.

For $\lambda_{1}$ positive-definite, the momentum $p$ is always positive.

  2. 2.

$\lambda_{1}$ negative and $|\lambda|<1/E$ provides a physical branch for the
particle.

  3. 3.

Whenever $\lambda_{1}>1/E$ for negative values of $\lambda_{1}$, then the
momentum is complex and the electromagnetic wave will decay in space.

Since the quantum gravity effects are introduced through modifications of the
relativistic dispersion relations, one could ask what are the consequences to
the statistical mechanics systems. In this sense, it is important to stress
the fact that the foundations of statistical mechanics are preserved in this
situation. However, there will be modifications in the density of accessible
states.

Let us then see how the density of states is modified in the present scenario.
To begin with, we use the grand-canonical ensemble to describe the photon gas.
For a system of $N$ massless bosons with energy spectrum $\beta\hbar{w}_{i}$,
we have that each state is labeled by $i$ $\left(i=1,2,...\right)$.

The connection between thermodynamics and statistical mechanics is given by
the gran-canonical potential,

$\displaystyle\phi=-\frac{1}{\beta}log\mathscr{Z},$ (3)

where the logarithmic of the partition function takes the following form:

$\displaystyle
log\mathscr{Z}=-\sum_{i}log\left(1-e^{-\beta\hbar{w}_{i}}\right).$ (4)

The average energy is then given by

$\displaystyle\langle{E}\rangle=\frac{1}{V}\sum_{i}\frac{\hbar{w}_{i}}{e^{\beta\hbar{w_{i}}}-1}.$
(5)

In the large-volume limit,

$\displaystyle\sum_{i}\rightarrow\frac{V}{\pi^{2}c^{3}}\int_{0}^{\infty}dww^{2}g\left(\lambda_{1};w\right),$
(6)

where the function $g\left(\lambda_{1};w\right)$ holds the quantum gravity
corrections and is given by

$\displaystyle
g\left(\lambda_{1};w\right)=\left(1+\lambda_{1}\hbar{w}\right)^{3/2}+\left(\frac{\lambda_{1}\hbar}{2}\right)w\sqrt{1+\lambda_{1}\hbar{w}}.$
(7)

From the above relation, one promptly notes that the density of accessible
states has changed from the standard result for photons.

According to the standard prescription, the number $N$ of modes in the momenta
interval $\left[k,k+dk\right]$ is

$\displaystyle g\left(k\right)dk=\frac{V}{2\pi^{2}}k^{2}dk.$ (8)

By performing a straightforward computation, and taking the modified
dispersion relation (2) into account, we promptly find in the frequency
$w-$space,

$\displaystyle
g\left(w\right)dw=\frac{Vw^{2}}{\pi^{2}c^{3}}g\left(\lambda_{1};w\right).$ (9)

The modified dispersion relation (8) depends only on the magnitude of the
momentum vector, i.e., $|\mathbf{p}|=p=\hbar{k}$. Therefore, the integration
in the phase space is trivial and yields $V$. In addition, the density of
momentum states remains the same since the periodic conditions are not
changed. From relation (6), becomes clear that the dispersion relation (8)
modifies the density of accessible states, which can be seen in relation (9).
As a consequence, there will be an increase in photons for each frequency
mode.

The classical equipartition theorem, in turn, states that the average energy
for photons in thermal equilibrium at temperature T is $k_{B}T$. Quantum
mechanics, on the other hand, assert that each mode of a photon gas satisfying
the Bose-Einstein statistics has an average energy given by
$\bar{E}=\hbar{w}/(e^{\beta\hbar{w}}-1)$. In the present framework, quantum
gravity effects modify the average energy, which generalizes to

$\displaystyle\bar{E}=\left(\frac{\hbar{w}}{e^{\beta\hbar{w}}-1}\right)g\left(\lambda_{1};w\right).$
(10)

##  III The deformed Planck spectrum

For a gas of noninteracting photons in this situation, the energy density in
thermal equilibrium at temperature $T$ in the frequency interval
$\left[w,w+dw\right]$ is given by

$\displaystyle u\left(w\right)dw$ $\displaystyle=$
$\displaystyle\frac{w^{2}}{\pi^{2}c^{3}}\left[\left(1+\lambda_{1}\hbar{w}\right)^{3/2}+\left(\frac{\lambda_{1}\hbar}{2}\right)w\sqrt{1+\lambda_{1}\hbar{w}}\right]$
(11) $\displaystyle\times\left(\frac{\hbar{w}}{e^{\beta\hbar{w}}-1}\right)dw.$

Note that whenever $\lambda_{1}\rightarrow 0$, one recovers the Planck
distribution. To get a better idea of the changes induced by the quantum
gravity effects, we plot, in Fig. (1), both the Planck distribution and the
deformed spectrum. To perform such analysis, we have considered $T\approx
1/\lambda_{1}$. There we can see that the quantum gravity effects are relevant
when we deal with high temperatures in the Universe, such as the inflationary
epoch temperature, for instance, where the temperature reaches $T=10^{13}GeV$.

Figure 1: Graph of the spectral energy density distribution for the quantum
gravity model under consideration for $T=10^{13}GeV$ and
$\lambda_{1}=10^{-13}GeV$. Here, we adopted $c=\hbar=k_{B}=1$. The red line
corresponds to the Planck spectrum, while the blue one is associated with the
quantum gravity model.

Usually, one would expect that quantum gravity corrections at the standard
thermodynamics results would be neglected at low temperatures. On the other
hand, these effects will depend on the value of $\lambda_{1}$. To evaluate if
low-temperature effects could be probed at this regime, let us consider the
Stefan-Boltzmann law for some values of $\lambda_{1}$.

Therefore, taking into account the Stefan-Boltzmann law

$\displaystyle u=\frac{4\sigma}{c}{T}^{4},$ (12)

where $\sigma=\left(2\pi^{5}k_{B}^{4}/15h^{3}c^{2}\right)$$\sigma$ is the
Stefan-Boltzmann constant, we obtain the following curves in Fig (2). For
$\lambda_{1}=10^{-10}eV^{-1}$ (purple line), one promptly notes that for
temperatures near $10^{8}eV$, small deviations from the standard Stefan-
Boltzman law can be observed. For values of $\lambda_{1}$ around
$10^{-10}eV^{-1}$, we can study the low-temperature regime and investigate how
the thermodynamics of the blackbody radiation is affected.

Figure 2: Graph of the energy density vs temperature. Here, we have adopted
$\lambda_{1}=10^{-10}eV^{-1}$ (purple line) and $\lambda_{1}=10^{-13}eV^{-1}$
(blue line) for the quantum gravity model under consideration, where the
curves were obtained through numerical integration. The standard Stefan-
Boltzman law is associated with the red line, which is overlapped by the blue
line. To evaluate the above graph, we have assumed $c=\hbar=k_{B}=1$.

##  IV Quantum gravity corrections at low temperature

Now, we would like to identify the leading quantum gravity correction to the
thermodynamic quantities at low temperatures, i.e., when $\beta\ll 1$. In this
regime, the factor $(e^{\beta\hbar{w}}-1)^{-1}$ decays exponentially, and then
one can expand $g(\lambda_{1},w)$ in Taylor series around zero. To perform
this approximation, the condition $\lambda_{1}k_{B}T\ll 1$ must be satisfied.
For compact objects such as newly formed neutron stars, which hold
temperatures around $10^{11}-10^{12}K$ inside its interior, this condition is
ensured if $\lambda_{1}\ll 10GeV^{-1}$.

Taking into account the mentioned assumptions, therefore, Eq. (11) reduces to

$\displaystyle
u\left(w\right)dw=\frac{w^{2}}{\pi^{2}c^{3}}\left(1+2\lambda_{1}\hbar{w}\right)\left(\frac{\hbar{w}}{e^{\beta\hbar{w}}-1}\right)dw.$
(13)

In the regime of long wavelength, the above frequency distribution takes the
following form:

$\displaystyle
u\left(w\right)=\frac{w^{2}}{\pi^{2}c^{3}}\left(1+2\lambda_{1}\hbar{w}\right)k_{B}T.$
(14)

The total energy, in turn, can now be computed, and give us

$\displaystyle
u\left(T\right)=\frac{\pi^{2}k_{B}^{4}T^{4}}{15c^{3}\hbar^{3}}+\frac{48\lambda_{1}\xi\left(5\right)T^{5}}{\pi^{2}c^{3}\hbar^{3}}.$
(15)

If one introduces the Stefan-Boltzmann constant
$\sigma=\left(2\pi^{5}k_{B}^{4}/15h^{3}c^{2}\right)$ in the above expression,
one promptly gets

$\displaystyle
u\left(T\right)=\frac{4\sigma}{c}{T^{4}}+\frac{48\lambda_{1}\xi\left(5\right)T^{5}}{\pi^{2}c^{3}\hbar^{3}}.$
(16)

For $\lambda_{1}\rightarrow 0$, one recovers the Stefan-Boltzmann law.
Equation (18) can also be written as

$\displaystyle
u\left(T\right)=\frac{4}{c}\sigma_{eff}\left(\lambda_{1};T\right){T^{4}},$
(17)

where

$\displaystyle\sigma_{eff}\left(\lambda_{1};T\right)=\sigma\left(1+\frac{720\lambda_{1}\xi\left(5\right)}{\pi^{4}}k_{B}T\right)$
(18)

is the effective Stefan-Boltzmann constant.

For high temperatures in dense stars, for instance, it is possible to observe
deviations from the Stefan-Boltzmann law.

Figure 3: Graph of the spectral energy density distribution for the quantum
gravity model under consideration for $T=86,17MeV$ and
$\lambda_{1}=10GeV^{-1}$. Here, we adopted $c=\hbar=k_{B}=1$. The red line
corresponds to the Planck spectrum, while the blue one is associated with the
quantum gravity model.

###  IV.1 Thermodynamics in the presence of modified dispersion relations

One of the most important results of modifying the relativistic dispersion
relation is the deformation of thermodynamical quantities. To investigate
this, we take into account the grand canonical potential, which can be used to
derive the pressure, entropy, and specific heat for the blackbody radiation.

For the pressure, we then find

$\displaystyle
p=\frac{\pi^{2}k_{B}^{4}T^{4}}{45c^{3}\hbar^{3}}+\frac{12\lambda_{1}\xi\left(5\right)T^{5}}{\pi^{2}c^{3}\hbar^{3}},$
(19)

while for the entropy and heat capacity densities, we have

$\displaystyle
s=\frac{4\pi^{2}k_{B}^{4}T^{3}}{45c^{3}\hbar^{3}}+\frac{60\lambda_{1}\xi\left(5\right)T^{4}}{\pi^{2}c^{3}\hbar^{3}},$
(20)

and

$\displaystyle
c_{V}=\frac{4\pi^{2}k_{B}^{4}T^{3}}{15c^{3}\hbar^{3}}+\frac{240\lambda_{1}\xi\left(5\right)T^{4}}{\pi^{2}c^{3}\hbar^{3}}.$
(21)

We would like to stress that depending on the sign of $\lambda_{1}$, the above
thermodynamic quantities can increase or decrease, which can influence the
dynamics of compact objects such as neutron stars, for instance. In addition,
for high enough temperatures, departures from the standard results of
thermodynamics should be measurable for these objects.

The equation of state, in turn, can be found by dividing the pressure (19) by
the energy (15), which gives us

$\displaystyle\frac{p}{u}=\frac{1}{3}-\frac{1620\lambda_{1}\xi(5)k_{B}T}{\pi^{4}}.$
(22)

In the limit $\lambda_{1}\rightarrow 0$, we have $p=3u$. The equation of state
is a function of the energy, i.e., $p=p(u)$. To eliminate the temperature
dependence in the above expression, one needs to expand the temperature as
$T=T_{0}+\lambda_{1}T_{1}$ and apply it recursively in the energy density
(18). Therefore, one promptly finds for the temperature the following form:

$\displaystyle
T=\left(\frac{15\hbar^{3}c^{3}}{\pi^{2}}\right)^{1/4}\frac{u^{1/4}}{k_{B}}+\frac{1620\xi(5)}{\pi^{5}}\left(15\hbar^{3}c^{3}\right)^{1/2}\frac{\lambda_{1}}{k_{B}}u^{1/2}.$

Inserting relation (IV.1) in (22), we have, up to the linear order in
$\lambda_{1}$,

$\displaystyle
p=\frac{u}{3}\left[1-\frac{540\lambda_{1}\xi(5)}{\pi^{4}}\left(\frac{15\hbar^{3}c^{3}}{\pi^{2}}\right)^{1/4}\lambda_{1}u^{1/4}\right].$
(24)

From the above equation, one promptly notes the emergence of quantum gravity
effects in the standard equation of state of the order
$\lambda_{1}(\hbar{c})^{3/4}u^{1/4}$.

##  V Final Remarks

In this work, we have considered the effect of quantum gravity effects in the
thermodynamics of blackbody radiation. Starting from a modified dispersion
relation in energy up to the leading order correction, we have explored the
behavior of the emission spectrum in both high and low temperatures. With
regards to the high-temperature regime, the deformed Planck distribution and
the generalized Stefan-Boltzmann law were derived. In the low-temperature
limit, on the other hand, the blackbody radiation laws and the thermodynamics
properties were obtained.

Some other theories have considered higher-order terms in relation (1). In
Ref. Nozari and Sefidgar (2008), for instance, based on arguments of black
hole thermodynamics, only even powers of the energy in the modified dispersion
relation were considered. Here, on the other hand, we have considered the
leading order term in the mentioned expansion. The presence of higher-order
terms would lead to different temperature dependence in the thermodynamics
quantities. Furthermore, the profile of the blackbody spectrum would also be
modified by these terms.

As a future prospect, we intend to extend our formalism here developed to
massive particles in the fermion sector and, then, study the consequences of
the hydrostatic equilibrium of compact objects such as neutron stars.

##  VI Acknowledgements

This work is a part of the project INCT-FNA proc. No. 464898/2014-5. RT
acknowledges financial support from the PCI program of the Brazilian agency
Conselho Nacional de Desenvolvimento Científico e Tecnológico – CNPq. SBD
thanks CNPq for partial financial support.

## References

  * Gross and Mende (1988) D. J. Gross and P. F. Mende, Nucl. Phys. B 303, 407 (1988).
  * Maggiore (1993) M. Maggiore, Phys. Lett. B 304, 65 (1993), arXiv:hep-th/9301067 .
  * Addazi  _et al._ (2022) A. Addazi _et al._ , Prog. Part. Nucl. Phys. 125, 103948 (2022), arXiv:2111.05659 [hep-ph] .
  * Amelino-Camelia (2000) G. Amelino-Camelia, Lect. Notes Phys. 541, 1 (2000), arXiv:gr-qc/9910089 .
  * Kostelecky (2004) V. A. Kostelecky, Phys. Rev. D 69, 105009 (2004), arXiv:hep-th/0312310 .
  * Amelino-Camelia  _et al._ (1998) G. Amelino-Camelia, J. R. Ellis, N. E. Mavromatos, D. V. Nanopoulos, and S. Sarkar, Nature 393, 763 (1998), arXiv:astro-ph/9712103 .
  * Ellis  _et al._ (1999) J. R. Ellis, N. E. Mavromatos, and D. V. Nanopoulos, Gen. Rel. Grav. 31, 1257 (1999), arXiv:gr-qc/9905048 .
  * Jacobson (1991) T. Jacobson, Phys. Rev. D 44, 1731 (1991).
  * Amelino-Camelia  _et al._ (2003) G. Amelino-Camelia, F. D’Andrea, and G. Mandanici, JCAP 09, 006 (2003), arXiv:hep-th/0211022 .
  * Alfaro and Palma (2003) J. Alfaro and G. Palma, Phys. Rev. D 67, 083003 (2003), arXiv:hep-th/0208193 .
  * Amelino-Camelia (2005) G. Amelino-Camelia, Lect. Notes Phys. 669, 59 (2005), arXiv:gr-qc/0412136 .
  * Alfaro  _et al._ (2004) J. Alfaro, M. Reyes, H. A. Morales-Tecotl, and L. F. Urrutia, Phys. Rev. D 70, 084002 (2004), arXiv:gr-qc/0404113 .
  * Magueijo and Smolin (2005) J. Magueijo and L. Smolin, Phys. Rev. D 71, 026010 (2005), arXiv:hep-th/0401087 .
  * Amelino-Camelia  _et al._ (2005) G. Amelino-Camelia, C. Lammerzahl, A. Macias, and H. Muller, AIP Conf. Proc. 758, 30 (2005), arXiv:gr-qc/0501053 .
  * Nozari and Sefidgar (2008) K. Nozari and A. S. Sefidgar, Chaos Solitons Fractals 38, 339 (2008).


# A Novel Vision Transformer based Load Profile Analysis using Load Images as
Inputs

Hyeonjin Kim, Yi Hu, Kai Ye, Ning Lu  North Carolina State University
Raleigh, NC 27606, USA
{hkim66, yhu28, kye3<EMAIL_ADDRESS>
###### Abstract

This paper introduces ViT4LPA, an innovative Vision Transformer (ViT) based
approach for Load Profile Analysis (LPA). We transform time-series load
profiles into load images. This allows us to leverage the ViT architecture,
originally designed for image processing, as a pre-trained image encoder to
uncover latent patterns within load data. ViT is pre-trained using an
extensive load image dataset, comprising 1M load images derived from smart
meter data collected over a two-year period from 2,000 residential users. The
training methodology is self-supervised, masked image modeling, wherein masked
load images are restored to reveal hidden relationships among image patches.
The pre-trained ViT encoder is then applied to various downstream tasks,
including the identification of electric vehicle (EV) charging loads and
behind-the-meter solar photovoltaic (PV) systems and load disaggregation.
Simulation results illustrate ViT4LPA’s superior performance compared to
existing neural network models in downstream tasks. Additionally, we conduct
an in-depth analysis of the attention weights within the ViT4LPA model to gain
insights into its information flow mechanisms.

###### Index Terms:

Image processing, Load analysis, Pre-trained model, Smart meter data, Vision
transformer

##  I Introduction

Pre-trained neural network models have been widely used in Natural Language
Processing (NLP) [1] and Computer Vision (CV) [2] tasks. Compared to training
a model from scratch, using pre-trained models reduces reliance on labeled
data and saves considerable time and computational resources when performing
downstream tasks [3]. In the field of natural language processing (NLP), the
adoption of large pre-trained language models, such as BERT and GPT, has
brought about a transformative shift in the domain of language comprehension
and generation. These models undergo extensive training on extensive text
corpora, endowing them with the capacity to apprehend intricate linguistic
structures and contextual nuances. Consequently, NLP has undergone remarkable
progress in the past couple of years, with pre-trained models emerging as
pivotal cornerstones for a diverse array of downstream applications, including
but not limited to sentiment analysis, machine translation, text
summarization, and question answering.

The transformer model exhibits exceptional scalability, surpassing traditional
machine learning based models such as Recurrent Neural Networks (RNNs) and
convolutional neural networks (CNNs). This led to the creation of the Vision
Transformer (ViT) [4], specifically designed for processing computer vision
(CV) tasks. In contrast to conventional CNNs, ViT uses a transformer
architecture by dividing images into fixed-size patches, linearly projecting
them, and employing self-attention mechanisms to capture long-range
dependencies among patches. ViT soon excels in conducting a wide variety of CV
tasks, including image classification, object detection, and segmentation. The
success of ViT in the realm of CV underscores the transformative potential of
applying transformer architectures across different domains.

In the field of power system analysis, load profile analysis (LPA) is becoming
increasingly vital for performing tasks like load disaggregation, model
parameterization, customer segmentation, load flexibility analysis, and
identifying behind-the-meter resources. However, a distinctive challenge in
this field is the scarcity of publicly accessible, non-sensitive power system
datasets. The sensitivity and proprietary nature of power grid operations make
it not only financially burdensome but often impractical for researchers and
developers to amass the substantial volume of training data required for
training resilient machine learning models.

Figure 1: An illustration of the ViT4LPA architecture. (a) Profile-to-image
conversion, (b) ViT4LPA workflow, and (c) Pre-training process.

Hence, leveraging extensive data for pre-training a model applicable to
subsequent LPA tasks, akin to the use of BERT in NLP and ViT in CV, can
effectively address the challenge of restricted data accessibility. To date,
the development of a versatile pre-trained model suitable for various
downstream tasks in the field of LPA has not received adequate attention from
researchers in the power system domain. Recent research endeavors have
prominently revolved around transfer learning, a technique wherein neural
networks, initially trained for specific supervised learning tasks, are
repurposed to perform related tasks by leveraging the knowledge they have
accrued. In [5], cross-domain transfer learning was introduced, aiming to
apply latent features learned from one appliance to another, particularly for
non-intrusive load monitoring tasks. In [6], the authors explored LPA through
representation learning using convolutional autoencoders. In our recent work
[7], we introduced a BERT-based load profile inpainting approach, serving as a
foundational model to streamline the restoration of missing data tasks.
Nonetheless, it is currently primarily focused on addressing a specific LPA
task.

Thus, in this paper, we present a pre-trained model framework tailored for
LPA. Our main contribution is the introduction of ViT4LPA, a novel approach
based on Vision Transformers (ViT) specifically designed for LPA tasks. This
approach marks a significant shift in the LPA domain by not only reducing the
dependence on labeled datasets but also enabling its seamless application
across a diverse spectrum of downstream tasks.

##  II Methodology

This section introduces ViT4LPA, an innovative Vision Transformer (ViT) based
approach for LPA.

###  II-A Profile-Image Conversion

One of our key contributions is to base the LPA study on load images derived
from these time-series load profiles. Traditional LPA studies typically use
time-series load profiles as inputs to capture temporal cyclic
characteristics, including daily, weekly, or monthly variation patterns.
However, in a more recent development, a few authors have proposed the
conversion of load profiles into color-coded load images for autoencoder-based
clustering, as seen in [6], and for BERT-based missing data recovery, as
demonstrated in [7]. Motivated by the accomplishments in those endeavors, we
use load images, as opposed to load profiles, to pre-train a ViT-based
foundational model that can be fine-tuned with a small amount of data when
conducting downstream LPA tasks.

As illustrated in Fig. 1(a), we generate load image comprising three distinct
channels. Channels 1 to 3 encompass data from smart meter loads, temperature
readings, and irradiance profiles, respectively. This approach allows the
encoder to capture not only hidden patterns within the load data but also the
correlation with temperature and irradiance variations. $x-$axis of the image
corresponds to the number of data ($N_{T}$) within one day and $y-$axis
represents the number of days ($N_{D}$). To convert a data point into a color
patch, we initially perform a linear projection of the data point onto the
[0,1] range. For instance, consider the case where the minimum load
consumption ($p^{-}$) is $-$4kW and the maximum load consumption ($p^{+}$) is
24kW. In this scenario, the load consumption at hour $t$, ${p_{t}}$, can be
normalized as follows: $(p_{t}-p^{-})/(p^{+}-p^{-})$.

Load images excel at condensing vast amounts of multi-modal data, encompassing
variables like load, temperature, and solar irradiance, facilitating the
streamlined consolidation of information. This proves especially advantageous
when compressing monthly or yearly data into load images, aligning profiles
based on their cyclical characteristics. Through the transformation of load
profiles into load images, we can seamlessly apply the ViT model [4],
initially designed for computer vision tasks, to effectively handle LPA tasks.

###  II-B Vision Transformer Encoder Inputs

In Fig. 1(b), we present the workflow of ViT4LPA. To directly apply the ViT
model architecture, we initially partition a load image comprising
$N_{D}\times N_{T}$ pixels into $\frac{{N_{D}}\times{N_{T}}}{{N_{P}}^{2}}$
image patches, ensuring each patch contains $N_{P}\times N_{P}$ pixels.
Subsequently, these image patches are flattened into a sequence of patches,
ready to be fed as inputs to the ViT model.

Note that the division of load patches is contingent on both the data
resolution ($N_{T}$) and data duration ($N_{D}$). In this paper, due to the
page limit, we will focus on introducing the ViT4LPA architecture and
workflow. The detailed discussion regarding the selection of $N_{P}$ in
relation to $N_{T}$ and $N_{D}$ will be presented in our follow-up journal
paper. Thus, we fix ${N_{T}}=24$ hours, ${N_{D}}=24$ days, and $N_{P}=4$. In
this setup, a load image will contain 24 rows (representing the number of
days) and 24 columns (representing the number of data points in a day). After
partitioning into $6\times 6$ image patches, each patch contains $4\times 4$
color patches. We then flatten the 36 image patches and send them directly to
the ViT encoder (the blue block in Fig. 1(b)) for processing. The descriptions
of the ViT encoder model can be found in [4].

The load embeddings produced by the pre-trained ViT encoder can be effectively
used for various downstream LPA tasks. In this paper, we use the behind-the-
meter load identification task and load disaggregation task as examples to
showcase the efficacy of ViT4LPA. ViT4LPA will be employed to identify the
following three types of distributed energy resources: electric vehicle (EV)
charging load, Photovoltaic (PV) generation, and Heating, Ventilation, and
Air-Conditioning (HVAC) load.

###  II-C Pre-training ViT4LPA using Masked Image Modeling Tasks

In Fig.1(c), we describe the architecture of the masked autoencoder network
used to pre-train the ViT4LPA encoder, following the masked autoencoder scheme
introduced in[2]. The training task involves restoring the masked image
patches using a reconstruction decoder. The input to this decoder is the load
image embedding produced by the ViT encoder using unmasked image patches as
inputs. In Fig. 1(c), since the ViT-generated embeddings are for the unmasked
image patches, it is necessary to incorporate the mask embeddings to ensure
that the input of the decoder matches the dimensionality of the original load
image.

As depicted in Fig. 1(c), 18 image patches are subjected to masking, leading
to the ViT encoder receiving only the visible 18 image patches as input.
Consequently, the output of the ViT encoder consists of only 18 embeddings,
requiring the insertion of 18 mask embeddings to achieve the necessary 36
embeddings for the reconstruction decoder. The 36 restored image patches by
the decoder will be rearranged back to a load image and compared with the
original load image to calculate losses.

Please note that various masking strategies, including grid-masking, random-
masking, span masking, and more, exist. However, due to page limitations, this
paper exclusively introduces the grid-masking approach. This choice aligns
with the primary focus of the downstream LPA tasks discussed in this paper,
i.e., the tasks of behind-the-meter load identification and load
disaggregation. In these tasks, grid masking has been found to demonstrate
superior training efficiency when compared to random masking. Consequently,
for this paper’s scope, we will exclusively showcase results obtained through
grid masking. A more in-depth examination of the impacts of different masking
strategies (e.g., random masking or variable grid masking) on downstream task
performances will be extensively addressed in our follow-up journal paper.

##  III Simulation Results

In this section, we showcase the efficacy of the proposed pre-trained ViT4LPA
model in three downstream load identification tasks, using smart meter data as
inputs.

###  III-A Simulation Setup

To construct the training dataset for the ViT encoder, we generate 4,000 sets
of 1-hour resolution yearly load profiles. The original datasets are collected
by the Pecan Street project [8]. The load profiles are metered from 150
households in Austin, Texas, including sub-metered PV, EV, and HVAC load
consumption. The dataset spans 2 years, with a 1-minute data resolution. To
streamline the load identification task, we downsample the 1-minute data to 1
hour. The data augmentation strategy introduced in [9] is used to generate
additional training examples for bolstering the model’s resilience. The
dataset is partitioned as follows: the training set comprises load profiles
from 100 users, while the testing dataset includes the remaining 50 customers.
This ensures a comprehensive evaluation of the proposed models.

###  III-B Performance on the Masked Image Modeling Tasks

The hyperparameters used for pre-training the ViT4LPA encoder model (see the
network architecture depicted in Fig. 1(c)) are listed in Table I. As shown in
Fig. 2(a), the reconstruction results demonstrate that the model can
reconstruct the load image from partial information (i.e., unmasked image
patches) with satisfactory performance. Interestingly, the masked autoencoder
demonstrates increased effectiveness in recovering missing patches during the
morning to noon period (8:00-12:00), which aligns with the typical peak
consumption hours for the specific customer. In Fig. 2(b), $n$MAE distribution
of the masked image reconstruction is represented. Error is concentrated in
the 1-2% range, where the mean and standard deviation of the error are 1.40%
and 0.515%, respectively.

TABLE I: Parameters Used in the ViT4LPA Pre-training Process Network | Layers | Heads | Proj. Dim. ($D$) | Parameters
---|---|---|---|---
Encoder | 3 | 4 | 128 | 1.0M
Decoder | 2 | 2 | 32 | 2.0M
Hyperparam. | Batch size | Dropout | Optimizer | Epochs
64 | 0.1 | Adam | 50

(a)

(b)

Figure 2:  ViT4LPA encoder performance evaluation. (a) An illustration of the
original, masked, and reconstructed load image, and (b) Reconstruction error
distribution.

Subsequently, the cosine similarity matrix is employed to assess the
similarity of position embeddings among patches. In Fig.3(a), we depict
heatmaps showcasing 36 similarity matrices. Each matrix corresponds to the
cosine similarity of position embeddings between an image patch at a specific
position in the load image, in relation to the 36 image patches, which include
itself. To illustrate, in Fig.3(b), the top-left similarity matrix represents
the cosine similarity of position embeddings between the initial image patch
(consisting of 48 data points from day 1-4 and hour 1-4) and the remaining 36
image patches.

Notably, when considering photo image patches, we observe significant
similarities between a patch and its neighboring patches, regardless of their
alignment in columns, diagonally, or rows, as discussed in [4]. In contrast,
for load image patches, high similarity is primarily found among patches
aligned in columns, indicating strong correlations among profiles at the same
time of day across different days. For instance, data from hours 1-4 exhibits
a strong correlation with data from hours 1-4 on all other days (i.e., in the
column-wise direction). However, this correlation may be less prominent when
comparing data from adjacent hours, such as between hours 1-4 and hours 5-8
(i.e., along the row-wise direction).

Figure 3: Analysis on position embeddings of the ViT4LPA encoder. (a) Heatmaps
of the 36 similarity matrices and (b) Similarity matrix of image patch 1.

In Figs. 4(a) to (c), the mean self-attention matrices, averaged across
multiple attention heads within each layer, are depicted. Similar to Fig.
3(a), the 1-by-36 attention weights are reshaped into a 6-by-6 matrix to
incorporate the positional information of patches. Analyzing the mean
attention matrices provides valuable insights into the local and global
information flow within each layer and the receptive field of the model. An
important insight obtained from the study is that in the first layer, the
attention weights are sharply focused on a small set of patches, indicating a
localized focus on specific details. As the network deepens, particularly in
layers 2 and 3, attention is dispersed across a larger number of patches
compared to the initial layer. This shift implies the model’s transition from
capturing localized features to embracing a more holistic perspective,
encompassing the entire load image to extract comprehensive global
information.

(a)  (b)  (c)

Figure 4: Heatmaps of the mean of the self-attention layers. (a) Layer 1, (b)
Layer 2, and (c) Layer 3.

###  III-C Downstream Task 1: Load Identification

To demonstrate the effectiveness of using pre-trained model on downstream
tasks, we first compare the performance of ViT4LPA on two load identification
tasks: behind-the-meter PV and EV identification. The identification models
have the same network structure for both the scenarios with and without
ViT4LPA. Note that the ViT4LPA encoder model can be further fine-tuned using a
small amount of labeled training data.

As illustrated in Fig. 5, a significant performance boost is evident,
particularly when the training dataset is limited, such as within the range of
10k to 150k examples, by incorporating the pre-trained ViT4LPA encoder. Even
with the utilization of the complete dataset comprising 500k examples, a
performance gain of 1-2% is attainable when employing ViT4LPA. This outcome
highlights that harnessing a pre-trained model not only effectively reduces
the reliance on extensive training datasets but also consistently yields
robust performance enhancements through comprehensive pre-training.

(a)

(b)

Figure 5: Comparison between identification models with or without pre-trained
encoder network. (a) PV generation identification, and (b) EV load
identification.

In Table II, we compare the proposed ViT-based identification model with
conventional CNN models. The CNN architectures used in our experiments are
based on the Inception model [10], with variations incorporating both 2D and
1D convolution layers, respectively. The comparison reveals a significant
performance gap between the proposed ViT-based model and the benchmark CNN
models. While the PV identification accuracy appears similar across the
models, a notable difference emerges in EV identification accuracy. For the PV
identification task, the abundance of available data makes pre-training less
impactful. However, when it comes to EV identification, labeled data is
scarce, underscoring the increased importance of leveraging the knowledge
acquired through pre-training for effective identification. Consequently, in
the EV identification task, using ViT4LPA results in a substantial performance
enhancement compared to using the two CNN models.

TABLE II: Performance Comparison for Conducting PV and EV Load Identification Tasks Model | Accuracy (PV) | Accuracy (EV) | Overall
---|---|---|---
ViT (proposed) | 0.984 | 0.943 | 0.964
Inception (2D) | 0.988 | 0.912 | 0.950
Inception (1D) | 0.987 | 0.911 | 0.949

###  III-D Downstream Task 2: HVAC Load Disaggregation

Next, we showcase the performance improvement when using the pre-trained
ViT4LPA encoder on load disaggregation tasks, using HVAC load disaggregation
as an example. In Table III, we compare ViT4LPA with two benchmark models: the
Inception model and the bidirectional Long Short-Term Memory (BiLSTM) model.
Note that BiLSTM is a RNN model widely employed for sequential data analysis.
To evaluate the models, we use the normalized mean absolute error ($n$MAE) to
assess point-to-point error and energy error (EE) to measure cumulative
estimation deviation. For more details on the performance metrics used for
load disaggregation, please refer to our prior work [9].

TABLE III: Performance Comparison (Load Disaggregation Task) Model |  $n$MAE (%) | EE (kWh) | std($n$MAE)
---|---|---|---
ViT (proposed) | 6.89 | 2.80 | 2.21
Inception | 8.89 | 4.11 | 2.59
BiLSTM | 8.01 | 3.53 | 2.26
Figure 6: Illustration of HVAC load disaggregation results over a three-day
period for a customer.

(a)

(b)

Figure 7: Empirical distribution of $n$MAE disaggregation methods. (a) ViT4LPA
model, and (b) BiLSTM model.

From the results, we made the following observations

  * •

ViT4LPA exhibits the lowest point-to-point error, energy error, and error
standard deviation.

  * •

As shown in Fig. 6, the HVAC load curve generated by ViT4LPA closely aligns
with the ground truth, while the estimates derived from the BiLSTM model
display noticeable errors, especially during peak HVAC load periods.

  * •

As shown in Fig. 7, ViT4LPA exhibits remarkable consistency, as evidenced by
the narrow standard deviation of $n$MAE. It’s important to note that each
error observation is derived from averaging the resultant errors across
multiple days for each customer, demonstrating the model’s robustness across
various customers. In the case of our proposed method, errors are densely
concentrated within the 4-8% range, indicating a high level of precision.
Conversely, the errors of the BiLSTM models exhibit a comparatively broader
distribution with a median value of 7.6%.

##  IV Conclusion

In this study, we introduce ViT4LPA, a pre-trained Vision Transformer
(ViT)-based load image encoder designed to generate load embeddings from load
images. By converting load, temperature, and solar irradiance profiles into
load images, we enable the direct adoption of ViT, a powerful image processing
model, for Load Profile Analysis (LPA). Pre-trained on masked image
restoration tasks, ViT4LPA captures correlations among load profile pixels to
enhance performance in downstream tasks. To gauge the performance of ViT4LPA,
we conducted tests on two popular downstream LPA tasks: load identification
and load disaggregation. Our results demonstrate significant improvements in
both tasks, with lower point-to-point errors, cumulative energy errors, and
reduced standard deviations of errors. This underscores the potential of pre-
trained models, particularly in machine learning and data analysis for power
system applications. The capability of ViT4LPA to fine-tune with limited
labeled data is especially valuable in applications constrained by data
scarcity or concerns related to data privacy.

Moving forward, our research will delve into analyzing suitable masking
strategies and patch sizes to further enhance the performance across multiple
downstream tasks.

## References

  * [1] J. Devlin, M.-W. Chang, K. Lee, and K. Toutanova, “Bert: Pre-training of deep bidirectional transformers for language understanding,” _arXiv preprint arXiv:1810.04805_ , 2018.
  * [2] K. He, X. Chen, S. Xie, Y. Li, P. Dollár, and R. Girshick, “Masked autoencoders are scalable vision learners,” in _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , 2022, pp. 16 000–16 009.
  * [3] A. Radford, K. Narasimhan, T. Salimans, I. Sutskever _et al._ , “Improving language understanding by generative pre-training,” 2018.
  * [4] A. Dosovitskiy, L. Beyer, A. Kolesnikov, D. Weissenborn, X. Zhai, T. Unterthiner, M. Dehghani, M. Minderer, G. Heigold, S. Gelly _et al._ , “An image is worth 16x16 words: Transformers for image recognition at scale,” _arXiv preprint arXiv:2010.11929_ , 2020.
  * [5] M. D’Incecco, S. Squartini, and M. Zhong, “Transfer learning for non-intrusive load monitoring,” _IEEE Transactions on Smart Grid_ , vol. 11, no. 2, pp. 1419–1429, 2019.
  * [6] S. Ryu, H. Choi, H. Lee, and H. Kim, “Convolutional autoencoder based feature extraction and clustering for customer load analysis,” _IEEE Transactions on Power Systems_ , vol. 35, no. 2, pp. 1048–1060, 2019.
  * [7] Y. Hu, K. Ye, H. Kim, and N. Lu, “Bert-pin: A bert-based framework for recovering missing data segments in time-series load profiles,” _arXiv preprint arXiv:2310.17742_ , 2023.
  * [8] “https://www.pecanstreet.org.”
  * [9] K. Ye, H. Kim, Y. Hu, N. Lu, D. Wu, and P. Rehm, “A modified sequence-to-point hvac load disaggregation algorithm,” in _2023 IEEE Power & Energy Society General Meeting (PESGM)_. IEEE, 2023, pp. 1–5.
  * [10] C. Szegedy, W. Liu, Y. Jia, P. Sermanet, S. Reed, D. Anguelov, D. Erhan, V. Vanhoucke, and A. Rabinovich, “Going deeper with convolutions,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2015, pp. 1–9.


# Wireless for Control: Over-the-Air Controller

Pangun Park, Piergiuseppe Di Marco, and Carlo Fischione P. Park is with the
Department of Radio and Information Communications Engineering, Chungnam
National University, Korea (e-mail: pgpark@cnu.ac.kr). P. Di Marco is with the
Department of Information Engineering, Computer Science and Mathematics,
University of L’Aquila, Italy (e-mail: piergiuseppe.dimarco@univaq.it). C.
Fischione is with the Department of Network and Systems Engineering, School of
Engineering, KTH Royal Institute of Technology, Sweden (e-mail:
carlofi@kth.se).

###### Abstract

In closed-loop wireless control systems, the state-of-the-art approach
prescribes that a controller receives by wireless communications the
individual sensor measurements, and then sends the computed control signal to
the actuators. We propose an over-the-air controller scheme where all sensors
attached to the plant simultaneously transmit scaled sensing signals directly
to the actuator; then the feedback control signal is computed partially over
the air and partially by a scaling operation at the actuator. Such over-the-
air controller essentially adopts the over-the-air computation concept to
compute the control signal for closed-loop wireless control systems. In
contrast to the state-of-the-art sensor-to-controller and controller-to-
actuator communication approach, the over-the-air controller exploits the
superposition properties of multiple-access wireless channels to complete the
communication and computation of a large number of sensing signals in a single
communication resource unit. Therefore, the proposed scheme can obtain
significant benefits in terms of low actuation delay and low wireless resource
utilization by a simple network architecture that does not require a dedicated
controller. Numerical results show that our proposed over-the-air controller
achieves a huge widening of the stability region in terms of sampling time and
delay, and a significant reduction of the computation error of the control
signal.

###### Index Terms:

Over-the-air computation, Wireless communication, Networked control systems.

##  I Introduction

Time-critical operations over a wireless network are the heart of essential
infrastructures for monitoring and control of systems for factory automation,
process control, and power distribution [1, 2]. In Wireless Networked Control
Systems (WNCS), distributed sensors, controllers, and actuators exchange
sensing and actuating signals over a wireless network to achieve a control
objective [1]. The WNCSs are fundamentally different from traditional
distributed systems since the network dynamics, such as time-varying capacity,
node faults, and stochastic delay and reliability, significantly affect the
physical dynamics of the control system.

In the state-of-the-art WNCS architecture, a number of sensor nodes sample the
plant states and then each sensor sends the samples to the controller; when
all sensor data arrive, the controller calculates the control signal and
transmit it to actuators; actuators operate the received control signal to
manipulate the plant [2]. The delay and losses of sensors-to-controller and
controller-to-actuators links are crucial for the stability of the system and
may degrade the control performance. Furthermore, as advanced systems using
microsensors and embedded computers determine an increased density and scale
of sensor networks, more bandwidth is required to collect the sensor
measurements. Since the design of WNCS requires the simultaneous interaction
of communication, computation, and control aspects, in addition to physical
phenomena, new computation and communication architectures and protocols are
needed for closed-loop wireless control systems [1].

Within the communication community, the over-the-air computation paradigm has
been recently proposed and investigated to efficiently compute linear
functions of sensor measurements by utilizing the superposition property of
multiple-access channels [3, 4]. Opening an innovative field of applications
for this paradigm, we recognize that a linear feedback controller, which is
one of the most practical controller design approaches, computes its control
signal in the form of a weighted sum of plant state measurements [5].

In this paper, we propose the concept of an over-the-air controller (AirCont)
by adopting the over-the-air computation concept to compute the control signal
of the closed-loop wireless control systems. Since AirCont uses the multiple-
access channel to compute the weighted sum of the sensor data, it
fundamentally overcomes the sensor-to-controller and controller-to-actuator
state-of-the-art architecture, by introducing a fundamentally new direct
sensor-to-actuator architecture. We first analyze the achievable benefits of
the proposed AirCont scheme in terms of control stability with the operational
feasibility constraints. We also investigate the computation error of existing
power control approaches to compute the control signal subject to the peak
power constraints of nodes. We demonstrate that AirCont yields significant
improvements in terms of the achievable stability region and reduces the
computation error for different power limits, channel noise, and a number of
nodes compared to the state-of-the-art. Finally, we illustrate the use of
AirCont in a case study for the control of a ball and beam system.

The rest of the paper is organized as follows. In Section II, we describe our
control system model and the controller design schemes over a wireless network
model. Following that, we analyze the control stability region with the
feasible operating constraints and the computation error using different
schemes in Section III. Then, we evaluate the performance benefits using
AirCont compared to the state-of-the-art control schemes in Section IV.
Finally, we present conclusions and discussion.

_Notations:_ Normal font $x$, boldface lowercase font $\mathbf{x}$, and
boldface uppercase font $\mathbf{X}$ denote scalar, vector, matrix,
respectively. $x_{i}$ (resp, $x_{i,j}$) shows element $i$ (resp. $(i,j)$) of
vector $\mathbf{x}$ (resp. matrix $\mathbf{X}$).

##  II System Model

This section describes the proposed AirCont system model and the differences
with respect to the state-of-the-art scheme for the feedback control system
over a wireless network.

###  II-A Control System Model

Figure 1: Over-the-air controller for closed-loop wireless control systems.

We consider the problem of stabilizing a closed-loop control system composed
of $N$ wireless sensors communicating with a single actuator where each sensor
measures the state of a plant as depicted in Fig. 1. The number of plant
states is equivalent to the number of sensors, $N$.

We assume a standard continuous-time linear plant model

$\displaystyle\dot{\mathbf{x}}(t)$
$\displaystyle=\mathbf{A}\mathbf{x}(t)+\mathbf{b}u(t)$ (1)

with plant state $\mathbf{x}\in\mathds{R}^{N}$, state matrix
$\mathbf{A}\in\mathds{R}^{N\times N}$, input matrix
$\mathbf{b}\in\mathds{R}^{N}$, and control signal $u\in\mathds{R}$ applied to
the plant. Each sensor attached to the plant periodically samples the plant’s
state and sends it by a sampling time of $\delta$.

The standard feedback controller receives the sampled plant state
$\mathbf{x}(t)$ after they are transmitted by the sensors and computes the
discrete-time control input $u(k\delta)=-\mathbf{k}^{\top}\mathbf{x}(k\delta)$
where control gain $\mathbf{k}\in\mathds{R}^{N}$ if there is no communication
delay. Due to the wireless communication delay $\tau$, which we assume to be
shorter than the sampling period $\tau\leq\delta$, two control signals
$u((k-1)\delta)$ and $u(k\delta)$ apply during $k$-th sampling period [2]. The
control feedback signal is

$\displaystyle u(t^{+})$
$\displaystyle=-\mathbf{k}^{\top}\mathbf{x}(t-\tau),t\in\left\\{k\delta+\tau,k=0,1,2,\ldots\right\\}$
(2)

where $u(t^{+})$ is a piecewise continuous control signal and changes only at
$k\delta+\tau$. We remark that the key interactive parameter between the
network and the plant is the control signal, a linear combination of the
sensors’ signals.

###  II-B Proposed Control Signal Computation

In AirCont, all sensors simultaneously transmit their data so that the
actuator receives the control signal directly over the wireless channel based
on the over-the-air computation concept. The scheme works as follows: each
sensor $i$ scales the measured signal $x_{i}\in\mathds{R}$, representing the
plant state, by the Tx scaling factor $\beta_{i}\in\mathds{R}$ and
concurrently transmits it to the actuator through a wireless channel with the
channel coefficient $h_{i}\in\mathds{R}$, as illustrated in Fig. 1. The
receiver adds additive white Gaussian noise (AWGN) $n\in N(0,\sigma^{2})$. The
actuator then scales the received signal with the Rx scaling factor $\alpha$
to compute the weighted sum of sensor data as

$\displaystyle u_{\rm
air}(t^{+})=-\alpha\left((\mathbf{h}\odot\bm{\beta})^{\top}\mathbf{x}(t-\tau_{\rm
air})+n\right)$ (3)

where $\tau_{\rm air}$ is the network delay, $n\in N(0,\sigma^{2})$ is the
receiver’s AWGN and $\odot$ denotes Hadamard product between the Tx scaling
factor and the channel coefficient for each sensor to the actuator. The Rx
scaling factor $\alpha$ applies to both the signal and the noise terms. We
consider a maximum transmit power constraint of each sensor,
$\beta_{i}^{2}\leq\overline{p}$, so, $\beta_{i}\in[0,\sqrt{\overline{p}}]$.
Furthermore, we make the natural assumption that the sensors know the channel
coefficients, and their transmissions are synchronized. The network delay of
AirCont is only one single time slot $\tau_{\rm air}=T_{s}$ where $T_{s}$
denotes the duration of the time slot, since it leverages simultaneous
coherent transmission from multiple sensors to the actuator.

As a benchmark, we consider a state-of-the-art scheme where each sensor
attached to the physical plant sends the scaled signal $\beta_{i}x_{i}$
separately to the controller. Through the time slot allocated by a static
scheduler, each sensor is only allowed to transmit once within a time frame
cycle $\delta$, which is equivalent to the sampling period. To compute the
control signal, the controller scales the received signal of sensor $i$ with
the Rx scaling factor $\alpha_{s,i}$. The wireless channel between sensor $i$
and controller has a channel coefficient $h_{i}$ with AWGN $n_{s,i}\in
N(0,\sigma_{s,i}^{2})$. The controller then transmits the computed control
signal to the actuator where we consider the channel coefficient $h_{a}$ with
AWGN $n_{a}\in N(0,\sigma_{a}^{2})$ for the controller-to-actuator link.
Eventually, the actuator scales the received control signal with $\alpha_{a}$
to compensate the channel $h_{a}$. Hence, the resulting control signal of the
state-of-the-art scheme is

$\displaystyle u_{\rm
sota}(t^{+})=-\alpha_{a}\left(h_{a}\bm{\alpha}_{s}^{\top}\mathbf{D}\mathbf{x}(t-\tau_{\rm
sota})+\mathbf{n}_{s}\right)+n_{a}$ (4)

where $\mathbf{D}$ is a diagonal matrix with $d_{i,i}=h_{i}\beta_{i}$ and
$\tau_{\rm sota}$ is the network delay for the state-of-the-art scheme.

There are two main sources of network delay for the state-of-the-art scheme,
namely the sensor-to-controller delay $\tau_{\rm sc}$ and the controller-to-
actuator delay $\tau_{\rm ca}$. The end-to-end delay between sampling instance
and actuating instance is the sum of the sensor-to-controller delay and
controller-to-actuator delay, $\tau_{\rm sota}=\tau_{\rm sc}+\tau_{\rm ca}$.
Since the state-of-the-art scheme requires a transmission from all sensors to
the controller then to the actuator, the minimum delay in actuating the
control signal in each sampling instance is $\tau_{\rm sota}=(N+1)T_{s}$. In
contrast, AirCont integrates the communication and computation of a large
number of sensor data in one time slot. It does not rely on the dedicated
controller since the actuator directly adapts the control signal as a linear
combination of the sensor data. Hence, it reduces 2-hop communication, namely,
sensor-to-controller and controller-to-actuator, to a single hop, namely,
sensor-to-actuator.

##  III Performance Analysis

In this section, we first derive the achievable stability region restricted by
the operational constraints of AirCont and the state-of-the-art scheme. We
then investigate the computation error of the control signal using different
scaling control policies.

###  III-A Stability Analysis

A higher sampling rate is generally desirable in the discrete-time system
since it approximates well the continuous-time system. However, a higher
sampling rate increases the network delay due to the possible congestion [1].
The WNCS design needs to find a sampling rate and delay that can both ensure
the control stability and be achievable by the wireless network.

Quantifying the stability boundary of the control system with respect to the
sampling period $\delta$ and delay $\tau$ is useful to understand the control
performance tradeoffs. We define two distinct notions, namely, maximum
stability region and achievable stability region.

######  Definition 1 (Maximum Stability Region).

It is the set of values of sampling period $\delta$ and delay $\tau$, which
guarantees the control stability of the closed-loop system.

######  Definition 2 (Achievable Stability Region).

It is the set of values of sampling period $\delta$ and delay $\tau$ that can
be supported by the wireless network, which guarantees the control stability
of the closed-loop system.

We extend the approach in [2] to analyze the stability of the AirCont-based
system and the benchmark state-of-the-art scheme. By considering the sampled
system with sampling period $\delta$, the expected linear difference equation
becomes

$\displaystyle\overline{\mathbf{x}}((k+1)\delta)=\Phi\overline{\mathbf{x}}(k\delta)+\Gamma_{0}(\tau)\overline{u}(k\delta)+\Gamma_{1}(\tau)\overline{u}((k-1)\delta)$

where $\overline{\mathbf{x}}$ and $\overline{u}$ are the expected value of
$\mathbf{x}$ and $u$ with respect to noise factors,
$\Phi=\mathrm{e}^{\mathbf{A}\delta},$

$\displaystyle\Gamma_{0}(\tau)=\int_{0}^{\delta-\tau}\mathrm{e}^{\mathbf{A}s}\,\mathbf{B}\,\mathrm{d}s\,,\,\,\,\,\Gamma_{1}(\tau)=\int_{\delta-\tau}^{\delta}\mathrm{e}^{\mathbf{A}s}\,\mathbf{B}\,\mathrm{d}s\,.$

$\Gamma_{0}(\tau)$ and $\Gamma_{1}(\tau)$ are the input matrix of two control
signals $\overline{u}(k\delta)$ and $\overline{u}((k-1)\delta)$ due to the
network delay $\tau$, respectively. Remind that the end-to-end network delay
from the sampling instance is shorter than the sampling period,
$\tau\leq\delta$.

By defining
$\mathbf{z}(k\delta)=[\overline{\mathbf{x}}^{\top}(k\delta)\;\overline{u}((k-1)\delta)]^{\top}$
as the augmented state vector, the augmented system becomes

$\displaystyle\mathbf{z}((k+1)\delta)=\tilde{\Phi}\mathbf{z}(k\delta)\,.$ (5)

The discrete-time linear system is asymptotically stable (in fact,
exponentially stable) if all eigenvalues of $\tilde{\Phi}$ have norm strictly
less than one, i.e., the spectral radius $\rho(\tilde{\Phi})<1$.

By using the AirCont scheme in Eq. (3), the augmented system matrix in (5)
becomes

$\displaystyle\tilde{\Phi}_{\rm air}=\begin{bmatrix}\Phi-\Gamma_{0}(\tau_{\rm
air})\alpha(\mathbf{h}\odot\bm{\beta})^{\top}&\Gamma_{1}(\tau_{\rm air})\\\
-\alpha(\mathbf{h}\odot\bm{\beta})^{\top}&0\end{bmatrix}\,.$ (6)

The matrix $\tilde{\Phi}_{\rm air}$ depend on both control aspects
(continuous-time plant dynamics $(\mathbf{A},\mathbf{b})$ and sampling period
$\delta$) and wireless communication aspects (channel coefficient
$\mathbf{h}$, delay $\tau_{\rm air}$ , and Tx-Rx scaling factors
$(\bm{\beta},\alpha)$).

In a similar way, the stability condition of the state-of-the-art scheme is
analyzed by considering Eq. (4), which gives the augmented system matrix

$\displaystyle\tilde{\Phi}_{\rm sota}=\begin{bmatrix}\Phi-\Gamma_{0}(\tau_{\rm
sota})\alpha_{a}h_{a}\bm{\alpha}_{s}^{\top}\mathbf{D}&\Gamma_{1}(\tau_{\rm
sota})\\\ -\alpha_{a}h_{a}\bm{\alpha}_{s}^{\top}\mathbf{D}&0\end{bmatrix}\,.$
(7)

The maximum stability region is determined by the spectral radius of the
extended closed-loop system matrices of Eqs. (6) and (7), whereas the
achievable stability region is a subset of the maximum stability region
restricted to the feasibility constraints $\tau_{\rm air}\leq\delta$ and
$\tau_{\rm sota}\leq\delta$.

###  III-B Mean Square Error Analysis

This section investigates the computation error (measured by the Mean Square
Error (MSE)) of different Tx-Rx scaling policies subject to individual Tx
scaling constraints of sensors. Some existing works [3, 6] have already
considered the Tx-Rx scaling optimization problem for the over-the-air
computation, where the objective is to minimize the computation error subject
to the Tx scaling limits of sensors. Inspired by [3, 6], we propose a
minimization of the average MSE with respect to the control signal
$\mathbf{k}^{\top}\mathbf{x}$. To simplify the analysis, we assume that
measured signals $x_{i}$ are independent and follow the standard normal
distribution $N(0,1)$.

We first analyze the computation distortion of the control signal using
AirCont. The MSE between Eq. (3) and the control signal
$\mathbf{k}^{\top}\mathbf{x}$ is

$\displaystyle\mathrm{MSE}_{\rm air}$
$\displaystyle=\mathds{E}\left[\left|\alpha\left((\mathbf{h}\odot\bm{\beta})^{\top}\mathbf{x}+n\right)-\mathbf{k}^{\top}\mathbf{x}\right|^{2}\right]$
$\displaystyle=\left(\alpha(\mathbf{h}\odot\bm{\beta})-\mathbf{k}\right)^{\top}\left(\alpha(\mathbf{h}\odot\bm{\beta})-\mathbf{k}\right)+\sigma^{2}\alpha^{2}$
(8)

where the expectation of MSE is calculated with respect to the distributions
of $\mathbf{x}$ and $n$. We obtain the optimal Tx-Rx scaling policy of the MSE
minimization without any constraints on the Tx scaling factor. The optimal Tx-
Rx scaling factors are $\beta_{i}\rightarrow k_{i}/(\alpha h_{i})$ and
$\alpha\rightarrow 0$ since these scaling factors minimize the first and
second terms of Eq. (8), respectively. By putting the optimal solutions of Tx-
Rx scaling factors to Eq. (8), the MSE using AirCont approaches $0$. However,
the optimization problem becomes non-convex with the Tx scaling limits. We
extend the existing Tx-Rx scaling policy [6] to the MSE minimization problem
of the weighted sum of sensor data and evaluate its performance.

As a benchmark, we investigate the optimal scaling factors of Tx scaling
factor $\bm{\beta}$ of sensors, Rx scaling factor $\bm{\alpha}_{s}$ of
controller, and Rx scaling factor $\alpha_{a}$ of actuator for the state-of-
the-art scheme. By considering Eq. (4), the computation distortion, MSE, of
the control signal $\mathbf{k}^{\top}\mathbf{x}$ is

$\displaystyle\mathrm{MSE}_{\rm sota}=$
$\displaystyle\,\mathds{E}\left[\left|\alpha_{a}\left(h_{a}\bm{\alpha}_{s}^{\top}\left(\mathbf{D}\mathbf{x}+\mathbf{n}_{s}\right)+n_{a}\right)-\mathbf{k}^{\top}\mathbf{x}\right|^{2}\right]$
$\displaystyle=$
$\displaystyle\,(\alpha_{a}h_{a}\mathbf{D}\bm{\alpha}_{s}-\mathbf{k})^{\top}(\alpha_{a}h_{a}\mathbf{D}\bm{\alpha}_{s}-\mathbf{k})$
$\displaystyle+\alpha_{a}^{2}h_{a}^{2}\sigma_{s}^{2}\bm{\alpha}_{s}^{\top}\bm{\alpha}_{s}+\alpha_{a}^{2}\sigma_{a}^{2}$
(9)

where the expectation of MSE is calculated with respect to $\mathbf{x}$,
$\mathbf{n}_{s}$, and $n_{a}$. Recall that the diagonal matrix $\mathbf{D}$
depends on the channel coefficient $\mathbf{h}$ and the Tx scaling factor
$\bm{\beta}$ as $d_{i,i}=h_{i}\beta_{i}$. The MSE of the state-of-the-art
scheme converges to $0$ when $\beta_{i}\rightarrow
k_{i}/(\alpha_{a}\alpha_{s,i}h_{i}),\alpha_{a}\rightarrow 0$ and
$\alpha_{s,i}\rightarrow 0$ if the Tx scaling factor is not constrained.

We separate the constrained optimization problem into two sub-problems,
namely, $\bm{\beta}^{\ast}$ and $\bm{\alpha}_{s}^{\ast}$ for sensor-to-
controller and $\alpha_{a}^{\ast}$ for controller-to-actuator. By considering
Eq. (9), $\mathbf{h},\mathbf{k}$ and $\sigma_{s}^{2}$, we reformulate the
optimization problem to optimize $\bm{\beta}$ and $\bm{\alpha}_{s}$ for
sensor-to-controller communications where the objective function is

$\displaystyle\mathrm{MSE}_{\rm sota}^{\rm sc}$
$\displaystyle=(\mathbf{D}\bm{\alpha}_{s}-\mathbf{k})^{\top}(\mathbf{D}\bm{\alpha}_{s}-\mathbf{k})+\sigma_{s}^{2}\bm{\alpha}_{s}^{\top}\bm{\alpha}_{s}$
(10)

subject to the Tx scaling limits. To minimize the first term of Eq. (10), we
optimize each pair of $\beta_{i}$ and $\alpha_{s,i}$ from sensor $i$ to the
controller to meet $\mathbf{D}\bm{\alpha}_{s}=\mathbf{k}$. Since the noise
term of Eq. (10) only depends on $\bm{\alpha}_{s}$, the optimal values are
$\beta_{i}^{\ast}=\sqrt{\overline{p}}$ and
$\alpha_{s,i}^{\ast}=(h_{i}k_{i}\beta_{i}^{\ast})/((h_{i}\beta_{i}^{\ast})^{2}+\sigma_{s}^{2})$.
Substituting $\beta_{i}^{\ast}$ and $\alpha_{s,i}^{\ast}$ to Eq. (9), we
obtain the optimal scaling factor of the actuator as
$\alpha_{a}^{\ast}=(h_{a}\bm{\alpha}_{s}^{{\ast}\top}\mathbf{D}^{\ast}\mathbf{k})/(h_{a}^{2}\bm{\alpha}_{s}^{{\ast}\top}(\mathbf{D}^{\ast})^{2}\bm{\alpha}_{s}+h_{a}^{2}\sigma_{s}^{2}\bm{\alpha}_{s}^{{\ast}\top}\bm{\alpha}_{s}^{\ast}+\sigma_{a}^{2})$
where the diagonal matrix $\mathbf{D}^{\ast}$ is calculated using
$\bm{\beta}^{\ast}$.

We note that the sensor-to-controller and controller-to-actuator scheme
significantly increases the complexity of the operation compared to AirCont
since the actuator requires the knowledge of the sensor-to-controller
communication including $\mathbf{h},\bm{\beta}^{\ast}$,
$\bm{\alpha}_{s}^{\ast}$ and $\sigma_{s}^{2}$ in addition to
$h_{a},\sigma_{a}^{2}$ and $\mathbf{k}$.

##  IV Performance Evaluation

This section evaluates the stability region and the average MSE of AirCont
compared to the state-of-the-art scheme.

(a) Stability region using AirCont.

(b) Stability region using state-of-the-art scheme.

Figure 2: Maximum stability region and achievable stability region using
AirCont and the state-of-the-art scheme with different sampling period
$\delta$ and delay $\tau$. All circular markers present the maximum stability
region for a given sampling period and delay. To quantify the network effects,
we represent the achievable stability region of WNCS using AirCont and the
state-of-the-art scheme with filled circular markers. The solid line shows the
lower bound of the feasible sampling period and delay $\tau\leq\delta$ using
each of the schemes.

###  IV-A Stability Performance

We first quantify the maximum stability region and achievable stability region
of the control system with respect to the sampling period $\delta$ and delay
$\tau$. We consider the position control for the ball and beam system where
the linear time-invariant model with $N=4$ sensors and $1$ actuator is used
[5]. Fig. 2 shows the maximum stability region and the achievable stability
region of both AirCont and state-of-the-art scheme with different sampling
period $\delta$ and delay $\tau$ for $\tau\leq\delta$. For a fair comparison,
the control gains for AirCont and state-of-the-art scheme are set to be equal,
with
$\alpha(\mathbf{h}\odot\bm{\beta})^{\top}=\alpha_{a}h_{a}\bm{\alpha}_{s}^{\top}\mathbf{D}=[6.67,11.09,41.15,11.27]$
where $\mathbf{D}$ is a function of $\mathbf{h}$ and $\bm{\beta}$ without any
noise.

Given sampling period $\delta$ and delay $\tau$, we plot the stability region
by evaluating the spectral radius of extended closed-loop system matrices of
Eqs. (6) and (7). In Fig. 2, a point is marked with circular markers to
represent the maximum stability region. As expected, the maximum stability
regions of AirCont and state-of-the-art scheme are identical since the network
effects are not explicitly considered. To quantify the network effects, the
achievable stability regions of both AirCont and state-of-the-art scheme are
reported with filled circular markers. The solid line shows the lower bound of
the feasible sampling period and delay $\tau\leq\delta$ using AirCont and the
state-of-the-art scheme in Fig. 2.

In Fig. 2, a lower sampling period and lower delay are generally desirable for
the maximum stability region since the control system ensures the stability
for a delay up to the full sampling period for $\delta\leq 0.07$s. As the
sampling period $\delta$ increases, the upper bound of $\tau/\delta$ of the
maximum stability region is considerably reduced. Note that the control system
becomes unstable even without delay for $\delta>0.26$s. However, it is not
trivial to quantify the maximum stability boundaries as we observe that the
lower delay $\tau/\delta\leq 0.11$ is worse for the control stability for
$0.18\mathrm{s}\leq\delta\leq 0.26$s, as shown in Fig. 2.

The network performance heavily affects the achievable stability region of the
control system. The lower delay is not achievable for the short sampling
period due to the fundamental congestion of the network performance. In Fig.
2(b), the state-of-the-art scheme is not able to provide the control stability
of the plant for $\delta<0.05$s due to the minimum network delay constraint
$\tau_{\rm sota}=0.05\mathrm{s}\leq\delta$. On the other hand, the minimum
sampling period of AirCont is $\delta=0.01$s, significantly lower than the one
of the state-of-the-art scheme in Fig. 2(a).

A large achievable stability region improves operating robustness against
uncertain losses and delays, and energy efficiency by reducing the sampling
rate. By comparing Figs. 2(a) and 2(b), the achievable stability region using
AirCont is $7.3$ times larger than that using the state-of-the-art scheme. The
main reason is that AirCont only takes a single time slot to compute and
communicate the control signal for the sensor-actuator, independently of the
number of sensors.

(a) Average control MSE versus different Tx scaling limits.

(b) Average control MSE versus different noise variance.

Figure 3: Average control MSE of both AirCont and state-of-the-art scheme with
$N=10,100$ for various Tx scaling limits and noise variance.

###  IV-B MSE Performance

Now, we investigate how AirCont and the state-of-the-art scheme behave with
different Tx scaling limits, noisy channels, and a number of sensors on the
average MSE. We set the peak transmit power $\overline{p}=2.5$ and the noise
variance $\sigma^{2}=\sigma_{s}^{2}=\sigma_{a}^{2}=0.5$ unless otherwise
stated. The channel coefficient is i.i.d Rayleigh fading with unit variance.
We set the control gain $\mathbf{k}$ with a uniformly generated random value
between $0$ and $100$. To normalize the effect of the control gain, we define
the average control MSE as the expected value of MSE dividing with the sum of
the squared control gain, namely,
$\mathds{E}\left[\mathrm{MSE}\right]/(\mathbf{k}^{\top}\mathbf{k})$. We
calculate the average control MSE using extensive Monte Carlo simulations with
random channel and control gain realizations.

Fig. 3(a) depicts the average control MSE of both AirCont and state-of-the-art
scheme with number of sensors $N=10,100$ while varying the peak power limits
$\overline{p}=0.1,\ldots,5$. The AirCont scheme outperforms consistently well
the state-of-the-art scheme with respect to varying Tx scaling constraints.
While the average control MSEs of both schemes decrease approximately
exponentially with the peak power limits, the gap between them increases as
the power constraint becomes strict. The AirCont scheme has control MSE
approaching $0$ as the peak power limits relaxes and the benefits of AirCont
drastically improve as the number of sensors increases, while this effect is
negligible for the state-of-the-art scheme.

Fig. 3(b) plots the average control MSEs of both AirCont and state-of-the-art
scheme with number of nodes $N=10,100$ while varying noise variances
$\sigma^{2}=\sigma_{s}^{2}=\sigma_{a}^{2}=0.1,\ldots,1$. The average control
MSEs of both schemes roughly increases linearly with the noise variance. The
AirCont scheme provides a remarkably lower average control MSE than that of
the state-of-the-art scheme throughout the whole considered range.

Comparing the MSE slopes between AirCont and state-of-the-art scheme, the
average control MSE of AirCont is less sensitive to the noise variance than
the state-of-the-art scheme. The state-of-the-art scheme possibly emphasizes
the noise effect of the multi-hop communication for sensor-to-controller-to-
actuator links due to the channel distortion and the Tx scaling limit, as
discussed in Section III-B. The MSE gain of AirCont increases as the noise
variance increases. Furthermore, while AirCont considerably improves the
average control MSE for a large number of nodes $N=100$, this effect is
negligible for the state-of-the-art scheme due to the presence of uncorrelated
noise for each sensor transmission. By analyzing Fig. 3, the AirCont scheme is
particularly attractive for WNCS relying on dense sensor networks with strict
Tx scaling limits and noisy channels.

Figure 4: Plant outputs of the ball and beam control system using an ideal
solution, state-of-the-art scheme, and AirCont

Finally, we illustrate how AirCont improves the control performance as a case
study using the ball and beam system [5]. In Fig. 4, we show $4$-th plant
output of the control system using an ideal solution, state-of-the-art scheme,
and AirCont. Note that the ideal solution means no delay and no channel noise
of the communication. The plant output of AirCont is very similar to that of
the ideal solution. The oscillations of the state-of-the-art scheme imply that
it does not essentially guarantee good control performance. The network
operating region of the state-of-the-art scheme is closer to the unstable
region, as shown in Fig. 2(b) due to the minimum sampling period
$\delta=0.05$s. Furthermore, the control stability using the state-of-the-art
scheme is significantly vulnerable to the noise since it increases the
computation error of the control signal.

##  V Conclusion

This paper presents AirCont, a novel paradigm that eliminates the control unit
and computes the control signal for closed-loop wireless control systems,
adopting the over-the-air computation concept. As opposed to the sensor-to-
controller and controller-to-actuator communication of state-of-the-art
schemes, AirCont effectively integrates communication and computation by
harnessing interference for computing the control signal as the weighted sum
of the sensor data. This approach simplifies the control system operation,
design, and analysis using direct sensor-to-actuator communications without
relying on a dedicated control unit. We demonstrated that AirCont can
dramatically improve the achievable stability region compared to the state-of-
the-art scheme. Moreover, numerical results confirmed that AirCont
substantially reduces the computation error of the control signal for various
power limits, channel conditions, and number of nodes.

Previous works on over-the-air function computation focus on optimizing the
Tx-Rx scaling factors to minimize the computation error. However, this
solution does not guarantee the optimal control cost. Inspired by this
observation, we aim at developing a Tx-Rx scaling policy to provide robust
control performance over uncertain channel distortion.

## References

  * [1] P. Park, S. C. Ergen, C. Fischione, C. Lu, and K. H. Johansson, “Wireless network design for control systems: A survey,” _IEEE Communications Surveys Tutorials_ , vol. 20, no. 2, pp. 978–1013, 2018.
  * [2] J. Lunze, _Control Theory of Digitally Networked Dynamic Systems_. pringer, Heidelberg, 2014,.
  * [3] O. Abari, H. Rahul, and D. Katabi, “Over-the-air function computation in sensor networks,” 2016, arXiv:1612.02307.
  * [4] H. Hellström, M. Barros da Silva Jr, M. Chen, M. Amiri, V. Fodor, V. Poor, and C. Fischione, “Wireless for machine learning,” 2020, arXiv:2008.13492.
  * [5] G. C. Goodwin, S. F. Graebe, and M. E. Salgado, _Control System Design_. Prentice Hall, 2000.
  * [6] W. Liu, X. Zang, Y. Li, and B. Vucetic, “Over-the-air computation systems: Optimization, analysis and scaling laws,” _IEEE Transactions on Wireless Communications_ , vol. 19, no. 8, pp. 5488–5502, 2020.