# Power-efficient Sampling Time and Resource Allocation in Cyber Physical
Systems with Industrial Application

Atefeh. Termehchi, and Mehdi. Rasti Atefeh. Termehchi and Mehdi. Rasti are
with the Department of Computer Engineering and Information Technology,
Amirkabir University of Technology, Tehran, Iran (e-mail: atefetermehchy,
rasti@aut.ac.ir).

###### Abstract

Cyber Physical Systems (CPSs) are the result of convergence of computation,
networking, and control of physical process. In this paper, we consider an
industrial CPS consisting of several control plants and Rate Constrained (RC)
users that communicate via a single cell OFDMA network. The problem of jointly
determining the sampling instant of each control plant and allocating power
and sub-carrier, in the CPS, is formulated. The problem is a multi objective
optimization with aims of determining the next maximum allowable sampling
instant of each control plant and minimizing the power consumption in uplink
and downlink, considering the dynamics and desired performance of control
plants, the quality of service requirement of RC users, power and sub-carrier
constraints. To solve the multi objective optimization problem, a novel
approach is proposed which decomposes the optimization problem into two
smaller loosely coupled problems. we show the effectiveness of the proposed
approach through simulation results.

###### Index Terms:

cyber physical system, power consumption minimization, sampling time, resource
allocation, cellular network.

##  I Introduction

Cyber Physical System is the result of integration of three subsystems
including computation, networking and control of physical process. CPSs are an
integral part of many applications such as healthcare industry, smart grid,
industrial automation, and avionic. With the introduction of CPS and Internet
of Things (IoT) concepts in the industrial automation, the fourth industrial
revolution, namely, Industry 4.0, is emerging [1]. The new technologies are
expected to make a tremendous change in future industrial systems [2, 3]. In
the direction of designing Industry 4.0, there are some key challenges that
should be taken into account, such as stability of control subsystem,
deterministic bounded delay, reliable communication, energy efficiency and
limited resources of wireless network [4, 5, 6]. In control subsystems, the
sampled data of sensors is sent to the controller(s), and then the control
messages are sent to the actuators from controller(s), through the network.
The stability and the performance of the control subsystem can be affected by
sampling rate and network based uncertainties like delay, jitter, packet loss,
and resource contention [7]. Therefore, significant researches have focused on
determining the sampling rate in control subsystem, guaranteeing deterministic
bounded delay and reliability in wireless network over the past decades [8, 9,
10, 11, 12].
The authors of [12] propose near optimal sampling rates of control plants
sharing a limited bandwidth in a WirelesHART network. The authors address the
problem of overall control cost minimization while all data flows meet their
end-to-end deadlines. In [8], a new routing protocol is proposed to minimize
the end-to-end delay for a cognitive radio network that supports real time
applications. On the other hand, due to energy crisis, carbon emission concern
and also given the fact that most wireless devices relay on batteries,
designing new energy-saving network is receiving substantial attention from
researchers [6, 13, 14]. However, there is a trade-off between the energy
saving and supplying the deterministic bounded delay [15]. To guarantee a
deterministic delay, a fixed delay bound have to be provided within a
connection even in a worst case. A method is proposed in [13] to schedule
sleep states of wireless computing systems to minimize the power consumption
such that a given tolerable delay constraint is not violated. It is assumed in
[13] that the sampling instants of control subsystem and its tolerable delay
are predefined. ̵َAlso, there is another trade-off between the energy saving
and the sampling rate. In the case of increasing the sampling rate, the number
of sensors sampled data transmissions and control messages will increase and
as a result, energy consumption will increase.

Generally, there are three methods to determine the sampling time in a control
subsystem: time-triggered (e.g periodic), event-triggered and self-triggered.
In event-triggered, instead of periodically updating the sampled data of
sensors and control input, a new sample is picked when a threshold, defined by
the state of the plant, is violated. On the other hand, under the self-
triggered implementations, the next sampling time is computed based on the
last state measurement and therefore a constant plant monitoring is not
required as in event-triggered method. The self-triggered is more appealing
scheme due to its irreplaceable advantage in reducing the number of
transmissions between sensors and controller, between controller and actuators
and, accordingly, computation load in the controller [16, 17]. A new dynamic
transmission scheduling is proposed in [18] to ensure the stability of several
processes controlled over a shared IEEE 802.15.4 wireless network and also to
reduce the number of sensors sampled data transmissions and increase the sleep
time of transmission nodes. A self-triggered sampler is proposed in [18] to
guarantee the stability of each control plant and reduce energy consumption.
In addition, in [18], it is shown that applying self-triggered sampling may
not provide any benefit if the MAC (Media Access Control) layer parameters are
not appropriately adjusted. In fact, in contrast to traditional MAC layer
designs, the transmission scheduling in CPSs should jointly consider sampling
instants, the stability of control subsystem, relevant deterministic delay
requirement and network specification [18, 19].

Several standards and relevant MAC protocol have been suggested for wireless
communication in Industrial CPSs such as, IEEE 802.15.4, IEEE 802.11ac and
WirelessHART so far. However, the fifth generation (5G) may provide an ideal
platform for these applications [20, 21]. It is predicted that 5G technologies
will be able to support ultra reliable low latency communications (URLLC) and
massive machine-type communications (mMTC), which are aligned with industrial
wireless network requirements. The 5G technologies may bring considerable
benefits to industrial network such as reliability, global connectivity and
security. However, a set of challenges need to be addressed before 5G
deployment in industrial networks. Two main challenges are scarce number of
radio resources and the problem of energy-efficient communication design while
providing deterministic delay guarantee [20, 21]. The problem of energy-
efficient resource allocation in an industrial wireless network based on 5G
technology is considered in [20]. The deterministic end to end delay
requirement is simplified as a minimum transmit rate constraint for each
sensor and actuator conservatively. However, this simplification is only
applicable if the sampling time in control subsystem is periodic. The number
of data transmissions in the case of using periodic sampling is usually more
than event triggered and self-triggered, consequently, it results in more
power consumption. Determining the sampling time in the control subsystem and
allocating the radio resources in the network subsystem, in a energy-efficient
manner, have been extensively studied in recent years. However, most of
existing works have considered these two issues separately or assumed simple
settings [12, 13, 18, 20]. To the best of our knowledge, the problem of
jointly sampling time determining in the control subsystem and power-efficient
resource allocation in network subsystem based on 5G technologies, in CPSs,
has not been investigated yet.

The main purpose of this paper is to minimize consumed power in industrial
CPSs. To this aim, we jointly determine the sampling time in control subsystem
and allocate the radio resources in network subsystem. We assume an industrial
CPS which is composed of several control plants and Rate Constrained (RC)
users that share a single cell Orthogonal Frequency-Division Multiple Access
(OFDMA) network. The major contributions of this paper are as follows:

  * •

The problem of jointly determining of sampling time in control subsystem and
allocating radio resources in network subsystem, in the CPS, is formulated.
The problem is a multi objective optimization with aims of determining the
next maximum allowable sampling time of control plants and minimizing the
aggregate power consumption in uplink and downlink, considering the dynamics
and desired performance of each control plant, the quality of service
requirement of RC users, power and sub-carrier constraints. The general goal
of the problem is power consumption minimization through decreasing sampling
rate and power-efficient resource allocation.

  * •

To solve the multi objective optimization problem, a novel approach is
proposed. The proposed approach decomposes the optimization problem into two
smaller loosely coupled problems. One of the two decomposed problems is
determining the next maximum allowable sampling time of control plants and
another one is the resource allocation problem. The approach is launched from
solving the control subsystem problem using the self-triggered method. Based
on using the self-triggered method, power saving is improved by minimizing the
number of sampling time and switching some sensors and actuators to sleep mode
when they are not working and wake them up when needed. The approach is
continued to solve power and sub-carrier allocation problem in the OFDMA
wireless network with the goal of power consumption minimization on the
condition of the deterministic delay which guarantees stability and desired
performance of control subsystem.

  * •

The power and sub-carrier allocation problem is formulated with aims of
minimizing aggregate consumed power (transmit and circuit power) in uplink and
downlink considering QoS (Quality of Service) requirement of RC users and
deterministic end to end delay requirement for each control plant. The QoS
requirement is considered in term of a minimum transmit rate. Thanks to self-
triggered method, the deterministic end to end delay requirement guarantees
the stability and desired performance of control subsystem. The problem with
aims of power consumption minimization in uplink and downlink is also a multi
objective optimization. The weighted sum technique is applied to convert two
objective functions into a single objective function. The resultant single
objective problem is mixed-integer nonlinear programing (MINLP). Directly
solving the MINLP problems suffers from high computational complexity. To deal
with it, two steps iterative approach is used which sequentially performs
power allocation and sub-carrier allocation in each iteration. These two steps
are proceeded iteratively until the predetermined convergence criterion is
satisfied.

  * •

We show the power saving effectiveness of using the proposed approach via
simulation results. Simulation results demonstrate that our approach
significantly improves total power saving in the CPS, while it is shown that
all control plants are properly controlled.

Notations: ”$||x||$” stands for Euclidean norm of the vector $x$. $\dot{x}(t)$
represents the derivative of variable x(t) with respect to the original time,
t.
The organization of the paper is as follows: The system model is described in
section II. The problem is formulated in Section III. In Section IV, the
proposed approach is presented. The proposed approach is evaluated through
simulation results in Section V. Finally, the paper is concluded in Section
VI.

##  II SYSTEM MODEL

Consider a model of CPSs comprising 1) $I$ numbers of independent plants, 2)
an OFDMA cellular network and 3) a central controller (Fig.1). Each plant $i$,
$i\in\\{1,2,...,I\\}$, has multiple sensors and multiple actuators that
communicate with the controller via Remote Transmitter Units (RTUs) through
the OFDMA network (Fig.1). RTUs of each plant are connected to the
corresponding sensors and actuators via a wired network in which the time
delay between the RTUs and the relative sensors and actuators is ignorable.
Given plant $i$, let $y_{i,k}$ be the $k^{th}$ sampled data measured by
sensors of plant $i$ which is sent to the controller at sampling time of
$t_{i,k}$ via the corresponding RTU, through the OFDMA network. The controller
calculates an appropriate control input ($u_{i,k}$) based on the measurements
($y_{i,k}=y_{i}(t_{i,k})$). The controller sends $u_{i,k}$ to the relevant RTU
(which is wired to actuators of plant $i$) through the OFDMA network. We
assume that the OFDMA network is responsible for all data exchanges in the
system model, including the communication between RTUs and the controller as
well as the communication between other users which send RC or Best Effort
(BE) packets.
Based on physical infrastructure of ITU’s (International Telecommunication
Union) 5G architecture, RAN (Radio Access Network) real time functions
including computing and data storing must be deployed on the central office
Data Center layer, located close to the base station and users [22]. So we
assume the central controller run on a dedicated platform or general platform
in the central office Data Center layer.
This system model is very appropriate for plenty of emerging industrial cyber
physical systems such as smart industry automation and smart grid. In the
following subsections, the CPS model are explained in details.

Figure 1: System model of the considered CPS

###  II-A Plants

Each plant $i$ has $n_{i}$ states, $m_{i}$ actuators and $q_{i}$ measurable
states ($q_{i}\leq n_{i}$) which measured by $q_{i}$ sensors. The dynamics of
each plant $i$ is linear in the form of:

$\begin{split}\dot{x}&=A_{i}x_{i}(t)+B_{i}u_{i}(t)+d_{i}(t),\\\
y_{i}(t)&=C_{i}x_{i}(t),\end{split}$ (1)

where $x_{i}(t)\in\mathbb{R}^{n_{i}}$, $u_{i}(t)\in\mathbb{R}^{m_{i}}$ and
$y_{i}(t)\in\mathbb{R}^{q_{i}}$ are the vector of the plant states, the
control input and plant output (measured by sensors) at time
$t\in\mathbb{R}^{+}$, respectively. $d_{i}(t)\in\mathbb{R}^{m_{i}}$ is used to
model measurement and actuation disturbances and other sources of uncertainty.
$A_{i}\in\mathbb{R}^{n_{i}\times n_{i}}$, $B_{i}\in\mathbb{R}^{n_{i}\times
m_{i}}$ and $C_{i}\in\mathbb{R}^{q_{i}\times n_{i}}$ are state, control and
output matrix, respectively. We suppose that the disturbance is non measurable
and bounded $\left\lVert d_{i}(t)\right\rVert\leq d$.

###  II-B Controller

The central controller receives measurements ($y_{i,k}$), sent by sensors, and
uses them to calculate the control input ($u_{i,k}$), afterward sends the
control input to the relevant RTU. The control input is calculated according
to the following control law:

$u_{i,k}=k_{i}x_{i}(t=t_{i,k}),$ (2)

