# Removing Speaker Information from Speech Representation using Variable-
Length Soft Pooling

###### Abstract

Recently, there have been efforts to encode the linguistic information of
speech using a self-supervised framework for speech synthesis. However,
predicting representations from surrounding representations can inadvertently
entangle speaker information in the speech representation. This paper aims to
remove speaker information by exploiting the structured nature of speech,
composed of discrete units like phonemes with clear boundaries. A neural
network predicts these boundaries, enabling variable-length pooling for event-
based representation extraction instead of fixed-rate methods. The boundary
predictor outputs a probability for the boundary between 0 and 1, making
pooling soft. The model is trained to minimize the difference with the pooled
representation of the data augmented by time-stretch and pitch-shift. To
confirm that the learned representation includes contents information but is
independent of speaker information, the model was evaluated with libri-light’s
phonetic ABX task and SUPERB’s speaker identification task.

Index Terms—  speech representation, self-supervised learning, contrastive
learning, soft pooling

##  1 Introduction

Numerous self-supervised speech representation learning methods, including
Contrastive Predictive Coding(CPC) [1], HuBERT [2], and wav2vec2.0 [3] aim to
encompass general information in speech representations, such as voice
characteristics, semantic content, and emotional nuiances. While these
representations offer distinct advantages in enabling diverse downstream
speech-related applications, it remains crucial to extract the specific
content information within them. This need is particularly prominent when
sharing the embedding space with other modals that have unrelated voice
characteristics. Examples of this include matching speech representations from
text for text-to-speech tasks or reconstucting speech from signals from EMG
sensors attached to the face in the field of silent speech [4]. Moreover, the
disentanglement of speech representations can alleviate the load associated
with spoken language modeling [5] and speech synthesis [6] by reducing the
non-linguistic information that is noisy to other modules.

To disentangle non-linguistic information, we will employ a prior assumption
that speech is constructed from discrete linguistic units, which are
discernible through distinct boundaries. [7] and [8] represent instances of a
variable-rate downsampling technique, which enables the model to leverage the
observation that the informational content of speech is distributed non-
uniformly. To employ this prior more clearly, there are also methods that
quantize the representation by k-means clustering and merge identical
codes.([9, 5]) While these methods are capable of acquiring highly effective
non-uniform representations, it is important to note that in order to predict
adjacent representations or reconstruct the original audio, the representation
should encompass speaker-related information.

Fig. 1: The concept of soft pooling

Our objective is to minimize the presence of speaker-related information by
augmenting the data through temporal misalignment and subsequently tasking the
model with predicting the altered boundaries. This approach falls under the
domain of contrastive learning through augmentation. What sets it apart is
that the positive pair is not predefined; instead, it varies based on the
model’s predictions of the boundary. The overarching objective of this
research is to enhance phonetic information while concurrently diminishing the
presence of speaker-related details within speech representations.

##  2 Related Works

Creating positive samples by augmentation is a common method used in
contrastive learning [10, 11, 12]. The most augmentation methods for speech
include pitch-shift, additive noise, applying room reverberation, or voice
conversion. Time-stretch is hardly used because the position of a frame with
the same information in augmented data is changed. However, as time-stretch
does not change the total number or order of phonemes, time-stretched samples
can be used as a positive sample to embed linguistic information. [10] used
speed perturbation not for linguistic information. In singing domain, time-
stretch was applied [13], but augmented data was used as a negative sample.

This paper is related to ContentVec [11] in that it calculates the contrastive
loss by creating positive samples by augmentation while maintaining the
contents. [11] proposed three methods for disentanglement of speaker
information in HuBERT framework: for student, teacher, and predictor. [12] is
an application of a similar idea in the wav2vec2.0 framework. [12] predicts a
quantized latent representation as a context representation like wav2vec2.0,
and yet the prediction target is a representation extracted from augmented
data. The concept of contrastive learning through augmentation shares
similarities with our approach, but our method deliberately introduces
temporal misalignment along the time axis to predict boundaries.

Also, this paper has similarities to segmental CPC (SCPC) [14] and
hierarchical CPC (hCPC) [7], which predict the boundary and conduct pooling.
In hCPC, the boundary is predicted to be 0 or 1 by sampling the action from
the policy function. The policy function is learned to minimize the objective
function of the high-level CPC. The differentiable pooling module in this
paper is distinguishable in that it performs a weighted sum through attention
rather than average pooling based on the boundary. Meanwhile, [15] proposed a
framework to detect phoneme boundaries in an unsupervised manner. Predicting a
boundary with similarity between neighboring representations based on the CPC
structure inspired this study. Our approach shares a fundamental idea with
these methods, as we aim to predict the boundaries within the speech
representation sequence. However, our method distinguishes itself by supplying
positive pairs from augmented data rather than employing the same utterance.
This unique approach allows us to determine boundaries based on content,
thereby obviating the need for the inclusion of speaker-related information in
the representation.

In this paper, we propose a self-supervised speech representation learning
framework for optimizing the contrastive loss on augmented data, which is
pooled based on the boundary. Note that there is no supervision of the
boundary and no text label. To the best of our knowledge, there was no self-
supervised learning framework that adaptively pools time-stretched data and
compares it with the original data.

Fig. 2: Overall architecture of the model. $g_{enc}$, $g_{ar}$ are feature
extractor and autoregressive network for extracting context vectors,
respectively. Contrastive loss is calculated from the pooled representations
of the original data and the augmented data. The blue solid line, the blue
dotted line, and the red dotted line represent the anchor, positive sample,
and negative sample, respectively.

##  3 Method

The overall architecture of proposed model and the soft pooling module are
depicted in Figure 3. The model is divided into two streams: a Contrastive
Predictive Coding(CPC) part and a contrastive loss with augmentation part.
Contrastive loss induces the boundary predictor to be based on linguistic
information so that speech pairs with the same linguistic information are
aligned on the time axis. The CPC loss is responsible for maintaining the
linguistic information that can be lost while the contrastive loss is
optimized.

###  3.1 Contrastive Predictive Coding

The procedure for obtaining the d-dimensional spectral representation $z$ from
speech using convolution layers and the context vector $c$ through Long Short-
Term Memory(LSTM) is identical to that employed in the CPC [1] framework. The
encoders are trained by minimizing the contrastive loss$L_{CPC}$ calculated
from $z$ and $c$

$\displaystyle
L_{CPC}=-\sum_{n}\sum_{k=1}^{K}log\frac{f_{k}(c_{n},z_{n+k})}{\Sigma_{z_{i}\in\mathcal{N}}f_{k}(c_{n},z_{i})},$
(1)

where $f_{k}$ is prediction head for k upcoming representation. we opted to
utilize CPC as a baseline and complement it with the incorporation of a soft
pooling module.

###  3.2 Augmentation

In our augmentation process, we incorporated pitch-shift and time-stretch
techniques that preserve the speech’s underlying content. When applying time-
stretch, it should be noted that the positive samples are not time-aligned in
the fixed-rate latent representation. If variable-rate pooling reduces it to
an event-based representation, the representation sequence will be aligned
regardless of the length at which a particular linguistic unit is pronounced.
In essence, time-stretching prompts the boundary predictor to discern
linguistic information boundaries. To create a more challenging scenario for
boundary prediction, we divided each speech sample into three segments,
subjected them to random time-stretching, and subsequently concatenated them.
Additionally, random pitch-shift was introduced to prevent boundary prediction
based solely on pitch alterations.

Fig. 3: Description of soft pooling module.

###  3.3 Soft Pooling Module

We used attention pooling [16] to pool the representation sequence with
predicted boundary probability. The probability is bounded from 0 to 1 using
the sigmoid function as the activation function.

$\displaystyle
s_{m}=\sum_{n=1}^{N}\alpha_{mn}z_{n}=\sum_{n=1}^{N}\frac{K_{m}(a_{n})}{\Sigma_{i=1}^{N}K_{m}(a_{i})}z_{n}$
(2) $\displaystyle K_{m}(a_{n})=\exp{(-0.5\sigma^{-2}(a_{n}-m)^{2}}),$ (3)

where $N$ is the number of sequence length, $K_{m}$ is a Gaussian kernel for
$m$-th attention head. As $a_{n}$ is cumulative sum of boundary, attention
becomes monotonic with $m$. The kernel causes $a_{n}$ to have more attention
as $a_{n}$ is closer to $m$. Ideally, if the boundaries only have values of 0
or 1 and the $\sigma$ is sufficiently reduced, attention would be equivalent
to averaging pooling representations between non-zero boundaries. In other
words, the clearer the boundary, the closer $a_{n}$ to an integer, and then
reducing the $\sigma$ helps to create attention that reduces side effects.
Because the boundary is not ideally trained, sufficient tolerance is required
for the kernel, and in this paper, the sigma is set to 0.5.

The number of attention heads M to be set smaller than N is also an important
hyperparameter. $M$ needs to be set large enough to prevent the attention of
$s_{M}$ from ending before reaching the end of the $z$ sequence. Since the
average number of phonemes per sec is about 10 and the frame rate of $g_{enc}$
is 100, $M=N/10$ can also be considered, but $M=N/4$ was conservatively set.
If the data is composed of few linguistic units (in the extreme case, silent
audio), $s_{m}$ will be close to 0-vector for large m, which will not have
much impact.

It is possible to substitute $c$ for $z$ used to calculate attention
weight(key, query) or $z$ to be pooled(value). However, when $c$ is used as a
key and query, the performance is slightly degraded, and when $c$ is used as a
value, it collapses to predict the boundary as 1 for all frames. Thus, only
$z$ is used in soft pooling module.

###  3.4 Contrastive Loss

Soft pooling module is trained by following objective function;

$\displaystyle
L_{contr}=-\sum_{m=1}^{M}\log\frac{f(s_{m},s^{\prime}_{m})}{\Sigma_{j=1}^{M}f(s_{m},s_{j})}$
(4) $\displaystyle
f(s^{(1)},s^{(2)})=\exp(\mathrm{sim}(s^{(1)},s^{(2)})/\tau),$ (5)

where $\tau$ is temperature and sim is similarity function such as cosine
similarity. Positive sample is representations of the same location in the
augmented data, and negative samples are representations of different
locations in the original data like disentanglement in student of ContentVec
[11]. This loss induces the downsampled representations $s$, $s’$ to be
similar frame-wise than different frames within the same $s$, making them
event-based representations.

##  4 Experiment

###  4.1 Experiment Setup

Proposed models and baseline CPC were trained on LibriSpeech 960hours dataset
[17]. Every model was trained with 2 RTX 3080ti GPUs and batch size is set to
32. We used RAdam optimizer with learning rate of 0.001. In order to see
feasibility with limited resources, the representation is set to 128
dimensions, and accordingly, the number of the model parameters is 600K. The
training took about 60 hours for 600,000 steps.

