# Energy control of HVAC units for provable ancillary service provision

Rupamathi Jaddivada, Marija D. Ilić, Mario Garcia-Sanz, Mirjana Marden, and
Sonja Glavaski R. Jaddivada and M.D. Ilic are with, Massachusetts Institute of
Technology, Cambridge, MA, 02139 USA e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>M.
Garcia-Sanz and M. Marden are with Advanced Research Project Agency (ARPA-E),
US Department of Energy, Washington, DC 20585 USA email:
<EMAIL_ADDRESS><EMAIL_ADDRESS>Glavaski is currently
with Pacific Northwest National Laboratory, 902 Battelle Boulevard, Richland,
WA USA<EMAIL_ADDRESS>
###### Abstract

In this paper, we consider the problem of controlling power consumption
dynamics of residential heating, ventilation and air conditioning (HVAC) units
so that they follow the grid-side power specifications. In order to do so, we
design a novel dynamical energy controller which ensures regulation of the
cumulative effects of power imbalances. For this, we derive a novel energy-
based model that relates the HVAC physics-based dynamics to both real and
reactive power balance at the point of interconnection with the grid. In
contrast to several other approaches in the literature, we show that a limited
number of HVAC units can meet the stringent performance metrics set by the
ARPA-E/NODES program on following the frequency regulation signal, while
maintaining consumer comfort. Theoretical and simulation-based model and
control validation is provided by making use of real-world HVAC data.

###### Index Terms:

Demand response, Load flexibility, Themostatically controlled loads,
Residential HVAC loads, Energy-based modeling and control

##  I Introduction

Renewable Portfolio Standards (RPS) being adopted by several states in the US
lead to large-scale deployment of solar panels and wind farms [1]. These
technologies are often not controlled by the grid operator and thus result in
hard-to-predict disturbances in the grid. Shown in Fig. 1 is the projected
hard-to-predict in near real-time net load change in California over the next
several years [2], leading to increased reliance on ancillary services.
Ancillary service includes frequency control, spinning reserves, and operating
reserves, all related to additional power adjustments needed to offset the
mismatch between the net demand and scheduled generation. Traditionally, large
conventional generation units have been utilized for providing ancillary
services. A major operational challenge is that these units have significant
inertia and thus incur increased wear and tear, as they follow the net load
patterns, such as the one indicated in Fig. 1.

Therefore, utilities are exploring the potential of the flexibility provided
by the residential consumer-end resources distributed throughout the grid,
such as the electric vehicles, batteries and thermostatically controlled loads
(TCLs) like electric water heaters, air conditioners and refrigeration units
[3, 4, 5, 6]. This paper’s particular interest is the flexibility offered by
the heating, ventilation, and air conditioning (HVAC) units. They comprise 40%
of the total residential power demand and thus present an enormous potential
for providing various types of ancillary services to the grid [7].
Furthermore, their inherent thermal storage can be leveraged to precisely
modify the power consumption and provide ancillary services while still
meeting the desired temperature requirements of the end-user.

Figure 1: The duck curve showing steep ramping needs and over generation of
wind farms in California [2]

Different control methods broadly categorized as direct and indirect load
control can enable HVACs to participate in ancillary services. In the direct
load control method, the utility or the aggregator models many HVAC units and
interrupts the power signals of several HVAC units to obtain desired aggregate
performance. This method requires accurate prediction of the thermal
flexibility of each HVAC unit to have a reliable supply of ancillary services
at an aggregate level [3, 8, 9]. Through the indirect load control method, the
consumer or the appliance automatically adjusts its switching cycles in
response to real-time electricity prices or frequency deviations in power
systems, thereby providing ancillary services indirectly [10, 11, 12]. The
relationship between the HVAC temperature setpoints and the real-time pricing
is often complicated. Thus indirect load control method results in deviations
from the power schedules, especially when the cluster of loads is not diverse
enough or if a significant fraction of them is operating at their limits [13,
14]. Some other control methods combining these different methods include
passivity-based distributed optimization [15, 16, 17], multi-stage robust
optimization considering detailed HVAC models [18], two layered-control
distributing the decision-making between local and global controllers [19].
However, these approaches do not consider the practical limitations of the
controllers or the availability of sensor measurements. For example, the HVAC
needs to remain in the OFF state for a specific length of time before
switching it back ON. These are called short-cycling constraints, which have
been considered to some extent in [8, 13, 20]. More importantly, to allow for
the participation of a cluster of loads, the NODES program of ARPA-E has
indicated the performance metrics in terms of the response time, ramp time,
and availability of units explained in Section 2 [21]. However, most methods
proposed in the literature can meet these metrics only upon considering a
large number of units that are diverse enough [13, 14, 22]. Thus in this
paper, we consider effects of a limited number of HVAC units and assess their
capability in following the scaled-down regulation signal, all while ensuring
provability. By provability we mean that the approach needs to result in
reproducible results with theoretical performance guarantees.

In Section II, we review conventional HVAC models, and identify sub-objectives
needed to pose the control problem. In Section III, we propose a novel energy-
space model for the HVAC system building upon our recent work [23, 24, 25,
26]. Utilizing this model, we next summarize our proposed multi-layered
energy-based control of HVAC in Section IV. This control comprises a device-
level sliding mode control tracking an output variable in energy-space to a
reference value, explained in Section V. Novel closed-loop droop relations
mapping the controlled output, reference signals and the regulation signals
have been derived next in Section VI. These are further utilized in a slower
time-scale MPC-based controller to maximize the efficiency and meet the
regulation signals. We finally conclude by summarizing contributions of this
paper and suggest directions for future research in Section VII.

##  II Problem formulation

In this section, we review residential HVAC modeling along with numerous
constraints imposed by the consumer, HVAC manufacturer and the utility
operator for it to qualify as a reliable source of regulation.

###  II-A Conventional models of residential HVAC system

HVAC system comprises the space or zone to be heated/cooled and the auxiliary
electrical equipment that injects hot/cool air into the space. As a result, we
can characterize two types of state variables, as shown in Fig. 2.

Figure 2: Schematic of a HVAC system in open loop: The space is modeled using
temperatures of multiple zones denoted using vector $x_{T}$. There are
additional auxillary states $x_{a}$ corresponding to the fan and
heater/cooling unit. The blower fan switch and the cooler/heater switch
$u_{c}$ are the control inputs

The first set of variables characterize the temperatures of multiple zones,
representing the different rooms, walls and are denoted as $x_{T}$ [27]. The
electrical equipment comprises a fan that circulates air into the zones, while
the heater/cooling unit maintains the temperature of the supply air. Let us
denote the respective state variables using $x_{a}$. The models of these sub-
systems are not exactly known.

Objective 1: The control design should be robust to model and parameter
uncertainties and be simple enough to utilize limited available sensor
measurements.

In residential HVAC systems, often single zone temperature is utilized to
characterize the thermal load [20]. Furthermore, the power rating of the fan
is orders of magnitude smaller than that of the compression system or arc
furnace for cooling or heating, respectively. Therefore, in the rest of the
paper, we only consider single control input, which is the switch position of
the cooling unit, denoted as $u$.

$\displaystyle C\dot{T}$ $\displaystyle=-\frac{1}{R}(T-T_{0})+{P^{rated}u}$
(1a) $\displaystyle P^{rated}$ $\displaystyle=\dot{m}_{a}C_{p}(T_{sup}-T)$
(1b)

$T$, $T_{0}$ and $T_{sup}$ respectively denote the zonal temperature, ambient
temperature and the supply temperature that is maintained by the compression
system in heater/cooler block. $C,R$ represent the thermal capacitance and
resistance, respectively. $P^{rated}u$ is the heating or cooling rate for with
$u$ being the switch position taking a positive or negative value
respectively. $\dot{m}_{a}$ is the airflow rate, $C_{p}$ is the specific heat
capacity of the air. $T_{sup},\dot{m}_{a}$ form the auxiliary state variables
in Fig. 2 and are typically assumed to be constant.

Objective 2: It is required to ensure that the temperature reaches the
setpoint $T^{ref}$ as set by the consumer with a tolerable temperature
deviations of $T_{db}$, i.e.
$T(t)\in\left[T^{ref}-T_{db},T^{ref}+T_{db}\right]\forall t$

The pre-programmed embedded automation $u_{em}$ in today’s HVACs is based on
PID controllers and satisfies Objective 2 under all circumstances. In context
of Eqn. (1a), $u_{em}$ is interpreted as a bang-bang control, with the
following logic assuming the heating mode of operation.

$\displaystyle u=u_{em}(t)$ $\displaystyle=1\quad\text{if}\qquad
T(t)<T^{ref}-T_{db}$ (2a) $\displaystyle=0\quad\text{if}\qquad
T(t)>T^{ref}+T_{db}$ (2b) $\displaystyle=u_{em}(t-\delta
t)\qquad\text{otherwise}$ (2c)

Here, $\delta t$ is the continuous time implementation timestep which is
dictated by the bandwidth of the available sensor measurements.

Objective 3: Several considerations from the perspective of the control
designer and manufacturer respectively must be met as follows:

  * •

The control needs to be implemented using the knowledge of temperature and
electrical power consumption alone. This is because sensors for other internal
variables are unavailable and are expensive.

  * •

The compressor needs to be left ON at least for 5 minutes and it can remain in
that state at most for 15 minutes.

###  II-B Ancillary services

The aggregator or service provider acts as the middle man between the
individual HVAC units and the grid operator. The aggregator’s responsibility
is to provide the regulation service as commanded by the grid operator over
minutes timescale as given by the discrete time samples $k$ shown in Figure 3.
However, there are thermal and electric power dynamic processes evolving at
finer time granularity of seconds and milliseconds, respectively. The
controller actions need to be taken at these finer timescales to have their
accumulated effect over minutes timescale track the commanded regulation
signal.

Figure 3: Multi-timescale evolution of control objectives: Aggregate
regulation signals arrive every $kT_{t}$ instant; Individual HVAC unit
commands need to arrive every $nT_{s}$ step while the power consumption
dynamics evolve at a much faster continuous millisecond timescale $(t)$

Objective 4: From the grid perspective, if single HVAC unit were to
participate in ancillary services, it should consume average power every
$kT_{t}$ instant equal to the commanded regulation signal $P^{reg}(kT_{t})$

$\sum_{t=(k-1)T_{t}}^{kT_{t}}{P(t)}={P^{reg}}(k{T_{t}})$ (3)

###  II-C ARPA-E performance metrics

The utility performs a few tests to assess if the units willing to participate
in ancillary services meet specific performance metrics before being deemed a
reliable provider. There is no clear rationale for setting these performance
metrics. These metrics are not standard even for current grid operations where
controllable generation units primarily provide ancillary services. Following
the same principles, ARPA-E has proposed performance metrics for different
types of ancillary services as needed to be met by a group of small-capacity
flexible units coordinated by and collectively referred to as Network
Optimized Distributed Energy System (NODES)[21]. We interpret these metrics in
the context of a single HVAC unit of interest in this paper.

Objective 5: The averaged power consumption over $T_{t}$ rate $P[k]$ of each
HVAC unit should observe the following specifications as it is commanded to
follow step change in regulation signal at time $t=0$, as shown in Fig. 4.

Figure 4: Performance metrics to be satisfied for participating in ancillary
services [21]

  * •

Response time ($t_{\text{re}}$): The time it takes for the actuators to
respond to the commanded signals - less than 5 seconds.

  * •

Reserve magnitude target (RMT): The maximum step change in regulation signal
that the HVAC unit can track - 7% of rated capacity 111It is stated 7% of the
peak inflexible demand in the original specifications.

  * •

Reserve magnitude variability tolerance (RMVT): Maximum deviation between the
actual power adjustments and the reserve magnitude target - 5% of ’RMT’

  * •

Ramp time ($t_{\text{ra}}$): The maximum time it takes for the HVAC unit to
reach the ’RMT’ and stay within ’RMVT’ deviations - 5 seconds

  * •

Availability time ($t_{\text{a}}$): The period of time for which the reserves
can be supplied - 3 hours

We ensure satisfaction of these performance metrics by explicitly considering
some of these metrics in the control design as will be explained in Section
IV.

##  III Energy-based modeling of HVAC

###  III-A Motivation for new modeling tools