where $k_{i}$ is calculated by the controller for plant $i$ such that the
dynamics of plant $i$ is stable ($||x(t)||<\infty:\forall t$). The control
input used by actuators is held constant by ZOH (Zero Order Hold) between two
successive control input updates. Therefore the control input used by
actuators is a piecewise constant signal:

$u_{i}(t)=u_{i,k},\quad t\in[t_{i,k}+\Delta^{total}_{i,k},\quad
t_{i,k+1}+\Delta^{total}_{i,k+1}),\\\ $ (3)

where $\Delta^{total}_{i,k}$ denotes the total delay (end to end delay) for
plant $i$ in the $k^{th}$ sampling. Considering (2) and (3), (1) may be
rewritten as bellow for
$t\in[t_{i,k}+\Delta^{total}_{i,k},t_{i,k+1}+\Delta^{total}_{i,k+1})$:

$\begin{split}\dot{x_{i}}&=A_{i}x_{i}(t)+B_{i}u_{i,k}+d_{i}(t),\\\
\dot{x_{i}}&=(A_{i}+B_{i})x_{i}(t)+B_{i}k_{i}e_{i,k}(t)+d_{i}(t),\\\
e_{i,k}(t)&=x_{i}(t)-x_{i}(t=t_{i,k}),\quad
t\in[t_{i,k}+\Delta^{total}_{i,k},t_{i,k+1}+\Delta^{total}_{i,k+1}),\end{split}$
(4)

where $e_{i,k}(t)$ denotes the deviation between $x_{i}(t)$ and
$x_{i}(t=t_{i,k})$ during
$t\in[t_{i,k}+\Delta^{total}_{i,k},t_{i,k+1}+\Delta^{total}_{i,k+1})$ which is
considered as an index of plant performance. $e_{i,k}(t)$ is bounded by
$\delta_{i}$, that is:

$||e_{i,k}(t)||\leqslant\delta_{i},$ (5)

where $\delta_{i}$ denotes the maximum deviation of $x_{i}(t)$ from
$x_{i}(t=t_{i,k})$ during
$t\in[t_{i,k}+\Delta^{total}_{i,k},t_{i,k+1}+\Delta^{total}_{i,k+1})$ which is
tolerated by plant $i$.

###  II-C Network

Consider a single cell OFDMA network which contains two types of users namely,
RC users and the control system users (Time Constrained-TC). RC users send
Rate Constrained or Best Effort packets and are irrelevant to the control
system. It is assumed that the network consists of M uplink users and N users
in downlink, forming the sets $\mathcal{U}_{u}=\\{1,2,...,M\\}$ and
$\mathcal{U}_{d}=\\{1,2,...,N\\}$, respectively. Uplink users include a set of
$M^{c}$ RTUs (control system users) denoted by
$\mathcal{U}^{c}_{u}=\\{1,2,...,M^{c}\\}$ and a set of $M^{u}$ RC users
denoted by $\mathcal{U}^{u}_{u}=\\{1,2,...,M^{u}\\}$. Users in downlink
include a set of $N^{c}$ RTUs (control system users) denoted by
$\mathcal{U}^{c}_{d}=\\{1,2,...,N^{c}\\}$ and a set of $N^{u}$ RC users
denoted by $\mathcal{U}^{u}_{d}=\\{1,2,...,N^{u}\\}$. It is assumed that each
plant enjoys two RTUs, one for sending its sensors data in uplink and another
one for receiving its control input data in downlink, hence $M^{c}=I$ and
$N^{c}=I$. The total bandwidth is divided in $L$ sub-carriers in set
$\mathcal{L}=\\{1,2,...,L\\}$. The channel gain for the uplink user $m$ and
the downlink user $n$ on the sub-carrier $l$ are represented by $g_{m,l}$ and
$g_{n,l}$, respectively. Let $p^{u}_{m,l}$ be the transmit power of the uplink
user $m$ on the sub-carrier $l$ and $p^{d}_{n,l}$ be the BS transmit power for
the downlink user $n$ on the sub-carrier $l$, consecutively. In addition
$\overline{P}_{m}$ and $\overline{P}_{BS}$ denote the peak transmit power of
the downlink user $m$ and the BS, respectively ($\forall
m:\sum_{l=1}^{L}p^{u}_{m,l}\leq\overline{P}_{m}$, and
$\sum_{n=1}^{N}\sum_{l=1}^{L}p^{d}_{n,l}\leq\overline{P}_{BS}$).

The uplink transmission rate for the uplink user $m$ is given by:

$R^{u}_{m}=\sum_{l=1}^{L}\alpha_{m,l}wlog(1+\frac{p^{u}_{m,l}g_{m,l}}{N_{0}}),$
(6)

where $\alpha_{m,l}$ is a binary variable and $\alpha_{m,l}=1$ if sub-carrier
$l$ is allocated to the uplink user $m$; otherwise, $\alpha_{m,l}=0$.
Additionally $N_{0}$ is the power of additive white Gausian noise and $w$
denotes the bandwidth of each sub-carrier. The downlink transmission rate for
the downlink user $n$ is given by:

$R^{d}_{n}=\sum_{l=1}^{L}\beta_{n,l}wlog(1+\frac{p^{d}_{n,l}g_{n,l}}{N_{0}}),$
(7)

where $\beta_{n,l}$ is a binary variable and $\beta_{n,l}=1$ if sub-carrier
$l$ is assigned to the downlink user $n$; otherwise, $\beta_{n,l}=0$.
The QoS requirement of RC users is stated in term of a minimum transmit rate.
The required QoS of RC users in uplink and downlink is given by:

$\displaystyle R^{u}_{m}$ $\displaystyle\geq\overline{R}^{u}_{m},\forall
m\in\mathcal{U}^{u}_{u},$ (8) $\displaystyle R^{d}_{n}$
$\displaystyle\geq\overline{R}^{d}_{n},\forall n\in\mathcal{U}^{u}_{d}.$

Given $k^{th}$ sampling of plant $i$, we consider three types of delays in
this system model including uplink transmission delay, downlink transmission
delay and computation delay, denoted by $\Delta^{u}_{i,k}$, $\Delta^{d}_{i,k}$
and $\Delta^{comp}_{i,k}$, respectively. Let $\Delta^{comp}_{max}$ be the
maximum computation delay that is
$\Delta^{comp}_{i,k}\leq\Delta^{comp}_{max}$. Conservatively, we assume
$\Delta^{total}_{i,k}$ is given by:

$\Delta^{total}_{i,k}=\Delta^{u}_{i,k}+\Delta^{d}_{i,k}+\Delta^{comp}_{max}.$
(9)

The RC and TC users send and receive data with size of $L_{RC}$ and $L_{TC}$
bit, respectively. At each $k^{th}$ sampling of plant $i$, uplink and downlink
transmission delay are given by:

$\displaystyle\Delta^{u}_{i,k}$ $\displaystyle=\frac{L_{TC}}{R^{u}_{i}},$ (10)
$\displaystyle\Delta^{d}_{i,k}$ $\displaystyle=\frac{L_{TC}}{R^{d}_{i}},$

where $R^{u}_{i}$ and $R^{d}_{i}$ are the uplink and downlink transmission
rate for the user of plant $i$ given by (6) and (7), respectively. It is
assumed that uplink users’ total power consumption of the CPS may be obtained
in terms of uplink power matrix
($\textbf{P}^{u}{\mathop{:}}=[{p^{u}_{m,l}}_{(m\in U_{u},l\in\textit{L})}]$)
and uplink sub-carrier assignment matrix
($\textbf{A}^{u}{\mathop{:}}=[{\alpha_{m,l}}_{(m\in U_{u},l\in\textit{L})}]$)
as the following:

$P^{u}_{total}(\textbf{A}^{u},\textbf{P}^{u})=\sum_{m=1}^{M}P^{cst}_{m}+\sum_{m=1}^{M}\sum_{l=1}^{L}\alpha_{m,l}p^{u}_{m,l},$
(11)

where $P^{cst}_{m}$ is the constant circuit power consumed by uplink user $m$.
Also the BS’s total power consumption at the modeled CPS is given by:

$P^{BS}_{total}(\textbf{B}^{d},\textbf{P}^{d})=P^{cst}_{BS}+\sum_{n=1}^{N}\sum_{l=1}^{L}\beta_{n,l}p^{d}_{n,l},$
(12)

where $P^{cst}_{BS}$ is the constant circuit power consumed by BS,
($\textbf{P}^{d}{\mathop{:}}=[{p^{d}_{n,l}}_{(n\in U_{d},l\in\textit{L})}]$)
and ($\textbf{B}^{d}{\mathop{:}}=[{\beta_{n,l}}_{(n\in
U_{d},l\in\textit{L})}]$) are downlink power matrix and downlink sub-carrier
assignment matrix, respectively.

##  III Problem Formulation

We consider the problem of consumed power minimization in the CPS. To realize
consumed power minimization in the CPS, the problem of jointly decreasing the
number of sampling time in control subsystem and allocating the radio
resources in network subsystem, in a power-efficient manner, is considered.
The objectives are determining the next maximum allowable sampling instant of
each plant, $t_{i,k}$, and minimizing aggregate power consumption in uplink
and downlink subject to the dynamics and desired performance requirement of
each plant, QoS of RC users, power and sub-carrier constraints. The multi
objective optimization problem is formulated as follows:

$\displaystyle{\text{maximize }}t_{i,k}\quad:\forall i\in\\{1,2,...,I\\}$ (13)
$\displaystyle\underset{\textbf{A}^{u},\textbf{P}^{u}}{\text{minimize
}}P^{u}_{total}$
$\displaystyle\underset{\textbf{B}^{d},\textbf{P}^{d}}{\text{minimize
}}P^{BS}_{total}$ $\displaystyle\text{subject to}:$ $\displaystyle
C_{1}:\dot{x_{i}}=(A_{i}+B_{i})x_{i}(t)+B_{i}k_{i}e_{i,k}(t)+d_{i}(t)$
$\displaystyle\quad:t\in[t_{i,k-1}+\Delta^{total}_{i,k-1},t_{i,k}+\Delta^{total}_{i,k})$
$\displaystyle C_{2}:||e_{i,k}(t)||\leqslant\delta_{i}$ $\displaystyle
C_{3}:\sum_{m=1}^{M}\alpha_{m,l}\leq 1,\forall l\in\mathcal{L}$ $\displaystyle
C_{4}:\sum_{n=1}^{N}\beta_{n,l}\leq 1,\forall l\in\mathcal{L}$ $\displaystyle
C_{5}:\alpha_{i,k}\in\\{0,1\\},\forall i,k$ $\displaystyle
C_{6}:\beta_{i,k}\in\\{0,1\\},\forall i,k$ $\displaystyle
C_{7}:\sum_{l=1}^{L}p^{u}_{m,l}\alpha_{m,l}\leq\overline{P}_{m},\forall
m\in\mathcal{U}_{u}$ $\displaystyle
C_{8}:\sum_{n=1}^{N}\sum_{l=1}^{L}p^{d}_{n,l}\beta_{n,l}\leq\overline{P}_{BS}$
$\displaystyle C_{9}:\frac{1}{2}R^{u}_{m}\geq\overline{R}^{u}_{m},\forall
m\in\mathcal{U}^{u}_{u}$ $\displaystyle
C_{10}:\frac{1}{2}R^{d}_{n}\geq\overline{R}^{d}_{n},\forall
n\in\mathcal{U}^{u}_{d},$

where the first two constraints denote the dynamics and the desired
performance requirement for each plant, respectively. The constraints $C_{3}$,
$C_{4}$, $C_{5}$ and $C_{6}$ correspond to the exclusive sub-carrier
allocation in uplink and downlink, respectively. It is noticeable that the
same subset of sub-carriers is used in uplink and downlink during the distinct
time intervals. The constraints $C_{7}$ and $C_{8}$ relate to the maximum
tolerated power consumption in each uplink user and BS. The last two
constraints guarantee the required QoS of RC users in uplink and downlink
consequently. Due to the similarity of uplink and downlink time length, we
share the total time length of transmission between uplink and downlink
equally. The factor of $\frac{1}{2}$ in $C_{7}$ and $C_{8}$ denotes this
sharing.
The optimization problem in (13) is a multi objective mixed integer nonlinear
programing. To address the multi objective optimization, we decompose the
original optimization problem into two smaller loosely coupled problems, as is
explained in the following section.

##  IV The Proposed Approach