As this paper introduces a framework that combines both the CPC loss and the
proposed contrastive loss, it is imperative to conduct a comparison with the
established CPC approach. For this purpose, we will consider CPC with a
reduced embedding dimension of 128 dimensions as our baseline. Additionally,
we will assess models trained solely with the proposed contrastive loss,
without the incorporation of CPC.

###  4.2 Evaluation Methods

Phonetic ABX task. As one of the black box and zero-shot metrics provided by
Zero-Resource Speech 2021 challenge[18], it evaluates whether the model
distinguishes phonemes of speech well. The evaluation data is a sample drawn
from libri-light and is given as a triphone. Given the speech samples of A, B,
and X, it determines whether X belongs to A or B based on the representation
sequence extracted by the model. A and B are speech of the same speaker, and
it is divided into within/across depending on whether X is the same speaker.
Through the phonetic ABX error rate, it is possible to evaluate the extent to
which the representation contains phonetic information.

Speaker Identification. It was used to evaluate speaker disentanglement in
ContentVec as a task to distinguish speakers in the SUPERB task. We used s3prl
toolkit for SID evaluation [19]. The lower the accuracy of speaker
classification, the less information about the speaker in the representation.

Phoneme Segmentation. TIMIT is often used as a benchmark for phoneme
segmentation because it provides labels for phonemes and their boundaries. The
precision, recall, and F1 scores are calculated by comparing the predicted
boundary with the ground truth. Through this, it is possible to check how well
the predicted boundary conforms to the phoneme boundary by learning with
contrastive learning without any labels. Since most of the speech
representation learning methods do not explicitly predict the boundary, only
the results of hCPC are reported for comparison.

Table 1: Phonetic ABX error rate and speaker identification accuracy(SID). Lower SID accuracy is desirable. Model | ABX $\downarrow$ | SID $\downarrow$
---|---|---
ContentVec [11] | 5.13 | 37.7
CPC baseline | 8.19 | 26.3
Proposed model | 7.61 | 17.7
Proposed model w/o CPC | 13.8 | 5.18
Table 2: Phoneme segmentation performance on TIMIT dataset. The metrics are
calculated with 20 ms tolerance

Model | Precision $\uparrow$ | Recall $\uparrow$ | F1 $\uparrow$ | R-val $\uparrow$
---|---|---|---|---
Proposed model | 73.31 | 75.01 | 74.15 | 77.79
Proposed model w/o CPC | 71.38 | 70.74 | 71.04 | 75.32

###  4.3 Experiment Results

The table 1 showed lower phonetic error rates and higher SID accuracy for the
proposed model compared to the baseline CPC model. This highlights the
effectiveness of the soft pooling module, optimized via contrastive loss, in
disentangling speaker-related information while enhancing linguistic content
embedding. However, when relying solely on contrastive loss (without CPC
loss), SID accuracy decreased, indicating the importance of minimizing
speaker-related data for optimizing the proposed contrastive loss. However, it
is worth noting that without CPC loss, some phonetic information was
sacrificed. In cases where both losses were jointly optimized, the performance
in phonetic tasks improved while concurrently achieving speaker information
disentanglement. Although the phonetic performance falls short of that
achieved by the HuBERT-based ContentVec model, a substantial reduction in
speaker-related information was accomplished.

The table 2 indicates that the boundaries predicted by the proposed model
align with phoneme boundaries. Remarkably, the segmentation F1 score reaches
74.15%, even in the absence of constraints explicitly dictating that the
predicted boundaries must coincide with phoneme boundaries. Similar to hCPC,
employing a prior that incorporates phoneme rates as a regularization
mechanism for boundary prediction has the potential to yield further
improvements.

The figure 4 displays the model’s predicted boundaries and attention patterns.
Although not all boundaries are precisely predicted, it’s notable that regions
with richer content engage more attention heads. For instance, the silent
segment (a) and the extended /ey/ phoneme section (c) primarily involve one or
two active attention heads, while the content-rich section (b) engages
multiple attention heads. Despite variations in content distribution over
time, both predicted boundaries and attention patterns share content-dependent
trends.

Fig. 4: The three figures above are the mel-spectrogram with phoneme
boundaries, predicted boundaries, and unnormalized attention weights of the
original data, respectively, and the three figure below shows the same
components for augmented data. Sample data is SI943 of FAKS0 from TIMIT.

##  5 Conclusion

Using the proposed soft pooling module has the effect of increasing phonetic
information while concurrently diminishing speaker information within the
speech representation for limited resources. In the context of self-supervised
speech representation learning, the boundary prediction has proven
advantageous in acquiring disentangled content-related information. An
intriguing observation arises from the fact that the boundaries identified as
optimal for the objective function closely align with the phoneme boundaries.
Our strategy entails enhancing the model’s capacity and implementing both
quantization techniques and a unit-language model.

## References

  * [1] Aaron van den Oord, Yazhe Li, and Oriol Vinyals,  “Representation learning with contrastive predictive coding,”  arXiv preprint arXiv:1807.03748, 2018.
  * [2] Wei-Ning Hsu, Yao-Hung Hubert Tsai, Benjamin Bolte, Ruslan Salakhutdinov, and Abdelrahman Mohamed,  “Hubert: How much can a bad teacher benefit asr pre-training?,”  in ICASSP 2021 - 2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2021, pp. 6533–6537.
  * [3] Alexei Baevski, Yuhao Zhou, Abdelrahman Mohamed, and Michael Auli,  “wav2vec 2.0: A framework for self-supervised learning of speech representations,”  in Advances in Neural Information Processing Systems, H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin, Eds. 2020, vol. 33, pp. 12449–12460, Curran Associates, Inc.
  * [4] Kevin Scheck and Tanja Schultz,  “Multi-speaker speech synthesis from electromyographic signals by soft speech unit prediction,”  in ICASSP 2023-2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2023, pp. 1–5.
  * [5] Kushal Lakhotia, Eugene Kharitonov, Wei-Ning Hsu, Yossi Adi, Adam Polyak, Benjamin Bolte, Tu-Anh Nguyen, Jade Copet, Alexei Baevski, Abdelrahman Mohamed, and Emmanuel Dupoux,  “On generative spoken language modeling from raw audio,”  Transactions of the Association for Computational Linguistics, vol. 9, pp. 1336–1354, 2021.
  * [6] Hyeong-Seok Choi, Juheon Lee, Wansoo Kim, Jie Lee, Hoon Heo, and Kyogu Lee,  “Neural analysis and synthesis: Reconstructing speech from self-supervised representations,”  Advances in Neural Information Processing Systems, vol. 34, pp. 16251–16265, 2021.
  * [7] Santiago Cuervo, Adrian Łańcucki, Ricard Marxer, Paweł Rychlikowski, and Jan Chorowski,  “Variable-rate hierarchical cpc leads to acoustic unit discovery in speech,”  arXiv preprint arXiv:2206.02211, 2022.
  * [8] Sander Dieleman, Charlie Nash, Jesse Engel, and Karen Simonyan,  “Variable-rate discrete representation learning,”  arXiv preprint arXiv:2103.06089, 2021.
  * [9] Tu Anh Nguyen, Maureen de Seyssel, Patricia Rozé, Morgane Rivière, Evgeny Kharitonov, Alexei Baevski, Ewan Dunbar, and Emmanuel Dupoux,  “The zero resource speech benchmark 2021: Metrics and baselines for unsupervised spoken language modeling,”  arXiv preprint arXiv:2011.11588, 2020.
  * [10] Dongwei Jiang, Wubo Li, Miao Cao, Wei Zou, and Xiangang Li,  “Speech simclr: Combining contrastive and reconstruction objective for self-supervised speech representation learning,”  arXiv preprint arXiv:2010.13991, 2020.
  * [11] Kaizhi Qian, Yang Zhang, Heting Gao, Junrui Ni, Cheng-I Lai, David Cox, Mark Hasegawa-Johnson, and Shiyu Chang,  “Contentvec: An improved self-supervised speech representation by disentangling speakers,”  in International Conference on Machine Learning. PMLR, 2022, pp. 18003–18017.
  * [12] Vasista Sai Lodagala, Sreyan Ghosh, and S Umesh,  “Ccc-wav2vec 2.0: Clustering aided cross contrastive self-supervised learning of speech representations,”  arXiv preprint arXiv:2210.02592, 2022.
  * [13] Hiromu Yakura, Kento Watanabe, and Masataka Goto,  “Self-supervised contrastive learning for singing voices,”  IEEE/ACM Transactions on Audio, Speech, and Language Processing, vol. 30, pp. 1614–1623, 2022.
  * [14] Saurabhchand Bhati, Jesús Villalba, Piotr Żelasko, Laureano Moro-Velazquez, and Najim Dehak,  “Segmental contrastive predictive coding for unsupervised word segmentation,”  arXiv preprint arXiv:2106.02170, 2021.
  * [15] Felix Kreuk, Joseph Keshet, and Yossi Adi,  “Self-supervised contrastive learning for unsupervised phoneme segmentation,”  arXiv preprint arXiv:2007.13465, 2020.
  * [16] Elizbar A Nadaraya,  “On estimating regression,”  Theory of Probability & Its Applications, vol. 9, no. 1, pp. 141–142, 1964.
  * [17] Vassil Panayotov, Guoguo Chen, Daniel Povey, and Sanjeev Khudanpur,  “Librispeech: An asr corpus based on public domain audio books,”  in 2015 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2015, pp. 5206–5210.
  * [18] Ewan Dunbar, Mathieu Bernard, Nicolas Hamilakis, Tu Anh Nguyen, Maureen de Seyssel, Patricia Rozé, Morgane Rivière, Eugene Kharitonov, and Emmanuel Dupoux,  “The Zero Resource Speech Challenge 2021: Spoken Language Modelling,”  in Proc. Interspeech 2021, 2021, pp. 1574–1578.
  * [19] Andy T. Liu, Shu-wen Yang, Po-Han Chi, Po-chun Hsu, and Hung-yi Lee,  “Mockingjay: Unsupervised speech representation learning with deep bidirectional transformer encoders,”  in ICASSP 2020 - 2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2020, pp. 6419–6423.


#  Statistical Mechanics of Thermostatically Controlled Multi-Zone Buildings