The thermal model in Eqn. (1a) is most commonly used in the literature. This
equation is based on the first law of energy conservation principle. However,
control designed based on the first law alone suffers major inconsistencies.
For instance, the power absorbed by the HVAC unit is given by

$P^{r,out}=\frac{1}{R}\left(T-T_{0}\right)+C\frac{dT}{dt}$ (4)

The first component is the useful work done by the HVAC unit and the second
term is the rate of change of stored energy. On the other hand, the injected
power into the HVAC unit for the embedded switching logic in Eqn. (2) is equal
to

$P^{r,in}=P^{rated}u_{em}$ (5)

If the control were to satisfy the energy conservation principles, the two
values of power interactions $P^{r,out}$ and $P^{r,in}$ should be equal.
However, this is not the case as evidenced by the the time trajectory plots in
Figure 5. It may appear that this imbalance stems from the discrete control
actions. However, the accumulated energy imbalance over time also doesn’t
average out to zero as shown in Figure 6.

Figure 5: The time trajectories of instantaneous power absorption and power
injections seen by the HVAC at the electrical terminal Figure 6: Accumulated
energy imbalance at the electrical terminal
$\int_{0}^{t}\left(P^{r,out}-P^{r,in}\right)dt$

This inconsistency stems fundamentally because the controllers based on first
law of energy conservation principle alone do not align the rate at which
power evolves thereby leaving a power imbalance offset that creates
accumulated energy imbalance over time. We thus introduce a novel energy
modeling approach to explicitly accommodate the rate of change of power
variables and control its cummulative effects on temperature and comfort.

###  III-B Energy modeling

The schematic of energy interaction for the detailed HVAC model is shown in
Fig. 7. The compressor system can operate either in heating/cooling mode that
dictates the supply temperature $T_{s}$ into the zone, used in the simplified
model in Eqn. (1a). The Air Handling Unit (AHU) injects air at desired
temperature into the zone, which circulates back through the fan. The zone is
associated with thermal energy $U$ while the AHU unit is associated with
stored energy $E_{a}$. The AHU comprises a blower fan, a and a compression
unit or chiller. These can both be individually controlled. We assume that
these units have internal automation. For enabling the HVAC units’
participation in ancillary services, we incorporate an additional switch shown
in Fig. 7 that modulates the rate at which power enters the HVAC unit. It is
important to note that we do not alter the existing power consumption as
needed by the chiller and fan while controlling the rate of reactive power
injections of the switching control action.

Figure 7: Schematic of HVAC system in energy space. The overall system is
associated with stored energy $E$ and its rate of change $p$ which can be
controlled by the switch shown in green. The switch itself injects non-zero
$\dot{Q}_{u}$, resulting in net absorbed HVAC power variables $P^{r,out}$,
$\dot{Q}^{r,out}$, while the grid power injections at instantaneous time are
$P^{r,in},\dot{Q}^{r,in}$

All the interactions between subsystems can get overly complex if modeled in
conventional state space. We instead assess the overall energy conversion
dynamics, and strive to design the control without exactly knowing the
internal physical models.

###  III-C Preliminary definitions of variables in energy space

In any energy domain, power variable are defined using the port variables
called effort flow variables. These pairs are (voltage, current) in electric
domain, (pressure, volume flow rate) in fluid domain and similarly based on
Eqn. (6), they are (Temperature, entropy flow) in thermal domain. The time
integral of this quantity is defined as the energy injected into the
component. We start by defining thermal energy.

######  Definition 1.

The internal or thermal energy of a thermodynamic fluid (air in the case of
HVAC) can be expressed as follows:

$U=\int_{0}^{t}TdS=\int_{0}^{t}T\frac{dS}{dt}dt=\int_{0}^{t}TS_{f}dt$ (6)

Here, $dS$ is the incremental entropy of the medium, $S_{f}=\dot{S}$ is the
net entropy flow into finite volume at a temperature $T$.

######  Definition 2.

The energy in tangent space is defined by replacing the effort-flow variables
in definition of stored energy with that of its derivatives.

$U_{t}=\int_{0}^{t}\frac{dS_{f}}{dt}\frac{dT}{dt}dt$ (7)

Irrespective of energy domain, each arrow in Fig. 7 models the power variables
defined next.

######  Definition 3.

Given the generalized effort $e$ and flow variables $f$ at a port, the
instantaneous power and rate of change of generalized reactive power absorbed
by the port is defined as in Eqn. (8) [24, 28].

$\displaystyle P$ $\displaystyle=ef$ (8a) $\displaystyle\dot{Q}$
$\displaystyle=e\frac{df}{dt}-\frac{de}{dt}f$ (8b)

We have defined here the thermal energy $U$ and thermal energy in tangent
space $U_{t}$. The similar characterization for electro-mechanical energy
domains is explained in [23, 24] and other energy domains in [25]. The total
energy of the HVAC system $E$ is the sum of the stored energy of the sub-
systems in AHU $E_{a}$ and thermal energy $U$. We have shown in [24] that the
dynamics of energy exchanges of a component with the rest of the system can be
captured by modeling the dynamics of aggregate variables energy $E$ and its
rate of change $p=\dot{E}$. The model was derived for electrical circuits and
was proven to hold for complex electromechanical systems by extending the
definitions of energy space variables through effort-flow analogy for multiple
energy domains. In the same way, extending the analogy to the hydraulic energy
domain with a continuum of space, we have derived models for turbo-machinery
in this energy space in [25]. We now extend the modeling approach to also
capture thermal interactions within the HVAC system.

###  III-D Stand-alone HVAC model in energy space

We now utilize the definitions of energy space variables to show that the
previously proposed energy space model in [24] also holds for thermal
interactions of HVAC.

######  Proposition 1.

The stand-alone HVAC component model shown in Figure 7 with an added switch at
the interface exhibits the following general dynamics in energy state space
comprising the states $\left[E,p\right]$

$\displaystyle\dot{E}$ $\displaystyle=-\frac{E}{\tau}+P^{r,out}=p$ (9a)
$\displaystyle\dot{p}$ $\displaystyle=4E_{t}-\dot{Q}_{u}-\dot{Q}^{r,out}$ (9b)

Here, $\frac{E}{\tau}$ represents the total damping losses. $E_{t}$ is the
total stored energy in tangent space. $\dot{Q}_{u}$ is the controllable
reactive power injected by the switch. The resulting internal state dynamics
result in instantaneous power and rate of reactive power absorption
$P^{r,out}$ and $\dot{Q}^{r,out}$ respectively.

###### Proof.

Starting from the model in Eqn. (1a), we show that the aggregate energy space
model holds in Appendix A. ∎

######  Remark 1.

Here ${p}=P^{r,out}$ for when the damping losses $\frac{E}{\tau}=0$. In such
cases $p$ is the supplied power. If there were no switching reactive power
injected, $\int_{0}^{t}4E_{t}dt$ is the power corresponding to the potential
to do useful work while the ${Q}^{r,out}$ represents component of power
corresponding to wasted work [24]. This wasted works is resulting from the
dynamic inefficiencies associated with the energy conversion processes at
instantaneous time. It must be noted that this inefficiency is different from
the damping losses from linear frictions, viscosity etc, which is captured in
the term $\frac{E}{\tau}$. This dynamic inefficiency arises due to the
difference in the rate at which thermal power gets absorbed by the zone and
the rate at which the power is absorbed from the grid. It is only when
$2\dot{Q}^{r,out}$ approaches $\dot{Q}_{T}$, maximum physical efficiency is
achieved.

###  III-E Validation of energy model using real data

In order to keep the analytical derivations simple, we have assumed in
Proposition 1 negligible stored energy in AHU when compared to that of the
zone. However, the energy space model postulated holds even if this assumption
were relaxed. To this end, we utilize real-world measurements of a detailed
HVAC system to validate our proposed energy model.

We performed a model validation experiment in [29] where the interface
variables of chiller and pumps in compression system and that of fans and
zones were used to conduct parameter estimation of the respective conventional
state space models. These parameters are the inertial and damping coefficients
of each of the components within a HVAC system. These parameters were used to
compute the overall stored energy $E$, damping $\frac{E}{\tau}$ and the stored
energy in tangent space, given the measurements of internal states
(temperature of zone and chiller in heater/cooler units, air flow rate in the
fan). The detailed definition of these energy variables is outside the scope
of this paper. The reader is referred to [29]. We further utilized the
measurements of voltages and currents at electrical terminals of compressors
and fans to obtain instantaneous and rate of reactive power using Eqn. (8).

Shown in Figure 8 is the perfect match in the comparison plots of left hand
side and right hand side of the proposed energy model of Eqn. (9) over longer
timescales. Shown in Figure 9 are the same plots zoomed-in over finer time
granularity.

(a) Zoomed-out Eqn. (9a)

(b) Zoomed-out Eqn. (9b)

Figure 8: Zoomed-out validation of unified energy model using real-world data

(a) Zoomed-in Eqn. (9a)

(b) Zoomed-in Eqn. (9b)

Figure 9: Zoomed-in validation of unified energy model using real-world data

###  III-F Interactive interconnected model

There exists a unique transformation of state and state derivatives called
interaction variable $z^{r,out}$ that has important structural properties. The
rate of change of interaction variable is zero when disconnected from the
grid. For details, see [30]. These outgoing interaction variables are defined
as:

$\dot{z}^{r,out}=\left[\begin{array}[]{l}{P^{r,out}}\\\
{{\dot{Q}}^{r,out}}\end{array}\right]=\left[\begin{array}[]{l}p+\frac{E}{\tau}\\\
4{E_{t}}-\dot{Q}_{u}-\dot{p}\end{array}\right]=\phi(\tilde{x},\dot{\tilde{x}})$
(10)

The definitions for $p$, $E_{t}$, $\dot{Q}_{T}$ for the simple dynamical model
in Eqn. (1a) are introduced in Appendix A. As a result the outgoing
interaction variable is abstracted as a function of extended states and its
derivatives which is defined as $\tilde{x}=\left[x,u\right]^{T}$. Here $x$ is
a vector of all HVAC state variables and $u$ is the switch control introduced
for ancillary service participation in Figure 7.

On the other hand, there is also incoming interactions into the HVAC from the
electric grid due to the energy conversion dynamics of rest of the components
connected to the grid i.e.

$\dot{z}^{r,in}=[P^{r,in},\dot{Q}^{r,in}]^{T}$ (11)

$\dot{z}^{r,out}$ is created through thermal processes, while $\dot{z}^{r,in}$
is the result of faster electric processes on the grid side. These two
quantities need to align with each other for system interconnection to be
stable and feasible [23, 30]. During transients this difference is not zero as
shown for the instantaneous power component with conventional control design
in Figure 5. Driving this imbalance to zero is set as the primary control
objective as will be explained in the next section.

Notably, the energy space model in Eqn. (9) is linear in terms of the energy
states $x_{z}=\left[E,p\right]^{T}$ where the control input in energy state
space is $u_{z}=\dot{Q}_{u}$ and the term
$\left(4E_{t}-\dot{Q}^{r,out}\right)$ is the state-dependent disturbance
entering the model. Setting bounds on this disturbance allows us to perform
linear control design in energy state space with provable design.

##  IV Proposed multi-layered energy-based control

We exploit the linear nature of energy model to first design a provable
primary control that addresses all the stated objectives. Shown in Figure 10
is the proposed multi-layered control architecture based on the timescales
illustrated in Figure 3.

Figure 10: Proposed multi-layered information exchange for control of HVAC.
Interface power measurements are utilized for sliding mode energy control at
continuous timescales to track a reference value $y_{z}^{ref}[n]$ obtained
over discrete $T_{s}$ timesteps such that the resulting HVAC power adjustments
could track a given regulation signals over even longer $T_{t}$ timescale.

In particular, at fastest timescales, we design a sliding mode control in
energy state space that relies on interface power measurements alone. The
sliding mode design inherently allows for the design to be robust to model and
parameter uncertainties. We further introduce an error threshold for the
sliding mode design to explicitly incorporate the compressor cycling
constraints. We can thus satisfy objectives 1 and 3.

In context of the timescales illustrated in Figure 3, we implement the the
sliding mode design to have a carefully selected output variable in energy
state space $y_{z}(t)$ track a consistence reference value $y_{z}^{ref}[n]$
varying at minutes timescale. Upon application of such design, it is possible
to characterize a quasi-static droop relations between output of interest in
energy state space, its reference value and the averaged power absorbed by the
HVAC unit over minutes timescale. This reference value is the secondary
control action designed in an interactive manner by accounting for the closed
loop HVAC dynamics with the primary control and the thermal comfort
constraints. As a result, we satisfy objective 2 over seconds timescale while
also satisfying objective 4 over minutes timescales.