We decompose the optimization problem into two smaller loosely coupled
problems. One of the two decomposed problems is determining the next maximum
allowable sampling time of control plants and another one is the power-
efficient resource allocation problem (Fig.2). First, the control subsystem
problem is addressed by using the self-triggered method. We use self-triggered
method to reduce the transmission number of control subsystem nodes and
increase their sleeping time. By determining the maximum allowable sampling
instant ($t_{i,k}$), it is possible to put some nodes in sleep mode when they
are not needed. Then, having determined the number of registered plants to
communicate at $k^{th}$ sampling time and their relevant end to end delay
requirement, the problem of resource allocation is addressed. This problem is
a joint power and sub-carrier allocation in uplink and downlink with aims of
uplink and downlink power minimization, considering the QoS requirement of RC
users and the deterministic end to end delay requirement for the plants which
was determined in previous step ($\Delta^{i}_{max}$). Thanks to self-triggered
method, the deterministic delay requirement guarantees stability and desired
performance of the plants.
It is possible to enlarge feasible region by increase in $\delta_{i}$
parameter such that $\delta_{i}\leq\delta_{i,max}$. By increasing
$\delta_{i}$, the inter-sampling time of the plant $i$ will increase and the
number of plants wanting to communicate will decrease. $R^{u}_{m}$ and
$R^{d}_{n}$ are determined through the solution of the resource allocation
problem and so $\Delta^{total}_{i,k}$ is revealed. $\Delta^{total}_{i,k}$ will
be used to determine the next sampling instant $t_{i,k+1}$.

Figure 2: Proposed approach

###  IV-A The Sampling Time Determining Problem

In the problem of determining the sampling time, we aim to determine the next
maximum allowable sampling time of control plants. The problem is formulated
as:

$\displaystyle{\text{maximize }}t_{i,k}\quad:\forall i\in\\{1,2,...,I\\}$ (14)
$\displaystyle\text{subject to}:$ $\displaystyle
C_{1}:\dot{x_{i}}=(A_{i}+B_{i})x_{i}(t)+B_{i}k_{i}e_{i,k}(t)+d_{i}(t)$
$\displaystyle\quad:t\in[t_{i,k-1}+\Delta^{total}_{i,k-1},t_{i,k}+\Delta^{total}_{i,k}),$
$\displaystyle C_{2}:||e_{i,k}(t)||\leqslant\delta_{i}.$

The self-triggered method is used to solve this problem [16]. We use the
proposed self-triggered method in [18]. According to self-triggered method,
the next maximum allowable sampling time is determined using the plant
dynamics model, its current state ($x_{i}(t=t_{i,k})$) and the total delay
(end to end delay) in the previous sampling time while the stability and
desired performance of the plant is guaranteed [18]. As stated, the proposed
self-triggered sampling method in [18] is adopted:

$\displaystyle t_{i,k+1}=$ $\displaystyle
t_{i,k}+min(\gamma(x_{i,k-1},x_{i,k},\hat{d}_{i,k-1},\hat{d}_{i,k},$ (15)
$\displaystyle\Delta^{total}_{i,k}),h_{max}),$

where $\hat{d}_{i,k}$ and $\hat{d}_{i,k-1}$ denote the estimation of
$d_{i}(t=t_{k})$ and $d_{i}(t=t_{k-1})$ respectively, and $h_{max}$ is the
maximum acceptable interval between two consecutive sampling times for plant
$i$ by which the stability of plant $i$ and
$||e_{i,k}(t)||\leqslant\delta_{i}$ are satisfied.
$\gamma(x_{i,k-1},x_{i,k},\hat{d}_{i,k-1},\hat{d}_{i,k},\Delta^{total}_{i,k})$
is defined as follows:

$\displaystyle\gamma(x$
$\displaystyle{}_{i,k-1},x_{i,k},\hat{d}_{i,k-1},\hat{d}_{i,k},\Delta^{total}_{i,k})$
(16) $\displaystyle\mathrel{\mathop{:}}=$
$\displaystyle\frac{1}{||A||}\ln\left(\frac{\phi(x_{i,k},\hat{d}_{i,k})}{\upsilon(x_{i,k-1},x_{i,k},\hat{d}_{i,k-1},\hat{d}_{i,k},\Delta^{total}_{i,k})}\right)$
$\displaystyle+$ $\displaystyle\Delta^{total}_{i,k}-\Delta^{i}_{max},$

where $\Delta^{i}_{max}$ is the maximum total delay (delay between sending
instant of measurements and receiving instant of input by actuators) tolerated
by plant $i$ for all $k$.
$\phi(x_{i,k-1},x_{i,k},\hat{d}_{i,k-1},\hat{d}_{i,k},\Delta^{total}_{i,k})\mathrel{\mathop{:}}=||A_{i}||\delta_{i}+||(A_{i}+B_{i}k_{i,k})x_{i,k}||+||\hat{d}_{i,k}||$
and $\upsilon$ is defined as follows:

$\displaystyle\upsilon(x_{i,k-1}$
$\displaystyle,x_{i,k},\hat{d}_{i,k-1},\hat{d}_{i,k},\Delta^{total}_{i,k})$
(17) $\displaystyle\mathrel{\mathop{:}}=$
$\displaystyle||A_{i}x_{i,k}-B_{i}k_{i,k}x_{k-1}||+||\hat{d}_{i,k}||*$
$\displaystyle(\exp(||A_{i}||\Delta^{total}_{i,k})-1)||(A_{i}+B_{i}k_{i,k})x_{i,k}||$
$\displaystyle+$ $\displaystyle||\hat{d}_{i,k}||.$

If the controller does not receive any new update from plant $i$ after
$h_{max}$ seconds, it will notify the relevant RTU to send the update.
In summing up, the controller uses (15) to determine the time of sampling
$t_{i,k}$. The stability of plant $i$ and the desired performance are
satisfied if $\Delta^{total}_{i,k}\leq\Delta^{i}_{max}$. The controller awakes
the RTU wired to sensors of plant $i$, at time of $t_{i,k}$, to pick $k^{th}$
measurement of plant $i$ ($y_{i,k}=y_{i}(t_{i,k})$) and calculates an
appropriate control input ($u_{i,k}$) based on the measurements
($y_{i,k}=y_{i}(t_{i,k})$). The controller sends $u_{i,k}$ to the relevant RTU
(which wired to plant $i$ actuators) through the OFDMA network.

###  IV-B Power-Efficient Resource Allocation Problem

The problem of power-efficient resource allocation is formulated as:

$\displaystyle\underset{\textbf{A}^{u},\textbf{P}^{u}}{\text{minimize
}}P^{u}_{total}$ (18)
$\displaystyle\underset{\textbf{B}^{d},\textbf{P}^{d}}{\text{minimize
}}P^{BS}_{total}$ $\displaystyle\text{subject to}:$ $\displaystyle
C_{3}:\sum_{m=1}^{M}\alpha_{m,l}\leq 1,\forall l\in\mathcal{L}$ $\displaystyle
C_{4}:\sum_{n=1}^{N}\beta_{n,l}\leq 1,\forall l\in\mathcal{L}$ $\displaystyle
C_{5}:\alpha_{i,k}\in\\{0,1\\},\forall i,k$ $\displaystyle
C_{6}:\beta_{i,k}\in\\{0,1\\},\forall i,k$ $\displaystyle
C_{7}:\sum_{l=1}^{L}p^{u}_{m,l}\alpha_{m,l}\leq\overline{P}_{m},\forall
m\in\mathcal{U}_{u}$ $\displaystyle
C_{8}:\sum_{n=1}^{N}\sum_{l=1}^{L}p^{d}_{n,l}\beta_{n,l}\leq\overline{P}_{BS}$
$\displaystyle C_{9}:\frac{1}{2}R^{u}_{m}\geq\overline{R}^{u}_{m},\forall
m\in\mathcal{U}^{u}_{u}$ $\displaystyle
C_{10}:\frac{1}{2}R^{d}_{n}\geq\overline{R}^{d}_{n},\forall
n\in\mathcal{U}^{u}_{d},$ $\displaystyle
C_{11}:\Delta^{comp}_{max}+\frac{L_{TC}}{\frac{1}{2}R^{u}_{m}}+\frac{L_{TC}}{\frac{1}{2}R^{d}_{n}}\leq\Delta^{i}_{max},$
$\displaystyle\forall m\in\mathcal{U}^{c}_{u},\forall
n\in\mathcal{U}^{c}_{d},\forall i\in\mathcal{I}^{r}$

where $C_{11}$ is the deterministic end to end delay requirement for the set
of plants ($\mathcal{I}^{r}$) determined by using (15) to solve
(LABEL:eqn:predict). $\mathcal{I}^{r}$ denotes the set of plants that their
RTUs should communicate with controller in the time of solving (18). Thanks to
self-triggered method, the stability and desired performance of the plants are
ensured provided that the $C_{11}$ constraint is satisfied. As stated before,
the factor of $\frac{1}{2}$ in $C_{9}$, $C_{10}$ and $C_{11}$ denotes the
equal time length sharing between uplink and downlink.
We use weighted sum technique to transform two objective functions into a
single objective function [23]. Therefore the transformed objective function
can be written as:

$\displaystyle\underset{\textbf{A}^{u},\textbf{B}^{d},\textbf{P}^{u},\textbf{P}^{d}}{\text{minimize}}aP^{u}_{total}+(1-a)P^{BS}_{total}$
(19) $\displaystyle\text{subject to}:C_{3}-C_{11},$

where $0\leq a\leq 1$ and denotes the priority of objective functions. The
optimization problem in (19) is non-convex due to coupled continuous and
integer variable in both the objective function as well as the constraints.
Notice that (19) is a Mixed Integer Nonlinear Programming (MINLP) which can be
solved using exhaustive search over all sub-carrer assignment choices. However
the complexity of MINLP problems raises fast when the number of users and sub-
carriers increases. In order to tackle the computational complexity of (19),
we transform the problem into two subproblems: i) sub-carrier assignment and
ii) power allocation. In other words, for a given power allocation, elements
of sub-carrier assignment matrices ($\textbf{A}^{u}$ and $\textbf{B}^{d}$) are
considered as variables in the first step. Thereafter, the power matrices of
users and BS ($\textbf{P}^{u}$ and $\textbf{P}^{d}$) are assumed variables
while using sub-carrier assignments given in previous step. The two steps are
proceeded iteratively until convergence is attained [24, 25]. This approach
has been widely used for resource allocation problems [24, 25, 26, 27].

Sub-carrier Assignment Subproblem: The sub-carrier assignment subproblem for a
given power allocation can be written as follows:

$\displaystyle\underset{\textbf{A}^{u},\textbf{B}^{d}}{\text{minimize
}}aP^{u}_{total}+(1-a)P^{BS}_{total}$ (20) $\displaystyle\text{subject
to}:C_{3}-C_{11}.$

The problem of (LABEL:eqn:opt3) is integer nonlinear programming. We solve
this problem using MATLAB/CVX with MOSEK solver [28].

Power Allocation subproblem: Given a sub-carrier assigned in (LABEL:eqn:opt3),
the power allocation subproblem is stated as:

$\displaystyle\underset{\textbf{P}^{u},\textbf{P}^{d}}{\text{minimize
}}aP^{u}_{total}+(1-a)P^{BS}_{total}$ (21) $\displaystyle\text{subject
to}:C_{7}-C_{11}.$

The objective function is linear with respect to the transmit power variables.
This subproblem is convex optimization problem and MATLAB/CVX can be applied
to numerically solve the optimal problem.

##  V Simulation Result and Discussion

The efficacy of the proposed approach is evaluated through numerical
simulations. We consider three independent plants that communicate over an
OFDMA cellular network and the simulation time of about 50s for each
simulation. The control plants and network specification are explained below:

OFDMA cellular network: We consider a single cell network where 5 RC users, 6
sensors and 6 actuators are randomly spread within the square cell. The
channel gain for each uplink (downlink) user is modeled as
$g_{m,l}=\overline{d}_{m,l}^{-3}\overline{h}$
($g_{n,l}=\overline{d}_{n,l}^{-3}\overline{h}$), where $\overline{d}_{m,l}$
($\overline{d}_{n,l}$) denote the distance between the transmitter and
receiver, $\overline{h}$ is the attenuation factor that indicates power
variations and $\overline{h}=0.09$. The simulation parameters of the
considered network are listed in Table I.

TABLE I: Values of Network parameters. Parameter | Value
---|---
Bandwidth of each sub carrier | 180KHz
Maximum transmit power of each user ($\overline{P}_{m}$) | 23dBm
Maximum transmit power of BS ($\overline{P}_{BS}$) | 43dBm
Constant Power of each user ($P^{cst}_{m}$) | 0.1dBm
Constant Power of BS ($P^{cst}_{BS}$) | 20dBm
Noise Power ($N_{0}$) | -62.24dBm
Distance between any user to BS | 10-50m
Minimum data rate for RC users in uplink ($\overline{R}^{u}_{m}$) | 50 bit/s
Minimum data rate for RC users in downlink ($\overline{R}^{d}_{n}$) | 100 bit/s
Maximum total delay ($\Delta^{i}_{max}$) | 1s
Length of TC users ($L_{TC}$) | 70 bit
Number of sub-carries (L) | 16