Lucas Fuentes Valenzuela  Department of Electrical Engineering, Stanford
University  Lindell Williams  Department of Physics, Cornell University
Program in Applied Mathematics & Department of Mathematics, University of
Arizona  Michael Chertkov  Program in Applied Mathematics & Department of
Mathematics, University of Arizona<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We study the collective phenomena and constraints associated with the
aggregation of individual cooling units from a statistical mechanics
perspective. These units are modelled as Thermostatically Controlled Loads
(TCLs) and represent zones in a large commercial or residential building.
Their energy input is centralized and controlled by a collective unit – the
Air Handling Unit (AHU) — delivering cool air to all TCLs, thereby coupling
them together. Aiming to identify representative qualitative features of the
AHU-to-TCL coupling, we build a realistic but also sufficiently simple model
and analyze it in two distinct regimes: the Constant Supply Temperature (CST)
and the Constant Power Input (CPI) regimes. In both cases, we center our
analysis on the relaxation dynamics of individual TCL temperatures to a
statistically steady state. We observe that while the dynamics are relatively
fast in the CST regime, resulting in all TCLs evolving around the control
setpoint, the CPI regime reveals emergence of a _bi-modal probability
distribution and two, possibly strongly separated, time scales_. We observe
that the two modes in the CPI regime are associated with all TCLs being in the
same low and high-temperature states, respectively, with occasional (and
therefore possibly rare) collective transition between the modes akin in the
Kramer’s phenomenon of statistical physics. To the best of our knowledge, this
phenomenon was overlooked in the context of the multi-zone energy building
engineering, even thought it has direct implications on the operations of
centralized cooling systems in buildings. It teaches us that a balance needs
to be struck between occupational comfort — related to variations in the
individual temperatures — and power output predictability — the main focus of
the DR schemes.

##  I Introduction

The essence of Demand Response (DR) lies in providing auxiliary services
helping power system operators manage uncertainty. The latter can emerge both
from variable generation, e.g. as wind and solar, or from electricity market
volatility [1]. The key idea behind DR consists in leveraging flexible and
inexpensive resources on the demand side of the power balance to ensure
stability. The essential source of flexibility lies in the fact that many
consumers of electricity can tolerate consumption delays, provided that some
constraints remain satisfied [2]. In addition to large and stable loads,
aggregations of many small loads, e.g. residential appliances, can also be
involved in DR services [3]. The heating and cooling system in residential
buildings is one such load that possesses inherent flexibility due to thermal
inertia, and thereby presents significant opportunities in a DR market.
Specifically, the potential for impacting total load via setpoint changes has
been investigated numerically [4, 5] and empirically [6, 7], as well as its
ramifications for occupant comfort [8].

###  I.1 Related Work

Theoretical studies on the matter have essentially anchored around the so-
called Thermostatically Controlled Loads (TCLs), denoting physical entities
whose temperature oscillates within a range or around a target value; examples
include rooms in building or refrigerators. Understanding the behavior of
aggregations of individual entities and the underlying potential for DR has
underpinned most of the interest in this discipline.

Initial studies focused on adapting existing and developing new approaches to
the statistics of TCLs [9], under different characteristics such as lifestyle
and weather [10], and proposed a methodology for the aggregation of individual
loads [11]. Large aggregates, also called ensembles, were studied with the
tools of statistical physics, such as Fokker-Planck equations [12, 13], and of
control and reinforcement learning, such as Markov Decision Processes [14,
15]. The main operational philosophies in the literature include individual
thermostat setpoints control [16, 17, 18] or randomization and automatic
feedback control at the individual TCL level, based on collective output [19,
20].

While standard TCL models do represent some installations like independent AC
units, or refrigerators, they do not directly capture the intricate dynamics
emerging from the coupling of those units via a district heating network [21],
or the grid services they provide [22, 23, 16]. In the context of heating and
cooling within a multi-zone building, individual zones are thermally regulated
by a small number of Air Handling Units (AHU), see Fig. 1 for a simplified
illustration. Each AHU is connected to some number of zones. Its role is to
cool and dehumidify a mix of outside and recirculated air to a given
temperature and relative humidity, and then to circulate this air — at a given
temperature, the supply temperature — throughout the building. The airflow is
distributed to each zone via Variable Air Volume (VAV) boxes. See [24] for an
accurate and detailed description of the operations of such systems.

Figure 1: Simplified illustration of the setup considered. The AHU is
represented by the fan, while individual units have their own temperatures.
The exhaust system from each unit is not represented for simplicity, and air
recirculation is neglected.

Standard TCL models fail to account for the specificities of such systems in
two ways. First, they usually assume units to be independent. However, they
are actually coupled through AHU – that is adjusting to the supply
temperature. This interaction and its implication for control has not been
studied extensively. Second, the dynamics of cooling in these models is very
different from the reality of residential and commercial buildings. Indeed, in
standard TCL models, cycling is inherent in the dynamics. A unit is cooled at
constant power until it crosses the lower limit of a deadband, and it is then
left to reheat until it crosses its upper limit, at which point the cycle
restarts. In AHU systems, such cycling is not enforced directly, and the
cooling power changes following changes in the supply temperature [25]. We
will see below that the abscence of a prescribed deadband has far-reaching
consequences on how fast the system can recover from a perturbation, e.g.
caused by the units following a DR request.

###  I.2 Contribution

In this manuscript, we suggest new models to represent TCLs coupled via an
AHU, and study the thermalization dynamics of the system subject to two-level
control – at the collective (AHU) level and at the individual (TCL) level.
Specifically, we consider two principally different models representing
distinct scenarios:

#### Constant Supply Temperature (CST)

The supply air temperature cooling each TCL is held constant over time,
thereby eliminating direct coupling between individual units. While
simplistic, this model enables the understanding of the fundamental dynamics
of TCL units.

#### Constant Power Input (CPI)

In order to reflect actual building operations, TCLs are coupled via the
supply temperature. The latter is adjusted in order to ensure constant power
input.

In both cases, we focus on analyzing the details of the stabilization of the
ensemble to a statistical steady state, established as a result of a balance
between exogenous fluctuations representing TCL-specific perturbations and
thermal relaxation. We combine analytical and numerical analyses to derive
insights into the dynamics of this system. In the CST model, we find that
thermalization is relatively fast, and results in TCLs stabilizing around the
prescribed temperature. However, in the CPI scenario, we observe the
appearance of more complex dynamics. When the level of stochastic fluctuations
is small, all units stabilize in either of the two states associated with low-
airflow and high-airflow. Transitions between both modes, and thus relaxation
to a statistically steady state, happens at a time scale which is much longer
that the natural time scale of the problem, associated with meta-stable
equilibration around one of the modes. Increasing the level of fluctuations
results in a phase transition to an "entropic" state distributed around the
control setpoint. We argue that this complex behavior of the CPI model is akin
to the phenomenon of the thermally-activated barrier crossing in statistical
mechanics, often referred to as the reaction-rate, or Kramers’ theory [26,
27].

These observations are directly relevant for the operations of those common
coupled systems. Under low uncertainty and perturbations, one can argue that
operating the system at fixed supply temperature (CST) will result in a
predictable outcome matching comfort requirements, with minor fluctuations in
the total power required. However, in the case of higher fluctuations, e.g.
associated with other uncertainties or driving forces such as sun irradiation,
the dynamics of such a system is non-trivial and implies a need to establish a
balance between comfort and requirements for utilization in DR.

###  I.3 Outline

The contents are presented as follows. In Section II, the problem is
formulated along with the two scenarios of interest. A statistical analysis of
the CST scenario is conducted in Section III. Analyses under Constant Power
Input are conducted in Section IV. We first consider the noiseless limit for 1
TCL, and extend its insights to the presence of noise. The general case with
an arbitrary number of TCLs is studied in Section IV.3. We conclude and
discuss the path forward in Section V.

##  II Problem Formulation

The system we consider consists of two levels. At the lower level, TCLs evolve
according to their individual dynamics. At the higher level, the Air Handling
Unit manages the energetic input to all TCLs, thereby coupling them to one
another (cf Fig. 1). In this Section, we lay out the models governing both
levels.

###  II.1 TCL dynamics

We model the dynamics of an individual TCL by the following Stochastic
Differential Equation (SDE):

$\displaystyle c\ dT_{i}=-f(T_{i}(t);s(T_{i}(t));T_{s})\ dt+\sqrt{2D}\
dW_{i}(t),$ (1)
$\displaystyle\begin{split}f(T_{i}(t);s(T_{i}(t));T_{s})&=\frac{T_{i}(t)-T_{o}}{r}\\\
&+\bar{\mu}s(T_{i}(t))c_{p}\left(T_{i}(t)-T_{s}\right),\\\ \end{split}$ (2)
$\displaystyle
s(T_{i}(t))=\left\\{\begin{array}[]{ll}s_{-},&T_{i}(t)\leq\bar{T}\\\
s_{+},&T_{i}(t)>\bar{T}\end{array}\right..$ (5)

where $t$ is time and $T_{i}(t)$ is the temperature of the $i$-th TCL unit
$i\in 1,...,N$. We assume that all units are identical, and characterized by
the following physical parameters: $r\ [K/kW]$ is the coefficient of thermal
resistance to the outside air; $c_{p}\ [kJ/kgK]$ is the specific heat capacity
of air; $c\ [kJ/K]$ is the capacitance coefficient representing the total air-
mass associated with an individual TCL; $W_{i}(t)$ is the Wiener process; and
$D\ [kJ^{2}]$ represents the amplitude of exogenous, stochastic, zero mean
thermal (white) noise affecting the TCLs (associated, e.g., with some
uncertainty around in-zone traffic and operations). $\bar{\mu}$ denotes the
maximum airflow that can get into the zone, while $s(\cdot)\in[0,1]$ indicates
how much air is actually flowing into the room at time $t$. In Eq. (2),
describing the temperature dynamics of a TCL unit,
$f(T_{i}(t);s(T_{i}(t));T_{s})$ denotes the thermal force associated with two
principally different terms. The first contribution to the thermal force
describes deterministic relaxation of the zone temperature, $T_{i}(t)$, with
the rate $r$ to the ambient temperature $T_{o}$, which is taken as constant in
this document. The second contribution to the thermal force represents the
injection of cool air delivered by VAV units at the mass flow rate
$\bar{\mu}s(T_{i}(t))$ and at a temperature equal to the supply temperature
$T_{s}$ into the $i$-th zone. In the present system, we assume the VAVs are
able to modify only the airflow entering each zone. While in general they are
also able to reheat the incoming air, we assume this is not applicable, as is
the case in warm regions where cooling is mostly needed. The dynamics of the
air flow at each VAV is a simplification of the dual-maximum control logic
[28]: we neglect the heating regime and we assume instantaneous switching as
the setpoint is crossed instead of progressive changes of the airflow (Eq.
(5)). The model proposed is a degenerate case of the more general one, with
the benefit of simplicity and tractability.

###  II.2 AHU dynamics

The model in Eqs (1)-(5) is incomplete until one provides a closure
relationship for $T_{s}$, controlling the second contribution to the thermal
force in Eq. (2). As detailed in the remainder of this Section, we consider
two scenarios which will be different in terms of the relation between $T_{s}$
and the vector of TCL temperatures, ${\bm{T}}=(T_{i}|i=1,\cdots,N)$. In this
paper, we consider simplified models of the AHU energy consumption. Indeed, we
assume that the units lie in areas where the energy attributed to cooling
dominates that associated to de-humidification. This means that the energy
consumption of the AHU can be modeled by the enthalpy change of the air as it
goes through the AHU [29, 30, 31, 32]. For the sake of simplicity, we also
neglect air recirculation within the building, and assume that only the
outside air is used as input into the AHU.

### (a) Constant Supply Temperature (CST)

The first scenario considered, indexed by $a$, consists in fixing the
temperature $T_{s}$ of the air delivered to the TCLs to a prescribed constant
value, $T_{s;a}$. The total power, i.e. energy injected into the ensemble per
unit time, then writes

$\displaystyle
P_{a}(t)=\mu(t)\left(T_{o}-T_{s;a}\right),\quad\mu(t)=\bar{\mu}\sum_{i}s(T_{i}(t)),$
(6)

where $\mu(t)$ is the aggregate airflow within the AHU. We assume that all the
airflows in the system are balanced at any moment of time. Notice that under
this CST scenario, where $T_{s}=\text{const}$, $s(T_{i}(t))$ changes with time
due to the evolution of $T_{i}$, thereby driving changes in $\mu(t)$. This
results in a change of the total power consumed by the ensemble, $P(t)$,
thereby inheriting some temporal dynamics.

### (b) Constant Power Input (CPI)

In this second scenario, indexed by $b$, we assume that the system operator
keeps the power input, $P$, constant at all times, i.e.
$P=P_{b}=\text{const}$. This is achieved by adjusting the supply temperature,
$T_{s}(t)$, according to the following modification of Eq. (6):

$\displaystyle
T_{s;b}({\bm{T}}(t))=T_{o}-\frac{{P}_{b}}{\bar{\mu}\sum_{i}s(T_{i}(t))}.$ (7)

We see that the supply temperature is continually adjusted based on the local
temperature dynamics.

##  III Statistical Analysis of the CST Scenario

The supply temperature $T_{s}$ is assumed constant, $T_{s}=T_{s;a}$. The
system of SDEs (1)-(5) is closed, thus translating into the so-called
Kolmogorov-Fokker-Planck (KPF) [33, 34] Partial Differential Equation (PDE)
for the joint probability distribution of the vector of temperatures
${\bm{T}}$ within the TCL ensemble:

$\displaystyle\partial_{t}{\cal P}_{a}({\bm{T}}|t)=\frac{1}{c}\sum_{i}\
\partial_{T_{i}}\left(f(T_{i};s(T_{i});T_{s;a})+D\partial_{T_{i}}\right){\cal
P}_{a}({\bm{T}}|t).$ (8)

We observe that the KFP Eq. (8) can be represented in the potential form,

$\displaystyle\partial_{t}{\cal P}_{a}({\bm{T}}|t)$ $\displaystyle=$
$\displaystyle-\sum_{i}\partial_{T_{i}}J_{i},$ (9) $\displaystyle J_{i}$
$\displaystyle=$
$\displaystyle-\frac{1}{c}\left(\partial_{T_{i}}U(\textbf{T}|T_{s;a})+D\partial_{T_{i}}\right){\cal
P}_{a}({\bm{T}}|t),$ (10) $\displaystyle U({\bm{T}}|T_{s;a})$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}U_{1}(T_{i}|T_{s;a}),$ (11) $\displaystyle
U_{1}(T_{i}|T_{s;a})$ $\displaystyle=$
$\displaystyle\frac{(T_{i}-T_{o})^{2}}{2r}$ $\displaystyle+$
$\displaystyle\frac{c_{p}\bar{\mu}s(T_{i})}{2}\left((T_{i}-T_{s;a})^{2}-(\bar{T}-T_{s;a})^{2}\right),$