The sliding mode design over fastest timescales together with model predictive
control over slower timescales allow us to ensure stability, regulations and
incorporation of ARPA-E performance metrics stated in Objective 5 in a
provable manner. Meeting the performance metrics, especially the RMT depends
on how large the target is, the thermal load of the HVAC unit and the number
of HVAC units available. With our proposed approach, it is possible to atleast
identify the extent to which the reserve magnitude target can be met apriori
while satisfying all other performance metrics. This allows for bringing in
additional HVAC units to thereby satisfy the reserve magnitude target as
needed. Each of these steps will be explained in more detail in the subsequent
sections.

##  V Provable primary control

###  V-A Control design in energy state space

HVAC systems already have an embedded proprietary digital control. We propose
to overwrite these signals through a robust primary control based on sliding
mode energy control design. Consider the energy space model in Eqn. (9) and
let $\dot{Q}_{u}$ be the available degree of control and let us select the
output variable of interest as in Eqn. (12a).

$y_{z}=P^{r,out}\left(E,p=0\right)=\frac{E}{\tau}=\frac{1}{R}\left(T-T_{0}\right)$
(12a) This quantity represents the device’s instantaneous power absorption
after the settling of stored energy dynamics. We set the reference value as
follows: $\displaystyle y_{z}^{ref}=P^{r,in,ref}=P^{r,out,ref}+\Delta
P^{r,out}[n]$ $\displaystyle=\frac{1}{R}\left(T^{ref}-T_{0}\right)+\Delta
P^{r,out}[n]\forall t\in[(n-1)T_{s},(nT_{s})]$ (12b) Here, the first term is
the desired baseline power consumption to satisfy the thermal comfort
requirements while the second term indicates the power consumption adjustments
for participation in ancillary services. $y_{z}^{ref}[n]$ is the feed-forward
signal that enters every $T_{s}$ timestep. The fast sliding mode design is
expected to ensure $y_{z}\rightarrow y_{z}^{ref}$ before $t=nT_{s}$. To this
end, we design a robust sliding mode control as in Eqn. (12d) that ensures
finite settling time. $\dot{Q}_{u}=\alpha
sign\left(\underbrace{y_{z}-y_{z}^{ref}}_{\sigma}\right)$ (12d) Here,
$\sigma=y_{z}-y_{z}^{ref}$ represents the sliding surface, and $\alpha$ is the
sliding mode gain.

###  V-B Physical control mapping

Let us now zoom-in to the switch in Figure 7 and characterize the reactive
power flows as shown in Figure 11. Let us denote left hand side and right-hand
side reactive power flows as $\dot{Q}_{sw,L},\dot{Q}_{sw,R}$. The switch on
left hand side interacts with the electrical grid while that on right side
interacts with the AHU and zone of the HVAC component.

Figure 11: Depiction of reactive power flows at the interfaces of the switch
considered in Figure 7.

Given the electrical effort flow variables $\left(v,i\right)$ modulated
through the switch $u$ at left and right hand side ports of the switch, we can
utilize the definition of reactive power in Eqn. 8 to obtain the following
expressions.

$\displaystyle\dot{Q}_{u}$ $\displaystyle=\dot{Q}_{sw,L}-\dot{Q}_{sw,R}$ (13a)
$\displaystyle=\left(v\frac{d}{dt}\left(ui\right)-\left(ui\right)\frac{dv}{dt}\right)-\left(\left(uv\right)\frac{di}{dt}-i\frac{d}{dt}\left(uv\right)\right)$
(13b) $\displaystyle=2vi\dot{u}=2P^{rated}\dot{u}$ (13c) Here, $P^{rated}$ is
the power flowing through the switch which is the same as the one defined in
Eqn. (1a).

Rearranging the terms in Eqn. (13), we obtain the dynamical switching control
expression in Eqn. (14).

$\dot{u}=\frac{1}{2P^{rated}}\dot{Q}_{u}$ (14)

######  Theorem 1.

The closed loop model (Eqn. (1a), (14)) exhibits following properties

  1. 1.

It is diffeomorphic to the energy state space model in Eqn. (9)

  2. 2.

The virtual control design in Eqn. (12d) ensures finite settling time if the
sliding mode gain is set equal to $\alpha=\overline{L}+K$ where
$\left|4E_{t}-\dot{Q}^{r,out}\right|\leq\overline{L}$, $K>0$

  3. 3.

The objective of $y_{z}\rightarrow y_{z}^{ref}$ is achieved within a finite
reaching time $t_{r}\leq\frac{K}{2|\sigma(0)|}$ where $\left|\sigma(0)\right|$
represents the distance of operating point from the sliding surface at initial
time.

###### Proof.

The proof is provided in the Appendix B. ∎

The finite reaching time result of the sliding mode control is an important
result. In particular, the sliding mode gain can be chosen to ensure the
reaching time $t_{r}$ is at least ten times smaller than secondary control
timescale $T_{s}$. We can further select $T_{s}$ to be equal to the ramp time
$t_{ra}$ stated in Objective 5.

For chosen values of $T^{ref}=75^{0}$ F, $\Delta P^{r,out}[n]=0$ and time
varying ambient temperature changes, we obtain finite settling time of the
temperature to the reference value as shown in Figure 12. The simulation here
has utilized a more detailed residential HVAC simulation model with AHU
dynamics simulated and with imperfect parameter corresponding to the damping
losses being utilized. Notice that the the energy control design utilizes
measurements of just the interface electrical power measurements and the zonal
temperature values which are typically available in the existing residential
smart meters and thermostats.

Figure 12: Temperature trajectories obtained with energy control applied to a
detailed HVAC simulation model

###  V-C Implementation-aware design

The energy control in Eqn. (12d) together with the design implementation in
Eqn. (14) may lead to frequent switching which may lead to violation of the
compressor cycling constraints. In order to accommodate the stated objective 3
with the compressor ON-time constraint of
$\tau^{min}\leq\tau_{ON}\leq\tau^{max}$, we incorporate a small threshold for
the sliding surface leading to the design as

$\begin{array}[]{l}\dot{Q}_{u}~{}=\alpha\qquad y_{z}\geq
y_{z}^{ref}+y_{z}^{+}\left(\tau^{max}\right)\\\ \quad\quad=-\alpha\qquad
y_{z}\leq y_{z}^{ref}+y_{z}^{-}\left(\tau^{min}\right)\\\
\qquad=0\quad\qquad{\rm{otherwise}}\end{array}$ (15)

The thresholds $y_{z}^{+}$ and $y_{z}^{-}$ are computed so that the
temperature excursion is between tolerable temperature bounds
$\left[T^{min},T^{max}\right]$ as permitted by the compressor cycling
constraints. The ON time based on the thermal dynamics are derived to be
related to the minimum and maximum temperature variations in [14] as

$\tau_{ON}=RClog\left(\frac{T^{max}-T_{0}+Ry_{z}}{T^{min}-T_{0}+Ry_{z}}\right)$
(16)

Rearranging terms and assuming $\frac{\tau_{ON}}{RC}<<1$ we obtain

$y_{z}=\frac{1}{R}\left(T^{min}-T0\right)\left(\frac{\tau_{ON}}{{RC}}\right)-\frac{1}{R}\left(T^{max}-T^{min}\right)$
(17)

With the limits of $\tau_{ON}$ given, we can obtain the permissible bounds on
$y_{z}$ as

$\displaystyle y_{z}^{+}$
$\displaystyle=\frac{1}{R}\left(T^{min}-T0\right)\left(\frac{\tau^{max}}{{RC}}\right)-\frac{1}{R}\left(T^{max}-T^{min}\right)$
(18a) $\displaystyle y_{z}^{-}$
$\displaystyle=\frac{1}{R}\left(T^{min}-T0\right)\left(\frac{\tau^{min}}{{RC}}\right)-\frac{1}{R}\left(T^{max}-T^{min}\right)$
(18b)

The responses for the same ambient temperature in Figure 12 with the sliding
mode energy control of Eqn. (15) with the thresholds defined in Eqn. (18) with
$\tau^{min}=5$ minutes and $\tau^{max}=15$ minutes is shown in Figure 13. If
the values of $R$ abd $C$ are not exactly known, we can further utilize
conservative bounds in Eqn. (18). The implementation with (15) on an average
corresponds to an equivalent continuous time signal of virtual control
$\dot{Q}_{u}$ in Eqn. (12d) over longer timescales.

Figure 13: Temperature trajectories obtained with implementation-aware energy
control applied to a detailed HVAC simulation model with compressor cycling
constraints considered

###  V-D Step-response for non-zero adjustment signal

Next, we assess if a non-zero step signal $\Delta P^{r,out}[n]=+0.2kW$ used in
Eqn. (LABEL:Eqn:yzRef) can be chased by our proposed sliding mode energy
control. We activate this step signal at $t=0.5$ hours.

The resulting plots are shown in Fig. 14(a) with perfect tracking performance
achieved. Notice that the response time of less than 5 seconds is achieved
because of the implementation timestep chosen as $0.001$ seconds. The ramp
time is clearly less than 5 minutes, which can further be improved with
selection of sliding mode gain $\alpha$ utilized in Eqn. (12d). The reserve
magnitude variability tolerance of less than 5% of the target is also achieved
as is evident from the figure. Furthermore, the duration of $30$ minutes of
reserves availability has also been validated.

The electrical power is also shown in Fig. 14(b) Here, the base signal
utilized for comparison is the one when the embedded automation does not
respond to the regulation signal, i.e., the step-change in $y_{z}^{ref}[n]$)
at 0.5 hours. After $0.5$ hours, notice that there is a difference of
approximately $0.2$ kW as required by the regulation reserve signal.

(a) Output of interest $y_{z}$

(b) Averaged power consumption of HVAC with and without demand response

(c) Temperature $T$ of the zone $y_{z}$

(d) Switch position $u$

Figure 14: HVAC primary control tracking output in energy space $y_{z}$, to
the value $y_{z}^{ref}$ provided by the secondary layer control

The temperature is also ensured to be within the permissible range of
$60-65^{0}F$ as shown in Figure 14(c). Finally, the imposition of the deadband
in Eqn. (15) of the control also ensures that the switching actions are not
too frequent, as shown in Figure 14(d).

We have thereby reassessed objective 5 in context of the unit testing of a
single device for step change in regulation signal. Here we have assumed that
consistent reference signal $y_{z}^{ref}$ is provided so that the comfort
metrics are not violated.

In order to explicitly accommodate the thermal comfort constraints, we pose a
secondary control problem that dictates the adjustment signal $\Delta
P^{r,out}[n]$ over faster timescale so that the step change in regulation
signal of $0.2kW$ can be better tracked over longer time horizons. In order to
ensure that, we pose a secondary control problem formulation over relatively
slower timescales.

##  VI  Efficient secondary layer control

With the primary control explained in Section V, we have ensured that the
tracking of $y_{z}$ to $y_{z}^{ref}$ happens within finite time. Once the
primary dynamics settle, we can derive a quasi-static droop relation as stated
next.

######  Proposition 2.

Given the closed loop model in Eqn. (12d) there exists a three-way incremental
droop relation between the output in energy state space $\Delta y_{z}[n]$, its
reference value $\Delta y_{z}^{ref}[n]$ and the the power consumption $\Delta
P^{r,out}[n]$ as shown in Eqn. (19) where $\sigma$ is called the droop
constant which is primary control dependent.

$\Delta y_{z}[n]=\left(1-\sigma\right)\Delta y_{z}^{ref}[n]-\sigma\Delta
P^{r,out}[n]$ (19)

These droop relations are utilized to solve an MPC problem posed over a
planning time horizon of 1 minute to ensure accumulated power consumption
adjustments $\Delta P^{r,out}[n]$ tack a regulation signal $P^{reg}[k]$
arriving every minute. The secondary layer control problem is thereby posed as
follows:

$\displaystyle\mathop{\min}\limits_{\Delta{y_{z}}^{ref}[n]}\left(\sum\limits_{n{T_{s}}=(k{T_{t}})}^{{n}{T_{s}}=(k+1){T_{t}}}{{\mu^{e}}{{P^{r,out}}[n]}}\right)+\qquad\qquad\qquad$
$\displaystyle{{\mu^{reg}}\left|\left(\sum\limits_{n{T_{S}}=(k{T_{t}})}^{{n}{T_{s}}=(k+1){T_{t}}}{{\Delta{P^{r,out}}[n]}}\right)-{P^{reg}}[k]\right|}$
(20a)
$\displaystyle{s.t.}\qquad\Delta{y_{z}}[n]=(1-\sigma)\Delta{y_{z}}^{ref}[n]-\sigma\Delta{P}^{r,out}[n]$
(20b)
$\displaystyle\qquad\qquad\qquad\qquad\qquad{y_{z}}\left(kT_{t}\right)={y_{z0}}$
$\displaystyle\begin{array}[]{l}\frac{1}{R}\left(T^{ref}-T^{db}-T_{0}\right)\leq{y_{z}}[n]\\\
\qquad\qquad\qquad\leq\frac{1}{R}\left(T^{ref}+T^{db}-T_{0}\right)\\\
\qquad\forall nT_{s}\in\left(kT_{T},(k+1)T_{t}\right)\end{array}$ (20f)

In the droop model in Eqn. (20b), we assume the droop $\sigma$ although is
operating conditions dependent, is changing much slower than $T_{s}$. The
droop equation forms the secondary control discrete time model with state
$y_{z}$, control $\Delta y_{z}^{ref}$ and the free variable $\Delta
P^{r,out}[n]$. $\mu^{reg}$ is the penalty of not following the regulation
signal, while $\mu^{e}$ is the fixed energy cost being paid by the device. As
a result, the objective of secondary control is to optimize the trade-offs
between energy consumption and provably supply of reserves. In addition we
incorporate thermal comfort thruogh the limits on $y_{z}$ in Eqn. (20f)
stemming from the definition of $y_{z}$ in Eqn. (12a) and to incorporate the
temperature bounds stated in Objective 2. The number of switching cycles can
also be accommodated as an additional constraint based on the bounds on
$y_{z}$ utilized for use in primary control in Eqn. (18). The result of this
optimization is the sequence of reference signals $y_{z}^{ref}[n]$ to which
the primary control in Eqns. (14), (15) responds to.

###  VI-A Case of single HVAC

The interactive primary and secondary controllers are simulated for a house
with HVAC with permissible temperature of with secondary control timestep of
$T_{s}=5$ minutes and a horizon length of 1 step to obtain secondary control
actions as the regulation signal arrives. The assumed costs are
$\mu^{reg}=100\$/kW$ and $\mu^{e}=10\$/kWh$. The trajectories obtained with
this approach is compared with that when the horizon length is $12$ steps,
amounting to 1 hour of planning.

The overlaid plots of temperature are for the cases with and without MPC are
shown in Fig. 16. Notice that the temperature crosses the limits at few time
instants because of the constraint on $y_{z}$ with upper and lower bounds
serving as proxies to temperature constraints were softened by plugging them
in objective function with a penalty factor of $1e2$.

Notice that the secondary control MPC problem computes the reference signals
that would result in the temperatures to be within pre-specified limits. As a
result, the primary control implementation would not lead to saturation.

(a) With MPC

(b) Without MPC

Figure 15: Tracking of regulation signal with interactive primary and
secondary control of HVAC system Figure 16: Temperature trajectories for the
cases with and without MPC-based secondary control interacting with primary
control of HVAC system $y_{z}$.

Fig. 15(a) and 15(b) respectively show the result of the control action for
the case with and without MPC respectively . In these plots, the primary
control action results in fast-changing values of $y_{z}$ shown in blue
perfectly chasing the reference signal $y_{z}^{ref}[n]$ in red computed every
$5$ minutes. Overlaid is also the regulation signal, which is approximately
same as the secondary control action $y_{z}^{ref}[n]$ since $\Delta y_{z}[n]$
is ensured to remain close to zero to minimize temperature deviations. Notice
from Figure 15(a) that the case with MPC results in perfect tracking of the
regulation signal while the case without MPC emphasizes more on minimizing
objectives at current time instant without any foresight of future temperature
trajectories.

###  VI-B Case of 50 HVACs

In the previous subsection, we softened the temperature constraint to prevent
infeasibility of the posed secondary control problem in Eqn. (20). However,
this results in unacceptable temperature deviations. To make this a hard
constraint, it becomes important to consider more HVAC units.

We thus consider next the problem of tracking the regulation signal by a group
of 50 HVAC units with the interactive primary and secondary control as shown
in Figure 10. A simple way to accommodate multiple HVAC units for secondary
control is to duplicate the constraints in Eqn. (19) and (20f) for each HVAC
unit and replace the objective for single HVAC power consumption
$P^{r,out}[n]$ with the power consumption of all HVAC units i.e. summation of
$P^{r,out}[n]$ of all HVAC units.

We utilize real data taken from Mueller community of Pecan Street Inc (PSI) to
obtain a regulation signal representative of the disturbances created by
renewable as seen by the Austin node of the transmission grid modeled by
ERCOT, the independent transmission system operator in Texas. We further
utilize the parameters of 50 HVAC units for simulations and the temperature
preferences of the respective units from the survey data as was made available
by PSI to carry out this experiment.

Shown in Fig. 17 are the individual power consumption trajectories of randomly
selected HVAC units. 50 HVACs together adjust their power consumption from
their baseline values as shown in red in Fig.18 overlaid on top of the
commanded regulation signal in blue. Notice that the perfect tracking was
achievable even with a small number of HVAC units.

Figure 17: Power consumption of randomly selected HVAC units  Figure 18:
Aggregate performance of 50 HVACs tracking the reserve signal

##  VII Conclusion

The contributions of this paper are multi-fold. First, a novel energy-based
approach is introduced to facilitate model-free control design utilizing
minimal sensor measurement data for implementation. Using this modeling, a
robust sliding mode control of energy dynamics is derived, along with
theoretical conditions under which proposed control design results in a
provable three-way droop relation between output variable of interest,
secondary control signals and the desired electrical power adjustments over
slower time scales. It is then shown how a secondary layer model predictive
controller accommodates provable droop relations, performance metrics and
implementation constraints. Finally, simulation-based evidence is presented
for backing up claims on efficiency and provability of the proposed design
method using real data and practical consideration of a HVAC unit in one
community in Texas. In particular, it is shown that the controlled real power
consumed by the HVAC comprises two components: one that does real work and the
other corresponding to wasted work. The former compensates for the effect of
energy imbalances created while participating in the provision of regulation
reserves. The latter corresponds to the wasted energy required to compensate
for the effect of thermal consumption created due to the temperature set-point
adjustment for tracking the regulation reserve signal. Notably, this first-of-
its-kind relation between electrical and thermal processes is viewed as being
fundamental to having sufficient granularity of the model. Further work is
needed toward deploying these controllers and fully utilizing their potential.

##  Appendix A Proof of Proposition 1

###### Proof.

We first derive the expressions for energy variables in context of the
simplified dynamical models in use in Eqn. (1a). In Definition 1, the entropy
is a complex function of temperature given by $S=K(T)$, which leads to the
simplification for stored energy as in Eqn. (21a) Since entropy and entropy
flow are abstract concepts and are not physically measurable, empirical
relations are often utilized. The term ${\left({T-T_{0}}\right)\frac{{\partial
K(T)}}{{dT}}}$ is assumed constant and is called the thermal capacitance
$C_{w}$ used in the simplified model in Eqn. (1a), thereby resulting in
following commonly used relation

$E\approx U=\int_{0}^{t}{T\frac{{\partial
K(T)}}{{dT}}\frac{{dT}}{{dt}}dt}=\underbrace{\frac{{\partial
K(T)}}{{dT}}{\left({T-{T_{0}}}\right)}}_{{C_{w}}}\left({T-{T_{0}}}\right)$
(21a)

Here, we have assumed initial temperature is the ambient temperature $T_{0}$.
Next, the stored energy in tangent space in Eqn. (7) can similarly be
simplified as

$\begin{array}[]{l}{E_{t}}=\int_{0}^{t}{\frac{{dT}}{{dt}}\frac{d}{{dt}}\left({\frac{{\partial
K(T)}}{{dT}}\frac{{dT}}{{dt}}}\right)dt}=\\\ =\int_{0}^{t}{\frac{{\partial
K(T)}}{{dT}}\frac{{dT}}{{dt}}\frac{{{d^{2}}T}}{{d{t^{2}}}}dt}+\int_{0}^{t}{\frac{d}{{dt}}\left({\frac{{\partial
K(T)}}{{dT}}}\right){{\left({\frac{{dT}}{{dt}}}\right)}^{2}}dt}\\\
=\frac{1}{2}\frac{{\partial
K(T)}}{{dT}}{\left({\frac{{dT}}{{dt}}}\right)^{2}}=\frac{{{C_{w}}}}{{2\left(T-T_{0}\right)}}\left({\frac{{dT}}{{dt}}}\right)^{2}\end{array}$
(21b)

In the second equation above, we have assumed the term ${\frac{{\partial
K(T)}}{{dT}}}$ is time independent. We have further assumed the time
derivative of temperature at initial time is zero. For the simplified model,
the injected power as viewed from the perspective of thermal energy domain is
given as in Eqn. (1b). However it also must satisfy the definition of
instantaneous power in Eqn. (8), thereby resulting in:

$\begin{array}[]{l}{S_{f}}\left({T-{T_{0}}}\right)={{\dot{m}}_{a}}{C_{p}}\left({{T_{sup}}-T}\right)\\\
\Rightarrow{S_{f}}=\left({{{\dot{m}}_{a}}\frac{{{C_{p}}}}{{T-{T_{0}}}}}\right)\left({{T_{sup}}-T}\right)\end{array}$
(21c)

Here, we assume the term $S_{f}$ is time-independent. With this effort flow
demarcation for the thermal processes involved, we can derive the rate of
reactive power based on Eqn. (8) as:

$\begin{array}[]{l}{{\dot{Q}}_{T}}=\left({T-{T_{0}}}\right)\frac{d}{{dt}}\left({\frac{{{{\dot{m}}_{a}}{C_{p}}}}{{T-{T_{0}}}}\left({{T_{\sup}}-T}\right)}\right)\\\
\qquad-\left({\frac{{{{\dot{m}}_{a}}{C_{p}}}}{{T-{T_{0}}}}\left({{T_{\sup}}-T}\right)}\right)\frac{d}{{dt}}\left({T-{T_{0}}}\right)\end{array}$
(21d)

Defining the damping losses $\frac{E}{\tau}$ with the time constant
$\tau=\frac{R}{C_{w}}$, we can derive the the first equation in Eqn. (9)
starting from the thermal model in Eqn. (1a) as follows
$\dot{E}=C_{w}\frac{dT}{dt}=-\frac{1}{R}\left(T-T_{0}\right)+P_{u}=-\frac{E}{\tau}+P_{u}$
(22a)

Now let us denote the integrand in Eqn. (21a) as rate of change of stored
energy $p$. Taking its time derivative, we obtain the following simplification

$\begin{array}[]{l}\dot{p}=\frac{{\partial
K(T)}}{{dT}}{\left({\frac{{dT}}{{dt}}}\right)^{2}}+\frac{{\partial
K(T)}}{{\partial T}}\left({T-{T_{0}}}\right)\frac{{{d^{2}}T}}{{d{t^{2}}}}\\\
=2{E_{t}}+{C_{w}}\frac{{{d^{2}}T}}{{d{t^{2}}}}\end{array}$ (22b)

Here we used the expression for stored energy in tangent space from Eqn.
(21b). Expanding the second term using the model in Eqn. (1a), we obtain

$\begin{array}[]{l}\dot{p}=2{E_{t}}+\frac{d}{{dt}}\left({{{\dot{m}}_{a}}{C_{p}}\left({{T_{\sup}}-T}\right)-\frac{1}{R}\left({T-{T_{0}}}\right)}\right)\\\
=2{E_{t}}+\frac{d}{{dt}}\left({{S_{f}}\left({T-{T_{0}}}\right)}\right)-\frac{1}{R}\frac{d}{{dt}}\left({T-{T_{0}}}\right)\\\
=2{E_{t}}+\left({\left({T-{T_{0}}}\right)\frac{{d{S_{f}}}}{{dt}}-{S_{f}}\frac{d}{{dt}}\left({T-{T_{0}}}\right)}\right)\\\
+{S_{f}}\frac{{dT}}{{dt}}-\frac{1}{R}\frac{{dT}}{{dt}}\\\
=2{E_{t}}+{{\dot{Q}}_{T}}+\frac{{{C_{w}}}}{{T-{T_{0}}}}{\left({\frac{{dT}}{{dt}}}\right)^{2}}=4{E_{t}}+{{\dot{Q}}_{T}}\end{array}$
(22c)