Control Plants: Similar to [18], we consider three control plants as bellow:

####  V-1 First Control plant

The state and control matrix of first control plant are set as bellow:

$\displaystyle A_{1}=\begin{bmatrix}-0.1&0.05\\\
0.2&0.1\end{bmatrix},B_{1}=\begin{bmatrix}0\\\ 1\end{bmatrix}$ (22)

The eigenvalues of first closed control loop are placed in
$\lambda_{1}^{1}(A_{1}+B_{1}K_{1})=-0.25$,
$\lambda_{2}^{1}(A_{1}+B_{1}K_{1})=-0.18$. The upper bound of external
disturbance is considered 0.6 and the initial values are set
$x_{1}^{1}(0)=-20$, $x_{2}^{1}(0)=15$.

####  V-2 Second Control plant

The state and control matrix of second control plant are set as bellow:

$\displaystyle A_{2}=\begin{bmatrix}0.01&0.2\\\
0.03&0\end{bmatrix},B_{2}=\begin{bmatrix}1\\\ 1\end{bmatrix}$ (23)

The eigenvalues of second closed control loop are placed in
$\lambda_{1}^{2}(A_{2}+B_{2}K_{2})=-0.15$,
$\lambda_{2}^{2}(A_{2}+B_{2}K_{2})=-0.3$. The upper bound of external
disturbance is considered 1.2 and the initial values are set
$x_{1}^{2}(0)=-12$, $x_{2}^{2}(0)=12$.

####  V-3 Third Control plant

The state and control matrix of third control plant are set as bellow:

$\displaystyle A_{3}=\begin{bmatrix}0.2&0.01\\\
0.3&-0.8\end{bmatrix},B_{2}=\begin{bmatrix}1\\\ 2\end{bmatrix}$ (24)

The eigenvalues of third closed control loop are placed in
$\lambda_{1}^{3}(A_{3}+B_{3}K_{3})=-0.4$,
$\lambda_{2}^{3}(A_{2}+B_{2}K_{2})=-0.6$. The upper bound of external
disturbance is considered 0.55 and the initial values are set
$x_{1}^{3}(0)=-5$, $x_{2}^{3}(0)=4$.
Each plant has two sensors and two actuators. We assume one central controller
in the central office layer which is located close to BS.
First, in table II, we compare the number of sensors transmission of each
plant when the proposed approach (using self-triggered method) and periodic
method are applied for a simulation time of 50s.

TABLE II: The number of sensors transmission of each plant when the proposed approach and periodic method are user. Plant | Proposed Approach | Periodic Method
---|---|---
First Plant | 37 | 555
Second Plant | 31 | 555
Third Plant | 25 | 555

We compare responses of the CPS when proposed approach and periodic method are
used in Fig.3. The results demonstrate that all control plants are properly
controlled. Although the number of transmissions is reduced considerably in
proposed approach, the control plants responses are similar to periodic method
implementation.

Figure 3: Analogy between CPS responses of proposed approach and periodic
method implementation, (a) First control plant response, (b) Second control
plant response, (c) Third control plant response

Figs.4 and 5 illustrate the total power consumption
($P^{u}_{total}+P^{BS}_{total}$), the power consumption in uplink
($P^{u}_{total}$) and downlink ($P^{u}_{total}$) for about 50s simulation
time, when periodic method and our proposed approach are applied,
respectively. In both Figs.4 and 5, the minimum point of the all power
consumption types is taken place at $a=0.5$, when the priority of objective
functions ($a$) is equal. Fig.6 clearly illustrates the comparison of the
power consumption in our proposed approach with Periodic method
Implementation. It is observed that the power consumption in our proposed
approach is significantly lower than periodic method implementation.

Figure 4: Power consumption in periodic method Implementation

Figure 5: Power consumption in our proposed approach

Figure 6: Comparison of the total power consumption in our proposed approach
with periodic method implementation

##  VI Conclusion and Future Work

We have studied the problem of power minimization in the industrial CPS. To
realize power minimization in the CPS, the problem of jointly decreasing
sampling times in control subsystem and allocating power-efficient radio
resources in network subsystem is formulated. The objectives of the problem
are determining the next maximum allowable sampling instant of each control
plant and power minimization in uplink and downlink, considering the dynamics
and desired performance of each control plant, the QoS requirement of RC
users, power and sub-carrier constraints. We proposed a novel approach to
address this multi objective problem. The proposed approach decomposes the
optimization problem into two smaller loosely coupled problems. We showed that
our proposed approach considerably decreases power consumption, whereas the
stability statuses and control performances of control plants responses are
satisfied. In this paper, we conservatively assumed a fixed maximum value for
computation delay variable. Considering an actual variable for computation
delay is a complicated scenario which will be studied in our future works.

## References

  * [1] M. Wollschlaeger, T. Sauter, and J. Jasperneite, “The future of industrial communication: Automation networks in the era of the internet of things and industry 4.0,” _IEEE Industrial Electronics Magazine_ , vol. 11, pp. 17–27, 2017.
  * [2] Z. Meng, Z. Wu, C. Muvianto, and J. Gray, “A data-oriented M2M messaging mechanism for industrial IOT applications.” _IEEE Internet of Things Journal_ , vol. 4, no. 1, pp. 236–246, 2017.
  * [3] C. Perera, C. H. Liu, and S. Jayawardena, “The emerging internet of things marketplace from an industrial perspective: A survey,” _IEEE Transactions on Emerging Topics in Computing_ , vol. 3, no. 4, pp. 585–598, 2015\.
  * [4] H. P. Breivold and K. Sandström, “Internet of things for industrial automation–challenges and technical solutions,” in _Data Science and Data Intensive Systems (DSDIS), 2015 IEEE International Conference on_. IEEE, 2015, pp. 532–539.
  * [5] L. Da Xu, W. He, and S. Li, “Internet of things in industries: A survey,” _IEEE Transactions on industrial informatics_ , vol. 10, no. 4, pp. 2233–2243, 2014.
  * [6] Z. Sheng, C. Mahapatra, C. Zhu, and V. Leung, “Recent advances in industrial wireless sensor networks towards efficient management in IOT,” _IEEE access_ , vol. 3, pp. 622–637, 2015.
  * [7] L. L. Bello, J. Åkerberg, M. Gidlund, and E. Uhlemann, “Guest editorial special section on new perspectives on wireless communications in automation: From industrial monitoring and control to cyber physical systems,” _IEEE Transactions on Industrial Informatics_ , vol. 13, no. 3, pp. 1393–1397, 2017.
  * [8] F. Tang, C. Tang, Y. Yang, L. T. Yang, T. Zhou, J. Li, and M. Guo, “Delay minimized routing in mobile cognitive networks for time critical applications,” _IEEE Transactions on Industrial Informatics_ , vol. 13, no. 3, pp. 1398–1409, 2017.
  * [9] C. Lu, A. Saifullah, B. Li, M. Sha, H. Gonzalez, D. Gunatilaka, C. Wu, L. Nie, and Y. Chen, “Real time wireless sensor actuator networks for industrial cyber physical systems,” _Proceedings of the IEEE_ , vol. 104, no. 5, pp. 1013–1024, 2016.
  * [10] J. Wang, W. Dong, Z. Cao, and Y. Liu, “On the delay performance in a large scale wireless sensor network: Measurement, analysis, and implications,” _IEEE/ACM Transactions on Networking_ , vol. 23, no. 1, pp. 186–197, 2015\.
  * [11] K. Yu, M. Gidlund, J. Åkerberg, and M. Björkman, “Performance evaluations and measurements of the realflow routing protocol in wireless industrial networks,” _IEEE Transactions on Industrial Informatics_ , vol. 13, no. 3, pp. 1410–1420, 2017.
  * [12] A. Saifullah, C. Wu, P. B. Tiwari, Y. Xu, Y. Fu, C. Lu, and Y. Chen, “Near optimal rate selection for wireless control systems,” _ACM Transactions on Embedded Computing Systems (TECS)_ , vol. 13, no. 4s, p. 128, 2014.
  * [13] K. Suto, H. Nishiyama, N. Kato, and C.-W. Huang, “An energy-efficient and delay aware wireless computing system for industrial wireless sensor networks,” _IEEE Access_ , vol. 3, pp. 1026–1035, 2015.
  * [14] N. Kaur and S. K. Sood, “An energy-efficient architecture for the internet of things (IOT),” _IEEE Systems Journal_ , vol. 11, no. 2, pp. 796–805, 2017.
  * [15] X. Chen, W. Chen, J. Lee, and N. B. Shroff, “Delay optimal buffer aware scheduling with adaptive transmission,” _IEEE Transactions on Communications_ , vol. 65, no. 7, pp. 2917–2930, 2017.
  * [16] Z. Heng, P. Chen, Z. Jin, and Z. Chu, “Event-triggered control in networked control systems: A survey,” in _Control and Decision Conference (CCDC), 2015 27th Chinese_. IEEE, 2015, pp. 3092–3097.
  * [17] W. Heemels, K. H. Johansson, and P. Tabuada, “An introduction to event-triggered and self-triggered control,” in _Decision and Control (CDC), 2012 IEEE 51st Annual Conference on_. IEEE, 2012, pp. 3270–3285.
  * [18] U. Tiberi, C. Fischione, K. H. Johansson, and M. D. Di Benedetto, “Energy-efficient sampling of networked control systems over IEEE 802.15. 4 wireless networks,” _Automatica_ , vol. 49, no. 3, pp. 712–724, 2013.
  * [19] H. Li, Z. Han, A. D. Dimitrovski, and Z. Zhang, “Data traffic scheduling for cyber physical systems with application in voltage control of distributed generations: A hybrid system framework,” _IEEE Systems Journal_ , vol. 8, no. 2, pp. 542–552, 2014.
  * [20] S. Li, Q. Ni, Y. Sun, G. Min, and S. Al-Rubaye, “Energy-efficient resource allocation for industrial cyber physical IOT systems in 5G era,” _IEEE Transactions on Industrial Informatics_ , vol. 14, no. 6, pp. 2618–2628, 2018.
  * [21] R. Atat, L. Liu, H. Chen, J. Wu, H. Li, and Y. Yi, “Enabling cyber physical communication in 5G cellular networks: challenges, spatial spectrum sensing, and cyber security,” _IET Cyber Physical Systems: Theory & Applications_, vol. 2, no. 1, pp. 49–54, 2017.
  * [22] “5G network architecture a high-level perspective,” _HUAWEI TECHNOLOGIES CO., LTD_ , 2016.
  * [23] J.-H. Cho, Y. Wang, R. Chen, K. S. Chan, and A. Swami, “A survey on modeling and optimizing multi-objective systems,” _IEEE Communications Surveys & Tutorials_, vol. 19, no. 3, pp. 1867–1901, 2017.
  * [24] L. Venturino, N. Prasad, and X. Wang, “Coordinated scheduling and power allocation in downlink multicell OFDMA networks,” _IEEE Transactions on Vehicular Technology_ , vol. 58, no. 6, pp. 2835–2848, 2009.
  * [25] N. Mokari, F. Alavi, S. Parsaeefard, and T. Le-Ngoc, “Limited feedback resource allocation in heterogeneous cellular networks.” _IEEE Trans. Vehicular Technology_ , vol. 65, no. 4, pp. 2509–2521, 2016.
  * [26] K. Son, S. Lee, Y. Yi, and S. Chong, “Refim: A practical interference management in heterogeneous wireless access networks,” _arXiv preprint arXiv:1105.0738_ , 2011.
  * [27] D. Tweed, M. Derakhshani, S. Parsaeefard, and T. Le-Ngoc, “Outage-constrained resource allocation in uplink NOMA for critical applications,” _IEEE Access_ , vol. 5, pp. 27 636–27 648, 2017.
  * [28] “Mosek optimization software,” https://mosek.com, May 2019.


# RydbergGPT