where $J_{i}$ is the probability current along $T_{i}$, and
$U({\bm{T}}|T_{s;a})$ and $U_{1}(T_{i}|T_{s;a})$ are the aggregated and
individual TCL thermal potentials, respectively. Notice that the KFP Eq. (8)
describes a Cauchy, i.e. an initial-value, problem. It should therefore be
equipped with an initial condition, e.g. that all TCLs are at the same
temperature $T_{0}$ at $t=0$. Then, ${\cal
P}_{a}(T|0)=\prod_{i}\delta(T_{i}-T_{0})$, where $\delta(\cdot)$ is the Dirac
$\delta$-function. More generally, if ${\cal P}_{a}({\bm{T}}|0)$ is factorized
into a product of marginal distributions of individual $T_{i}$, the solution
to the KFP equations at all $t$ is also factorized into the product of the
corresponding marginal probability distributions, i.e. ${\cal
P}_{a}({\bm{T}}|t)=\prod_{i}{\cal P}_{a;1}(T_{i}|t)$.

###  III.1 Stationary Distribution

Due to the factorized structure of the differential operator on the right hand
side of Eq. (8), the steady-state solution to Eq. (8), ${\cal
P}_{a;st}({\bm{T}})$, i.e. one acquired at $t\to\infty$, is also factorized
into the product of the respective marginals, ${\cal
P}_{a;st}({\bm{T}})=\prod_{i}{\cal P}_{a;1;st}(T_{i}|T_{s;a})$. Moreover,
${\cal P}_{a;1;st}(T_{i}|T_{s;a})$ satisfies a second order ODE which can be
solved for any values of $T_{s;a}$, considered as a parameter:

$\displaystyle{\cal P}_{a;1;st}(T_{i}|T_{s;a})$
$\displaystyle=\frac{1}{Z}\exp\left(-\frac{cU_{1}(T_{i}|T_{s;a})}{D}\right),$
(13) $\displaystyle Z$
$\displaystyle=\int\limits_{-\infty}^{\infty}\\!\\!dT\exp\left(-\frac{cU_{1}(T|T_{s;a})}{D}\right).$
(14)

The single-TCL thermal potential, $U_{1}(T_{i}|T_{s;a})$, defined in Eq. (11),
is continuous in $T_{i}$. It is however not smooth as the derivative jumps at
$T_{i}=\bar{T}$, and attains a single minimum at