Also for a switching device, considering the left hand side and right-hand
side reactive power flows $\dot{Q}_{sw,L},\dot{Q}_{sw,R}$ as depicted in
Figure 11, we have

$\dot{Q}_{sw,R}=\left(\dot{Q}_{sw,L}+\dot{Q}_{u}\right)$ (22d)

Also, since $\dot{Q}_{sw,L}=\dot{Q}^{r,out}$ and $\dot{Q}_{sw,R}=-\dot{Q}_{T}$
in this case, we can further simplify Eqn. (22c) as

$\displaystyle\dot{p}$
$\displaystyle=4E_{t}-\dot{Q}_{sw,R}=4E_{t}-\left(\dot{Q}_{sw,L}+\dot{Q}_{u}\right)$
(22e) $\displaystyle=4E_{t}-\dot{Q}^{r,out}-\dot{Q}_{u}$

∎

##  Appendix B Proof of Theorem 1

###### Proof.

(1): From the expression for stored energy in Eqn. (21a),
$E=C_{w}\left(T-T_{0}\right)$ and
$p=\dot{E}=-\frac{1}{R}\left(T-T_{0}\right)+P^{rated}u=-\frac{1}{R}\left(\frac{E}{C_{w}}-T_{0}\right)++P^{rated}u$.
Clearly, there is a one-one mapping from $\left[T,u\right]$ to
$\left[E,p\right]$. We can thereby utilize energy space model in place of the
conventional state space model for the rest of analysis and control design.
Even if the AHU dynamics were modeled, we could establish diffeomporshism
between states appearing at the interfaces and the energy space variables. For
details, refer to [30].

(2)&(3): Let $\sigma=y_{z}-y_{z}^{ref}$, Starting from Eqn. (12a) and
(LABEL:Eqn:yzRef) and taking the time derivative we have

$\begin{array}[]{l}\frac{d\sigma}{dt}=\frac{d}{{dt}}\left({{y_{z}}-{y_{z}}^{ref}}\right)=-\frac{1}{R}\frac{{dT}}{{dt}}-{{\dot{P}}_{u}}\\\
=-\frac{d}{{dt}}\left({\frac{1}{R}\left({T-{T_{0}}}\right)+{P_{u}}}\right)\\\
=-{C_{w}}\dot{T}=-\dot{p}=-\left({4{E_{t}}+2{{\dot{Q}}_{T}}}\right)+{{\dot{Q}}_{u}}\end{array}$
(23)

Plugging in the control design in Eqn. (12d), we obtain

${\sigma}=\leq\left(\overline{L}-\alpha\right)sign(\sigma)\leq-Ksign(\sigma)$
(24)

Consider next a storage function $V=\frac{1}{2}\sigma^{2}$. We obtain

$\dot{V}=\sigma\dot{\sigma}\leq-K\left|\sigma\right|=-K(2V)^{1/2}$ (25)

From this expression, clearly $\dot{V}\leq 0$ for positive gain $K$ and by
LaSalle’s invariance theorem, we obtain the desired result. Furthermore,
taking the time integral of Eqn. (25), we obtain

$\begin{array}[]{l}{V^{1/2}}(t)-{V^{1/2}}(0)\leq-\sqrt{2}Kt\\\
\Rightarrow|\sigma(t)|-|\sigma(0)|\leq-\sqrt{2}Kt\end{array}$ (26)

For $\sigma(t)\rightarrow 0$, rearranging terms, we see that the maximum
reaching time is $t_{r}\leq\frac{\sigma(0)|}{\sqrt{2}K}$ ∎

##  Appendix C Proof of Proposition 2

Consider the closed loop dynamical model in Eqn. (1a), (12d) and (14). Let us
denote $P^{\prime}=P^{rated}u$ in Eqn. (1a) and let us characterize its
dynamics by reexpressing rate of reactive power in terms of rate of
instantaneous power following the definitions in Eqn. (8). We thereby obtain
$\begin{array}[]{l}\dot{T}=-\frac{1}{{R{C_{w}}}}(T-{T_{0}})+\frac{1}{{{C_{w}}}}{P^{\prime}}\\\
{{\dot{P}}^{\prime}}=2\frac{{{P^{\prime}}}}{{T-{T_{0}}}}\frac{{dT}}{{dt}}-\dot{Q}_{sw,L}\\\
=2\frac{{{P^{\prime}}}}{{T-{T_{0}}}}\frac{{dT}}{{dt}}+\alpha
sign(\sigma)-\dot{Q}_{T}\end{array}$ (27a) The second equation simplification
is a result of Eqn. (22d) and further substuting of the expression of
$\dot{Q}_{u}$ from Eqn. (12d). Simplifying the resulting equation further by
substituting the temperature dynamics, we obtain
${{\dot{P}}^{\prime}}=\underbrace{\left({\frac{{-2}}{{R{C_{w}}}}}\right)}_{a}{P^{\prime}}+\underbrace{\frac{\alpha}{{|\sigma|}}}_{b(\tilde{x})}\sigma+\underbrace{2\frac{1}{{{C_{w}}(T-{T_{0}})}}{P^{\prime}}^{2}+\dot{Q}_{T}}_{c(\tilde{x})}$
(27b) At time t = $nT_{s}$ taking the time derivative to zero. Repeating the
same at $t=(n-1)T_{s}$ and then taking the difference between the two
equations, we obtain $a\Delta{P^{\prime}}[n]=b(\tilde{x})\Delta\sigma[n]$
(27c) Here,
$\Delta\left(.\right)[n]=\left(.\right)(nT_{s})-\left(.\right)((n-1)T_{s})$.
We have ignored the incremental effects of $c(\tilde{x})$ in Eqn. (27b) to
obtain Eqn. (27c). We have also assumed the coefficient $b(\tilde{x})$ is
almost constant between subsequent $nT_{s}$ timesteps. Since
$\sigma=y_{z}-y_{z}^{ref}$, we can establish the following relation
$\Delta\sigma[n]=\Delta y_{z}[n]-\Delta y_{z}^{ref}[n]$ (27d) Furthermore,
from Eqn. (LABEL:Eqn:yzRef), we can write an incremental relation as $\Delta
P^{\prime}[n]=\Delta y_{z}^{ref}[n]-\Delta P^{r,out}[n]$ (27e) Substituting
Eqn. (27e) and (27d) into Eqn. (27c), we obtain
$a(\Delta{y_{z}}^{ref}[n]-\Delta{P^{reg}}[n])=b(\tilde{x})\left({\Delta{y_{z}}[n]-\Delta{y_{z}}^{ref}[n]}\right)$
(27f) Rearranging the terms, we obtain
$\Delta{y_{z}}[n]=\underbrace{\left({1-b{{(\tilde{x})}^{-1}}a}\right)}_{1-\sigma(\tilde{x})}\Delta{y_{z}}^{ref}[n]-\underbrace{\left({b{{(\tilde{x})}^{-1}}a}\right)}_{\sigma(\tilde{x})}\Delta{P^{reg}}[n]$
(27g)

######  Remark 2.

Notice the similarity of HVAC droop relation in Eqn. (27g) to the one
typically utilized for generation resources [31, 32]

$\Delta\omega[n]=\left(1-\sigma_{G}D_{G}\right)\Delta\omega^{ref}[n]-\sigma_{G}\Delta
P^{reg}[n]$ (28)

Here $\sigma_{G}$ is the generator droop and $D$ is the rotor damping.

## Acknowledgment

DISTRIBUTION STATEMENT A. Approved for public release. Distribution is
unlimited.

This material is based upon work supported by the Department of Energy under
Air Force Contract No. FA8702-15-D-0001. Any opinions, findings, conclusions
or recommendations expressed in this material are those of the author(s) and
do not necessarily reflect the views of the Department of Energy.

©2021 Massachusetts Institute of Technology

Delivered to the U.S. Government with Unlimited Rights, as defined in DFARS
Part 252.227-7013 or 7014 (Feb 2014). Notwithstanding any copyright notice,
U.S. Government rights in this work are defined by DFARS 252.227-7013 or DFARS
252.227-7014 as detailed above. Use of this work other than as specifically
authorized by the U.S. Government may violate any copyrights that exist in
this work.

## References

  * [1] Ryan Wiser, Ryan Wiser, Galen Barbose, Lori Bird, Susannah Churchill, Jeff Deyette, and Ed Holt.  Renewable portfolio standards in the united states-a status report with data through 2007.  Technical report, Lawrence Berkeley National Lab.(LBNL), Berkeley, CA (United States), 2008.
  * [2] John Goodin.  California independent system operator demand response & proxy demand resources.  In 2012 IEEE PES Innovative Smart Grid Technologies (ISGT), pages 1–3. IEEE, 2012.
  * [3] Duncan 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(5):1389–1400, 2009.
  * [4] Stephan Koch, Johanna L Mathieu, and Duncan S Callaway.  Modeling and control of aggregated heterogeneous thermostatically controlled loads for ancillary services.  In Proc. PSCC, pages 1–7. Citeseer, 2011.
  * [5] Junjie Qin, Yinlam Chow, Jiyan Yang, and Ram Rajagopal.  Modeling and online control of generalized energy storage networks.  In Proceedings of the 5th international conference on Future energy systems, pages 27–38, 2014.
  * [6] Mahnoosh Alizadeh, Anna Scaglione, Andy Applebaum, George Kesidis, and Karl Levitt.  Scalable and anonymous modeling of large populations of flexible appliances.  arXiv preprint arXiv:1404.1958, 2014.
  * [7] US DOE.  Quadrennial technology review: An assessment of energy technologies and research opportunities, 2015.
  * [8] Wei Zhang, Jianming Lian, Chin-Yao Chang, and Karanjit Kalsi.  Aggregated modeling and control of air conditioning loads for demand response.  IEEE transactions on power systems, 28(4):4655–4664, 2013.
  * [9] Mads Almassalkhi, Jeff Frolik, and Paul Hines.  Packetized energy management: asynchronous and anonymous coordination of thermostatically controlled loads.  In 2017 American Control Conference (ACC), pages 1431–1437. IEEE, 2017.
  * [10] Ning Lu and David P Chassin.  A state-queueing model of thermostatically controlled appliances.  IEEE Transactions on Power Systems, 19(3):1666–1673, 2004.
  * [11] Krishnappa Subbarao, Jason C Fuller, Karanjit Kalsi, Abhishek Somani, Robert G Pratt, Steven E Widergren, and David P Chassin.  Transactive control and coordination of distributed assets for ancillary services.  Technical report, Pacific Northwest National Lab.(PNNL), Richland, WA (United States), 2013.
  * [12] Haiyan Shu, Rongshan Yu, and Susanto Rahardja.  Dynamic incentive strategy for voluntary demand response based on tdp scheme.  In Proceedings of The 2012 Asia Pacific Signal and Information Processing Association Annual Summit and Conference, pages 1–6. IEEE, 2012.
  * [13] He Hao, Borhan M Sanandaji, Kameshwar Poolla, and Tyrone L Vincent.  Aggregate flexibility of thermostatically controlled loads.  IEEE Transactions on Power Systems, 30(1):189–198, 2014.
  * [14] Ning Lu.  An evaluation of the hvac load potential for providing load balancing service.  IEEE Transactions on Smart Grid, 3(3):1263–1270, 2012.
  * [15] Venkatesh Chinde, Krishna Chaitanya Kosaraju, Atul Kelkar, Ramkrishna Pasumarthy, S Sarkar, and Navdeep M Singh.  A passivity-based power-shaping control of building hvac systems.  Journal of Dynamic Systems, Measurement, and Control, 139(11), 2017\.
  * [16] Sumit Mukherjee, Sandipan Mishra, and John T Wen.  Building temperature control: A passivity-based approach.  In 2012 IEEE 51st IEEE Conference on Decision and Control (CDC), pages 6902–6907. IEEE, 2012.
  * [17] Takeshi Hatanaka, Xuan Zhang, Wenbo Shi, Minghui Zhu, and Na Li.  An integrated design of optimization and physical dynamics for energy efficient buildings: A passivity approach.  In 2017 IEEE Conference on Control Technology and Applications (CCTA), pages 1050–1057. IEEE, 2017.
  * [18] Faran A Qureshi and Colin N Jones.  Hierarchical control of building hvac system for ancillary services provision.  Energy and Buildings, 169:216–227, 2018.
  * [19] Thomas Navidi, Abbas El Gamal, and Ram Rajagopal.  A two-layer decentralized control architecture for der coordination.  In 2018 IEEE Conference on Decision and Control (CDC), pages 6019–6024. IEEE, 2018.
  * [20] Borhan M Sanandaji, Tyrone L Vincent, and Kameshwar Poolla.  Ramping rate flexibility of residential hvac loads.  IEEE Transactions on Sustainable Energy, 7(2):865–874, 2015.
  * [21] Network optimized distributed energy systems (nodes) funding opportunity no. de-foa-0001289.  Technical report, US DOE Advanced Research Projects Agency-Energy (ARPA-E), 2015.
  * [22] Andrey Bernstein.  Final report for arpa-e nodes” real-time optimization and control of next-generation distribution infrastructure” project.  Technical report, National Renewable Energy Lab.(NREL), Golden, CO (United States), 2021.
  * [23] Marija D Ilic and Rupamathi Jaddivada.  Fundamental modeling and conditions for realizable and efficient energy systems.  In 2018 IEEE Conference on Decision and Control (CDC), pages 5694–5701. IEEE, 2018.
  * [24] Marija D Ilić and Rupamathi Jaddivada.  Multi-layered interactive energy space modeling for near-optimal electrification of terrestrial, shipboard and aircraft systems.  Annual Reviews in Control, 45:52–75, 2018.
  * [25] Marija Ilic and Rupamathi Jaddivada.  Exergy/energy dynamics-based integrative modeling and control for difficult hybrid aircraft missions.  In AIAA Propulsion and Energy 2019 Forum, page 4501, 2019.
  * [26] Marija Ilic and Rupamathi Jaddivada.  Unified value-based feedback, optimization and risk management in complex electric energy systems.  Optimization and Engineering, pages 1–57, 2020.
  * [27] Kun Deng, Prabir Barooah, Prashant G Mehta, and Sean P Meyn.  Building thermal model reduction via aggregation of states.  In Proceedings of the 2010 American Control Conference, pages 5118–5123. IEEE, 2010.
  * [28] JL Wyatt and M Ilic.  Time-domain reactive power concepts for nonlinear, nonsinusoidal or nonperiodic networks.  In IEEE international symposium on circuits and systems, pages 387–390. IEEE, 1990.
  * [29] Dan Wu, Janak Agrawal, Pallavi Bharadwaj, Le Li, Jinyi Zhang, and Marija Ilic.  On the validity of decomposition for distributed parameter estimation in complex dynamical systems: The case of cooling systems.  In 2020 52nd North American Power Symposium (NAPS), pages 1–6. IEEE, 2021.
  * [30] Rupamathi Jaddivada.  A unified modeling for control of reactive power dynamics in electrical energy systems.  PhD thesis, Massachusetts Institute of Technology, 2020.
  * [31] Marija D Ilić, Nipun Popli, Jhi-Young Joo, and Yunhe Hou.  A possible engineering and economic framework for implementing demand side participation in frequency regulation at value.  In 2011 IEEE Power and Energy Society General Meeting, pages 1–7. IEEE, 2011.
  * [32] Marija D Ilic and John Zaborszky.  Dynamics and control of large electric power systems.  Wiley New York, 2000.