David Fitzek  Department of Microtechnology and Nanoscience, Chalmers
University of Technology, 412 96 Gothenburg, Sweden  Volvo Group Trucks
Technology, 405 08 Gothenburg, Sweden  Yi Hong Teoh  Department of Physics and
Astronomy, University of Waterloo, 200 University Ave. West, Waterloo, Ontario
N2L 3G1, Canada  Hin Pok Fung  Department of Physics and Astronomy, University
of Waterloo, 200 University Ave. West, Waterloo, Ontario N2L 3G1, Canada
Gebremedhin A. Dagnew  Department of Physics and Astronomy, University of
Waterloo, 200 University Ave. West, Waterloo, Ontario N2L 3G1, Canada  Ejaaz
Merali  Department of Physics and Astronomy, University of Waterloo, 200
University Ave. West, Waterloo, Ontario N2L 3G1, Canada  M. Schuyler Moss
Department of Physics and Astronomy, University of Waterloo, 200 University
Ave. West, Waterloo, Ontario N2L 3G1, Canada  Benjamin MacLellan  Department
of Physics and Astronomy, University of Waterloo, 200 University Ave. West,
Waterloo, Ontario N2L 3G1, Canada  Roger G. Melko  Department of Physics and
Astronomy, University of Waterloo, 200 University Ave. West, Waterloo, Ontario
N2L 3G1, Canada  Perimeter Institute for Theoretical Physics, Waterloo,
Ontario N2L 2Y5, Canada

(August 28, 2024)

###### Abstract

We introduce a generative pretained transformer (GPT) designed to learn the
measurement outcomes of a neutral atom array quantum computer. Based on a
vanilla transformer, our encoder-decoder architecture takes as input the
interacting Hamiltonian, and outputs an autoregressive sequence of qubit
measurement probabilities. Its performance is studied in the vicinity of a
quantum phase transition in Rydberg atoms in a square lattice array. We
explore the ability of the architecture to generalize, by producing
groundstate measurements for Hamiltonian parameters not seen in the training
set. We focus on examples of physical observables obtained from inference on
three different models, trained in fixed compute time on a single NVIDIA A100
GPU. These can act as benchmarks for the scaling of larger RydbergGPT models
in the future. Finally, we provide RydbergGPT open source, to aid in the
development of foundation models based off of a wide variety of quantum
computer interactions and data sets in the future.

##  I Introduction

Generative models have emerged as a key technology for predicting the
probabilistic behavior of a quantum system. Their most basic task is to
produce a target sequence representing qubit measurement outcomes; e.g.
projective measurements distributed according to the Born rule of quantum
mechanics. A wide variety of architectures have been explored for this
purpose, from restricted Boltzmann machines [1, 2, 3, 4, 5, 6, 7, 8, 9], to
recurrent neural networks [10, 11, 12, 5, 13, 14, 15, 16] and transformers
[17, 18, 19, 20, 21, 22, 23, 24]. While the above models have been trained to
output target sequences drawn from a single parameterized distribution, the
utility of generative models can be expanded by adopting modern encoder-
decoder architectures. In this case, the encoder takes a source sequence and
maps it to a context vector. This is then passed to the decoder, which
combines it with other inputs and maps it to the target sequence.

The most powerful architecture to emerge recently is the attention-based
transformer [25], variations of which underlie most of the large language
models (LLMs) being developed by academia and industry today [26, 27, 28, 29,
30, 31]. A key challenge in physics that can be tackled by the vanilla
encoder-decoder transformer structure is the task of predicting the output of
a quantum computer given its experimental settings. In the simplest example,
the source sequence is a description of an interacting Hamiltonian, and the
target sequence is the conditional probability of each qubit measurement
outcome, selected from a dictionary of size two. Viewing the task as a
language problem, the self-attention mechanism encodes context, or
correlations, between qubit measurement outcomes. Self-attention has been
shown to be a suitable mechanism of encoding the possible long-range
correlations – including entanglement – found in quantum systems without the
need to pass hidden or latent vectors along the sequence [19, 17, 20]. The
impressive demonstrations of scaling in large language models gives hope that
the transformer architecture will benefit quantum computers as they continue
to grow in size and capability [32].

In this work, we train an attention-based transformer to learn the
distribution of qubit measurement outcomes corresponding to a Hamiltonian
describing an array of interacting Rydberg atoms. Using the encoder-decoder
architecture, we use simulated Rydberg occupation data drawn from a two-
dimensional qubit array to train three different models. Once trained, each
model can produce new data via autoregressive sampling of the decoder. In the
case of a Rydberg atom array in its groundstate, the autoregressive output can
be interpreted as the complete quantum wavefunction, allowing for the
calculation of many physical observables of interest. We report results for
estimators of the energy and other physical observables. By varying the
Hamiltonian that is input to the encoder, we find that the transformer is able
to accurately predict qubit measurement outcomes in regions outside of the
training regime. Our results give optimism that attention-based transformer
models are well-positioned to provide a powerful approach for predicting the
properties of quantum many-body systems.

##  II Rydberg atom array physics

Rydberg atoms arrays are a powerful platform for quantum information
processing where interacting qubits are encoded with the electronic ground
state $|g\rangle$ and excited (Rydberg) state $|r\rangle$ for each atom [33,
34, 35, 36, 37]. We consider a system of $N=L\times L$ atoms arranged on a
square lattice. The governing Hamiltonian defining the Rydberg atom array
interactions has the following form:

$\displaystyle\hat{H}$
$\displaystyle=\sum_{i<j}\frac{C_{6}}{\lVert\mathbf{r}_{i}-\mathbf{r}_{j}\rVert^{6}}\hat{n}_{i}\hat{n}_{j}-\delta\sum_{i=1}^{N}\hat{n}_{i}-\frac{\Omega}{2}\sum_{i=1}^{N}\hat{\sigma}^{x}_{i},$
(1) $\displaystyle C_{6}$
$\displaystyle=\Omega\left(\frac{R_{b}}{a}\right)^{6},\quad
V_{ij}=\frac{a^{6}}{\lVert\mathbf{r}_{i}-\mathbf{r}_{j}\rVert^{6}},$ (2)

where $\hat{\sigma}^{x}_{i}=|g\rangle_{i}\langle r|_{i}+|r\rangle_{i}\langle
g|_{i}$, the occupation number operator
$\hat{n}_{i}=\frac{1}{2}\left(\hat{\sigma}_{i}+\mathbb{1}\right)=|r\rangle_{i}\langle
r|_{i}$ and $\hat{\sigma}_{i}=|r\rangle_{i}\langle r|_{i}-|g\rangle_{i}\langle
g|_{i}$. The experimental settings of a Rydberg atom array are controlled by
the detuning from resonance $\delta$, Rabi frequency $\Omega$, lattice length
scale $a$ and the positions of the atoms $\\{\mathbf{r}_{i}\\}_{i}^{N}$. From
Eq. (2), we obtain a symmetric matrix $\mathbf{V}$, that encapsulates the
relevant information about the lattice geometry, and derive the Rydberg
blockade radius $R_{b}$, within which simultaneous excitations are penalized
[33, 34]. Finally, for the purposes of our study, the atom array is considered
to be affected by thermal noise, in equilibrium at a temperature $T=1/\beta$.
The experimental settings are thus captured by the set of parameters
$\mathbf{x}=(\Omega,\delta/\Omega,R_{b}/a,\mathbf{V},\beta\Omega)$. By
adjusting the free parameters in the Rydberg Hamiltonian, the system can be
prepared in various phases of matter, separated by distinct phase transitions
[37, 38, 39, 40, 41, 42].

Rydberg atom array experiments are able to perform projective measurements of
the qubit states in the occupation basis. These measurements are
informationally complete at $T=0$ due to the positive real-valued ground state
wavefunction [43, 37, 44, 45]. In the next section, we consider these
projective measurements as training vectors for an unsupervised learning
strategy for our generative model. The goal of the generative model after
training will be to predict, given a Hamiltonian, the distribution of
projective measurements in the quantum system.

##  III Transformer architecture

Figure 1: Overview of RydbergGPT. The bottom section represents the data from
the physical system used in the architecture. On the left is a diagram
representing the experimental settings $\mathbf{x}$ of the system and on the
right is a corresponding projective measurement $\bm{\sigma}$. The
experimental settings are inputs to the encoder portion of the model, composed
of a graph neural network and the transformer encoder, and is processed to
form a context vector. This context vector is combined with the sample to
obtain a chain of conditional probabilities from the transformer decoder, the
product of which is the probability of obtaining the sample from the system.
The sequence of conditional probabilities is mapped onto the 2D lattice with a
“snake” path as illustrated.

In order to learn the qubit measurement distribution of the Rydberg atom
array, we employ the transformer encoder-decoder architecture [25]. These
generative models are universal sequence-to-sequence function approximators
and, granted sufficient tunable parameters $\theta$, can encode arbitrarily
complex functions [46]. Transformers belong to the class of generative models
known as autoregressive models. Autoregressive models learn conditional
probability distributions over the $N$ individual variables, such as qubit
states $\bm{\sigma}=\\{\sigma_{i}\\}_{i}^{N}$,
$p_{\theta}(\sigma_{i}|\sigma_{i-1},\ldots,\sigma_{1})$. The joint
distribution that arises from these conditionals can be used to represent the
likelihood of measurement outcomes,

$\displaystyle
p_{\theta}(\bm{\sigma})=\prod_{i=1}^{n}p_{\theta}\left(\sigma_{i}\mid\sigma_{i-1},\ldots,\sigma_{1}\right).$
(3)

This is the decoder output, which can be sampled to infer any number of qubit
measurements. This output can be related to a pure quantum wavefunction for
example through the Born rule
$p_{\theta}(\bm{\sigma})=|\Psi_{\theta}(\bm{\sigma})|^{2}$. Additionally,
under the assumption of a positive real-valued wavefunction, the transformer’s
output captures the full quantum state, i.e.
$\Psi_{\theta}(\bm{\sigma})=\sqrt{p_{\theta}(\bm{\sigma})}$ [10, 33], where
$\theta$ denotes the model parameters.

An encoder that maps the experimental Hamiltonian settings to a context vector
can be used to condition this decoder. The encoder-decoder architecture that
we use is illustrated in Fig. 1 and follows the overall design introduced by
Vaswani et al. [25].

Note that the source sequence for our encoder, the experimental Hamiltonian
settings $\mathbf{x}$, have a natural graph representation. The interactions
$\mathbf{V}$ can be considered as the edge attributes and the remaining
variables node attributes. The field of machine learning has proposed a
multitude of models for processing such graphs [47]. In our architecture, we
choose to process the graph with a convolutional graph neural network [47,
48]. This allow us to map the graph, which scales quadratically in the number
of atoms $N$, to a sequence representation that scales linearly. To further
capture correlations, we feed the sequence to the vanilla transformer encoder,
consisting of multiple self multi-head attention and feedforward blocks, and
obtain the context vector. The decoder autoregressively parameterizes the
equilibrium state of the system with respect to the context via an attention
mechanism between the context and the intermediate state within the decoder.
The choice of hyperparameters for the architecture is outlined in Table 1.

Table 1: The transformer encoder-decoder architecture and training parameters. Including the dimension of the feedforward network layer ($d_{\mathrm{ff}}$), the model dimension ($d_{\mathrm{model}}$), the graph dimension ($d_{\mathrm{graph}}$) as well as the total number of trainable parameters. Parameter | Value
---|---
Neural network architecture
$d_{\mathrm{ff}}$ | 128
$d_{\mathrm{model}}$ | 32
$d_{\mathrm{graph}}$ | 64
num heads | 8
num blocks encoder | 1
num blocks decoder | 3
num graph layers | 2
trainable params | 66562
Training hyperparameters
batch size | 1024
optimizer | AdamW
dropout | 0.1
learning rate | 0.001
learning rate schedule | Cosine annealing warm start
$T_{0}$ | 1
$T_{\mathrm{mult}}$ | 2
$\eta_{\mathrm{min}}$ | 0.00001
dataset buffer | 50

A powerful property of autoregressive models is the ability to perform
independent and identically distributed sampling without the use of Markov
chain Monte Carlo (MCMC). To perform sampling, we input the experimental
settings of interest to our encoder and obtain the context vector as output.
Following that, the decoder, conditioned on this context, takes as input a
starting token and samples the first element of the sequence. Then, we
iteratively sample the $i^{\mathrm{th}}$ token conditioned on all tokens
preceding it. We traverse the lattice of atoms systematically with the snake
path convention, shown in the bottom right of Fig. 1.

##  IV Training transformer models

We follow a data-driven training procedure, where we train transformers on a
synthetic dataset $\mathcal{D}$ composed of pairs of source and target
sequences. The former corresponds to experimental settings $\mathbf{x}$, input
to the encoder, and the latter corresponds to binary qubit measurement data,
input to the decoder. Our target sequences are simulated Rydberg occupation
measurements, produced via Stochastic Series Expansion quantum Monte Carlo
(QMC) method described in Ref. [44].