$\displaystyle
T_{a;min}(T_{s;a})=\left\\{\begin{array}[]{ll}\frac{T_{o}+T_{s;a}\mu^{(-)}cr}{1+\mu^{(-)}cr},&T_{s;a}\leq\beta^{(-)};\\\
\bar{T},&\beta^{(-)}\leq T_{s;a}\leq\beta^{(+)};\\\
\frac{T_{o}+T_{s;a}\mu^{(+)}cr}{1+\mu^{(+)}cr},&\beta^{(+)}\leq
T_{s;a}.\end{array}\right.$ (18)

where $\beta^{(\pm)}=\bar{T}-(T_{o}-\bar{T})/(\mu^{(\pm)}cr)$. The most
probable value of the TCL temperature in the steady state does not depend on
the amplitude of the thermal noise, $D$.

We report in Fig. 2 the dependence of the steady-state individual TCL
probability distribution ${\cal P}_{a;1;st}(T_{i}|T_{s;a})$ on $T_{i}$ at
different values of $T_{s;a}$ 111Unless noted otherwise, parameters for
numerical results are: $T_{o}=30^{\circ}C$, $s_{-}=0.2$, $s_{+}=1$,
$\bar{\mu}=1\ kg/s$, $r=2\ K/kW$, $c=15\ kJ/K$, $\bar{T}=25^{\circ}C$. . Note
that there is a range of supply temperature values where the most probable
value is achieved at $T_{i}=\bar{T}$, e.g. for $T_{s}=20^{\circ}C$ in Fig. 2.
This maximum appears because fluctuations drive the temperature across
$\bar{T}$, thereby implying repetitive switching of the airflow $s(T_{i})$,
and on average stabilizing the temperature around the setpoint.

Figure 2: Temperature distributions for an individual TCL under the CST
scenario after initial transients, for different values of the supply
temperature $T_{s}$. $c=20.5\ kJ/K$. Filled: Experimental results. Dotted
lines: Theoretical distributions obtained according to Eq. (13).

###  III.2 Relaxation to the Steady Distribution

Even in the case when the initial probability distribution cannot be
factorized into a product of independent terms, independence will appear
dynamically at sufficiently large $t\gg\tau_{a}$, where $\tau_{a}$ is the so-
called mixing time222$\tau_{a}$ is defined as a time when ${\cal
P}_{a;1}(T|t;T_{s;a})$ gets sufficiently close, wrt norm $L_{2}$ or Kullback-
Leibler divergence, to its steady asymptotics, i.e. $\|{\cal
P}_{a;1;st}(T|T_{s;a})-{\cal
P}_{a;1}(T|t;T_{s;a})\|\propto\exp(-t/\tau_{a})\ll 1$. One can extract
$\tau_{a}$ and, specifically, its dependence on the supply temperature
$T_{s;a}$ from the spectral analysis of the the $i$-th component of the
differential operator entering the KFP Eq. (8). The resulting dependence does
not hold any specific and unexpected features, and is therefore not reported
here.

##  IV Statistical Analysis of the CPI Scenario

The most significant difference between the CST and the CPI scenarios lies in
the stochastic coupling between individual TCLs. According to Eq. (7), the
supply temperature $T_{s}$ directly depends on all TCL temperatures. In this
Section, we discuss the significance of this modification.

Informally, we expect that this case is special because the stochastic
dynamics of individual TCLs are no longer decoupled. General theory of KFP
equations suggests [33, 34] that at $N>1$ in the CPI scenario (a) the thermal
force can no longer be represented as a gradient of a potential; (b) detailed
balance is broken; (c) the steady distribution is no longer a Gibbs
distribution and cannot be factorized into a product of components, each
representing an individual TCL.

In this Section, we provide a quantitative analysis of the aforementioned
expected qualitative behavior of the system. In Section IV.1, we analyze the
system in the noiseless limit, i.e. $D=0$, and show that the system is
strongly sensitive to initial conditions via the identification of two
different fixed points governing the long-time behavior of the system. In
Section IV.2, we reintroduce noise in the context of a single TCL as it
enables some analytical treatment. We leverage the latter in Section IV.3
where the general case of $N>1$ TCLs with noise is considered and where we
discover that several features from the $N=1$ case surprisingly hold.

###  IV.1 Noiseless Limit

In order to extract generic statements about the behavior of this coupled
system, we analyze the fixed points of Eq. (1) with $T_{s}$ substituted by
$T_{s;b}({\bm{T}}(t))$ from Eq. (7) in the noiseless case, i.e. in the
deterministic regime. As TCLs are indistinguishable, a fixed point solution is
determined by the number of TCLs in the high-airflow regime,
$N_{+}(t)=\sum_{i}\theta(T_{i}(t)>\bar{T})$, where $\theta(\cdot)$ is the
indicator function. We observe that there are potentially $N+1$ fixed point
solutions to Eq. (1) complemented with the power balance in the CPI scenario,
i.e. $N_{+}(t)\in\\{0,\dots,N\\}$. According to Eq. (7), $T_{s;b}$ is directly
determined by $N_{+}(t)$:
$T_{s;b}(N_{+}(t))=T_{o}-P_{b}/\bar{\mu}(N_{+}s_{+}+(N-N_{+})s_{-})$. For each
of the $N+1$ possible states, the $N$ TCLs are split into two groups: a group
of $N_{+}$ units in the high-flow regime and another $N-N_{+}$ in the low-flow
regime. In addition, in the deterministic case, all TCLs in a given group are
at the same temperature at steady-state, given by

$\displaystyle
T_{\pm}(N_{+}(t))=\frac{T_{o}+T_{s;b}(N_{+}(t))\mu_{\pm}cr}{1+\mu_{\pm}cr}.$
(19)

An immediate consequence of this equation is that
$T_{+}(N_{+}(t))<T_{-}(N_{+}(t))$, providing that $\mu_{+}>\mu_{-}$. However,
this is in contradiction with the requirement that
$T_{+}(N_{+}(t))>\bar{T}>T_{-}(N_{+}(t))$. In other words, there can be no
stable fixed points, i.e. states realized dynamically at $t\to\infty$ in the
noiseless regime, with a nonzero number of TCLs in both states (high-flow and
low-flow). Therefore two options are left for stable fixed points:
$N_{+}(t)=0$ or $N_{+}(t)=N$. Which of the two stable fixed points is attained
at long times will depend on the value of $\bar{T}$. Indeed, the $N_{+}(t)=0$
state is valid (self-consistent) if $T_{-}(N_{+}(t)=0)<\bar{T}$ and the
$N_{+}(t)=N$ state is valid if $T_{+}(N_{+}(t))>\bar{T}$. In addition,
$T_{-}(N_{+}(t)=0)<T_{+}(N_{+}(t)=N)$. Therefore, the steady state that will
be attained will consist of (i) $N_{+}(t)=N$ if
$\bar{T}<T_{-}(N_{+}(t)=0)<T_{+}(N_{+}(t)=N)$, (ii) $N_{+}(t)=0$, if
$T_{-}(N_{+}(t)=0)<T_{+}(N_{+}(t)=N)<\bar{T}$, or (iii) either of the two if
$T_{-}(N_{+}(t)=0)<\bar{T}<T_{+}(N_{+}(t)=N)$, in which case the final state
is directly dictated by the initial conditions $T(0)$. Fig. 3 illustrates the
dependence of the realized fixed points on the values of $\bar{T}$ for
different initial conditions in the case of $N=2$ TCLs. $T_{\pm}$ do indeed
act as stable fixed points in the noiseless limit. Depending on $\bar{T}$,
either one or both of them are accessible to the system.

Figure 3: Location of the two fixed points for $N=2$ in the CPI scenario. The
identity is represented by the gray diagonal line. Green line, dashed:
Expected position of the two fixed points according to Eq. (19). Blue dots:
Final temperatures of both TCLs attained after $100\mathrm{s}$ under multiple
initializations.

###  IV.2 Case of $N=1$, $D>0$.

As analysis in the general case with $N>1,\ D>0$ is difficult, we here focus
on the stochastic case for 1 TCL. Eqs (III), (13) still hold in this context,
with the difference that the supply temperature is dictated by Eq. (7). A
similar analysis as that conducted in Section III is hence possible, and
results are reported in Fig. 4.

(a) Temperature distribution of 1 TCL under the CPI scenario after initial
transients, for different values of the noise amplitude $D_{i}$. Filled:
Experimental results. Dashed line: Theoretical distributions obtained
according to Eq. (13).

(b) Time series of 1 TCL under the CPI scenario, with $D_{i}=15\ kJ^{2}$. Red
lines, dashed: Locations of the theoretical fixed points $T_{\pm}$ (cf. Eq.
(19)). Filled background: A change in background color indicates a change in
the value of $T_{min}(T_{s})$.

Figure 4: Dynamics and steady state distribution of 1 TCL under the CPI
scenario.

In contrast to the single-peak distributions observed in the CST scenario (cf.
Fig. 2), one observes the emergence of two maxima corresponding with the two
fixed points identified in Section IV.1. The position of these stable
temperatures is directly determined by the physical parameters of the system,
while their relative importance can be tuned by the setpoint $\bar{T}$. We
attribute the mismatch between analytical and empirical distributions in Fig.
4(a) to a finite computational budget and step size. The stochastic nature of
the system actually enables switching between both modes. As illustrated in
Fig. 4(b), a single TCL driven by noise successively switches between the two
equilibrium points — thereby implying an adaptation of the supply air
temperature $T_{s}$, not shown here.

As the potential $U_{1}(T)$ has two distinct minima in some regimes, Kramers’
theory [26, 27] suggests the emergence of two distinct temporal scales. One
short time scale $\tau_{s}$ should be associated to equilibrium fluctuations
around either of the two minima, while the transition between either minimum
as seen in Fig. 4(b) is characterized by a longer one $\tau_{l}$. These are
related via $\tau_{l}\sim\tau_{s}\exp(\Delta/D)$, where $\Delta$ is the
potential barrier at $T_{1}=\bar{T}$. Strong separation of these timescales,
i.e. $\tau_{l}\gg\tau_{s}$, emerges when $D\ll\Delta$.

###  IV.3 Numerical Experiments for $N>1,D>0.$

Even if, strictly speaking, insights from Kramers’ theory only apply to the
$N=1$ case discussed above, most qualitative features seem to translate to the
general $N>1$ case.

#### Bi-modal dynamics and emergence of two timescales

Up to moderate values of the noise amplitude $D$ a similar bi-modal
distribution of accessed temperatures is observed (cf. Fig. 5(a), with $D=10\
kJ^{2}$). Oscillations of temperatures around two central values are clearly
visible in Fig. 5(b). Note that both TCLs illustrated exhibit coordinated
switching from one fixed point to the other, resulting from the AHU coupling.
Both observations are qualitatively analogous to the $N=1$ case, cf. Fig. 4.
Specifically, one can also expect switching and the natural appearance of two
characteristic timescales in the system. As shown in Fig. 6(a), the
autocorrelation $\rho(T(t)-T_{min}(t))$ of a given TCL temperature decays
exponentially fast. However, the decay of $\rho(T_{min}(t))$, associated with
the $N$ TCLs switching from one stable point to the next, is significantly
slower in Fig. 6(b). This suggests that memory is indeed present in the system
and that switching is a rarer event with increasing values of $N$. Fig. 7
provides another illustration of this strong $N$-dependence. The switching
time over trials increases with $N$, confirming the intuition provided by Fig.
6(b). While the dependence of the switching time with $N$ seems sub-
exponential, those results were obtained with a finite computational budget
and might not reflect the true scaling. Further investigation is needed, via
alternative theoretical and computational methods, to resolve it.

(a) Empirical temperature distributions of TCLs under the CPI scenario with
$N=5$, after initial transients and for different values of the noise
amplitude $D_{i}$.

(b) Time series of 2 TCLs under the CPI scenario, with $N=2$. Red lines,
dashed: Locations of the theoretical fixed points $T_{\pm}$. Filled
background: A change in background color indicates a switch in the value of
$T_{min}(T_{s})$.

Figure 5: Temperature distribution and dynamics in the general case of
$N>1,D>0$ CPI case.

(a) Autocorrelation of $T(t)-T_{min}(t)$. The exponential decay is indicative
of a short associated timescale.

(b) Autocorrelation of $T_{min}$, defined in (18). The decay is much slower
than for Fig. 6(a), and the characteristic timescale increases with the number
of TCLs in the system.

Figure 6: Statistical analysis and observation of two different timescales in
the system. Figure 7: Scaling of the switching time with the number $N$ of
TCLs in the CPI scenario. Dots indicate the average switching time. Errorbars
indicate the 75th and 25th percentile of the distributions.

#### Emergence of an entropic state

Large noise amplitudes will prevent the system from settling in either of the
stable points. Indeed, the two modes of the temperature distribution disappear
as $D$ increases (see Fig. 5(a)). Instead, they are slowly replaced by an
entropic peak at the location of the setpoint $\bar{T}$. From the perspective
of a single TCL, this entropic state represents a balance between stochastic
and deterministic forces resulting in the TCL switching back and forth between
low- and high-flow regimes, and hence meandering around $\bar{T}$. In high
noise regimes, the effect of the potential barrier disappears and the
resulting distribution is akin in nature to the CST case in some regimes as
shown in Fig. 2. However, the noise intensity required to reach such regimes
is such that the distribution broadens significantly with potential impact on
occupants.

##  V Conclusions and Path Forward

In this manuscript, we formulate a novel model describing the stochastic
dynamics of coupled TCLs. It differs from standard TCL models in that it more
accurately replicates cooling via a centralized Air Handling Unit. The main
consequence of this modification is the resulting coupling between individual
units, which is usually disregarded in other models.

We analyzed the aggregate dynamics in two different scenarios. First, we
considered the uncoupled dynamics of each unit, by fixing the supply air
temperature (CST regime). While the CST regime is highly stable and
predictable, fluctuations of the aggregated power consumption of the building,
associated with airflow adjustments in each unit, might be undesirable. In the
second CPI scenario, the power input is set constant, leading to significant
coupling between individual TCLs and revealing interesting collective
phenomena, with ramifications for the use of the ensemble in DR. The system’s
response to a perturbation may result in the co-existence of two quasi-steady
modes associated with stable fixed points of the dynamics in the low-noise
regime, which transition to an entropic state when the level of the noise
increases. Dynamic consequence of this complex behavior is the emergence of
long, and generally parasitic, transients associated with the transitions
between the two modes in the low-noise regime.

This new model and associated results set the stage for several future
challenges:

  1. 1.

Investigating the effects associated with inhomogeneity (disorder) in physical
parameters, thereby describing populations of non-identical units, would
provide a more accurate picture of real-life systems.

  2. 2.

Developing methods leveraging both unit-level and system-level control to
strike a balance between individual comfort and collective benefits would
further our understanding of the inherent flexibility in such systems, and
potentially result in innovative solutions for AHU control — e.g. to remove
the anomalously long relaxation times observed in the low-noise CPI regime.
Ideas in recent papers [19, 20] could be translated to the present case.
Specifically, understanding the non-equilibrium nature of the TCL ensemble in
the CPI regime and its practical significance for operations and control will
be interesting. While similar to typical Kramers phenomena, the dynamics in
the CPI regime are not at statistical equilibrium. Detailed balance is broken
and fluxes emerge. The implications of such complex dynamics on
controllability remain to be explored.

  3. 3.

Confronting this new model with real zone-level time series and experiments
[37, 38] will provide crucial insight for its validity and highlight where
improvements, modifications or relaxing assumptions are needed.

  4. 4.

Analyzing the integration of such systems into larger-scale DR schemes, e.g.
in the spirit of [39, 40], and involving both multiple buildings as well as
other appliances, will be critical in the deployment of such solutions.

##  VI Acknowledgments

LFV acknowledges Jacques de Chalendar for inspiring and helpful discussions,
as well as for proof-reading the manuscript. Research Experience for
Undegraduates work of LW at UArizona during the summer of 2021 was funded by
NSF # 1937229 "RTG: Applied Mathematics and Statistics for Data-Driven
Discovery".

## References

  * Palensky and Dietrich [2011] P. Palensky and D. Dietrich, “Demand side management: Demand response, intelligent energy systems, and smart loads,” IEEE Transactions on Industrial Informatics 7, 381–388 (2011).
  * Oconnell  _et al._ [2014] N. Oconnell, P. Pinson, H. Madsen, and M. Omalley, “Benefits and challenges of electrical demand response: A critical review,” Renewable and Sustainable Energy Reviews 39, 686–699 (2014).
  * D’hulst  _et al._ [2015] R. D’hulst, W. Labeeuw, B. Beusen, S. Claessens, G. Deconinck, and K. Vanthournout, “Demand response flexibility and flexibility potential of residential smart appliances: Experiences from large pilot test in Belgium,” Applied Energy 155, 79–90 (2015).
  * Yin  _et al._ [2016] R. Yin, E. C. Kara, Y. Li, N. DeForest, K. Wang, T. Yong, and M. Stadler, “Quantifying flexibility of commercial and residential loads for demand response using setpoint changes,” Applied Energy 177, 149–164 (2016).
  * Finck  _et al._ [2018] C. Finck, R. Li, R. Kramer, and W. Zeiler, “Quantifying demand flexibility of power-to-heat and thermal energy storage in the control of building heating systems,” Applied Energy 209, 409–425 (2018).
  * Xu, Haves, and Piette [2010] P. Xu, P. Haves, and M. A. Piette, “Demand Shifting With Thermal Mass in Large Commercial Buildings: Field Tests, Simulation and Audits,” , 35–43 (2010).
  * Kara  _et al._ [2014] E. C. Kara, M. D. Tabone, J. S. MacDonald, D. S. Callaway, and S. Kiliccote, “Quantifying flexibility of residential thermostatically controlled loads for demand response: A data-driven approach,” BuildSys 2014 - Proceedings of the 1st ACM Conference on Embedded Systems for Energy-Efficient Buildings , 140–147 (2014).
  * Aghniaey and Lawrence [2018] S. Aghniaey and T. M. Lawrence, “The impact of increased cooling setpoint temperature during demand response events on occupant thermal comfort in commercial buildings: A review,” Energy and Buildings 173, 19–27 (2018).
  * Chong and Debs [1979] C. Y. Chong and A. S. Debs, “Statistical synthesis of power system function load models,” Proceedings of the IEEE Conference on Decision and Control 1, 264–269 (1979).
  * Ihara and Schweppe [1981] S. Ihara and F. C. Schweppe, “Physically based modeling of cold load pickup,” IEEE Transactions on Power Apparatus and Systems PAS-100, 4142–4150 (1981).
  * Malhamé and Chong [1984] R. Malhamé and C. Y. Chong, “Statistical synthesis of physically based load models with applications to cold load pickup,” IEEE Transactions on Power Apparatus and Systems PAS-103, 1621–1628 (1984).
  * Malhamé and Chong [1985] R. Malhamé and C. Y. Chong, “Electric Load Model Synthesis by Diffusion Approximation of a High-Order Hybrid-State Stochastic System,” IEEE Transactions on Automatic Control 30, 854–860 (1985).
  * Malhamé and Chongt [1988] R. Malhamé and C.-Y. Chongt, “On the statistical properties of a cyclic diffusion process arising in the modeling of thermostat-controlled electric power system loads,” SIAM J. APPL. MATH 48 (1988).
  * Chertkov and Chernyak [2017] M. Chertkov and V. Chernyak, “Ensemble of Thermostatically Controlled Loads: Statistical Physics Approach,” Scientific Reports 7, 8673 (2017).
  * Chertkov, Deka, and Dvorkin [2018] M. Chertkov, D. Deka, and Y. Dvorkin, “Optimal ensemble control of loads in distribution grids with network constraints,” in  _2018 Power Systems Computation Conference (PSCC)_ (2018) pp. 1–7.
  * Callaway [2009] D. S. Callaway, “Tapping the energy storage potential in electric loads to deliver load following and regulation, with application to wind energy,” Energy Conversion and Management 50, 1389–1400 (2009).
  * Callaway and Hiskens [2011] D. S. Callaway and I. A. Hiskens, “Achieving controllability of electric loads,” Proceedings of the IEEE 99, 184–199 (2011).
  * Bashash and Fathy [2011] S. Bashash and H. K. Fathy, “Modeling and control insights into demand-side energy management through setpoint control of thermostatic loads,” Proceedings of the American Control Conference , 4546–4553 (2011).
  * Métivier, Luchnikov, and Chertkov [2019] D. Métivier, I. Luchnikov, and M. Chertkov, “Power of Ensemble Diversity and Randomization for Energy Aggregation,” Scientific Reports 9, 5910 (2019).
  * Métivier and Chertkov [2020] D. Métivier and M. Chertkov, “Mean-field control for efficient mixing of energy loads,” Physical Review E 101, 022115 (2020).
  * Claessens  _et al._ [2018] B. J. Claessens, D. Vanhoudt, J. Desmedt, and F. Ruelens, “Model-free control of thermostatically controlled loads connected to a district heating network,” Energy and Buildings 159 (2018), 10.1016/j.enbuild.2017.08.052.
  * Webborn and Mackay [2017] E. Webborn and R. S. Mackay, “A stability analysis of thermostatically controlled loads for power system frequency control,” Complexity 2017 (2017), 10.1155/2017/5031505.
  * Bashash and Fathy [2013] S. Bashash and H. K. Fathy, “Modeling and control of aggregate air conditioning loads for robust renewable power management,” IEEE Transactions on Control Systems Technology 21 (2013), 10.1109/TCST.2012.2204261.
  * Ma  _et al._ [2012] Y. Ma, A. Kelman, A. Daly, and F. Borrelli, “Predictive Control for Energy Efficient Buildings with Thermal Storage,” IEEE Control Systems 32, 44–64 (2012).
  * MacDonald, Vrettos, and Callaway [2020] J. S. MacDonald, E. Vrettos, and D. S. Callaway, “A Critical Exploration of the Efficiency Impacts of Demand Response from HVAC in Commercial Buildings,” Proceedings of the IEEE 108, 1623–1639 (2020).
  * Kramers [1940] H. Kramers, “Brownian motion in a field of force and the diffusion model of chemical reactions,” Physica 7, 284–304 (1940).
  * Hänggi, Talkner, and Borkovec [1990] P. Hänggi, P. Talkner, and M. Borkovec, “Reaction-rate theory: fifty years after kramers,” Rev. Mod. Phys. 62, 251–341 (1990).
  * Taylor  _et al._ [2012] S. T. Taylor, J. Stein, G. Paliaga, and H. Cheng, “Dual maximum VAV box control logic,” ASHRAE Journal 54, 16–24 (2012).
  * Sterling  _et al._ [2014] R. Sterling, G. Provan, J. Febres, D. O’Sullivan, P. Struss, and M. M. Keane, “Model-based fault detection and diagnosis of air handling units: A comparison of methodologies,” in  _Energy Procedia_, Vol. 62 (2014).
  * Razban  _et al._ [2019] A. Razban, A. Khatib, D. Goodman, and J. Chen, “Modelling of air handling unit subsystem in a commercial building,” Thermal Science and Engineering Progress 11 (2019), 10.1016/j.tsep.2019.03.019.
  * Zajic  _et al._ [2011] I. Zajic, T. Larkowski, M. Sumislawska, K. J. Burnham, and D. Hill, “Modelling of an air handling unit: A Hammerstein-bilinear model identification approach,” in  _Proceedings - ICSEng 2011: International Conference on Systems Engineering_ (2011).
  * Afram and Janabi-Sharifi [2014] A. Afram and F. Janabi-Sharifi, “Review of modeling methods for HVAC systems,” Applied Thermal Engineering 67, 507–519 (2014).
  * Gardiner [2004] C. W. Gardiner,  _Handbook of stochastic methods for physics, chemistry and the natural sciences_ , Vol. 13 (Berlin: Springer-Verlag, 2004).
  * Kampen [2007] N. V. Kampen,  _Stochastic processes in physics and chemistry_ (North Holland, 2007).
  * Note [1] Unless noted otherwise, parameters for numerical results are: $T_{o}=30^{\circ}C$, $s_{-}=0.2$, $s_{+}=1$, $\bar{\mu}=1\ kg/s$, $r=2\ K/kW$, $c=15\ kJ/K$, $\bar{T}=25^{\circ}C$.
  * Note [2] $\tau_{a}$ is defined as a time when ${\cal P}_{a;1}(T|t;T_{s;a})$ gets sufficiently close, wrt norm $L_{2}$ or Kullback-Leibler divergence, to its steady asymptotics, i.e. $\|{\cal P}_{a;1;st}(T|T_{s;a})-{\cal P}_{a;1}(T|t;T_{s;a})\|\propto\mathop{exp}\nolimits(-t/\tau_{a})\ll 1$.
  * Keskar  _et al._ [2018] A. Keskar, S. Lei, T. Webb, S. Nagy, H. Lee, I. Hiskens, J. Mathieu, and J. Johnson, “Stay cool and be flexible : Energy-efficient grid services using commercial buildings HVAC systems,” (2018).
  * Li  _et al._ [2015] P. Li, D. Vrabie, D. Li, S. C. Bengea, S. Mijanovic, and Z. D. O’Neill, “Simulation and experimental demonstration of model predictive control in a building hvac system,” Science and Technology for the Built Environment 21, 721–733 (2015).
  * Hassan  _et al._ [2019] A. Hassan, R. Mieth, M. Chertkov, D. Deka, and Y. Dvorkin, “Optimal load ensemble control in chance-constrained optimal power flow,” IEEE Transactions on Smart Grid 10, 5186–5195 (2019).
  * Hassan  _et al._ [2020] A. Hassan, S. Acharya, M. Chertkov, D. Deka, and Y. Dvorkin, “A hierarchical approach to multienergy demand response: From electricity to multienergy applications,” Proceedings of the IEEE 108, 1457–1474 (2020).


# Arbitrary Ground State Observables from Quantum Computed Moments

Harish J. Vallury1 and Lloyd C. L. Hollenberg2  0009-0007-3355-5458
0000-0001-7672-6965 School of Physics, The University of Melbourne
Parkville, VIC 3010, Australia
Email<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The determination of ground state properties of quantum systems is a
fundamental problem in physics and chemistry, and is considered a key
application of quantum computers. A common approach is to prepare a trial
ground state on the quantum computer and measure observables such as energy,
but this is often limited by hardware constraints that prevent an accurate
description of the target ground state. The quantum computed moments (QCM)
method has proven to be remarkably useful in estimating the ground state
energy of a system by computing Hamiltonian moments with respect to a
suboptimal or noisy trial state. In this paper, we extend the QCM method to
estimate arbitrary ground state observables of quantum systems. We present
preliminary results of using QCM to determine the ground state magnetisation
and spin-spin correlations of the Heisenberg model in its various forms. Our
findings validate the well-established advantage of QCM over existing methods
in handling suboptimal trial states and noise, extend its applicability to the
estimation of more general ground state properties, and demonstrate its
practical potential for solving a wide range of problems on near-term quantum
hardware.

##  I Introduction

Many difficult problems with practical application in physics and chemistry
involve determining the ground state properties of a quantum system. This has
been a popular road to follow towards the development of near-term
applications of quantum computers, with much of the recent focus primarily on
variational quantum algorithms [1, 2, 3]. These hybrid approaches involve the
preparation of complicated states designed to represent the ground state of
the desired quantum system, via the measurement and minimisation of energy [4,
5, 6]. Most interesting problems at scale require high-depth circuits to
prepare such states, and are thus subject to noise which is often severely
detrimental to the quality of results [7, 8].

Some time ago, an analytic cluster expansion of the Lanczos recursion was
uncovered [9]. Following from this, it was found that for a given quantum
system Hamiltonian $H$, one could compute Hamiltonian moments $\langle
H^{k}\rangle$ with respect to a trial state with some overlap with the ground
state to produce a powerful approximation to the ground state energy [10, 11],
shown below to fourth order in the moments:

$E_{\rm
0}^{\rm{L}(4)}=c_{1}-\frac{c_{2}^{2}}{c_{3}^{2}-c_{2}c_{4}}\left[\sqrt{3c_{3}^{2}-2c_{2}c_{4}}-c_{3}\right],$
(1)

where the cumulants $c_{n}$ are defined as:

$c_{n}=\langle H^{n}\rangle-\sum_{p=0}^{n-2}\binom{n-1}{p}c_{p+1}\langle
H^{n-1-p}\rangle.$ (2)

The rapid development of quantum computer technology has provided a fitting
application for this result, allowing for advanced state preparation and
measurement of these Hamiltonian moments in order to solve interesting
problems. We refer to this as the quantum computed moments (QCM) method, and
recent work applied to problems in quantum magnetism [12, 13] and chemistry
[14] has demonstrated the effectiveness of QCM in estimating the ground state
energy of a system, even when moments are measured with respect to a
suboptimal or noisy trial state. This technique can correct for a low
complexity ansatz and a high degree of noise, proving to be remarkably useful
on present-day quantum hardware for a variety of different problem instances.

In recent times, these QCM studies have been accompanied elsewhere in the
literature by a general renewed interest in Hamiltonian moments in a quantum
computing context [15, 16, 17, 18], with the primary focus being on the ground
state energy problem. However, a problem which is much more practically useful
and indeed more computationally complex is that of determining arbitrary
ground state observables [19]. It turns out that Hamiltonian moments can also
be used to estimate general ground state properties via a Hellman-Feynman
approach [20, 21]. Here we develop this technique in the QCM context. When
applied to models of interest, our results show that the crucial robustness
property of QCM for energy problems indeed carries over to this more general
ground state observable procedure.

##  II Method

Suppose we wish to find the expectation value of some general observable
(given by Hermitian operator $A$) with respect to the ground state
$\ket{\psi_{0}}$ of a Hamiltonian $H$. Consider the following parameterised
Hamiltonian:

$\displaystyle H_{\lambda}=H+\lambda A$ (3)

The parameterised Hamiltonian $H_{\lambda}$ has ground state
$\ket{\psi_{\lambda,0}}$ with energy $E_{\lambda,0}$ which we can approximate
with:

$E_{\lambda,0}\approx E^{{\rm L}(4)}_{\lambda,0}|_{\psi_{t}}.$ (4)

The approximation in Equation (4) is expressed entirely in terms of
Hamiltonian moments $\langle\psi_{t}|(H_{\lambda})^{k}|\psi_{t}\rangle$, up to
$k=4$, with respect to some trial state $\ket{\psi_{t}}$ which has nonzero
overlap with $\ket{\psi_{\lambda,0}}$.

Starting with the Hellmann-Feynman theorem,

$\langle\psi_{\lambda,0}|\,\diff{H_{\lambda}}{\lambda}\,|\psi_{\lambda,0}\rangle=\diff{E_{\lambda,0}}{\lambda},$
(5)

one can substitute the derivative of the Hamiltonian in Equation (3) and use
the approximation in Equation (4) to obtain:

$\langle\psi_{\lambda,0}|\,A\,|\psi_{\lambda,0}\rangle\approx\frac{\rm
d}{\mathrm{d}\lambda}\left(E^{{\rm L}(4)}_{\lambda,0}|_{\psi_{t}}\right).$ (6)

Thus if we want to compute the observable $A$ with respect to the ground state
$\ket{\psi_{0}}$ of $H$, we can choose $\ket{\psi_{t}}$ to have some overlap
with $\ket{\psi_{0}}$ and evaluate the derivative in Equation (6) around
$\lambda=0$:

$\langle\psi_{0}|\,A\,|\psi_{0}\rangle\approx\frac{\rm
d}{\mathrm{d}\lambda}\left(E^{{\rm
L}(4)}_{\lambda,0}|_{\psi_{t}}\right)_{\lambda=0}.$ (7)

This leads to a practically useful final expression for the ground state
observable:

$\langle A\rangle_{0}^{\rm
L(4)}\equiv\frac{1}{2\varepsilon}\left(E_{+\varepsilon,0}^{\rm
L(4)}-{E_{-\varepsilon,0}^{\rm
L(4)}}\right)\approx\langle\psi_{0}|\,A\,|\psi_{0}\rangle,$ (8)

where $\varepsilon>0$ is sufficiently small as to find the $\lambda=0$
derivative, and the $E_{\pm\varepsilon,0}^{\rm L(4)}$ are calculated via
Equation (1) from the moments $\langle\,(H\pm\varepsilon A)^{k}\,\rangle$ with
respect to some trial state $\ket{\psi_{t}}$ up to $k=4$.

We expand upon the analytic studies of this method, first introduced in [20,
21], by now considering the approach in a quantum computing context. Note
that, in practice, the $E_{+\varepsilon,0}^{\rm L(4)}$ and
${E_{-\varepsilon,0}^{\rm L(4)}}$ in Equation 8 are not computed separately,
but rather come from the same set of measurements on the quantum computer. The
derivative is taken with respect to a parameter only introduced in the
classical post-processing stage, so remains unaffected by the stochastic
nature of the quantum measurement. We can thus expect to see similar
robustness behaviour as has been observed in the QCM energy estimate, and in
the following section, we test the effectiveness of this arbitrary observables
technique by performing noisy simulations on a variety of trial states and
Hamiltonians.

##  III Results

Figure 1: (a) Ground state energy estimates for the Heisenberg XXZ model as a
function of $x$ in the Hamiltonian. Green shows the exact ground state energy,
$E_{0}$. The energy estimates, $\langle H\rangle$ (blue) and $E_{0}^{\rm
L(4)}$ (orange), are evaluated with respect to three trial states: the $x=-1$
(ferromagnetic) ground state for $-1<x<-\frac{1}{2}$, the $x=0$ (Néel) ground
state for $-\frac{1}{2}<x<\frac{1}{2}$, and the $x=1$ (antiferromagnetic)
ground state for $\frac{1}{2}<x<1$. (b) Estimates of ground state $ZZ$
correlation between two next-nearest-neighbour spins (marked in red) versus
$x$, using the same three trial states. Like before, the exact correlation is
shown in green, the direct expectation value in blue, and the moments-based
estimate in orange. Figure 2: Staggered magnetisation for a 6-site 1D
Heisenberg antiferromagnet with staggered external field, as a function of
field parameter $g$ in the Hamiltonian (Eq. 11). Comparison of methods for
obtaining ground state staggered magnetisation – direct $\langle M\rangle$
(blue) and QCM $\langle M\rangle_{0}^{\rm L(4)}$ (orange) estimates vs. exact
value (green). To compute the estimates, trial states $\ket{\psi_{t}^{g}}$
were generated at each value of $g$ by starting at the ground state
$\ket{\psi_{0}^{g}}$ and applying a random rotation as
$\ket{\psi_{t}^{g}}=e^{-i\theta M}\ket{\psi_{0}^{g}}$, where $M$ is a random
Hermitian matrix sampled from the Gaussian unitary ensemble. By tuning the
parameter $\theta$, trial states of varying degrees of overlap with the ground
state $F=\left|\langle\psi_{t}^{g}|\psi_{0}^{g}\rangle\right|^{2}$ were
generated. These trial states ${\rho}=\ket{\psi_{t}^{g}}\\!\bra{\psi_{t}^{g}}$
were then passed through a depolarising error channel
${\rho}\mapsto(1-\frac{3p}{4}){\rho}+\frac{p}{4}(X{\rho}X+Y{\rho}Y+Z{\rho}Z)$
with noise parameter $p$, which was varied from 1% up to 50%.

Of particular interest in the study of quantum magnetism is the Heisenberg XXZ
model, given by the Hamiltonian:

$H=\frac{1}{4q}\sum_{\langle
i\,j\rangle}\left(Z_{i}Z_{j}+x(X_{i}X_{j}+Y_{i}Y_{j})\right).$ (9)

Here, $q$ is the number of qubits, the $S_{i}$ are the Pauli operators
$S\in\\{X,Y,Z\\}$ acting on the $i^{\rm th}$ qubit, the $\langle i\,j\rangle$
are the nearest neighbour couplings in the lattice. The parameter $x$ controls
the strength of the transverse spin components relative to the Ising-like
interaction. As a first test case of the moments-based ground state observable
estimate (Eq. 8), we take the XXZ model over a 2D square lattice, and consider
the $ZZ$ correlation operator,

$C_{ij}=Z_{i}Z_{j},$ (10)

acting on a pair of next-nearest-neighbour qubits $i$ and $j$. Focusing on
these longer range interactions provides additional insights into the system
behaviour beyond what is captured explicitly by the Hamiltonian.

On a quantum computer, to simulate and study the ground state behaviour (e.g.
the energy and spin-spin correlations) of this system as the parameter $x$
varies, one might need to prepare the ground state using a different circuit
at each value of $x$. However, one can perform far fewer computations by
leveraging the fact that the QCM approach has reduced sensitivity to the
choice of trial state. A well-chosen state that accurately describes the
solution at a given value of $x$ will also produce an accurate QCM estimate
for nearby values of $x$. Figure 1(a) shows that the ground state energy of
the Heisenberg XXZ model can be fully characterised for all $-1<x<1$ using the
QCM estimate $E_{0}^{\rm L(4)}$ with respect to just three different trial
states (ground states at $x\in\\{-1,0,1\\}$), in stark contrast to $\langle
H\rangle$ from the standard variational approach.

In Figure 1(b), the difference between the two methods is made even clearer
when these states are used to find estimates of other ground state properties.
QCM correctly estimates the correlation in the ground state between the corner
site and its next-nearest-neighbour in the middle of the $4\times 3$ square
lattice, remaining remarkably close to the exact value as $x$ varies,
diverging only slightly at $x=\pm\frac{1}{2}$ where the ground state overlap
is at a minimum.

In practice on a real device, measuring the moments up to $H^{4}$ requires one
to multiply out the terms (or Pauli strings) in the Hamiltonian and determine
their expectation values. These Pauli strings can be grouped into
simultaneously measurable tensor product basis (TPB) sets, the total number of
which scales logarithmically with the number of Pauli strings [12]. Each TPB
corresponds to one measurement on the quantum computer – thus when determining
a general ground state observable using QCM, if the operator describes
interactions present explicitly in the Hamiltonian, no additional measurements
need to be made over the energy calculation. In our Heisenberg XXZ case, to go
from the energy calculation in Figure 1(a) to the next-nearest-neighbour
correlation in Figure 1(b), the total number of Pauli strings increases from
66 343 to 68 960, and the total number of TPB measurements goes from 1 906 to
1 973. Note that these measurements need to only be made once for each of the
three trial states, after which the entire curve for $-1<x<1$ can be filled
out via simple classical post-processing.

So far, we have shown that the QCM estimate of a ground state observable other
than energy shows the same resilience to an imperfect choice of trial state,
however the presence of the noise robustness property still remains to be
verified. For this, we turn our attention to the antiferromagnetic Heisenberg
model with an external (staggered) field:

$H=\frac{1}{4}\sum_{\langle
i\,j\rangle}(X_{i}X_{j}+Y_{i}Y_{j}+Z_{i}Z_{j})+\frac{g}{2}\sum_{k}(-1)^{k}Z_{k},$
(11)

where $g$ is the external field parameter. An operator of interest in this
case is the staggered magnetisation $M$, given by:

$M=\frac{1}{2}\sum_{k}(-1)^{k}Z_{k}.$ (12)

We look at this model over 6 qubits in 1D. To simulate the results one might
typically see in a variational quantum algorithm, given that preparing ansaẗze
that can fully represent the target state is often challenging, we approximate
an intermediate step of such an algorithm by starting with the exact ground
state vector $\ket{\psi_{0}^{g}}$ for each $g$ and applying random rotations
to obtain a “trial state” $\ket{\psi_{t}^{g}}$. Both the direct and moments-
based estimates of the staggered magnetisation are then calculated with
respect to this trial state after applying various levels of depolarising
noise. Figure 2 shows the interplay between these two phenomena in our
simulations for varying average ground state fidelity,
$F=\left|\langle\psi_{t}^{g}|\psi_{0}^{g}\rangle\right|^{2}$, and depolarising
noise,
${\rho}\mapsto(1-\frac{3p}{4}){\rho}+\frac{p}{4}(X{\rho}X+Y{\rho}Y+Z{\rho}Z)$
with parameter $p$. We observe that the robustness of the QCM technique to low
fidelity trial states and noise, previously only seen in the context of energy
problems [12, 13, 14], persists when considering arbitrary ground state
observables.

We see even for trial states with $\sim$40% ground state fidelity and 50%
depolarising error, at each value of $g$, the moments based estimate $\langle
M\rangle_{0}^{\rm L(4)}$ provides a powerful correction to the direct
measurement of $\langle M\rangle$, remaining remarkably close to the exact
value for the ground state staggered magnetisation. This is more or less the
level of noise one would expect to see in a typical variational circuit on
real hardware with 1% CNOT error, so the result bodes well for the prospect of
applying this technique on present-day devices.

##  IV Conclusion

We have presented an extension of the QCM method to estimate arbitrary ground
state observables of quantum systems, with a focus on the Heisenberg model. We
observe the robustness of the QCM technique to low fidelity trial states and
noise, previously only seen in the context of ground-state energy problems. In
particular, the results presented in this paper demonstrate the ability of QCM
to estimate ground state properties such as magnetisation and spin-spin
correlations of the Heisenberg model to a high precision, even with suboptimal
trial states and levels of noise present in near-term quantum hardware.

Further advancement of this technique will lead to improved understanding of
the conditions under which various QCM observable estimates hold greater or
lesser validity. We reiterate that the results presented here are from
simplified noisy simulations, and the key next steps will involve performing
the technique on real quantum hardware and under improved noise models –
previous analysis of the QCM energy estimate in [13] in addition to the
preliminary results presented here for the generalisation to arbitrary
observables suggests that we should also expect this technique to work quite
well on a real noisy device. Furthermore, the examples we have studied here
have only scratched the surface of possible Hamiltonians and observables that
could be suitable for this approach. Many intriguing problems invite future
application of this technique, such as calculating ground state observables
for the Fermi-Hubbard model [22] and detecting phase transitions in more
unconventional systems by measuring nonlocal order parameters [23].

QCM has so far proven to be an exciting and potentially useful technique for
obtaining the ground-state of many body systems on quantum computers, even in
the presence of noise and with suboptimal trial states. Our findings here
generalise the method to other observables of interest and further highlight
its practical potential for solving a wide range of problems on near-term
quantum hardware.

## Acknowledgements

This research was supported by the University of Melbourne through the
establishment of the IBM Quantum Network Hub at the University. HJV is
supported by an Australian Government Research Training Program Scholarship.

## References

  * [1] Yudong Cao et al.  “Quantum chemistry in the age of quantum computing”  In _Chemical reviews_ 119.19 ACS Publications, 2019, pp. 10856–10915  DOI: 10.1021/acs.chemrev.8b00803
  * [2] Sam McArdle et al.  “Quantum computational chemistry”  In _Reviews of Modern Physics_ 92.1 APS, 2020, pp. 015003  DOI: 10.1103/RevModPhys.92.015003
  * [3] Marco Cerezo et al.  “Variational quantum algorithms”  In _Nature Reviews Physics_ 3.9 Nature Publishing Group UK London, 2021, pp. 625–644  DOI: 10.1038/s42254-021-00348-9
  * [4] Alberto Peruzzo et al.  “A variational eigenvalue solver on a photonic quantum processor”  In _Nature communications_ 5.1 Nature Publishing Group, 2014, pp. 1–7  DOI: 10.1038/ncomms5213
  * [5] Abhinav Kandala et al.  “Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets”  In _Nature_ 549.7671 Nature Publishing Group, 2017, pp. 242–246  DOI: 10.1038/nature23879
  * [6] Yunseong Nam et al.  “Ground-state energy estimation of the water molecule on a trapped-ion quantum computer”  In _npj Quantum Information_ 6.1 Nature Publishing Group UK London, 2020, pp. 33  DOI: 10.1038/s41534-020-0259-3
  * [7] Samson Wang et al.  “Noise-induced barren plateaus in variational quantum algorithms”  In _Nature communications_ 12.1 Nature Publishing Group, 2021, pp. 1–11  DOI: 10.1038/s41467-021-27045-6
  * [8] Daniel Stilck França and Raul Garcia-Patron  “Limitations of optimization algorithms on noisy quantum devices”  In _Nature Physics_ 17.11 Nature Publishing Group UK London, 2021, pp. 1221–1227  DOI: 10.1038/s41567-021-01356-3
  * [9] Lloyd C L Hollenberg  “Plaquette expansion in lattice Hamiltonian models”  In _Physical Review D_ 47.4 APS, 1993, pp. 1640  DOI: 10.1103/PhysRevD.47.1640
  * [10] Lloyd C L Hollenberg and N S Witte  “General nonperturbative estimate of the energy density of lattice Hamiltonians”  In _Physical Review D_ 50.5 APS, 1994, pp. 3382  DOI: 10.1103/PhysRevD.50.3382
  * [11] Lloyd C L Hollenberg and N S Witte  “Analytic solution for the ground-state energy of the extensive many-body problem”  In _Physical Review B_ 54.23 APS, 1996, pp. 16309  DOI: 10.1103/PhysRevB.54.16309
  * [12] Harish J Vallury et al.  “Quantum computed moments correction to variational estimates”  In _Quantum_ 4 Verein zur Förderung des Open Access Publizierens in den Quantenwissenschaften, 2020, pp. 373  DOI: 10.22331/q-2020-12-15-373
  * [13] Harish J. Vallury et al.  “Noise-robust ground state energy estimates from deep quantum circuits”  In _Quantum_ 7 Verein zur Förderung des Open Access Publizierens in den Quantenwissenschaften, 2023, pp. 1109  DOI: 10.22331/q-2023-09-11-1109
  * [14] Michael A Jones et al.  “Chemistry beyond the Hartree–Fock energy via quantum computed moments”  In _Scientific Reports_ 12.1 Nature Publishing Group, 2022, pp. 1–9  DOI: 10.1038/s41598-022-12324-z
  * [15] Karol Kowalski and Bo Peng  “Quantum simulations employing connected moments expansions”  In _The Journal of Chemical Physics_ 153.20 AIP Publishing LLC, 2020, pp. 201102  DOI: 10.1063/5.0030688
  * [16] Kazuhiro Seki and Seiji Yunoki  “Quantum power method by a superposition of time-evolved states”  In _PRX Quantum_ 2.1 APS, 2021, pp. 010333  DOI: 10.1103/PRXQuantum.2.010333
  * [17] Philippe Suchsland et al.  “Algorithmic error mitigation scheme for current quantum processors”  In _Quantum_ 5 Verein zur Förderung des Open Access Publizierens in den Quantenwissenschaften, 2021, pp. 492  DOI: 10.22331/q-2021-07-01-492
  * [18] Joseph C Aulicino, Trevor Keen and Bo Peng  “State preparation and evolution in quantum computing: A perspective from Hamiltonian moments”  In _International Journal of Quantum Chemistry_ 122.5 Wiley Online Library, 2022, pp. e26853  DOI: 10.1002/qua.26853
  * [19] Ruizhe Zhang, Guoming Wang and Peter Johnson  “Computing ground state properties with early fault-tolerant quantum computers”  In _Quantum_ 6 Verein zur Förderung des Open Access Publizierens in den Quantenwissenschaften, 2022, pp. 761  DOI: 10.22331/q-2022-07-11-761
  * [20] Lloyd C L Hollenberg and Michael J Tomlinson  “Staggered magnetisation in the Heisenberg antiferromagnet”  In _Australian journal of physics_ 47.2 CSIRO Publishing, 1994, pp. 137–144  DOI: 10.1071/PH940137
  * [21] N S Witte, Lloyd C L Hollenberg and Zheng Weihong  “Two-dimensional XXZ model ground-state properties using an analytic Lanczos expansion”  In _Physical Review B_ 55.16 APS, 1997, pp. 10412  DOI: 10.1103/PhysRevB.55.10412
  * [22] Stasja Stanisic et al.  “Observing ground-state properties of the Fermi-Hubbard model using a scalable algorithm on a quantum computer”  In _Nature Communications_ 13.1 Nature Publishing Group UK London, 2022, pp. 5743  DOI: 10.1038/s41467-022-33335-4
  * [23] Ken N Okada et al.  “Identification of topological phases using classically-optimized variational quantum eigensolver”  In _arXiv preprint arXiv:2202.02909_ , 2022  DOI: 10.48550/arXiv.2202.02909