capbtabboxtable[][]

# Simulating open quantum systems using noise models and NISQ devices with
error mitigation

Mainak Roy  Jessica John Britto  Indian Institute of Technology Kharagpur,
Kolkata, India  Ryan Hill  qBraid, 60615, Chicago, United States  Victor
Onofre<EMAIL_ADDRESS>Quantum Open Source Foundation

(August 28, 2024)

###### Abstract

###

In this work, we present simulations of two Open Quantum System models,
Collisional and Markovian Reservoir, with noise simulations, the IBM devices
(ibmq_kyoto, ibmq_osaka) and with the OQC device Lucy. Extending the results
of García-Pérez, et al. [npj Quantum Information 6.1 (2020): 1]. Using the
Mitiq toolkit, we apply Zero-Noise extrapolation (ZNE), an error mitigation
technique and analyze their deviation from the theoretical results for the
models under study. For both models, by applying ZNE, we were able to reduce
the error and overlap it with the theoretical results. All our simulations and
experiments were done in the qBraid environment.

Keywords—Open Quantum Systems, NISQ devices, Quantum Error Mitigation

##  1 Introduction

In the field of quantum computing, all real quantum systems are open and can’t
be perfectly isolated. As a result, there is always an external environment
affecting the system, and the extent of this impact is influenced by various
factors that are difficult to control. To understand the behavior of these
open quantum systems and control it, accurate simulation is crucial,
especially in the gate-based approach.

In this work, we simulated the correlated collisional model and Markovian
reservoir engineering (MRE) based on the results presented in García-Pérez
_et al._ (2020). The models were simulated using Oxford Quantum Computers
(OQC)’s and the latest IBM quantum systems (ibmq_kyoto and ibmq_osaka, IBM
Quantum Eagle type processors), both based on superconducting qubits . To
update and extend the results from the original paperGarcía-Pérez  _et al._
(2020) using IBM Hardware 2020.

We also did an analysis of the effect of noise on the quantum circuits with
depolarising noise models. In the majority of the cases that we considered,
the effects of noise were overcome with Zero Noise Extrapolation (ZNE)
Giurgica-Tiron  _et al._ (2020) using the Mitiq toolkit LaRose  _et al._
(2022). We run all the experiments in the qBraid-SDK environment (a Python
toolkit for cross-framework abstraction, transpilation, and execution of
quantum programs on hardware and simulators Hill  _et al._ (2023)). The
circuits were build with Cirq Developers (2023), using Qiskit Qiskit
contributors (2023) to run on the IBM devices with the IBM Quantum Services
IBM Quantum (2023), and AWS braket Amazon Web Services (2020) for OQC Lucy
Oxford Quantum Circuits (2023).

##  2 Open quantum systems circuit models

###  2.1 Markovian reservoir engineering

This model was developed by García-Pérez  _et al._ (2020), with the initial
ideas for ion traps presented in Barreiro  _et al._ (2011). A carefully
designed environment can give rise to interactions that drive the system to
prepare a maximally entangled state as a result of the dissipative open
system. The four Bell states can be identified as the eigenstates of the
operators $\sigma_{x}^{(1)}\otimes\sigma_{x}^{(2)}$ and
$\sigma_{z}^{(1)}\otimes\sigma_{z}^{(2)}$ as shown in table 1. We can design
two channels, the XX pump and the ZZ pump. The action of each is to pump
states from the +1 eigenspace to the $-1$ eigenspace of the corresponding
operator. The only state belonging to the $-1$ eigenspace of both
$\sigma_{x}^{(1)}\otimes\sigma_{x}^{(2)}$ and
$\sigma_{z}^{(1)}\otimes\sigma_{z}^{(2)}$ is $\ket{\psi^{-}}$. If we compose
the two channels and apply them together, we can pump everything to a single
state, which is $\ket{\psi^{-}}$.

Bell state |  Eigenvalue $\sigma_{x}^{(1)}\otimes\sigma_{x}^{(2)}$ |  Eigenvalue $\sigma_{z}^{(1)}\otimes\sigma_{z}^{(2)}$
---|---|---
$\ket{\phi^{+}}=\frac{1}{\sqrt{2}}\left(\ket{00}+\ket{11}\right)$ |  +1 |  +1
$\ket{\phi^{-}}=\frac{1}{\sqrt{2}}\left(\ket{00}-\ket{11}\right)$ |  $-1$ |  +1
$\ket{\psi^{+}}=\frac{1}{\sqrt{2}}\left(\ket{01}+\ket{10}\right)$ |  +1 |  $-1$
$\ket{\psi^{-}}=\frac{1}{\sqrt{2}}\left(\ket{01}-\ket{10}\right)$ |  $-1$ |  $-1$
Table 1: Eigenvalues of the bell state with operators
$\sigma_{x}^{(1)}\otimes\sigma_{x}^{(2)}$ and
$\sigma_{z}^{(1)}\otimes\sigma_{z}^{(2)}$

Applying any one channel by itself will obviously pump states to a mixture of
the states in the corresponding $-1$ eigenspace (which is
{$\ket{\phi^{-}}$,$\ket{\psi^{-}}$} for the XX pump and
{$\ket{\psi^{+}}$,$\ket{\psi^{-}}$} for the ZZ pump). The forms of the XX pump
and ZZ pump are given (following the notation and results of García-Pérez  _et
al._ (2020)) as
$\Phi_{xx}\rho_{s}=E_{1x}\rho_{s}E_{1x}^{\dagger}+E_{2x}\rho_{s}E_{2x}^{\dagger}$
and
$\Phi_{zz}\rho_{s}=E_{1z}\rho_{s}E_{1z}^{\dagger}+E_{2z}\rho_{s}E_{2z}^{\dagger}$
with,

$E_{1z}=\sqrt{p}\mathbb{I}^{(1)}\otimes\sigma_{x}^{(2)}\frac{1}{2}(\mathbb{I}+\sigma_{z}^{(1)}\otimes\sigma_{z}^{(2)})$
(1)
$E_{2z}=\frac{1}{2}(\mathbb{I}-\sigma_{z}^{(1)}\otimes\sigma_{z}^{(2)})+\frac{\sqrt{1-p}}{2}(\mathbb{I}+\sigma_{z}^{(1)}\otimes\sigma_{z}^{(2)})$
(2)

with $0\leq p\leq 1$, and $E_{1x}$ and $E_{2x}$ are the same with the
replacements $\sigma_{z}^{(2)}\rightarrow\sigma_{x}^{(2)}$ and
$\sigma_{z}^{(1)}\otimes\sigma_{z}^{(2)}\rightarrow\sigma_{x}^{(1)}\otimes\sigma_{x}^{(2)}$.

Given a Bell state, we can map it to a factorised state where one qubit is an
eigenstate of $\sigma_{z}$ and the other is an eigenstate of $\sigma_{x}$ by
applying a CNOT between them. Applying the CNOT on the first system qubit with
the second one as control gives us the table 2.

Bell state | Mapped state
---|---
$\ket{\phi^{+}}=\frac{1}{\sqrt{2}}\left(\ket{00}+\ket{11}\right)$ | $\ket{0}\ket{+}$
$\ket{\phi^{-}}=\frac{1}{\sqrt{2}}\left(\ket{00}-\ket{11}\right)$ | $\ket{0}\ket{-}$
$\ket{\psi^{+}}=\frac{1}{\sqrt{2}}\left(\ket{01}+\ket{10}\right)$ | $\ket{1}\ket{+}$
$\ket{\psi^{-}}=\frac{1}{\sqrt{2}}\left(\ket{01}-\ket{10}\right)$ | $\ket{1}\ket{-}$
Table 2: Mapping of the bell state

Notice how we can now clearly tell which states need to be pumped. We’ll use
this to our advantage.

####  2.1.1 ZZ pump

Figure 1: Quantum circuit for the ZZ pump with inital state $\ket{10}$. With
$[q1,q2]$ as the system qubits and $q0$ as the environment ancillae. With the
angle of the control RY rotation depending on the parameter p.

On operating the ZZ pump on the two qubits, we want the populations of
$\ket{\psi^{+}}$ and $\ket{\psi^{-}}$ to increase since they belong to the
$-1$ eigenspace of the $\sigma_{z}^{(1)}\otimes\sigma_{z}^{(2)}$ operator.

The pump works by applying an X rotation on the qubit which is an eigenstate
of the $\sigma_{z}$ operator. If the initial Bell state belonged to the $+1$
eigenspace, we will rotate the state to the $-1$ eigenspace with some
probability $p$ (the parameter mentioned in Eq. 1 and 2) using a controlled Y
rotation with the angle $\theta=\arccos{(1-2p)}$. If the initial Bell state
belonged to the $-1$ eigenspace, we will leave it as it is.