We employ our QMC to produce uncorrelated samples, which we collect across a
wide range of settings. Specifically, we choose values of the inverse
temperature $\beta\Omega$ within the set $\\{0.5,1,2,4,8,16\\}$, the detuning
parameter $\delta/\Omega$ over the range $[-0.364,3.173]$, the Rydberg
blockade $R_{b}/a$ within $\\{1.05,1.15,1.3\\}$, and the linear system size
$L$ for the values $\\{5,6,11,12,15,16\\}$. Without loss of generality, we
take $\Omega=1$ henceforth. For each configuration, we generate a dataset
containing $\lvert\mathcal{D}\rvert=10^{5}$ entries, each corresponding to a
binary measurement sampled from the QMC simulation. The datasets are selected
to explore $\delta/\Omega$ values near a quantum phase transition in the
square-lattice Rydberg atom arrays, which occurs at approximately
$\delta/\Omega=1.1$ [44].

A transformer architecture with parameters $\theta$ is trained in an
unsupervised manner with a dataset $\mathcal{D}$ by minimizing the Kullback-
Leibler (KL) divergence, resulting in the loss function:

$\mathcal{L}(\theta)\approx-\frac{1}{|\mathcal{D}|}\sum_{\bm{\sigma}\in\mathcal{D}}\ln
p_{\theta}(\bm{\sigma}).$ (4)

This loss is minimized when the decoder output and the target distributions
are equal. We initialize each model using the method introduced by Glorot et
al. [49] We use the AdamW optimizer [50], which is an extension to stochastic
gradient descent with an adaptive learning rate and weight regularization, in
combination with a cosine annealing warm restart learning rate scheduler [51].
Details for the training hyperparameters can be found in Table 1.

In this manuscript, we train three models with different subsets of the
simulated data. Specifically, the models are all trained on data at
$R_{b}/a=1.15$ and across the values for $\delta/\Omega$ and $\beta\Omega$.
However, each model is trained with different sets of linear size $L$. Model
$M_{1}$ is trained with data for systems of size $L=5,6$, model $M_{2}$ with
$L=5,6,11,12$ and model $M_{3}$ with $L=5,6,11,12,15,16$, refer to Table 2 for
the total number of tokens, i.e. measurement outcome of a single qubit, for
each lattice size. As we generate $10^{5}$ samples for each lattice size, the
larger lattice sizes have more tokens compared to the smaller ones.

Table 2:  Number of tokens, in units of $10^{8}$, contained in our dataset for various sizes $L$. Each token corresponds to the measurement outcome of a single Rydberg qubit. For each configuration, i.e. size and experimental settings, we have a total of $10^{5}$ samples, with each sample containing a total of $N=L\times L$ measurements.  Lattice size, $L$ | $5$ | $6$ | $11$ | $12$ | $15$ | $16$
---|---|---|---|---|---|---
Tokens ($10^{8}$) | $1.5$ | $2.16$ | $7.26$ | $8.64$ | $13.5$ | $15.4$
Table 3:  Details on the training of each model, including epochs trained, the amount of tokens in the training dataset, in units of $10^{8}$, and the total number of tokens processed during training, in units of $10^{10}$. The models were each trained for a constant amount of time, 85 hours, on an Nvidia A100.  Model | $M_{1}$ | $M_{2}$ | $M_{3}$
---|---|---|---
Epochs | 116 | 40 | 17
Tokens in dataset ($10^{8}$) | 3.66 | 19.6 | 48.4
Tokens processed ($10^{10}$) | 4.25 | 7.82 | 8.23

We conducted training on a single A100 GPU for 85 hours. The training datasets
for all three models are of the order of $10^{8}$, while during training, all
models processed tokens on the order of $10^{10}$. We summarize the relevant
training information in Table 3.

##  V Inference with trained models

In this section, we perform inference using the trained model and study the
disordered-to-checkerboard phase transition of the Rydberg atom array,
governed by the dimensionless detuning $\delta/\Omega$. First, samples are
drawn from the model in accordance to the method described at the end of
Section III. The samples are then used to obtain predictions of physical
observables. We compare our model’s predictions with estimates obtained from
QMC simulations [44]. For Rydberg atom arrays, QMC is able to obtain an
estimate of the observables that can be considered exact, within statistical
errors. In showing that our model’s predictions of these observables are in
good agreement with the QMC values, we illustrate the efficacy of data-driven
training as a technique for predicting Rydberg atom array measurement outcomes
for ground states. Further, we explore our trained models’ ability to
generalize away from the ground state and also to system sizes not represented
in the training data.

In order to construct our estimators below, we draw samples from our model
decoder. In our implementation, the complexity of sampling with respect to the
system size $N=L\times L$ is of order $\mathcal{O}(N^{3})$. However, this
complexity can theoretically be reduced to $\mathcal{O}(N^{2})$ by caching the
appropriate intermediate values in the attention calculation.

###  V.1 Physical observables

Physical quantities that are diagonal in the occupation basis can be directly
computed from the samples drawn from the model. One such quantity is the
staggered magnetization for the square-lattice Rydberg atom array. The
staggered magnetization is significant as it is the order parameter for the
disorder-to-checkerboard quantum phase transition [52], and it can be
calculated with:

$\displaystyle\langle\hat{\sigma}^{stag}\rangle\approx\frac{1}{N_{s}}\sum_{\bm{\sigma}\sim
p_{\theta}(\bm{\sigma})}\frac{1}{N}\left|\sum_{i=1}^{N}(-1)^{i}\left(n_{i}(\bm{\sigma})-1/2\right)\right|,$
(5)

where $i$ runs over all $N=L\times L$ atoms using the snake path convention of
Fig. 1. Here, $n_{i}(\bm{\sigma})=\langle\bm{\sigma}|r_{i}\rangle\langle
r_{i}|\bm{\sigma}\rangle$ is the occupation number operator acting on atom $i$
in a given configuration $\bm{\sigma}$.

In contrast, we consider an off-diagonal observable, where we must make use of
the ground state wave function amplitudes of the inferred samples
$\Psi(\bm{\sigma})=\sqrt{p_{\theta}(\bm{\sigma})}$. As an example, we examine
the spatially averaged expectation value of $\hat{\sigma}^{x}$, which is
defined as

$\displaystyle\langle\hat{\sigma}^{x}\rangle$
$\displaystyle\approx\frac{1}{N_{s}}\sum_{\bm{\sigma}\sim
p_{\theta}(\bm{\sigma})}\frac{1}{N}\sum_{\bm{\sigma}^{\prime}\in\mathrm{SSF}(\bm{\sigma})}\frac{\Psi_{\theta}(\bm{\sigma}^{\prime})}{\Psi_{\theta}(\bm{\sigma})},$
(6)

where $\mathrm{SSF}(\bm{\sigma})$ is the set of configurations that are
connected to $\bm{\sigma}$ by a single spin flip (SSF).

Finally, we consider an estimate of the ground state energy $\langle
E\rangle$, which is defined as

$\displaystyle\langle E\rangle$
$\displaystyle\approx\frac{1}{N_{s}}\sum_{\bm{\sigma}\sim
p_{\theta}(\bm{\sigma})}\frac{\langle\bm{\sigma}|\hat{H}|\Psi_{\theta}\rangle}{\langle\bm{\sigma}|\Psi_{\theta}\rangle}.$
(7)

This observable can be efficiently evaluated for local non-diagonal operators
[3, 53, 54, 55]. More detailed derivations of Eq. (6) and Eq. (7) can be found
in Appendix A.

It is important to emphasize that these estimates of off-diagonal observables
are appropriate only when the output of the decoder can be interpreted as the
normalized wavefunction amplitude,
$\Psi_{\theta}(\bm{\sigma})=\sqrt{p_{\theta}(\bm{\sigma})}$. In particular,
this correspondence is valid when the system is a pure state and the
wavefunction can be assumed to be real and positive in the occupation basis,
as is the case for the Rydberg Hamiltonian in its $T=0$ ground state.
Importantly, the use of the decoder output as a wave function amplitude is no
longer valid for finite temperature states. While the above estimates can
still be computed, $\Psi_{\theta}(\bm{\sigma})$ is not able to capture the
physics of the thermal state, as is demonstrated in our results which follow.

Figure 2:  Model predictions of multiple observables across the disordered-to-
checkerboard phase transition, which is governed by the dimensionless
detuning, $\delta/\Omega$. The observables include, (a) energy, (b) site-
averaged $x$-magnetization and (c) the staggered magnetization. Here, blue
points correspond to observables estimated from projective measurements
sampled from the trained model (i.e., out-of-distribution), red points
indicate parameters contained within the training dataset (i.e., in-
distribution), and dotted black lines correspond to the exact observables.

###  V.2 Results

We now demonstrate the results of the trained models across parameter regimes
and observables. In Fig. 2 we illustrate the performance of model $M_{1}$ on
the three observables introduced in Section V.1 as a function of
$\delta/\Omega$, for $R_{b}/a=1.15,\beta\Omega=16.0$. The estimators are
calculated over $N_{s}=10^{5}$ samples. The predictions are benchmarked
against the QMC observables calculated from $10^{5}$ decorrelated samples,
which we take as the ground truth. At $\beta\Omega=16$, the temperature is
sufficiently low for the Rydberg atom array to be in its groundstate to a good
approximation. The model shows good agreement with the QMC estimates for the
energy and staggered magnetization. As for the $x$-magnetization, we observe a
slight deviation from the QMC estimates in the tails. We conclude that the
trained model is able to generalize well to dimensionless detunings,
$\delta/\Omega$, that lie outside of the training data set, when the system is
restricted to the groundstate.

Next, we investigate how well the model generalizes away from the groundstate.
In Fig. 3 we illustrate predictions made by model $M_{1}$ with the same
strategy, except temperature $T=1/\beta$ is varied instead of the detuning. We
see a strong disagreement between the exact observables and the model-obtained
estimates at high temperatures. In our architecture, the decoder represents a
probability distribution over configurations in the Rydberg occupation basis,
which can be mapped one-to-one onto a positive wavefunction of a pure state.
As mentioned, the $T=0$ ground state of the Rydberg Hamiltonian Eq. (1) is
known to take such a form, and thus the model performs well for the ground
state, e.g. low temperatures. However, as the temperature is increased, the
system becomes a well-mixed Gibbs state, which can no longer be accurately
represented with a positive wavefunction. As such, the output of the decoder
is incapable of accurately representing this state. The model is merely
approximating a probability distribution that results in identical occupation-
basis measurement likelihoods. This results in the incorrect estimate of off-
diagonal observables, which require the wave function interpretation of the
decoder output $\Psi_{\theta}(\bm{\sigma})=\sqrt{p_{\theta}(\bm{\sigma})}$, as
in Eq. (6). To appropriately capture such high temperature states, one might
employ more sophisticated methods such as a complex density matrix [56] or
minimally entangled typical thermal states (METTS) [57, 58].

This is particularly evident in Fig. 3 at the limit of large temperature.
Here, the system tends to the fully-mixed state where all states are equi-
probable. The positive wavefunction that shares the same computational basis
probability as this state is the equi-superposition of all states. This state
corresponds to the maximum eigenstate of the x-magnetization, leading to the
divergence in Fig. 3 between the exact QMC and the trained model estimators at
high temperatures, where $\langle\hat{\sigma}^{x}\rangle$ becomes
significantly over-estimated. This is further supported by the fact that the
high-temperature estimates of the staggered magnetization are somewhat
accurate. Again, this behavior is due to the fact that this diagonal estimator
relies only on the computational basis probabilities with which certain
configurations are sampled.

Figure 3:  Model predictions of multiple observables as the dimensionless
temperature $T/\Omega$ is increased. The observables include, (a) energy, (b)
site-averaged $x$-magnetization and (c) the staggered magnetization. Again,
blue points correspond to observables estimated from projective measurements
sampled from the trained model (i.e., out-of-distribution), red points
indicate parameters contained within the training dataset (i.e., in-
distribution), and dotted black lines correspond to the exact observables.

Finally, we assess the ability of our transformer to generalize to larger
system sizes beyond those included in the training set. For this we use all
three trained models, $M_{1}$, $M_{2}$ and $M_{3}$. Figure 4 presents the
energy evaluations from each of our three trained models, averaged over
$10^{4}$ samples. These models were trained on consecutively larger datasets
as previously described in Table 3. We first observe that the model $M_{1}$
trained on the smallest dataset, containing only data from lattice sizes
$L=5,6$, captures the $L=5$ energy accurately but struggles to generalize to
larger system sizes. Interestingly, the models trained on more extensive
datasets do not accurately capture the properties of their respective training
sets, with the energy estimates deviating significantly from expected values.

These results illustrate how physical estimators are affected by bias in the
training datasets. As shown in Table 2, we always accumulate $10^{5}$ samples
for each grid size. Therefore, the number of training tokens or measurements
varies for each lattice size; more precisely, the dataset is biased towards
larger lattice sizes. This explains why models $M_{2}$ and $M_{3}$ perform
worse than $M_{1}$ for the lattice sizes $L=5,6$, since the fraction of
measurements for smaller sizes is less represented in their training dataset.
However, they perform better for larger grid sizes, where a larger relative
number of training tokens are available.

Figure 4: Extrapolation over lattice size $L$. The ground-state energy, as
estimated by quantum Monte Carlo simulations and the models trained on three
different datasets, demonstrates how training data affect the ability to
extrapolate to system sizes not included in the training.

##  VI Conclusion

In this paper we have explored the training of attention-based transformer
models for the problem of mapping a quantum Hamiltonian to the qubit
measurement outcomes of the corresponding quantum state. We construct an
encoder-decoder architecture suitable for Rydberg atom arrays, where the
encoder input is a representation of the interacting Hamiltonian, and the
decoder output represents measurements of the Rydberg occupation.

Since transformers have been developed as autoregressive models for natural
language processing, they naturally represent probability distributions, such
as Born-rule-distributed measurement outcomes. In the case where a quantum
state can be assumed to be pure, real and positive, this probability
distribution has a one-to-one mapping to the wavefunction, giving the output
of the decoder a powerful interpretation. In the groundstate of Rydberg atom
arrays governed by Hamiltonian Eq. (1), this assumption is valid, and we use
it to demonstrate that a trained transformer model can accurately reproduce
physical estimators in this case.

We have further explored the behavior of the transformer model when the
quantum state is no longer real and positive, i.e. in a thermally mixed state
at higher temperature. We clearly see the breakdown in our ability to
interpret the output of the decoder as the corresponding quantum state, as the
physical estimators significantly deviate from their exact values for large
$T$. However, we propose that our decoder architecture could be modified to
represent a density matrix [56] or METTS [57, 58], which would be able to
represent noisy, thermal and dynamical states, extending the capabilities of
the model.

Finally, we attempt to demonstrate the scaling potential of the encoder-
decoder transformer. Namely, whether the transformer model, when trained on a
dataset composed of smaller many-body system, is able to make predictions for
larger system sizes. Our results shows signatures of generalization, however,
it also suggest that we have not attained the necessary scale in our
transformer to exploit its full potential. Hopefully, by training models with
more data on a diverse set of lattice sizes, systematic improvements in model
performance will be observed. Indeed, such systematic scaling is the hallmark
of large transformers used in industry, giving reason for optimism [32].

More generally, just like in modern LLMs, our transformer architecture should
allow for vastly improved performance with increases in the number of
trainable parameters $\theta$, larger data sets, and more computational
resources. In addition, since these models are being designed to represent
physical qubit systems, a number of architectural improvements should be
explored. For example, there may be more efficient ways to use tokens than to
represent a dictionary of only two individual qubit measurements. In
particular, transformers designed for natural language processing typically
manage dictionaries comprising tens to hundreds of thousands of tokens [29].
An alternative approach involves using patches of qubits as tokens, which can
reduce some computational demands [20]. However, this method introduces
constraints on the system geometries, potentially limiting the applicability
to arbitrary quantum array configurations. This balance between computational
efficiency and system flexibility remains a critical area for future research.

Notably, all three models discussed in this study were trained using only a
single NVIDIA A100 GPU over a duration of 85 hours. Looking forward, we can
use this as a benchmark as the size and capability of such transformer models
is increased towards the goal of predicting the behavior of larger quantum
systems. Inspired by the scaling observed in large language models within
industry today, we imagine that significantly larger transformers trained on
vastly more data may be able to make highly useful predictions, e.g. about the
behavior of a quantum computer in parameter regimes for which no data is
currently available.

Finally, while our models in this study learned from synthetic data, our
transformer architecture could instead be trained on qubit measurement data
obtained from real experimental devices. Such would be the beginning of a new
type of foundation model, trained on real quantum computer output to generate
predictions on a wide variety of unseen input settings. Like today’s LLMs, the
capabilities of such large quantum models could increase dramatically with
increased number of parameters, data set size, and compute, potentially
ushering in a new era of scaling in AI and quantum computer co-design.

The RydbergGPT code and model weights for the models discussed in this study
are available at [59].

###### Acknowledgements.

We acknowledge support from the Natural Sciences and Engineering Research
Council of Canada (NSERC) and the Perimeter Institute for Theoretical Physics.
This research was supported in part by grant NSF PHY-2309135 to the Kavli
Institute for Theoretical Physics (KITP). Computational support was provided
by the facilities of the Shared Hierarchical Academic Research Computing
Network (SHARCNET:www.sharcnet.ca) and Compute/Calcul Canada. Research at
Perimeter Institute is supported in part by the Government of Canada through
the Department of Innovation, Science and Economic Development Canada and by
the Province of Ontario through the Ministry of Economic Development, Job
Creation and Trade. D.F. acknowledges the Knut and Alice Wallenberg (KAW)
Foundation for funding through the Wallenberg Centre for Quantum Technology
(WACQT). The computations were enabled by resources provided by Chalmers
e-Commons at Chalmers University of Technology.