The simplest way to do that is to initialise an enviroment qubit to $\ket{1}$,
then apply a CNOT on it controlled by the relevant system qubit, and then
apply the pumping X rotation controlled by the environment qubit. The final
circuit for is shown in figure 1

The XX pump follows the same principles as the ZZ pump, but we are dealing
with $\ket{+}$ and $\ket{-}$ here instead of $\ket{0}$ and $\ket{1}$. So we
first apply a Hadamard gate and then we can do the same thing as before.

####  2.1.2 ZZ-XX pump

Figure 2: Quantum circuit for the ZZ-XX pump with inital state $\ket{00}$.
With $[q1,q2]$ as the system qubits and $[q0,q3]$ as the environment
ancillae.With the angle of the control RY rotation depending on the parameter
p.

We will concatenate the ZZ and XX pumps and remove the final CNOT of the ZZ
pump and the first CNOT of the XX pump. The final circuit for is shown in
figure 2

Figure 3: Transpiled quantum circuit for the ZZ-XX pump with inital state
$\ket{10}$. With gate basis $[cx,id,rz,sx]$

In figure 3 we can observe the transpiled circuit for the ZZ-XX pump with the
gate basis used in the IBM devices, $[cx,id,rz,sx]$. We can observe the
increase in depth and in the number of CNOTs used.

###  2.2 Collisional model

This model was developed as the Markovian reservoir by García-Pérez  _et al._
(2020), following the work presented in Filippov  _et al._ (2017). In this
model, a collision consists of applying the unitary operator
$U=e^{ig\tau\sigma_{z}}\otimes\ket{0}_{k}\bra{0}+e^{-ig\tau\sigma_{z}}\otimes\ket{1}_{k}\bra{1}$
between the system qubit and the $k^{th}$ environmental qubit. Here, $g$ is
the coupling strength and $\tau$ is the time duration over which one collision
takes place and $n$ collisions occurs in the time $t=n\tau$. Now we consider
two cases.

####  2.2.1 Correlated case

We’ll take our environment qubits prepared in the classically correlated state
$\rho=\frac{1}{2}(\ket{0}^{\otimes n}\bra{0}^{\otimes n}+\ket{1}^{\otimes
n}\bra{1}^{\otimes n})$. If we consider
$\rho_{S}=a\outerproduct{0}{0}+b\outerproduct{0}{1}+c\outerproduct{1}{0}+d\outerproduct{1}{1}$,
then, $a=b=c=d=\frac{1}{2}$.

While density matrix elements can be complex in general, in this case it is
enough to only get the real part of $\rho_{12}$.

The dynamical map after $n$ collisions (following the results of García-Pérez
_et al._ (2020)) is given by

$\displaystyle\Phi_{t}\rho_{S}$ $\displaystyle=Tr_{E}[U_{n}\cdots
U_{2}U_{1}(\rho_{S}\otimes\rho_{corr})U_{1}^{\dagger}U_{2}^{\dagger}\cdots
U_{n}^{\dagger}]$ (3)
$\displaystyle=\cos^{2}{(ng\tau)}\rho_{S}+\sin^{2}{(ng\tau)}\sigma_{z}\rho_{S}\sigma_{z}$

(since there is only one system qubit), then applying the dynamical map above
gives us

$\displaystyle\Phi_{t}\rho_{S}$
$\displaystyle=a\outerproduct{0}{0}+(\cos^{2}{(ng\tau)}-\sin^{2}{(ng\tau)})b\outerproduct{0}{1}+$
(4)
$\displaystyle(\cos^{2}{(ng\tau)}-\sin^{2}{(ng\tau)})c\outerproduct{1}{0}+d\outerproduct{1}{1}$

Only affecting the coherences, and it multiplies the same factor to both.
Given this, only need to check $\rho_{12}$ in order to quantify the effect of
the dynamical map.

In the general structure of each circuit, first, the system qubit is
initialised to $\ket{+}$ using a Hadamard gate. Then we prepare the
environment qubits, a series of $i$ collisions follow where $i\in[1,c]$.
Finally, we apply a Hadamard gate on the system qubit, with the following
effect,

$\displaystyle H\rho_{S}H^{\dagger}$
$\displaystyle=\frac{1}{\sqrt{2}}\begin{bmatrix}1&1\\\
1&-1\end{bmatrix}\begin{bmatrix}\frac{1}{2}&\frac{x}{2}\\\
\frac{x}{2}&-\frac{1}{2}\end{bmatrix}\frac{1}{\sqrt{2}}\begin{bmatrix}1&1\\\
1&-1\end{bmatrix}$ (5) $\displaystyle=\begin{bmatrix}\frac{1+x}{2}&0\\\
0&\frac{1-x}{2}\end{bmatrix}$

We can then measure the populations of $\ket{0}$ and $\ket{1}$, subtract them,
and then divide by number of shots to find out $x$ (where $x$ is the factor
the coherences are multiplied by in each case). Then we know that
$\rho_{12}=\frac{x}{2}$.

A collision is performed in the following way. We apply a CNOT gate on the
system qubit controlled by the relevant environment qubit, then apply a Z
rotation by the required factor, then apply the same CNOT again.

Figure 4: Circuit diagram for 2 collisions in the correlated case.

In order to prevent the compiler from optimizing the circuit in an undesired
way, we have instructed it to run the circuit without any further gate
optimization through the transpiler. We used the verbatim feature in AWS
Braket and the Qiskit tranpiler options. This is necessary because the
compiler may remove one of the CNOT gates that are involved in each collision,
and combine them into a single rotation gate, which is not what we want to
achieve.

####  2.2.2 Uncorrelated case

In this case, the environment qubits are prepared in the uncorrelated state
$\ket{+}^{\otimes n}$. The dynamical map after $n$ collisions (following the
results of García-Pérez  _et al._ (2020)) is given by

$\Phi_{t}\rho_{S}=\frac{1}{2}(1+\cos^{n}{(2g\tau)})\rho_{S}+\frac{1}{2}(1-\cos^{n}{(2g\tau)})\sigma_{z}\rho_{S}\sigma_{z}$
(6)

Like before, considering
$\rho_{S}=a\outerproduct{0}{0}+b\outerproduct{0}{1}+c\outerproduct{1}{0}+d\outerproduct{1}{1}$
gives us

$\displaystyle\Phi_{t}\rho_{S}$
$\displaystyle=a\outerproduct{0}{0}+\cos^{n}{(2g\tau)}b\outerproduct{0}{1}+$
(7) $\displaystyle\cos^{n}{(2g\tau)}c\outerproduct{1}{0}+d\outerproduct{1}{1}$

The form of the result is similar as the correlated case, so the same
conclusions follow.

Figure 5: Circuit diagram for 3 collisions in the uncorrelated case.

In this case, we simply apply a Hadamard gate on each environment qubit to get
them to the state $\ket{+}^{\otimes n}$. Then we perform collisions with each
of the environment qubits sequentially. The circuit is shown in figure 5.

##  3 Noise models

A method to study the behavior of Quantum Circuits before running on real
quantum computers is to create a noise model. It can be consider that some of
the main errors in a QPU are decoherence and gate errors Georgopoulos  _et
al._ (2021). In this work, we focused on a noise model with different
fidelities of the CNOT gate and one qubit gates. Decreasing the fidelity will
affect the results considerably. Finding a reasonable fidelity in the
simulations will give us an estimate of the resources we will need for a
perfect outcome of the circuits in real QPU. For QPUs from IBM, taking the
information from the IBM quantum platform for the ibmq_peekskill at 13 of
october of 2023. The single-qubit gates’ error rates are approximately
$10^{-4}$ ( $\approx 99.98\%$) and the two-qubit CNOT gate’s error rate is
$10^{-3}$ ( $\approx 99.41\%$). This errors could change depending on the day
and the device.

In the noise models used a simple model with errors from 1 and 2 qubit gate,
which was modeled by applying depolarising noise. For the real hardware we
used the latest devices from IBM with 127 qubits, ibmq_osaka and ibmq_tokyo,
and the OQC device available trough the qBraid environment in AWS.

##  4 Quantum Error Mitigation

We have applied Quantum Error Mitigation using the Zero-Noise Extrapolation
(ZNE)Giurgica-Tiron  _et al._ (2020); Temme  _et al._ (2017) technique from
the Mitiq toolkit LaRose  _et al._ (2022). ZNE extrapolate the expectation
value with zero noise of an observable from a list of expectation values
computed at different noise scales. In this case, the extrapolation is done
over the measure state overlap, different extrapolations methods were tested
with a minimum of 4 scale factors and a maximum of 8. To scaled the noise the
unitary fold gates at random method has been used, which performs the
transformation $G\rightarrow GG^{\dagger}G$ where $G$ is some gate (or to a
random subset of gates). The unitary folding approach has a major advantage
that distinguishes it from other methods. It is entirely digital, which means
that the noise can be scaled using a high level of abstraction from the
physical hardware. Furthermore, this approach does not require knowledge of
the intricate details of the underlying noise model. This makes it an
attractive option for a reliable solution that is easy to implement.

Readout Error Mitigation (REM) Bravyi  _et al._ (2021) was implemented as
well. Errors that occur during the measurement of qubits in a quantum device
can negatively impact quantum computation. These errors are known as readout
errors. By thoroughly analyzing and understanding the readout errors, it is
possible to create a readout error mitigator. This tool can be used to obtain
more accurate output distributions, as well as more precise measurements of
expectation values for measurables. In this work we used the REM technique
with the confusion matrix from Mitiq, a square matrix that encodes the
measurement error for each pair of basis states.

Comparison of the results for the mitigated and unmitigated cases are shown
for each model. The code to reproduce this results are in the Github
repository 111https://github.com/JessicaJohnBritto/QOSF_Cohort6

##  5 Results

###  5.1 Simulation Markovian reservoir

In figure 6 we can see the results of the Markovian reservoir for the ZZ-XX
pump with noisy simulations. Also, we can observe the improvements at the
moment of using ZNE. Each plot compares the overlap between the state of the
system and the Bell states (dots in figure 6 and 7).

Figure 6: (Top) Results with noise simulation for the ZZ-XX pump with fidelity
of CNOT gates of 97 % and the one qubit gates of 99.97 %. (Bottom) Results
with ZNE using a linear extrapolation with 4 scaling factors. The dashed line
is the analytical prediction.In both cases 1024 shots is used.

Figure 7: Results from ibmq_kyoto for the ZZ-XX pump (Top) Results with ZNE
with quadratic extrapolation and 1024 shots. (Bottom) Results with ZNE and REM
using the confusion matrix method. The points correspond to the experimental
results, while the dashed lines show the theoretical prediction.

In figure 7 we can observe the results of running the simulation in a real
quantum device, ibmq_kyoto, for the ZZ-XX pump where the population of states
$\ket{\psi^{-}}$ are clearly increasing with $p$. To compare the results of
the OQC device Lucy, we run the same circuits, taking into account Lucy’s
topology, particularly the direction of two-qubit gates. In the case of the
states $\ket{\phi^{+}}$ and $\ket{\phi^{-}}$ are decreasing, but slowly,
because their populations are increased by one pump and decreased by the
other. The population of state $\ket{\phi^{+}}$ decreases the fastest because
its population is decreased by both pumps. ZNE is applied with a quadratic
extrapolations and scale factors of [1, 3, 5, 7].

Figure 8: Results from ibm_kyoto (Left) Mitigated Results using ZNE with
Mitiq. Different Extrapolations - Richardson, Linear and Quadratic for
$\ket{\psi^{-}}$.

In figure 8, we used different extrapolations that are available in Mitiq -
Richardson, Linear and Quadratic, in an attempt to find out which
extrapolation gives mitigated results that are closest to the theoretical
results. This was particularly done for the case of $\ket{\psi^{-}}$. And, for
the same case, we used scale factors of [1, 3, 5, 7]. The points in figure 8
correspond to the experimental results with 1024 shots, while the dashed lines
show the theoretical prediction. From figure 8, it is clear that Quadratic
Extrapolation gives the results closest to the theoretical results. In the
case of the noisy simulations, the best extrapolation method was the linear.
This difference can be attribute at the complexity of the noise in a real
quantum device compare with the depolarising noise simulations used.

It is important to highlight the improvements on previous works, in García-
Pérez  _et al._ (2020) for the Markovian reservoir was necessary to use 8192
shots to obtain accurate results. In our case we observe that with 1024 shots
is enough for satisfactory results, even though, at the moment of using ZNE
the circuits are being run multiple times. In Abu-Nada and Salhab (2023), the
number of shots and error mitigation technique was not specified.

Figure 9: Results from OQC Lucy for the ZZ-XX pump. The points correspond to
the experimental results, while the dashed lines show the theoretical
prediction. Given the noise of the system the results are inconsistent.

In figure 9 we can observe the results of the Markovian model in the Lucy
device from OQC. Given the noise of the system the results shown an
inconsistent simulation of the model, even though we chose a qubit mapping
with the best quality of connections. Applying ZNE to this high noisy results
will not improve the outcome.

###  5.2 Simulation of the Collisional Model

####  5.2.1 Correlated Case

Figure 10: Results with different fidelity associated with the two-qubit gates
with a constant fidelity on the one qubit gates of 99.97 %. With 1024 shots
and 20 correlated collisions

In figure 10 we can see that with a decrease in fidelity, the results deviate
from the theoretical curve, such that the system undergoes damping. The region
surrounding the maximas and minimas have maximum deviation from the
theoretical results.

Figure 11: Simulated noisy results (coherences vs time $t$) with 99.14 %
fidelity in the CNOT gates for the collisional model with the correlated case.
For the theoretical, mitigated and unmitigated cases using ZNE with 1024 shots
and 20 collisions.

In figure 11 we can observe the results with noisy simulations using the ZNE
technique for error mitigation with a linear extrapolation and 4 scale
factors. We can observe a corrections almost perfect with the mitigated
results.

Running the collisional model in real quantum hardware is more complicated
than the Markovian Reservoir model. The increase in depth of the circuits is
linear with the number of collisions. To understand how this affects the
circuits we can see in figure 12, the depth with different number of
collisions for the case of ibm_kyoto. If the qubit mapping is not specified
the depth of the circuit becomes unstable, in our case we chose a mapping
given the best qubits in the calibrations data of IBM Quantum. For the case of
OQC_Lucy, we used qubits 4 and 5 from the hardware based on their connectivity
(as shown in fig.21) and fidelity rates. With IBM we ran the same circuits,
taking into account Lucy’s topology, as with the Markovian reservoir model.

Figure 12: Circuit Depth w.r.t number of collisions (Left) After Qubit
Mapping. (Right) Without using Qubit Mapping

In figure 13, we observe a comparison between different extrapolation methods
for scale factors of [1, 3, 5, 7]. The experimental results are from ibm_osaka
and with 1024 number of shots. The linear extrapolations shows the best
results. In this case the best extrapolations method with the noise models and
the quantum device are the same. Comparing with the results in García-Pérez
_et al._ (2020), we can observe an improvement in the unmitigated values with
only 1024 shots, which is improve even more in certain cases with ZNE.

Figure 13: Results from ibmq_kyoto for the collisional model in the correlated
case. The results are with different extrapolations methods using ZNE. The
grey points correspond to the experimental,while the dashed lines show the
theoretical prediction.

In the figure 15, we can see the simulation results of using the circuit shown
in Figure 4, where we utilized the same qubit for all collisions. To avoid any
undesired optimization, we enabled the verbatim feature in AWS Braket, which
runs the circuit without any gate optimizations. We also experimented with the
option of not using the verbatim feature, and in this case, we changed the
circuit to use two different qubits for each collision to prevent any unwanted
optimization. The results are shown in the figure 14. However, in both cases,
just like with the Markovian model simulations, the system’s noise was too
high, affecting the simulation output. Each simulation produced inconsistent
results, despite selecting an appropriate qubit mapping for the hardware. It
is worth noting that when we used four qubits without the verbatim feature,
the results began oscillating between similar values. This behavior can be
attributed to a uniform type of noise that the quantum device produces.

Figure 14: Results from OQC Lucy for the correlated collisional model. In this
case we iterate trough two different qubits for each collision. The red points
in the graph represent the experimental results, while the dashed lines
indicate the theoretical prediction. However, due to the noise present in the
system, the results obtained are not consistent.  Figure 15: Results from OQC
Lucy for correlated collisional model. In this case we using the same circuit
shown in figure 4, using the same qubit for the different collisions. To avoid
unwanted optimization we used the verbatim option in AWS Braket.

####  5.2.2 Uncorrelated case

Figure 16: Results of the noisy simulation of the collisional model for the
uncorrelated case with different fidelities of CNOT gates and one fidelity of
one qubit gates of 99.97 % using 1024 shots with 12 collisions.

In figure 16 we can observe the results of the noisy simulations of the
uncorrelated collisional model by varying the fidelity associated with the
two-qubit gates, the CNOT. In all cases we use 1024 shots with 12 collisions.

Figure 17: Simulated noisy results (coherences vs time $t$) with 99.14 %
fidelity in the CNOT gates for the collisional model with the uncorrelated
case. For the theoretical, mitigated and unmitigated cases using ZNE with 8
scale factors with 1024 shots and 12 collisions.

Figure 18: Results from ibmq_osaka for the uncorrelated collisional model with
ZNE using different extrapolation methods with (Top) 1024 and (Bottom) 4096
shots.The grey points correspond to the experimental, while the dashed lines
show the theoretical prediction.

In figure 17 we can observe the results with noisy simulations using ZNE with
a linear extrapolation using 8 scale factors. In this case the results using
the mitigation technique doesn’t show an improvement compared with the
unmitigated one. This results may be because of the number of qubits necessary
to simulated each collision in the uncorrelated case.

In figure 18, we observe a comparison between different extrapolation methods
for scale factors of [1, 3, 5, 7]. The experimental results are from ibm_osaka
with 1024 and 4096 number of shots. The Richardson extrapolations shows the
best results. As with the Markovian model simulations, the difference in
parameters using ZNE with the noise simulations and the real quantum hardware
are because of the complexity of real quantum noise. Also, using REM shows not
significant change in the results. We can observe an improvement in the error
mitigation when increasing the number of shots but still the results are not
stable, given the number of qubits needed for each collision. Comparing with
the results in García-Pérez  _et al._ (2020), we observe an an overall
improvement in the values with less shots. Given the number of qubits and the
levels of noise of OQC Lucy, this model was not run on their device.

##  6 Conclusion

We have updated and expanded the results of García-Pérez, et al. García-Pérez
_et al._ (2020), incorporating error mitigation techniques into the
experiments and using the latest quantum devices. The results show that it is
not only possible to conduct small experiments on Open Quantum Systems using
gate-based NISQ devices, but also the exponential improvements in the
performance of these devices and their cloud services by comparing with
previous results. The noise simulations, which take into account the gate
errors in both the CNOT and one-qubit gates, provide a simple way to study the
circuit’s behavior in a real QPU before sending it to the cloud service.
However, the parameters required for error mitigation may vary depending on
the actual noise present in the quantum hardware. Our results indicate a
decrease in the error rate caused by noisy qubits when Error Mitigation
techniques are employed. In this case, we used the Zero-Noise Extrapolation
method for the Markovian Reservoir and collisional models. For the
uncorrelated collisional model, the results indicate instability when the
number of collisions (qubits) increases, which is consistent with the current
NISQ devices.

The experiments reveal that IBM’s hardware provides highly accurate results
even without error mitigation, whereas OQC hardware produces unstable results
due to its noise, which is consistent with previous quantum volume
benchmarking of the system Pelofske  _et al._ (2022). We expect that as cloud
access to quantum devices expands, we will have greater opportunities to
explore larger Open Quantum Systems with different types of hardware.

## Acknowledgments

This work was supported by qbraid and resulted from the Quantum Open Source
Foundation (QOSF) mentorship program.

## References

  * García-Pérez  _et al._ (2020) G. García-Pérez, M. A. Rossi, and S. Maniscalco, npj Quantum Information 6, 1 (2020).
  * Giurgica-Tiron  _et al._ (2020) T. Giurgica-Tiron, Y. Hindy, R. LaRose, A. Mari, and W. J. Zeng, in  _2020 IEEE International Conference on Quantum Computing and Engineering (QCE)_ (IEEE, 2020) pp. 306–316.
  * LaRose  _et al._ (2022) R. LaRose, A. Mari, S. Kaiser, P. J. Karalekas, A. A. Alves, P. Czarnik, M. El Mandouh, M. H. Gordon, Y. Hindy, A. Robertson, _et al._ , Quantum 6, 774 (2022).
  * Hill  _et al._ (2023) R. J. Hill, R. Young, E. Weis, T. Jun Liang, G. Jacobson, H. Gupta, M. M. Louamri, C. McIrvin, S. Purohit, J. Necaise, E. A. Vara, P. Chakraborty, J. Liu, A. W. Coladangelo, P. Kakhandiki, H. Makhanov, and K. Setia, “qBraid-SDK: Python toolkit for cross-framework abstraction of quantum programs.”  (2023).
  * Developers (2023) C. Developers, “Cirq,”  (2023).
  * Qiskit contributors (2023) Qiskit contributors, “Qiskit: An open-source framework for quantum computing,”  (2023).
  * IBM Quantum (2023) IBM Quantum, “IBM Quantum,”  (2023).
  * Amazon Web Services (2020) Amazon Web Services, “Amazon Braket,”  (2020).
  * Oxford Quantum Circuits (2023) Oxford Quantum Circuits, “Oxford Quantum Circuits,”  (2023).
  * Barreiro  _et al._ (2011) J. T. Barreiro, M. Müller, P. Schindler, D. Nigg, T. Monz, M. Chwalla, M. Hennrich, C. F. Roos, P. Zoller, and R. Blatt, Nature 470, 486 (2011).
  * Filippov  _et al._ (2017) S. N. Filippov, J. Piilo, S. Maniscalco, and M. Ziman, Physical Review A 96, 032111 (2017).
  * Georgopoulos  _et al._ (2021) K. Georgopoulos, C. Emary, and P. Zuliani, Physical Review A 104, 062432 (2021).
  * Temme  _et al._ (2017) K. Temme, S. Bravyi, and J. M. Gambetta, Physical review letters 119, 180509 (2017).
  * Bravyi  _et al._ (2021) S. Bravyi, S. Sheldon, A. Kandala, D. C. Mckay, and J. M. Gambetta, Physical Review A 103, 042605 (2021).
  * Note (1) https://github.com/JessicaJohnBritto/QOSF_Cohort6.
  * Abu-Nada and Salhab (2023) A. A. Abu-Nada and M. A. Salhab, International Journal of Quantum Information , 2340008 (2023).
  * Pelofske  _et al._ (2022) E. Pelofske, A. Bärtschi, and S. Eidenbenz, IEEE Transactions on Quantum Engineering 3, 1 (2022).

##  7 Appendices

##  8 Architecture of Quantum Hardware

In figure 19 and 20, we can see the topology of the IBM hardware used for the
experiments. For OQC Lucy, we can see in figure 21 the topology of the device.
In the case of Lucy, the two qubit connections are different than with the IBM
devices. The calibration data for the IBM hadware is shown in table 3 and 4,
for OQC Lucy in table 5

Property | Median
---|---
T1 ($\mu s$) | 270.36
T2 ($\mu s$) | 170.43
SX error | 1.800e-2
ECR error | 8.257e-3
Readout error | 1.8e-2
Table 3: Calibration data for ibm_osaka from IBM quantum IBM Quantum (2023) with gate basis [ECR, ID, RZ, SX, X].  Figure 19:  Architecture of ibm_osaka Property | Median
---|---
T1 ($\mu s$) | 223.49
T2 ($\mu s$) | 118.92
SX error | 2.484e-2
ECR error | 8.12e-3
Readout error | 1.53e-2
Table 4: Calibration data for ibm_kyoto from IBM quantum IBM Quantum (2023) with gate basis [ECR, ID, RZ, SX, X].  Figure 20:  Architecture of ibm_kyoto Figure 21: Architecture of OQC Lucy Property | Average | Median
---|---|---
T1 ($\mu s$) | 36.460 | 34.469
T2 ($\mu s$) | 33.401 | 37.972
Fidelity (RB) (%) | 99.886 | 99.912
Readout Fidelity (%) | 88.675 | 89.250
ECR gate fidelity (%) | 94.672 | 95.589
Table 5: Calibration data for OQC Lucy from AWS services Amazon Web Services
(2020) with gate basis [ect, i, rz, v, x].