## References

  * Nomura  _et al._ [2017] Y. Nomura, A. S. Darmawan, Y. Yamaji, and M. Imada, Restricted boltzmann machine learning for solving strongly correlated quantum systems, Phys. Rev. B 96, 205152 (2017).
  * Gao and Duan [2017] X. Gao and L.-M. Duan, Efficient representation of quantum many-body states with deep neural networks, Nature Communications 8, 662 (2017).
  * Carleo and Troyer [2017] G. Carleo and M. Troyer, Solving the quantum many-body problem with artificial neural networks, Science 355, 602 (2017).
  * Torlai and Melko [2018] G. Torlai and R. G. Melko, Latent space purification via neural density operators, Phys. Rev. Lett. 120, 240503 (2018).
  * Carrasquilla and Torlai [2021] J. Carrasquilla and G. Torlai, How to use neural networks to investigate quantum many-body physics, PRX Quantum 2, 040201 (2021).
  * Viteritti  _et al._ [2022] L. L. Viteritti, F. Ferrari, and F. Becca, Accuracy of restricted Boltzmann machines for the one-dimensional $J_{1}-J_{2}$ Heisenberg model, SciPost Phys. 12, 166 (2022).
  * Torlai  _et al._ [2019] G. Torlai, B. Timar, E. P. L. van Nieuwenburg, H. Levine, A. Omran, A. Keesling, H. Bernien, M. Greiner, V. Vuletić, M. D. Lukin, R. G. Melko, and M. Endres, Integrating neural networks with a quantum simulator for state reconstruction, Phys. Rev. Lett. 123, 230504 (2019).
  * Carrasquilla  _et al._ [2019] J. Carrasquilla, G. Torlai, R. G. Melko, and L. Aolita, Reconstructing quantum states with generative models, Nat. Mach. Intell. 1, 155 (2019).
  * Czischek  _et al._ [2018] S. Czischek, M. Gärttner, and T. Gasenzer, Quenches near ising quantum criticality as a challenge for artificial neural networks, Phys. Rev. B 98, 024311 (2018).
  * Moss  _et al._ [2024] M. S. Moss, S. Ebadi, T. T. Wang, G. Semeghini, A. Bohrdt, M. D. Lukin, and R. G. Melko, Enhancing variational monte carlo simulations using a programmable quantum simulator, Phys. Rev. A 109, 032410 (2024).
  * Czischek  _et al._ [2022] S. Czischek, M. S. Moss, M. Radzihovsky, E. Merali, and R. G. Melko, Data-Enhanced Variational Monte Carlo Simulations for Rydberg Atom Arrays, Physical Review B 105, 205108 (2022), arxiv:2203.04988 .
  * Carrasquilla and Melko [2017] J. Carrasquilla and R. G. Melko, Machine learning phases of matter, Nature Physics 13, 431 (2017).
  * Morawetz  _et al._ [2021] S. Morawetz, I. J. S. De Vlugt, J. Carrasquilla, and R. G. Melko, U(1)-symmetric recurrent neural networks for quantum state reconstruction, Phys. Rev. A 104, 012401 (2021).
  * Hibat-Allah  _et al._ [2023] M. Hibat-Allah, R. G. Melko, and J. Carrasquilla, Investigating topological order using recurrent neural networks, Phys. Rev. B 108, 075152 (2023).
  * Hibat-Allah  _et al._ [2021] M. Hibat-Allah, E. M. Inack, R. Wiersema, R. G. Melko, and J. Carrasquilla, Variational neural annealing, Nature Machine Intelligence 3, 952–961 (2021).
  * Hibat-Allah  _et al._ [2024] M. Hibat-Allah, R. G. Melko, and J. Carrasquilla, Supplementing recurrent neural network wave functions with symmetry and annealing to improve accuracy (2024), arXiv:2207.14314 [cond-mat.dis-nn] .
  * Wang  _et al._ [2022] H. Wang, M. Weber, J. Izaac, and C. Y.-Y. Lin, Predicting Properties of Quantum Systems with Conditional Generative Models (2022), arxiv:2211.16943 [quant-ph] .
  * Carrasquilla  _et al._ [2021] J. Carrasquilla, D. Luo, F. Pérez, A. Milsted, B. K. Clark, M. Volkovs, and L. Aolita, Probabilistic simulation of quantum circuits using a deep-learning architecture, Physical Review A 104, 032610 (2021).
  * Zhang and Di Ventra [2023] Y.-H. Zhang and M. Di Ventra, Transformer quantum state: A multipurpose model for quantum many-body problems, Physical Review B 107 (2023).
  * Sprague and Czischek [2024] K. Sprague and S. Czischek, Variational Monte Carlo with large patched transformers, Communications Physics 7, 90 (2024).
  * Cha  _et al._ [2022] P. Cha, P. Ginsparg, F. Wu, J. Carrasquilla, P. L. McMahon, and E.-A. Kim, Attention-based Quantum Tomography, Machine Learning: Science and Technology 3, 01LT01 (2022).
  * Viteritti  _et al._ [2023] L. L. Viteritti, R. Rende, and F. Becca, Transformer variational wave functions for frustrated quantum spin systems, Phys. Rev. Lett. 130, 236401 (2023).
  * Viteritti  _et al._ [2024] L. L. Viteritti, R. Rende, A. Parola, S. Goldt, and F. Becca, Transformer wave function for the shastry-sutherland model: emergence of a spin-liquid phase (2024), arXiv:2311.16889 [cond-mat.str-el] .
  * Lange  _et al._ [2024] H. Lange, G. Bornet, G. Emperauger, C. Chen, T. Lahaye, S. Kienle, A. Browaeys, and A. Bohrdt, Transformer neural networks and quantum simulators: a hybrid approach for simulating strongly correlated systems, (2024).
  * Vaswani  _et al._ [2017] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, Attention is all you need, in  _Advances in Neural Information Processing Systems_ , Vol. 2017-Decem (2017) pp. 5999–6009, arxiv:1706.03762 .
  * Zhang  _et al._ [2023] C. Zhang, C. Zhang, S. Zheng, Y. Qiao, C. Li, M. Zhang, S. K. Dam, C. M. Thwal, Y. L. Tun, L. L. Huy, D. kim, S.-H. Bae, L.-H. Lee, Y. Yang, H. T. Shen, I. S. Kweon, and C. S. Hong, A complete survey on generative ai (aigc): Is chatgpt from gpt-4 to gpt-5 all you need? (2023), arXiv:2303.11717 [cs.AI] .
  * Gozalo-Brizuela and Garrido-Merchan [2023] R. Gozalo-Brizuela and E. C. Garrido-Merchan, Chatgpt is not all you need. a state of the art review of large generative ai models (2023), arXiv:2301.04655 [cs.LG] .
  * Wu  _et al._ [2023] S. Wu, O. Irsoy, S. Lu, V. Dabravolski, M. Dredze, S. Gehrmann, P. Kambadur, D. Rosenberg, and G. Mann, BloombergGPT: A Large Language Model for Finance (2023), arxiv:2303.17564 [cs, q-fin] .
  * Brown  _et al._ [2020] T. Brown, B. Mann, N. Ryder, M. Subbiah, J. D. Kaplan, P. Dhariwal, A. Neelakantan, P. Shyam, G. Sastry, A. Askell, S. Agarwal, A. Herbert-Voss, G. Krueger, T. Henighan, R. Child, A. Ramesh, D. Ziegler, J. Wu, C. Winter, C. Hesse, M. Chen, E. Sigler, M. Litwin, S. Gray, B. Chess, J. Clark, C. Berner, S. McCandlish, A. Radford, I. Sutskever, and D. Amodei, Language models are few-shot learners, in  _Advances in Neural Information Processing Systems_ , Vol. 33, edited by H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin (Curran Associates, Inc., 2020) pp. 1877–1901.
  * Touvron  _et al._ [2023] H. Touvron, T. Lavril, G. Izacard, X. Martinet, M.-A. Lachaux, T. Lacroix, B. Rozière, N. Goyal, E. Hambro, F. Azhar, A. Rodriguez, A. Joulin, E. Grave, and G. Lample, LLaMA: Open and Efficient Foundation Language Models (2023), arxiv:2302.13971 [cs] .
  * Nakaji  _et al._ [2024] K. Nakaji, L. B. Kristensen, J. A. Campos-Gonzalez-Angulo, M. G. Vakili, H. Huang, M. Bagherimehrab, C. Gorgulla, F. Wong, A. McCaskey, J.-S. Kim, T. Nguyen, P. Rao, and A. Aspuru-Guzik, The generative quantum eigensolver (gqe) and its application for ground state search (2024), arXiv:2401.09253 [quant-ph] .
  * Kaplan  _et al._ [2020] J. Kaplan, S. McCandlish, T. Henighan, T. B. Brown, B. Chess, R. Child, S. Gray, A. Radford, J. Wu, and D. Amodei, Scaling laws for neural language models (2020), arXiv:2001.08361 [cs.LG] .
  * Wu  _et al._ [2021] X. Wu, X. Liang, Y. Tian, F. Yang, C. Chen, Y.-C. Liu, M. K. Tey, and L. You, A concise review of Rydberg atom based quantum computation and quantum simulation, Chinese Physics B 30, 020305 (2021).
  * Lukin  _et al._ [2001] M. D. Lukin, M. Fleischhauer, R. Cote, L. M. Duan, D. Jaksch, J. I. Cirac, and P. Zoller, Dipole Blockade and Quantum Information Processing in Mesoscopic Atomic Ensembles, Physical Review Letters 87, 037901 (2001).
  * Browaeys and Lahaye [2020] A. Browaeys and T. Lahaye, Many-body physics with individually controlled Rydberg atoms, Nature Physics 16, 132 (2020).
  * Jaksch  _et al._ [2000] D. Jaksch, J. I. Cirac, P. Zoller, S. L. Rolston, R. Côté, and M. D. Lukin, Fast Quantum Gates for Neutral Atoms, Physical Review Letters 85, 2208 (2000).
  * Endres  _et al._ [2016] M. Endres, H. Bernien, A. Keesling, H. Levine, E. R. Anschuetz, A. Krajenbrink, C. Senko, V. Vuletic, M. Greiner, and M. D. Lukin, Atom-by-atom assembly of defect-free one-dimensional cold atom arrays, Science 354, 1024 (2016).
  * Ebadi  _et al._ [2021] S. Ebadi, T. T. Wang, H. Levine, A. Keesling, G. Semeghini, A. Omran, D. Bluvstein, R. Samajdar, H. Pichler, W. W. Ho, S. Choi, S. Sachdev, M. Greiner, V. Vuletić, and M. D. Lukin, Quantum phases of matter on a 256-atom programmable quantum simulator, Nature 595, 227 (2021).
  * Xu  _et al._ [2021a] W. Xu, A. V. Venkatramani, S. H. Cantú, T. Šumarac, V. Klüsener, M. D. Lukin, and V. Vuletić, Fast Preparation and Detection of a Rydberg Qubit Using Atomic Ensembles, Physical Review Letters 127, 050501 (2021a).
  * Samajdar  _et al._ [2021] R. Samajdar, W. W. Ho, H. Pichler, M. D. Lukin, and S. Sachdev, Quantum phases of rydberg atoms on a kagome lattice, Proceedings of the National Academy of Sciences 118 (2021).
  * Miles  _et al._ [2023] C. Miles, R. Samajdar, S. Ebadi, T. T. Wang, H. Pichler, S. Sachdev, M. D. Lukin, M. Greiner, K. Q. Weinberger, and E.-A. Kim, Machine learning discovery of new phases in programmable quantum simulator snapshots, Phys. Rev. Res. 5, 013026 (2023).
  * Kalinowski  _et al._ [2022] M. Kalinowski, R. Samajdar, R. G. Melko, M. D. Lukin, S. Sachdev, and S. Choi, Bulk and boundary quantum phase transitions in a square rydberg atom array, Phys. Rev. B 105, 174417 (2022).
  * Bravyi  _et al._ [2007] S. Bravyi, D. P. DiVincenzo, R. I. Oliveira, and B. M. Terhal, The complexity of stoquastic local hamiltonian problems (2007), arXiv:quant-ph/0606140 [quant-ph] .
  * Merali  _et al._ [2023] E. Merali, I. J. S. De Vlugt, and R. G. Melko, Stochastic Series Expansion Quantum Monte Carlo for Rydberg Arrays (2023), arxiv:2107.00766 .
  * Xu  _et al._ [2021b] W. Xu, A. V. Venkatramani, S. H. Cantú, T. Šumarac, V. Klüsener, M. D. Lukin, and V. Vuletić, Fast preparation and detection of a rydberg qubit using atomic ensembles, Phys. Rev. Lett. 127, 050501 (2021b).
  * Yun  _et al._ [2020] C. Yun, S. Bhojanapalli, A. S. Rawat, S. J. Reddi, and S. Kumar, Are Transformers universal approximators of sequence-to-sequence functions? (2020), arxiv:1912.10077 [cs, stat] .
  * Kipf and Welling [2017] T. N. Kipf and M. Welling, Semi-supervised classification with graph convolutional networks, in  _5th International Conference on Learning Representations, ICLR 2017 - Conference Track Proceedings_ (2017) arxiv:1609.02907 .
  * Wang  _et al._ [2019] Y. Wang, Y. Sun, Z. Liu, S. E. Sarma, M. M. Bronstein, and J. M. Solomon, Dynamic graph cnn for learning on point clouds, ACM Trans. Graph. 38 (2019).
  * Glorot and Bengio [2010] X. Glorot and Y. Bengio, Understanding the difficulty of training deep feedforward neural networks, in  _Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics_, Proceedings of Machine Learning Research, Vol. 9, edited by Y. W. Teh and M. Titterington (PMLR, Chia Laguna Resort, Sardinia, Italy, 2010) pp. 249–256.
  * Loshchilov and Hutter [2019] I. Loshchilov and F. Hutter, Decoupled weight decay regularization, in  _International Conference on Learning Representations_ (2019).
  * Loshchilov and Hutter [2017] I. Loshchilov and F. Hutter, SGDR: Stochastic gradient descent with warm restarts, in  _International Conference on Learning Representations_ (2017).
  * Samajdar  _et al._ [2020] R. Samajdar, W. W. Ho, H. Pichler, M. D. Lukin, and S. Sachdev, Complex density wave orders and quantum phase transitions in a model of square-lattice rydberg atom arrays, Phys. Rev. Lett. 124, 103601 (2020).
  * Torlai  _et al._ [2018] G. Torlai, G. Mazzola, J. Carrasquilla, M. Troyer, R. Melko, and G. Carleo, Neural-network quantum state tomography, Nature Physics 14, 447 (2018).
  * Hibat-Allah  _et al._ [2020] M. Hibat-Allah, M. Ganahl, L. E. Hayward, R. G. Melko, and J. Carrasquilla, Recurrent neural network wave functions, Physical Review Research 2, 023358 (2020).
  * Melko  _et al._ [2019] R. G. Melko, G. Carleo, J. Carrasquilla, and J. I. Cirac, Restricted Boltzmann machines in quantum physics, Nature Physics 15, 887 (2019).
  * Kothe and Kirton [2023] S. Kothe and P. Kirton, Liouville space neural network representation of density matrices (2023), arXiv:2305.13992 [quant-ph] .
  * White [2009] S. R. White, Minimally entangled typical quantum states at finite temperature, Phys. Rev. Lett. 102, 190601 (2009).
  * Hendry  _et al._ [2022] D. Hendry, H. Chen, and A. Feiguin, Neural network representation for minimally entangled typical thermal states, Phys. Rev. B 106, 165111 (2022).
  * Fitzek and Teoh [2024] D. Fitzek and Y. H. Teoh, Rydberggpt (2024), code repository available at: https://github.com/PIQuIL/RydbergGPT/tree/main.

##  Appendix A Estimating expectation values of off-diagonal operators

In this work, we estimate the ground state expectation values of two off-
diagonal operators, $\langle\hat{\sigma}^{x}\rangle$ and $\langle E\rangle$,
using samples inferred from our trained GPT model. Here, we provide more
detail for how to arrive at equations Eq. (6) and Eq. (7).

Considering a general, off-diagonal operator $\hat{\mathcal{O}}$, we can write
down its expectation value with respect to the probability distribution
encoded in the trained GPT model,

$\langle\hat{\mathcal{O}}\rangle=\frac{\langle\Psi_{\theta}|\hat{\mathcal{O}}|\Psi_{\theta}\rangle}{\langle\Psi_{\theta}|\Psi_{\theta}\rangle},$

where $\Psi_{\theta}(\bm{\sigma})=\sqrt{p_{\theta}(\bm{\sigma})}$ is the wave
function interpretation of the joint probability distribution. Inserting the
identity in the numerator and the denominator, and multiplying by 1, the above
expectation value can be equivalently written as

$\langle\hat{\mathcal{O}}\rangle=\frac{\sum_{\bm{\sigma}}|\langle\Psi_{\theta}|\bm{\sigma}\rangle|^{2}\times\frac{\langle\bm{\sigma}|\hat{\mathcal{O}}|\Psi_{\theta}\rangle}{\langle\bm{\sigma}|\Psi_{\theta}\rangle}}{\sum_{\bm{\sigma}\prime}|\langle\Psi_{\theta}|\bm{\sigma}\prime\rangle|^{2}}.$

One can recognize
$\frac{|\langle\Psi_{\theta}|\bm{\sigma}\rangle|^{2}}{\sum_{\bm{\sigma}\prime}|\langle\Psi_{\theta}|\bm{\sigma}\prime\rangle|^{2}}$
as some normalized probability $P(\bm{\sigma})$ which can be approximated
using importance sampling. Here we also define the local estimator of our
operator $\hat{\mathcal{O}}$,

$\hat{\mathcal{O}}_{\text{loc}}(\bm{\sigma})=\frac{\langle\bm{\sigma}|\hat{\mathcal{O}}|\Psi_{\theta}\rangle}{\langle\bm{\sigma}|\Psi_{\theta}\rangle},$

such that the expectation value of $\hat{\mathcal{O}}$ can be estimated as

$\langle\hat{\mathcal{O}}\rangle\approx\frac{1}{N_{s}}\sum_{\bm{\sigma}\sim
p_{\theta}(\bm{\sigma})}\hat{\mathcal{O}}_{\text{loc}}(\bm{\sigma}).$

If we take $\hat{\mathcal{O}}=\hat{H}$, then we have exactly recovered
equation Eq. (7).

In order to arrive at equation Eq. (6), we will consider what happens when the
off-diagonal operator $\hat{\mathcal{O}}$ acts on the bra
$\langle\bm{\sigma}|$ in the definition of
$\hat{\mathcal{O}}_{\text{loc}}(\bm{\sigma})$.

If our off-diagonal operator $\hat{\mathcal{O}}$ is purely off-diagonal, it
will map a given configuration configuration $\bm{\sigma}$ to a linear
combination of some set of other configurations $\\{\bm{\sigma}\prime\\}$,
that is
$\hat{\mathcal{O}}|\bm{\sigma}\rangle=\sum_{\\{\bm{\sigma\prime}\\}}\text{C}_{\mathcal{O}}(\bm{\sigma\prime})|\bm{\sigma}\prime\rangle$,
and we can rewrite our expression for
$\hat{\mathcal{O}}_{\text{loc}}(\bm{\sigma})$ as,

$\hat{\mathcal{O}}_{\text{loc}}(\bm{\sigma})=\sum_{\bm{\sigma}\prime\in\mathrm{SSF}(\bm{\sigma})}\text{C}_{\mathcal{O}}(\bm{\sigma}\prime)\frac{\Psi_{\theta}(\bm{\sigma}\prime)}{\Psi_{\theta}(\bm{\sigma})},$

where we have used the fact that
$\Psi_{\theta}(\bm{\sigma})=\langle\bm{\sigma}|\Psi_{\theta}\rangle$.

We can note that if $\hat{\mathcal{O}}=\hat{\sigma}^{x}$, then
$\mathrm{SSF}(\bm{\sigma})$ is the set of configurations connected to
$\bm{\sigma}$ by a single spin flip (SSF) and
$\text{C}_{\mathcal{O}}(\bm{\sigma}\prime)=1$. Then it is clear that the inner
sum in equation Eq. (6) takes the form of the equation above.