# A unified framework for coordination of thermostatically controlled loads

Austin Coffman<EMAIL_ADDRESS>Ana Bušić<EMAIL_ADDRESS>Prabir Barooah
<EMAIL_ADDRESS>University of Florida, Gainesville, FL, USA  Inria, Paris,
France  University of Florida, Gainesville, FL, USA

###### Abstract

A collection of thermostatically controlled loads (TCLs) – such as air
conditioners and water heaters – can vary their power consumption within
limits to help the balancing authority of a power grid maintain demand supply
balance. Doing so requires loads to coordinate their on/off decisions so that
the aggregate power consumption profile tracks a grid-supplied reference. At
the same time, each consumer’s quality of service (QoS) must be maintained.
While there is a large body of work on TCL coordination, there are several
limitations. One is that they do not provide guarantees on the reference
tracking performance and QoS maintenance. A second limitation of past work is
that they do not provide a means to compute a suitable reference signal for
power demand of a collection of TCLs. In this work we provide a framework that
addresses these weaknesses. The framework enables coordination of an arbitrary
number of TCLs that: (i) is computationally efficient, (ii) is implementable
at the TCLs with local feedback and low communication, and (iii) enables
reference tracking by the collection while ensuring that temperature and
cycling constraints are satisfied at every TCL at all times. The framework is
based on a Markov model obtained by discretizing a pair of Fokker-Planck
equations derived in earlier work by Malhame and Chong [21]. We then use this
model to design randomized policies for TCLs. The balancing authority
broadcasts the same policy to all TCLs, and each TCL implements this policy
which requires only local measurement to make on/off decisions. Simulation
results are provided to support these claims.

###### keywords:

Distributed control, Grid support, Randomized control, Thermostatically
controlled loads.

††thanks: This paper was not presented at any IFAC meeting. Corresponding
author A. Coffman. The research reported here has been partially supported by
the NSF through awards 1646229 (CPS-ECCS) and 1934322 (CPS-ECCS), and the
French National Research Agency grant ANR-16-CE05-0008.

, ,

###### Contents

  1. 1 Introduction
    1. 1.1 Literature review and contribution
    2. 1.2 Notation
  2. 2 Modeling: Individual TCL
    1. 2.1 Temperature dynamics of TCLs
      1. 2.1.1 Policy (at the TCL)
    2. 2.2 PDE model
  3. 3 Markov model from PDE Discretization
    1. 3.1 Spatial discretization
    2. 3.2 Temporal discretization
  4. 4 Discrete space model of a TCL: structure and grid friendly policies
    1. 4.1 Discrete state space
    2. 4.2 Conditional independence in $P_{k}$
      1. 4.2.1 Constructing the factorization
    3. 4.3 BA control command $=$ policy
  5. 5 Proposed framework
    1. 5.1 Individual TCL model with cycling
    2. 5.2 Aggregate model of a collection of TCLs
      1. 5.2.1 Evaluating the aggregate model
    3. 5.3 Grid support Policy design
      1. 5.3.1 Convex control synthesis
      2. 5.3.2 Computational considerations
      3. 5.3.3 Communication burden
  6. 6 Numerical experiments
    1. 6.1 Planning
    2. 6.2 Real time control
  7. 7 Conclusion
  8. A Proofs
    1. A.1 Proof of Lemma 1
    2. A.2 Proof of Lemma 3
    3. A.3 Proof of Lemma 4
    4. A.4 Proof of Theorem 1
      1. A.4.1 $\eta^{*}_{\text{CVX}}\leq\eta^{*}_{\text{NCVX}}$
      2. A.4.2 $\eta^{*}_{\text{NCVX}}\leq\eta^{*}_{\text{CVX}}$
  9. B PDE discretization
    1. B.1 Internal CV’s
    2. B.2 Boundary CV’s
      1. B.2.1 Additional conditions
    3. B.3 Overall system

##  1 Introduction

Many loads are flexible in their power demand: they can vary their demand
around a baseline without adversely affecting consumers’ quality of service
(QoS). The flexibility can be used by a balancing authority (BA) to balance
supply and demand in a power grid. The baseline demand refers to the power
demand under normal operation, when each load operates only to meet its
consumer’s QoS without any interference from the BA. Since the rated power of
each load is small, it is necessary to use a collection of loads. To provide
grid support, the collection has to vary its demand from its baseline. It is
envisioned that the BA would supply a reference signal for power demand and
the actions of the loads in a collection would be coordinated so that their
total demand tracks this reference.

Thermostatically controlled loads (TCLs) - such as residential air
conditioners, heat pumps, and water heaters - are recognized to be valuable
sources of flexible demand [4, 6, 23, 18]. For an air conditioner or a heat
pump, baseline demand is largely dictated by ambient weather conditions. There
are at least two QoS requirements: the indoor temperature must be maintained
within a prespecified range and compressor short-cycling must be avoided,
meaning, once the compressor turns on it cannot turn off until a prespecified
time period elapses, and vice versa. Coordination of TCLs involves two
conflicting requirements: (i) the TCLs collectively need to track the
reference power demand signal, and (ii) every TCL’s QoS need to be maintained.

The actuation at each TCL is discrete: it can either be on or off. Direct load
control [7], in which a centralized controller at the BA directly commands
on/off status of each TCL is not scalable to large populations. A more
scalable idea, that subsequent works on TCL coordination use, is for the BA to
broadcast a low dimensional control command to all TCLs, which is translated
by each TCL into its actuation command with a local policy. To avoid confusion
between the decision making at the BA and a TCL, we use the word “policy” to
mean the algorithm at a TCL that makes on/off decisions. The literature on
decentralized coordination of TCLs differ in their choice of the broadcast
signal (i.e., BA’s control command) and the policy at the TCL that translates
this broadcast to on/off decisions. Coordination architectures can be divided
into two broad categories based on these choices: (i) thermostat set point
change [4, 1] and (ii) probabilistic control [23, 20, 6, 9]. These are
discussed in more detail in Section 1.1.

A framework for coordinating TCLs needs two parts. The coordination scheme is
one part. The other part is reference computation: the framework must provide
the BA with a method to determine a suitable reference signal for the TCLs.
That is, the reference must be such that the TCLs can collectively track the
signal while each TCL maintains its QoS. Otherwise, even the best coordination
scheme will fail to meet either the BA’s need, which is reference tracking, or
the consumers’ need, which is maintaining indoor temperature etc., or both.

This work presents a unified framework for coordination of a collection of
TCLs for providing grid support services. The framework enables both of the
above mentioned components, i.e., (i) planning a suitable reference for a
collection of TCLs and (ii) designing a randomized policy for coordination of
the individual TCLs, so that both the BA’s requirement and consumers’ QoS are
satisfied. In the proposed framework, the BA computes randomized control
policies for the TCLs and broadcasts them to all the TCLs. Each TCL receives
the same policy and implements it using locally measurable information. The
framework is computationally tractable for an arbitrary number of TCLs. The
communication burden is low: only a few numbers need to be broadcast by the BA
at every sampling instant. Feedback from TCLs to the BA can be infrequent.

Underlying the framework is: (i) a Markov chain model that is derived from
partial differential equations developed in the early work of Malhame and
Chong [21], (ii) state augmentation to incorporate cycling constraints, and
(iii) convexification of the non-convex problem that appears in the design of
the randomized control policy for the individual TCL. Additionally, we show
that the assumption made about the effect of weather in earlier work [3] on
randomized control, under certain conditions, is in fact true.

###  1.1 Literature review and contribution

Before reviewing coordination methods, we discuss two interrelated modeling
approaches that underpin many of the ideas in the TCL control architectures.
These are the Markov chain and partial differential equation (PDE) models [16,
21, 26, 17, 28, 24], which stem from the early work of Malhame and Chong [21].
In [21] a pair of coupled Fokker-Planck equations are developed to model a
collection of TCLs under thermostat control. The Fokker-Planck equations are
PDEs that describe the time evolution of a certain probability density
functions (pdf) over the state space of temperature and on/off mode. The PDEs
can be used to model the entire collection or a single TCL: the probability
that a single TCL is “on” is approximately the fraction of TCLs that are “on”.
Discretizing the PDE yields a Markov chain model, though some works have
obtained Markov models without using the PDEs. Hence, _one_ set of PDEs can
model a collection of TCLs. Thus, methods that base control design on the PDE
or Markov chain framework scales well with the number of TCLs.

Due to the lack of scalability of direct load control, we limit our attention
to the two broad classes mentioned earlier: (i) thermostat set point changes,
(ii) probabilistic policy. There are many forms of probabilistic policy, which
can be roughly subdivided into two sub categories: (ii-A) “bin switching” and
(ii-B) “randomized policy”. We discuss these in detail below.

In the thermostat setpoint change coordination architecture, a time-varying
thermostat set point is broadcast to all TCLs, and each TCL makes on/off
decisions based on this new setpoint [4, 1]. This approach may ask for an
extremely small change in thermostat setpoint, far below the resolution of the
temperature sensor at each TCL. Or it may ask for large changes in thermostat
setpoint which will violate occupant comfort.

In a probabilistic policy architecture, the TCL policy - the mapping from BA’s
broadcast command to a TCL’s on/off decision - is a non-deterministic mapping.
Works in this category typically first model the population of TCLs under
thermostat control, which is a deterministic policy, as a Markov chain. The
continuous temperature range is divided into a number of discrete bins. A
finite dimensional state vector, a probability mass function, is then defined.
Each entry of the state vector represents “the fraction of TCLs that are on
(or off) and has temperature in a certain range.”

Since the basic Markov model is derived for the thermostat policy,
introduction of the BA’s control to manipulate TCLs’ on/off state is somewhat
ad-hoc. In the the bin switching literature, the control command from the BA
is chosen so as to affect the fraction of TCLs in the temperature bins
directly. In [23], the BA’s control command is chosen to be another vector,
whose $i^{\text{th}}$ entry represents “the fraction of TCLs in bin $i$ to
increase/decrease”. A policy is then proposed to translate this command to
on/off action at each TCL, which requires knowledge of the state of the Markov
model. In [20], BA’s control command is chosen to be a scalar. The probability
of a TCL turning on or off is proportional to this scalar. Subsequent works
have proposed various refinements, such as BA’s command affecting the rate of
fractions to switch instead of fraction to switch [26]. Providing performance
guarantees with bin switching architecture has proved challenging, either on
reference tracking or on QoS maintenance for individual TCLs.

An alternative to bin switching that still uses probabilistic on/off decision
making is randomized policy [3, 6]. A randomized policy is a specification of
the conditional probability of turning on or off given the current state.
On/off decisions are computed with the help of a random number generator and
the policy. In this architecture it is envisioned that the thermostat policy
at the TCL is replaced with a randomized policy. In [3, 6], the policy is
parameterized by a scalar $\zeta(t)$. Coordination of the population is then
achieved by appropriate design of $\zeta(t)$, which is computed and broadcast
by the BA. This architecture also uses a Markov model of the evolution of
binned temperature, but assumes a certain factorization: the next values of
the temperature and mode are conditionally independent given the current joint
pair of temperature and mode values under the effects of the randomized policy
and exogenous disturbances, especially weather. That is, the transition matrix
of the state process is a point wise product of two controlled transition
matrices. In an optimal control setting, computation of the BA’s control
command, $\zeta(t)$, for reference tracking is a non-convex optimization
problem [11]. The probability of turning on when temperature exceeds the upper
limit, or off when temperature dips below the lower limit, is set to 1 by
design. This will ensure the temperature constraint is maintained. Attempts
have been made to maintain the cycling constraint [9]. But a formal design
method to incorporate the cycling constraint has been lacking.

A complete framework for coordination of TCL needs not only a control
algorithm to make decisions at TCLs, but also a method to compute a _feasible_
reference signal for the collection’s power demand. Feasible means that no TCL
needs to violate local constraints in order for the collection to track the
reference. The topic is sometimes described as “flexibility capacity” and has
been examined in many recent works, with various definitions of flexibility
[15, 25, 13, 10]. A unified treatment of reference design and coordination
algorithm design that would provide a complete framework is lacking.

In short, existing work on TCL coordination has a number of scattered
disadvantages. Direct load control suffers from scalability/privacy issues and
thermostat set-point methods have implementation issues. Bin switching does
not provide guarantees on reference tracking and often requires solving a
challenging state estimation problem. Prior work on randomized control
requires non-convex optimization and is based on an assumed conditional
independence. Finally, there is a lack of unified treatment of the reference
design and policy design problems.

In this work we develop a unified framework for coordination of TCLs that
addresses the weaknesses of prior work described above. Our major
contributions are as follows.

  1. 1.

We provide a complete framework that allows the BA to compute (a) an optimal
reference signal that is feasible for the collection and (b) optimal
randomized policies for the TCLs. When the TCLs implement these policies,
their total power demand collectively tracks the reference signal and the
policies guarantee that temperature and cycling QoS requirements at each TCL
are satisfied. Optimal reference means it is closest to what the BA wants
while being feasible for the TCLs. Implementation of the policy at a TCL is
easy; it requires only local measurements. The communication burden for
coordination is also low. At each sampling time, a randomized control policy -
parameterized by a few numbers - is broadcast to all TCLs. Feedback from TCLs
to the BA can be infrequent.

  2. 2.

Our framework is based on a careful discretization of the partial differential
equation (PDE) model described in [21]. This discretization shows that a
certain “conditional independence” that was assumed in [3] indeed holds. The
conditional independence separates the effects of the policy at the TCL
(control) and weather (disturbance) on the transition matrix, and greatly
facilitates computation of policies.

  3. 3.

Numerical experiments are provided to illustrate the efficacy of the
framework. Simulations show that TCLs are able to track the optimal reference
collectively while each TCL is able to maintain both temperature and cycling
constraints. Matlab implementation is made publicly available at [8].

Figure 1 illustrates the two parts of the proposed framework.

Figure 1: Coordination architecture with the proposed framework.

The Markov model obtained by discretizing a PDE was presented in [12]. For
completeness, we include the discretization in this paper as an Appendix.

###  1.2 Notation

The symbol $\mathbb{1}$ denotes the vector of all ones, $\mathbf{e}_{i}$
denotes the ith canonical basis vector, and $\mathbf{0}$ denotes the zero
matrix or vector, all of appropriate dimension. For a vector $v$,
$\text{diag}(v)$ denotes the diagonal matrix with entries of $v$, i.e.,
$\text{diag}(v)\mathbb{1}=v$. Further, $\otimes$ denotes matrix Kronecker
product and $\mathbf{I}_{A}(\cdot)$ the indicator function of the set $A$.

##  2 Modeling: Individual TCL

A thermostatically controlled load (TCL) is an on/off device that ensures the
temperature of a given environment remains within a specified region, e.g., an
air conditioner. During its operation, the TCL must adhere to certain
operational requirements (QoS constraints). We consider two: the temperature
constraint and the cycling constraint. The temperature constraint is that the
TCL’s temperature must remain within a prespecified deadband,
$[\lambda^{\min},\lambda^{\max}]$. This is achieved by switching the TCL on or
off when it is too hot or cold. The cycling constraint is that the TCL can
only change from “on” to “off” or vice versa once every $\tau$ (discrete) time
instants, where $\tau$ is a prespecified constant. The cycling constraint is
to ensure the mechanical hardware is not damaged. In both cases, ensuring the
two constraints amounts to appropriately deciding when to switch the TCL on or
off.

###  2.1 Temperature dynamics of TCLs

The typical model for the TCL’s temperature $\theta(t)$ in the literature is
the following ordinary differential equation (ODE),

$\displaystyle\begin{split}\frac{d}{dt}\theta(t)&=f_{m}(\theta,t),\quad\text{with}\\\
f_{m}(\theta,t)&=-\frac{1}{RC}\left(\theta-\theta^{a}(t)\right)-m(t)\frac{\eta
P_{0}}{C}.\end{split}$ (1)

The rated electrical power consumption is denoted $P_{0}$ with coefficient of
performance (COP) $\eta$. The parameters $R$ and $C$ denote thermal resistance
and capacitance, respectively. The signal $\theta^{a}(t)$ is the ambient
temperature. The quantity $m(t)$ is the on/off mode, and in the following we
identify $m(t)=1$ and $m(t)=$ on, as well as $m(t)=0$ and $m(t)=$ off. We
denote arbitrary temperature values through the variable $\lambda$, and the
thermostat setpoint as $\lambda^{\text{set}}$. The values
$\lambda^{\text{max}}$ and $\lambda^{\text{min}}$ set the upper and lower
limit for the temperature deadband.

A model for the temperature state that accounts for modeling errors in (1) and
will be crucial in developing the content in Section 2.2 is the following Itô
stochastic differential equation (SDE),

$\displaystyle d\theta(t)=f_{m}(\theta,t)dt+\sigma dB(t).$ (2)

The term $B(t)$ is Brownian motion with parameter $\sigma>0$, and the quantity
$\sigma dB(t)$ captures modeling errors in (1). In either model, the baseline
power for the TCL is the value of $P$ so that
$f_{1}(\lambda^{\text{set}},t)=0$, solving yields:

$\displaystyle\text{Baseline
Power:}\quad\bar{P}^{\text{ind}}(t)=\frac{\theta^{a}(t)-\lambda^{\text{set}}}{\eta
R}.$ (3)

For ${\sf{N_{tcl}}}$ TCLs the baseline power $\bar{P}(t)$ and maximum power
$P_{\text{\footnotesize{agg}}}$ are,

$\displaystyle\bar{P}(t)\triangleq{\sf{N_{tcl}}}\bar{P}^{\text{ind}}(t),\quad\text{and}\quad
P_{\text{\footnotesize{agg}}}\triangleq{\sf{N_{tcl}}}P_{0}.$ (4)

The total electrical power consumption of the collection, whether with
thermostat policy or some other policy, is denoted by $y(t)$:

$\displaystyle y(t)\triangleq P_{0}\sum_{\ell=1}^{{\sf{N_{tcl}}}}m^{\ell}(t)$
(5)

where $m^{\ell}(t)$ is the on/off state of the $\ell$-th TCL.

####  2.1.1 Policy (at the TCL)

The mode state of a TCL evolves according to a policy. The following policy,
which we denote as the _thermostat policy_ , ensures the temperature
constraint:

$\displaystyle\lim_{\epsilon\rightarrow 0}\
m(t+\epsilon)=\begin{cases}1,&\theta(t)\geq\lambda^{\text{max}}.\\\
0,&\theta(t)\leq\lambda^{\text{min}}.\\\ m(t),&\text{o.w.}\end{cases}$ (6)

We add the following set of assumptions about the individual TCL discussed so
far.

  1. A.1

The thermostat policy does not violate the cycling constraint.

  2. A.2

For all $t\geq 0$ and $\theta\in[\lambda^{\text{min}},\lambda^{\text{max}}]$,
$f_{\text{\footnotesize{on}}}(\theta,t)\leq 0$ and
$f_{\text{\footnotesize{off}}}(\theta,t)\geq 0$.

  3. A.3

The TCL’s cycling and temperature constraint are both simultaneously feasible.

The sizing/design of the TCL is most likely to ensure that A.1 holds. With A.1
, we depart from discussing the cycling constraint until Section 5 since up to
that point the mode state is assumed to evolve according to (6).

Assumption A.2 states that when the TCL is on, the temperature does not
increase and when the TCL is off the temperature does not decrease. All prior
works focusing on cooling TCLs (e.g., air conditioners) implicitly make this
assumption. Every result that is to follow is also valid for heating TCLs
(e.g., a water heater or a heat pump) with a sign reversal.

Like A.2, assumption A.3 is also implicit in any work that considers both the
TCLs temperature and cycling constraint.

###  2.2 PDE model

We now describe a PDE model of a TCL’s temperature with thermostat policy
originally derived in [21]. Consider the following marginal pdfs
$\mu_{\text{\footnotesize{on}}},\mu_{\text{\footnotesize{off}}}$:

$\displaystyle\mu_{\text{\footnotesize{on}}}(\lambda,t)d\lambda$
$\displaystyle={\sf P}\left((\lambda<\theta(t)\leq\lambda+d\lambda),\
m(t)=\text{on}\right),$ (7)
$\displaystyle\mu_{\text{\footnotesize{off}}}(\lambda,t)d\lambda$
$\displaystyle={\sf P}\left((\lambda<\theta(t)\leq\lambda+d\lambda),\
m(t)=\text{off}\right),$ (8)

where ${\sf P}(\cdot)$ denotes probability, $\theta(t)$ evolves according to
(2) and for now $m(t)$ evolves according to (6). It was shown in [21] that the
densities $\mu_{\text{\footnotesize{on}}}$ and
$\mu_{\text{\footnotesize{off}}}$ satisfy the Fokker-Planck equations,

$\displaystyle\frac{\partial}{\partial
t}\mu_{\text{\footnotesize{on}}}(\lambda,t)$
$\displaystyle=\frac{\sigma^{2}}{2}\nabla^{2}_{\lambda}\mu_{\text{\footnotesize{on}}}(\lambda,t)-\nabla_{\lambda}\Big{(}f_{\text{\footnotesize{on}}}(\lambda,t)\mu_{\text{\footnotesize{on}}}(\lambda,t)\Big{)}$
(9) $\displaystyle\frac{\partial}{\partial
t}\mu_{\text{\footnotesize{off}}}(\lambda,t)$
$\displaystyle=\frac{\sigma^{2}}{2}\nabla^{2}_{\lambda}\mu_{\text{\footnotesize{off}}}(\lambda,t)-\nabla_{\lambda}\big{(}f_{\text{\footnotesize{off}}}(\lambda,t)\mu_{\text{\footnotesize{off}}}(\lambda,t)\big{)}$
(10)

that are coupled through their boundary conditions [21]. The boundary
conditions are listed in Appendix B.2.

######  Remark 1.

The coupled equations (9)-(10) can be used to model either: (i) a _single_ TCL
or (ii) a _collection_ of TCLs. For (i) the quantities (7)-(8) represent the
_probability_ that a single TCLs temperature and on/off mode reside in the
respective region. For (ii) the quantities (7)-(8) represent the _fraction_ of
TCLs whose temperature and on/off mode reside in the respective region. How
the equations (9)-(10) (specifically their discretized form) can be used to
model an ensemble is discussed further in Section 5.2.

Figure 2: The control volumes (CVs). The colors correspond to the colors found
in Figure 3. The values in each CV represent the nodal temperature for the CV.
The arrows describe the sign of the convection of the TCL through the CVs. The
values are such that $N=m+q$. The terms involving $\alpha$ model rate of
transfer between the corresponding CVs due to the thermostat policy, where
$\alpha=\gamma+\frac{\sigma^{2}}{(\Delta\lambda)^{2}}$. The parameter
$\gamma>0$ is a design parameter; see Remark 4.3.

##  3 Markov model from PDE Discretization

We use the finite volume method (FVM) to discretize the PDEs (9) and (10). The
discretization of (9) and (10) yields a finite dimensional probabilistic model
for a single TCL (equation (16)). We discretize the PDEs (9) and (10) in a way
that a control input for the BA can then be identified. More on this point
will be discussed in Section 4, however the discretization here will play a
role.

###  3.1 Spatial discretization

The FVM bins the continuous temperature into $N$ control volumes (CV). The
layout of the CVs is shown in Figure 2. The $N$ CVs for both the on and off
mode state, as shown in Figure 2, are defined through the nodal temperature
values ($\lambda_{\text{\footnotesize{on}}}$ and
$\lambda_{\text{\footnotesize{off}}}$) and their boundaries
($\lambda^{+}_{\text{\footnotesize{on}}}$ and
$\lambda^{+}_{\text{\footnotesize{off}}}$) and
($\lambda^{-}_{\text{\footnotesize{on}}}$ and
$\lambda^{-}_{\text{\footnotesize{off}}}$):

$\displaystyle\lambda_{\text{\footnotesize{on}}}=(\lambda^{i}_{\text{\footnotesize{on}}})_{i=1}^{N},\quad\lambda_{\text{\footnotesize{on}}}^{+}=\lambda_{\text{\footnotesize{on}}}+\frac{\Delta\lambda}{2},\quad\lambda_{\text{\footnotesize{on}}}^{-}=\lambda_{\text{\footnotesize{on}}}-\frac{\Delta\lambda}{2},$
$\displaystyle\lambda_{\text{\footnotesize{off}}}=(\lambda^{i}_{\text{\footnotesize{off}}})_{i=1}^{N},\quad\lambda_{\text{\footnotesize{off}}}^{+}=\lambda_{\text{\footnotesize{off}}}+\frac{\Delta\lambda}{2},\quad\lambda_{\text{\footnotesize{off}}}^{-}=\lambda_{\text{\footnotesize{off}}}-\frac{\Delta\lambda}{2},$

where $\Delta\lambda$ is the CV width. All intermediate values of
$\lambda_{\text{\footnotesize{on}}}$ and $\lambda_{\text{\footnotesize{off}}}$
are separated from each other by $\Delta\lambda$. The values in
$\lambda^{+}_{\text{\footnotesize{on}}}$ (respectively,
$\lambda^{+}_{\text{\footnotesize{off}}}$) are the right edges of the CVs and
the values $\lambda^{-}_{\text{\footnotesize{on}}}$ (respectively,
$\lambda^{-}_{\text{\footnotesize{off}}}$) are the left edges of the CVs, for
example, $\lambda^{1,-}_{\text{\footnotesize{off}}}=\lambda^{\text{low}}$. The
quantities $\lambda^{\text{min}}$ and $\lambda^{\text{max}}$ specify the
thermostat deadband, and are _different_ from the quantities
$\lambda^{\text{high}}$ and $\lambda^{\text{low}}$ (see Figure 2).

Figure 3: Sparsity pattern of the matrix $A(t)$ for $N=51$ CVs for both the on
and off state. The colors correspond to the colors found in Figure 2.

The steps taken to obtain the spatially discretized PDEs is detailed in
Appendix B. To give an overview, the discretization is done in two parts: (i)
for the internal CV’s (Appendix B.1) and (ii) for the boundary CV’s (Appendix
B.2). We describe here the end result of the derivation in Appendix B. First,
define the following quantities

$\displaystyle\nu_{\text{\footnotesize{off}}}(\lambda^{i},t)$
$\displaystyle\triangleq\mu_{\text{\footnotesize{off}}}(\lambda^{i},t)\Delta\lambda,\quad\text{and}$
(11) $\displaystyle\nu_{\text{\footnotesize{on}}}(\lambda^{i},t)$
$\displaystyle\triangleq\mu_{\text{\footnotesize{on}}}(\lambda^{i},t)\Delta\lambda,$
(12)

then construct the row vector,
$\nu(t)=[\nu_{\text{\footnotesize{off}}}(t),\nu_{\text{\footnotesize{on}}}(t)]$.
with

$\displaystyle\nu_{\text{\footnotesize{off}}}(t)$
$\displaystyle\triangleq[\nu_{\text{\footnotesize{off}}}(\lambda^{1},t),\dots,\nu_{\text{\footnotesize{off}}}(\lambda^{N},t)],\quad\text{and}$
(13) $\displaystyle\nu_{\text{\footnotesize{on}}}(t)$
$\displaystyle\triangleq[\nu_{\text{\footnotesize{on}}}(\lambda^{1},t),\dots,\nu_{\text{\footnotesize{on}}}(\lambda^{N},t)].$
(14)

By combining all the ordinary differential equations (ODEs) for the
$\nu_{\text{\footnotesize{off}}}(\lambda^{i},t),\nu_{\text{\footnotesize{on}}}(\lambda^{i},t)$
for all the $i$’s, we obtain the linear time varying system

$\displaystyle\frac{d}{dt}\nu(t)=\nu(t)A(t).$ (15)

The sparsity pattern of $A(t)$ is shown in Figure 3. The system (15) is the
spatially discretized version of the PDEs (9)-(10). The matrix $A(t)$ also
satisfies the properties of a transition rate matrix, described in the
following lemma.

######  Lemma 1.

For all $t$, the matrix $A(t)$ is a transition rate matrix. That is, for all
$t$

(i): $\displaystyle\quad A(t)\mathbb{1}=\mathbf{0}.$ (ii):
$\displaystyle\quad\text{for all}\ i,\ A_{i,i}(t)\leq 0,\ \text{and}\
\text{for all}\ j\neq i\ A_{i,j}(t)\geq 0.$

###### Proof.

See Appendix A.1. ∎

######  Remark 2.

The choice of the FVM and how we discretize the convection and diffusion terms
appearing in (9)-(10) is important for $A(t)$ to satisfy the conditions in
Lemma 1. This issue is well known in the CFD literature, and also recognized
in the related work [2]. If a finite difference method had been used with
central differences for both diffusion and convection terms, the resulting
$A(t)$ would require restrictive conditions on both $\sigma^{2}$ and
$\Delta\lambda$ to satisfy the properties in Lemma 1 [27].

###  3.2 Temporal discretization

To temporally integrate the dynamics (15) we use a first order Euler
approximation with time step $\Delta t>0$. Making the identifications
$\nu_{k}\triangleq\nu(t_{k})$ and $A_{k}\triangleq A(t_{k})$ we have

$\displaystyle\nu_{k+1}$ $\displaystyle=\nu_{k}P_{k},\quad\text{with}\quad
P_{k}=I+\Delta tA_{k}.$ (16)

In the continuous time setting elements of the vector $\nu(t)$ were referred
to as, for example, $\nu_{\text{\footnotesize{on}}}(\lambda^{i},t)$. The
counterpart to this, in the discrete time setting, is referring to elements of
$\nu_{k}$ as, for example, $\nu_{\text{\footnotesize{on}}}[\lambda^{i},k]$. We
further have the following.

######  Lemma 2.

The matrix $P_{k}$ is a Markov transition probability matrix if

$\displaystyle\forall\ i,\ \text{and}\ \forall\ k,\quad 0<\Delta
t\leq\left|[A_{k}]_{i,i}\right|^{-1}.$

where $[A_{k}]_{i,i}$ is the $i^{th}$ diagonal element of the matrix $A_{k}$.

###### Proof.

From Lemma 1 we have that $P_{k}\mathbb{1}=I\mathbb{1}+\Delta
tA_{k}\mathbb{1}=\mathbb{1}$ since $A_{k}\mathbb{1}=0$. Also from Lemma 1,
every element of $A_{k}$ is non-negative, save for the diagonal elements.
Under the hypothesis on $A_{k}$, then every diagonal element of $I+\Delta
tA_{k}$ will be in $[0,1]$. ∎

######  Remark 3.

The bound on the time step $\Delta t$ given in Lemma 2 is $O(\Delta\lambda)$,
which follows from the PDE discretization; see Appendix B . Since
$\Delta\lambda=\frac{\lambda^{\text{high}}-\lambda^{\text{min}}}{N}$, as the
temperature resolution $\Delta\lambda$ becomes finer the time resolution
$\Delta t$ must also become finer at the same rate. See also Remark 4.3 for a
related comment.

##  4 Discrete space model of a TCL: structure and grid friendly policies

Recall that the dynamics (16) derived in the previous section was for the
thermostat policy. We now delve into the structure of these dynamics so to
introduce a BA control input. We first formalize a discrete state space for
the dynamics (16). We will then show that the transition matrix in (16) can be
written as $P_{k}=\Phi G_{k}$ where $\Phi$ depends on the thermostat policy
and $G_{k}$ on the TCL temperature dynamics and weather. The isolation of the
policy then indicates how a BA could introduce grid friendly policies in place
of the thermostat policy $\Phi$.

###  4.1 Discrete state space

When the conditions of Lemma 2 are met $P_{k}$ is a transition matrix and
hence each $\nu_{k}$ is a marginal pmf if $\nu_{0}$ is a pmf. The structure of
this marginal is given from (7) for the on state (a similar interpretation
holds for the off state) as,

$\displaystyle\nu_{\text{\footnotesize{on}}}[\lambda^{i},k]$
$\displaystyle={\sf P}\left(\theta(t_{k})\in\text{CV}(i),\
m(t_{k})=\text{on}\right),$ (17)

where $\theta(t_{k})$ is the temperature. Now denote,
$\theta_{k}\triangleq\theta(t_{k})$, $m_{k}\triangleq m(t_{k})$, and

$\displaystyle
I_{k}\triangleq\sum_{i=1}^{N}i\mathbf{I}_{\text{CV}(i)}(\theta_{k},m_{k}).$
(18)

The quantity $I_{k}$ indicates which CV the TCLs temperature resides in at
time $k$. It also is a function of $m_{k}$ since the CV index for the on mode
is different from the index for the off mode. We then define the following
discrete state space:

$\displaystyle{\sf
Z}\triangleq\\{m\in\\{\text{\footnotesize{on}},\text{\footnotesize{off}}\\},\
I\in\\{1,\dots,N\\}\\},$ (19)

with cardinality $\left|{\sf Z}\right|=2N$. Using the newly defined quantity
$I_{k}$ we rewrite the marginals
$\nu_{\text{\footnotesize{on}}}[\lambda^{i},k]$ and
$\nu_{\text{\footnotesize{off}}}[\lambda^{i},k]$ as functions on ${\sf Z}$,

$\displaystyle\nu_{\text{\footnotesize{on}}}[\lambda^{i},k]$
$\displaystyle={\sf P}\left(I_{k}=i,\ m_{k}=\text{on}\right),\quad\text{and}$
(20) $\displaystyle\nu_{\text{\footnotesize{off}}}[\lambda^{i},k]$
$\displaystyle={\sf P}\left(I_{k}=i,\ m_{k}=\text{off}\right).$ (21)

From the above, the matrix $P_{k}$ (with the conditions of Lemma 2 satisfied)
is the transition matrix for the joint process $(I_{k},m_{k})$ on the state
space ${\sf Z}$. The dynamic equation $\nu_{k+1}=\nu_{k}P_{k}$ is then a
probabilistic model for a TCL with state space ${\sf Z}$ and operating under
the thermostat policy.

###  4.2 Conditional independence in $P_{k}$

In the following, we refer to the values of $I_{k}$ with $i$ and $j$ and the
values of $m_{k}$ with $u$ and $v$. We introduce the following notation to
refer to the elements of the transition matrix $P_{k}$:

$\displaystyle P_{k}((i,u),(j,v))\triangleq$ (22) $\displaystyle{\sf
P}\Big{(}I_{k+1}=j,\ m_{k+1}=v\ \Big{|}\ I_{k}=i,\ m_{k}=u,\
\theta^{a}_{k}=w_{k}\Big{)}.$

Recall, the matrix $P_{k}$ is derived for the thermostat policy. We will now
show that the matrix $P_{k}$ can be written as the product of two matrices.
One depends only on the thermostat policy (control) and the other depends only
on weather and TCL temperature dynamics. That is, we show that each entry of
$P_{k}$ factors as

$\displaystyle P_{k}((i,u),(j,v))=\phi^{\text{TS}}_{u}(v\ |\ i)P^{u}_{k}(i,j)$
(23)

where, for each given values of $\theta^{a}_{k}$, $P^{u}_{k}(i,j)$ is a
_controlled transition matrix_ on Z:

$\displaystyle P^{u}_{k}(i,j)\triangleq{\sf P}\left(I_{k+1}=j\ |\ I_{k}=i,\
m_{k}=u,\ \theta^{a}_{k}=w_{k}\right)$ (24)

and $\phi^{\text{TS}}_{u}(v\ |\ i)$ is an instance of a _randomized policy_
$\phi_{u}(v\ |\ i)$ on Z:

$\displaystyle\phi_{u}(v\ |\ i)\triangleq{\sf P}\left(m_{k+1}=v\ |\ I_{k}=i,\
m_{k}=u\right).$ (25)

We show the factorization (23) through construction next.

####  4.2.1 Constructing the factorization

The quantity $\phi^{\text{TS}}_{u}(v\ |\ i)$ in (25) is the thermostat policy
on ${\sf Z}$, which is formally defined as follows.

######  Definition 1.

The thermostat policy on ${\sf Z}$ is specified by the two vectors,
$\phi^{\text{TS}}_{\text{\footnotesize{off}}},\phi^{\text{TS}}_{\text{\footnotesize{on}}}\in\mathbb{R}^{N}$,
where
$\phi^{\text{TS}}_{\text{\footnotesize{off}}}\triangleq\phi^{\text{TS}}_{\text{\footnotesize{off}}}(\text{\footnotesize{on}}\
|\ \cdot)=\mathbf{e}_{N}$,
$\phi^{\text{TS}}_{\text{\footnotesize{on}}}\triangleq\phi^{\text{TS}}_{\text{\footnotesize{on}}}(\text{\footnotesize{off}}\
|\ \cdot)=\mathbf{e}_{1}$, and
$\phi^{\text{TS}}_{\text{\footnotesize{off}}}(\text{\footnotesize{off}}\ |\
\cdot)\triangleq 1-\phi^{\text{TS}}_{\text{\footnotesize{off}}}$,
$\phi^{\text{TS}}_{\text{\footnotesize{on}}}(\text{\footnotesize{on}}\ |\
\cdot)\triangleq 1-\phi^{\text{TS}}_{\text{\footnotesize{on}}}$.

The quantity $P^{u}_{k}(i,j)$ in (24) represents the policy-free (open loop)
evolution of the TCL on ${\sf Z}$. That is, it describes how the TCLs
temperature evolves under a fixed mode. We define matrices with entries
$P^{u}_{k}(i,j)$ next.

######  Definition 2.

Let
$P_{k}^{\text{\footnotesize{off}}},P_{k}^{\text{\footnotesize{on}}}\in\mathbb{R}^{N\times
N}$ have $(i,j)$ entries given by,

$\displaystyle P_{k}^{\text{\footnotesize{off}}}(i,j)$
$\displaystyle=P_{w_{k}}((i,\text{\footnotesize{off}}),(j,\text{\footnotesize{off}})),\quad
i\neq N\ \text{and}\ j\neq N,$ $\displaystyle
P_{k}^{\text{\footnotesize{on}}}(i,j)$
$\displaystyle=P_{w_{k}}((i,\text{\footnotesize{on}}),(j,\text{\footnotesize{on}})),\quad
i\neq 1\ \text{and}\ j\neq 1,$

with $P_{k}^{\text{\footnotesize{off}}}(N,N)=1$ and
$P_{k}^{\text{\footnotesize{on}}}(1,1)=1$.

The quantities defined in Definition 1 and 2 correspond to entries of $P_{k}$.
To construct the promised factorization, from these definitions, the idea is
to construct its four sub-matrices that correspond to all possible
combinations of
$u,v\in\\{\text{\footnotesize{on}},\text{\footnotesize{off}}\\}$ (see Figure
3). For example, the $\text{\footnotesize{off}}-\text{\footnotesize{off}}$
quadrant of $P_{k}$ is given by the matrix product

$\displaystyle\big{(}I-\text{diag}(\phi^{\text{TS}}_{\text{\footnotesize{off}}})\big{)}P^{\text{\footnotesize{off}}}_{k}.$

However, since the temperature associated with the $i^{th}$ CV for the on mode
is not the same temperature associated with the $i^{th}$ CV for the off state
(see Figure 2) it is _not_ true that the
$\text{\footnotesize{off}}-\text{\footnotesize{on}}$ quadrant of $P_{k}$ is
given as
$\text{diag}(\phi^{\text{TS}}_{\text{\footnotesize{off}}})P^{\text{\footnotesize{off}}}_{k}$.
The entries of the matrix $P^{\text{\footnotesize{off}}}_{k}$ need to be re-
arranged so to correctly account for the difference in CV index between the
on/off mode. We define such correctly re-arranged matrices next.

######  Definition 3.

Let $I^{\text{\footnotesize{off}}}=\\{m,\dots,N\\}$,
$I^{\text{\footnotesize{on}}}=\\{1,\dots,q\\}$, $m^{-}=m-1$, and
$S_{k}^{\text{\footnotesize{off}}},S_{k}^{\text{\footnotesize{on}}}\in\mathbb{R}^{N\times
N}$ with $(i,j)$ entries

$\displaystyle S_{k}^{\text{\footnotesize{off}}}(i,j-m^{-})$
$\displaystyle=\begin{cases}P_{k}^{\text{\footnotesize{off}}}(i,j)&i,j\in
I^{\text{\footnotesize{off}}}\\\ 0&\text{otherwise}.\end{cases}$ (26)
$\displaystyle S_{k}^{\text{\footnotesize{on}}}(i,j+m^{-})$
$\displaystyle=\begin{cases}P_{k}^{\text{\footnotesize{on}}}(i,j)&i,j\in
I^{\text{\footnotesize{on}}}\\\ 0&\text{otherwise}.\end{cases}$ (27)

The above definition is based on the construction that $N=q+m$. The quantities
in Definition 3 let us construct, e.g., the
$\text{\footnotesize{off}}-\text{\footnotesize{on}}$ quadrant of $P_{k}$ as
$\text{diag}(\phi_{\text{\footnotesize{off}}}^{\text{TS}})S_{k}^{\text{\footnotesize{off}}}$.

The next result shows that $P_{k}=\Phi^{\text{\footnotesize{TS}}}G_{k}$ under
certain conditions and for appropriate choices of the matrices
$\Phi^{\text{\footnotesize{TS}}}$ and $G_{k}$.

######  Lemma 4.1.

Let the time discretization period $\Delta t$ and the parameter $\alpha$ that
appears as a design choice in discretizing the PDEs to ODEs be chosen to
satisfy $\alpha=(\Delta t)^{-1}$. Let
$\Phi^{\text{TS}}_{\text{\footnotesize{off}}}=\text{diag}(\Phi^{\text{TS}}_{\text{\footnotesize{off}}})$
and
$\Phi^{\text{TS}}_{\text{\footnotesize{on}}}=\text{diag}(\Phi^{\text{TS}}_{\text{\footnotesize{on}}})$,
and

$\displaystyle\Phi^{\text{\footnotesize{TS}}}$
$\displaystyle\triangleq\begin{bmatrix}I-\Phi^{\text{TS}}_{\text{\footnotesize{off}}}&\Phi^{\text{TS}}_{\text{\footnotesize{off}}}&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&\Phi^{\text{TS}}_{\text{\footnotesize{on}}}&I-\Phi^{\text{TS}}_{\text{\footnotesize{on}}}\end{bmatrix}\quad\text{and}$
(28) $\displaystyle G_{k}$
$\displaystyle\triangleq\begin{bmatrix}\mathbf{0}&S_{k}^{\text{\footnotesize{off}}}&\mathbf{0}&P_{k}^{\text{\footnotesize{on}}}\\\
P_{k}^{\text{\footnotesize{off}}}&\mathbf{0}&S_{k}^{\text{\footnotesize{on}}}&\mathbf{0}\end{bmatrix}^{T},$
(29)

then

$\displaystyle P_{k}=\Phi^{\text{\footnotesize{TS}}}G_{k}.$ (30)

######  Proof 4.2.

See Appendix A.2.

######  Remark 4.3.

The condition $\alpha=1/\Delta t$ can be satisfied as long as time and
temperature discretization intervals are chosen to satisfy $\Delta
t<(\Delta\lambda)^{2}/\sigma^{2}$. To understand how, recall that in the
discretizing the PDE to the coupled ODEs, a design parameter $\gamma>0$
appears: some rate of density is transferred out of the control volume
$\lambda^{N}_{\text{\footnotesize{off}}}$ and into the CV
$\lambda^{q}_{\text{\footnotesize{on}}}$ (as depicted in Figure 2) due to
thermostatic control. The rate of the density transfer is then given as
$-\gamma\nu_{\text{\footnotesize{off}}}(\lambda^{N},t)$, where $\gamma>0$ is a
modeling choice and a constant of appropriate units that describes the
discharge rate. We then define $\alpha\triangleq D+\gamma$ where
$D=\frac{\sigma^{2}}{(\Delta\lambda)^{2}}$. Recall that $\sigma^{2}$ is the
variance in the Fokker-Planck equation (9)-(10) and $\Delta\lambda$ is the
temperature discretization interval. Thus, as long as $1/\Delta t>D$, a
positive $\gamma$ can be chosen while meeting the condition $\alpha=1/\Delta
t$. The condition $1/\Delta t>D$ is equivalent to $\Delta
t<(\Delta\lambda)^{2}/\sigma^{2}$.

######  Remark 4.4.

The conditional independence factorization (30) has been a useful assumption
in the design of algorithms in [3]. In the present it is a byproduct of our
spatial and temporal discretization of the PDEs (9)-(10). There are other
works [2, 1, 25] that develop Markov models for TCLs through discretization of
PDEs. However, to our knowledge, our work is the first to uncover this
factorization.

Lemma 4.1 informs us how to define the dynamics of the marginals (20) under a
different policy than the thermostat policy, which is described next.

###  4.3 BA control command $=$ policy

In light of the previous section, an arbitrary randomized policy can replace
the thermostat policy to control the state process on ${\sf Z}$. That is
equivalent to replacing $\Phi^{\text{\footnotesize{TS}}}$ in (30) with a new
matrix $\Phi$ that corresponds to a policy designed for grid support. From the
viewpoint of the BA this randomized policy _is_ the control input that it must
design and broadcast to a TCL. The TCL now implements this policy to make
on/off decisions instead of using the thermostat policy. As we shall soon see,
if the BA appropriately designs and sends the randomized policy to multiple
TCLs it can achieve coordination of the TCLs for grid support.

To distinguish from thermostat policy
$\phi^{\text{TS}}_{\text{\footnotesize{off}}}$ and
$\phi^{\text{TS}}_{\text{\footnotesize{on}}}$ in the prior section that only
maintains temperature, we denote the newly introduced policies for providing
grid support with the superscript ‘GS’. We require the policies,
$\phi^{\text{GS}}_{\text{\footnotesize{on}}}$ and
$\phi^{\text{GS}}_{\text{\footnotesize{off}}}$, to have the following
structure

$\displaystyle\phi^{\text{GS}}_{\text{\footnotesize{off}}}(\text{on}\ |\
j)=\begin{cases}\kappa^{\text{\footnotesize{on}}}_{j},&(m+1)\leq
j\leq(N-1).\\\ 1,&j=N.\\\ 0,&\text{o.w.}\end{cases}$ (31)
$\displaystyle\phi^{\text{GS}}_{\text{\footnotesize{on}}}(\text{off}\ |\
j)=\begin{cases}\kappa^{\text{\footnotesize{off}}}_{j},&2\leq j\leq(q-1).\\\
1,&j=1.\\\ 0,&\text{o.w.}\end{cases}$ (32)

with $\phi^{\text{GS}}_{\text{\footnotesize{off}}}(\text{off}\ |\
\cdot)=1-\phi^{\text{GS}}_{\text{\footnotesize{off}}}(\text{on}\ |\ \cdot)$
and $\phi^{\text{GS}}_{\text{\footnotesize{on}}}(\text{on}\ |\
\cdot)=1-\phi^{\text{GS}}_{\text{\footnotesize{on}}}(\text{off}\ |\ \cdot)$
and
$\kappa^{\text{\footnotesize{on}}}_{j},\kappa^{\text{\footnotesize{off}}}_{j}\in[0,1]$
for all $j$. The policies could also be time varying, for example:
$\kappa^{\text{\footnotesize{off}}}_{j}[k]$ and
$\kappa^{\text{\footnotesize{on}}}_{j}[k]$. The dependence of the policies on
time is denoted as $\phi^{\text{GS}}_{\text{\footnotesize{off}}}[k]$ and
$\phi^{\text{GS}}_{\text{\footnotesize{on}}}[k]$. Designing the grid support
control policies is then equivalent to choosing the values of
$\kappa^{\text{\footnotesize{on}}}_{j}[k]$ and
$\kappa^{\text{\footnotesize{off}}}_{j}[k]$ for all $j$ and $k$.

We have required $\phi^{\text{GS}}_{\text{\footnotesize{off}}}(\text{on}\ |\
j)=0$ for $1\leq j\leq m$ since the temperatures corresponding to these
indices are below the permitted deadband temperature, $\lambda^{\text{min}}$.
Hence, turning on at these temperature does not make physical sense. The
arguments for the zero elements in
$\phi^{\text{GS}}_{\text{\footnotesize{on}}}$ are symmetric.

######  Remark 4.5.

From the individual TCL’s perspective, implementing grid support randomized
policies of the form (31)-(32) is straightforward: (i) the TCL measures its
current temperature and on/off status, (ii) the TCL “bins” this temperature
value according to (18) and (iii) the TCL flips a coin to decide its next
on/off state according to the probabilities given in (31)-(32). Note that the
thermostat policy is a special case of the grid support control policy, and
both policies enforce the temperature constraint.

##  5 Proposed framework

We are now in a position to present our unified framework for coordination of
TCLs. We first expand the state of the model (16) so to incorporate cycling,
following [19, 26]. We then shift the viewpoint from a single TCL to that of a
collection of TCLs (recall Remark 1) to develop our control oriented aggregate
model. Using this model we develop a method for designing both reference and
policy through convex optimization.

###  5.1 Individual TCL model with cycling

We now augment the model for a TCL’s temperature evolution with cycling
dynamics. Recall the cycling constraint: as soon as a TCL switches its mode,
the TCL becomes stuck in that mode for $\tau$ time instances. This constraint
can be represented as the evolution of a state, specifically, a counter
variable. First defining the binary variable $s_{k}$ as $s_{k}=1$ if the TCL
is stuck in the current mode at time $k$ and $0$ if it is not stuck. The
counter variable is defined as follows

$\displaystyle{\sf{L}}_{k+1}\triangleq\begin{cases}{\sf{L}}_{k}+1,&s_{k}=1.\\\
0,&s_{k}=0.\end{cases}$ (33)

This variable denotes the time spent in the “stuck” mode ($s_{k}=1$). A TCL
has flexibility to help the grid only when ${\sf{L}}_{k}=0$, which means it is
not stuck in either the on or off mode. If ${\sf{L}}_{k}>0$, it is stuck in
either the on or off mode, and switching the mode to help the grid will
violate the cycling constraint.

Recall, the discrete state space ${\sf Z}$ for a TCL included binned
temperature and on/off mode. The space ${\sf Z}$, the policies
$\phi^{\text{GS}}_{\text{\footnotesize{on}}}$ and
$\phi^{\text{GS}}_{\text{\footnotesize{off}}}$, the marginal pmf $\nu_{k}$,
and the transition matrix $P_{k}$ (and consequently its factors $\Phi$ and
$G_{k}$) now all need to be expanded to be defined over a state space
consisting of $(I_{k},m_{k},{\sf{L}}_{k})$. This expansion is described next.

We denote this newly expanded state space as the set of values: ${\sf
X}\triangleq$

$\displaystyle\Big{\\{}m\in\\{\text{\footnotesize{on}},\text{\footnotesize{off}}\\},\
I\in\\{1,\dots,N\\},\ {\sf{L}}\in\\{0,\dots,\tau\\}\Big{\\}},$ (34)

with cardinality $|{\sf X}|=2N(\tau+1)$. The policies on the expanded state
space are:

$\displaystyle\phi^{\text{E}}_{\text{\footnotesize{off}}}$
$\displaystyle=\mathbf{I}_{\\{0\\}}({\sf{L}})\phi^{\text{GS}}_{\text{\footnotesize{off}}}+(1-\mathbf{I}_{\\{0\\}}({\sf{L}}))\phi^{\text{TS}}_{\text{\footnotesize{off}}},\quad\text{and}$
(35) $\displaystyle\phi^{\text{E}}_{\text{\footnotesize{on}}}$
$\displaystyle=\mathbf{I}_{\\{0\\}}({\sf{L}})\phi^{\text{GS}}_{\text{\footnotesize{on}}}+(1-\mathbf{I}_{\\{0\\}}({\sf{L}}))\phi^{\text{TS}}_{\text{\footnotesize{on}}}.$

To ensure that expanded policy (35) will enforce the cycling constraint, we
impose the following restriction at the design stage: a TCL with
${\sf{L}}_{k}>0$ will only implement the thermostat policy, and a TCL with
${\sf{L}}_{k}=0$ will make on/off decisions based on the grid support policy.
The construction in this way ensures a TCL will not violate its cycling and
temperature constraints under the conditions in Assumption A.2 and A.3.

Each entry of the expanded policy is denoted as
$\phi^{\text{E}}_{\text{\footnotesize{off}}}(u\ |\ j,\ l)$ and
$\phi^{\text{E}}_{\text{\footnotesize{on}}}(u\ |\ j,\ l)$. The expanded
marginals are $\nu_{\text{\footnotesize{off}}}[\lambda^{j},l,k]$ and
$\nu_{\text{\footnotesize{on}}}[\lambda^{j},l,k]$, and
$\nu_{\text{\footnotesize{off}},l}$ (resp.,
$\nu_{\text{\footnotesize{on}},l}$) is shorthand for
$\nu_{\text{\footnotesize{off}}}[\cdot,l,k]$ (resp.,
$\nu_{\text{\footnotesize{on}}}[\cdot,l,k]$). In vectorized form, the expanded
marginal is
$\nu^{\text{E}}=[\nu_{\text{\footnotesize{off}}}^{\text{E}},\nu_{\text{\footnotesize{on}}}^{\text{E}}]$
where
$\nu_{\text{\footnotesize{off}}}^{\text{E}}=[\nu_{\text{\footnotesize{off}},0},\dots,\nu_{\text{\footnotesize{off}},\tau}]$
and
$\nu_{\text{\footnotesize{on}}}^{\text{E}}=[\nu_{\text{\footnotesize{on}},0},\dots,\nu_{\text{\footnotesize{on}},\tau}]$.
Define

$\displaystyle G^{\text{E}}_{k}$
$\displaystyle\triangleq\begin{bmatrix}\mathbf{0}&D_{\tau}\otimes
S_{k}^{\text{\footnotesize{on}}}&\mathbf{0}&C_{\tau}\otimes
P_{k}^{\text{\footnotesize{off}}}\\\ C_{\tau}\otimes
P_{k}^{\text{\footnotesize{on}}}&\mathbf{0}&D_{\tau}\otimes
S_{k}^{\text{\footnotesize{off}}}&\mathbf{0}\end{bmatrix}^{T},$ (36)

where
$D_{\tau}\triangleq\mathbb{1}^{T}\otimes\mathbf{e}_{2}\in\mathbb{R}^{\tau+1\times\tau+1}$
and

$\displaystyle C_{\tau}\triangleq\begin{bmatrix}1&0&\mathbf{0}_{\tau-1}^{T}\\\
\mathbf{0}_{\tau-1}&\mathbf{0}_{\tau-1}&I_{\tau-1}\\\
1&0&\mathbf{0}_{\tau-1}^{T}\end{bmatrix}\in\mathbb{R}^{(\tau+1)\times(\tau+1)}.$
(37)

We define the matrix $\Phi_{k}^{\text{E}}$ as having the same structure as
(28), but with the expanded policies
$\phi^{\text{E}}_{\text{\footnotesize{off}}}$ and
$\phi^{\text{E}}_{\text{\footnotesize{on}}}$, i.e.,

$\displaystyle\Phi^{\text{E}}_{k}\triangleq\begin{bmatrix}I-\Phi^{\text{E}}_{\text{\footnotesize{off}}}[k]&\Phi^{\text{E}}_{\text{\footnotesize{off}}}[k]&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&\Phi^{\text{E}}_{\text{\footnotesize{on}}}[k]&I-\Phi^{\text{E}}_{\text{\footnotesize{on}}}[k]\end{bmatrix},$
(38)

where
$\Phi^{\text{E}}_{\text{\footnotesize{off}}}[k]\triangleq\text{diag}(\phi^{\text{E}}_{\text{\footnotesize{off}}}[k])$
and
$\Phi^{\text{E}}_{\text{\footnotesize{on}}}[k]\triangleq\text{diag}(\phi^{\text{E}}_{\text{\footnotesize{on}}}[k])$.
The model of a TCL with cycling dynamics and grid support policy becomes

$\displaystyle\nu^{\text{E}}_{k+1}=\nu^{\text{E}}_{k}\Phi_{k}^{\text{E}}G_{k}^{\text{E}}.$
(39)

The structure of the transition matrix $\Phi_{k}^{\text{E}}G_{k}^{\text{E}}$
is shown in Figure 4. For comparison, the transition matrix with policy
$\phi^{\text{GS}}$ and _without_ the cycle counter variable would simply be
the four red shaded blocks appearing in their respective quadrant. In the
expanded system, an on to off mode switch forces probability mass from the red
shaded region ($l=0$ and $m=\text{\footnotesize{on}}$) to the green shaded
region ($l=1$ and $m=\text{\footnotesize{off}}$). Mass must then transition
through the chain of $\tau$ green blocks until it reaches the red block again,
so to respect the cycling constraint.

Figure 4: The sparsity pattern of the expanded transition matrix (the dots
represent non-zero entries in the matrix) with $\tau=5$. Each shaded block is
over the entire range of temperature values.

###  5.2 Aggregate model of a collection of TCLs

We now transition from the viewpoint of a single TCL to that of a collection
of ${\sf{N_{tcl}}}$ TCLs: $\ell=1,\dots,{\sf{N_{tcl}}}$. For example,
$m_{k}^{\ell}$ and $I_{k}^{\ell}$ are the mode and binned temperature of the
$\ell^{th}$ TCL at time $k$. Recall Remark 1, the model (39) also describes an
entire collection of TCLs. For a single TCL, we view the state
$\nu_{k}^{\text{E}}$ as a marginal but for a collection of TCLs we expect the
marginal pmf $\nu_{k}^{\text{E}}$ to approximate the histogram

$\displaystyle h_{k}[u,i,l]$
$\displaystyle\triangleq\frac{1}{{\sf{N_{tcl}}}}\sum_{\ell=1}^{{\sf{N_{tcl}}}}\Big{(}\mathbf{I}_{\\{i\\}}(I^{\ell}_{k})\mathbf{I}_{\\{u\\}}(m^{\ell}_{k})\mathbf{I}_{\\{l\\}}({\sf{L}}^{\ell}_{k})\Big{)},$
(40)

for each state $(u,i,l)\in{\sf X}$ as ${\sf{N_{tcl}}}\rightarrow\infty$. In
the same regard, we define

$\displaystyle\gamma_{k}^{\text{E}}\triangleq\nu_{k}^{\text{E}}C^{\text{E}},\quad\text{where}\quad
C^{\text{E}}\triangleq[\mathbf{0}^{T},P_{\text{\footnotesize{agg}}}\mathbb{1}^{T}]^{T},$
(41)

where $P_{\text{\footnotesize{agg}}}$ is the maximum possible power of the
collection, defined in (4). We expect $\gamma_{k}^{\text{E}}$ to approximate
the total power consumption $y_{k}$ of the collection of ${\sf{N_{tcl}}}$
TCLs:

$\displaystyle y_{k}$ $\displaystyle\triangleq
P\sum_{\ell=1}^{{\sf{N_{tcl}}}}m_{k}^{\ell}.$ (42)

which is the discrete-time equivalent of $y(t)$ defined in (5). That is, we
expect $\gamma_{k}^{\text{E}}\approx y_{k}$ for large ${\sf{N_{tcl}}}$, based
on a law of large numbers argument [5]. The _control oriented aggregate model
of a TCL collection_ is the dynamics (39) together with the output (41):

$\displaystyle\nu^{\text{E}}_{k+1}=\nu^{\text{E}}_{k}\Phi_{k}^{\text{E}}G_{k}^{\text{E}}.\qquad\text{and}\qquad\gamma_{k}^{\text{E}}=\nu_{k}^{\text{E}}C^{\text{E}}.$
(43)

####  5.2.1 Evaluating the aggregate model

Before proceeding to policy design with our developed model (43), we first
show that it is effective in modeling a population of TCLs. We do this by
comparing the state of the model to (40) and (42) obtained from a simulation
of ${\sf{N_{tcl}}}$= 50,000 air conditioning TCLs.

The comparison results are shown in Figure 5 and Figure 6. The mode state of
each TCL evolves according to a control policy, where the
$\phi^{\text{GS}}_{\text{\footnotesize{off}}}$ and
$\phi^{\text{GS}}_{\text{\footnotesize{on}}}$ portion are shown in Figure 6
(bottom). The policy is arbitrary, designed merely to be an example of a non-
thermostat policy. This policy satisfies the structure in (31) and (32) so
that both temperature and cycling constraints are satisfied at each TCL. The
temperature evolution evolves according to (2). We see the state
$\nu^{\text{E}}_{k}$ matches the histogram $h_{k}$ of the collection for the
devices that are not stuck (Figure 5 (top)) and for the devices that are stuck
(Figure 5 (bottom)). Additionally, the output of the aggregate model,
$\gamma_{k}^{\text{E}}$, matches it’s empirical counterpart $y_{k}$ (shown in
Figure 6 (top)).

Figure 5: (Top): Histogram of the collection for the devices that are on and
not stuck. (Bottom): Histogram of the collection for the devices that are on
and are stuck.

Figure 6: (Top): Comparison of the output of the expanded aggregate model
$\gamma_{k}^{\text{E}}$ and the ensembles power consumption $y_{k}$. (Bottom):
The policies $\phi^{\text{GS}}_{\text{\footnotesize{off}}}$ and
$\phi^{\text{GS}}_{\text{\footnotesize{on}}}$ used for the numerical
experiment in Section 5.2.

###  5.3 Grid support Policy design

The goal of coordinating TCLs is to help the BA balance supply and demand of
electricity in the grid. We denote $r^{\text{BA}}_{k}$ as the desired demand
from all flexible loads and batteries that will reduce the imbalance to 0. It
is unreasonable to expect any collection of TCLs to meet the entire desired
demand $r^{\text{BA}}_{k}$ while maintaining their QoS. Only a portion of
$r^{\text{BA}}_{k}$ can be supplied by TCLs, and we denote this portion by
$r_{k}$. Determining $r_{k}$ becomes an optimal control problem due to the
time coupling produced by the TCL dynamics. We consider a planning horizon of
$T_{\text{plan}}$. To simultaneously design grid support control policies
$\phi^{\text{GS}}_{\text{\footnotesize{off}}}[k]$ and
$\phi^{\text{GS}}_{\text{\footnotesize{on}}}[k]$ and determine a suitable
reference signal $r_{k}$ over $T_{\text{plan}}$ the BA solves the following
optimization problem,

$\displaystyle\eta^{*}=\min_{\nu^{\text{E}}_{k},\Phi^{\text{E}}_{k}}\ $
$\displaystyle\eta(\hat{\nu})=\sum_{k\in{\sf
T}}\Big{(}r^{\text{BA}}_{k}-\gamma^{\text{E}}_{k}\Big{)}^{2}$ (44) s.t.
$\displaystyle\nu^{\text{E}}_{k+1}=\nu^{\text{E}}_{k}\Phi_{k}^{\text{E}}G_{k}^{\text{E}},\quad\nu^{\text{E}}_{{\sf
T}(0)}=\hat{\nu},$ (45)
$\displaystyle\gamma^{\text{E}}_{k}=\nu_{k}^{\text{E}}C^{\text{E}},\quad\nu^{\text{E}}_{k}\in[0,1],\quad\Phi^{\text{E}}_{k}\in\varPhi.$
(46)

The solution at time $k$ is denoted $r_{k}\triangleq\gamma_{k}^{\text{E},*}$,
$\phi^{\text{GS},*}_{\text{\footnotesize{off}}}[k]$, and
$\phi^{\text{GS},*}_{\text{\footnotesize{on}}}[k]$. We have ${\sf
T}\triangleq\\{{\sf T}(0),\dots,{\sf T}(0)+T_{\text{plan}}-1\\}$ is the index
set of times, ${\sf T}(0)$ denotes the initial time index, $\hat{\nu}$ is the
initial condition, and $\nu_{k}^{\text{E}}\in[0,1]$ holds elementwise. The set
$\varPhi$ collects all of the constraints on the policy. This includes the
equality constraints set by the structural requirements in (31)-(32) and (35)
as well as the structural requirement in (38). These constraints require
certain elements of the policy to be either zero or one. The policy should
also be a valid conditional pmf and its elements in $[0,1]$. Hence, the set
$\varPhi$ is the following convex set

$\displaystyle\varPhi\triangleq\Big{\\{}\Phi\in\mathbb{R}^{|{\sf X}|\times
2|{\sf X}|}_{[0,1]}\ \big{|}\ $ $\displaystyle\Phi\
\text{satisfies}~{}\eqref{eq:fullBAcontPol},\ \mathbb{1}=\Phi\mathbb{1},$
$\displaystyle\phi^{\text{GS}}_{\text{\footnotesize{off}}}\
\text{satisfies}~{}\eqref{eq:randPolOff2On},$
$\displaystyle\phi^{\text{GS}}_{\text{\footnotesize{on}}}\
\text{satisfies}~{}\eqref{eq:randPolOn2Off},\ \text{and}$
$\displaystyle\phi^{\text{E}}_{\text{\footnotesize{off}}}\ \text{and}\
\phi^{\text{E}}_{\text{\footnotesize{on}}}\
\text{satisfy}~{}\eqref{eq:expPolStruct}\Big{\\}}.$ (47)

Where, e.g., $\mathbb{R}^{|{\sf X}|\times|{\sf X}|}_{[0,1]}$ is the set of
$|{\sf X}|\times|{\sf X}|$ matrices with elements in $[0,1]$.

##### QoS + Solution of (44)

  1. 1.

The equality constraints in $\Phi$ are present to ensure the individual TCL’s
QoS constraints: the structure (35) ensures the cycling constraint and the
structure (31)-(32) ensures the temperature constraint. Recall that this
structure guarantees QoS by requiring the policy to place zero probability on
state transitions that would violate QoS.

  2. 2.

A solution to (44) yields, for $k\in{\sf T}$, two things: (i) the optimal
randomized policies $\phi^{\text{GS},*}_{\text{\footnotesize{off}}}[k]$ and
$\phi^{\text{GS},*}_{\text{\footnotesize{on}}}[k]$ and (ii) an optimal
reference for the power demand of the TCL collection
$r_{k}(=\gamma^{\text{E},*}_{k})$. The reference is optimal in the following
sense: among all power demand signals the collection can track without
requiring any TCL to violate its local QoS constraints in so doing, it is the
closest to the BA’s desired demand $r^{\text{BA}}$ in 2-norm. The reference is
also the predicted power consumption of the TCLs whilst using the policies
$\phi^{\text{GS},*}_{\text{\footnotesize{off}}}[k]$ and
$\phi^{\text{GS},*}_{\text{\footnotesize{on}}}[k]$.

######  Remark 5.1.

Since the reference $r_{k}(=\gamma^{\text{E},*}_{k})$ from (44) is the best
the TCLs can do to help the BA without any TCL having to violate its QoS,
Problem (44) therefore also provides an answer to the “aggregate flexibility”
question: how much can a collection of TCLs vary their demand while
maintaining their local QoS constraints. This question has been investigated
by many works [25, 13, 10, 15].

####  5.3.1 Convex control synthesis

The problem (44) is non-convex due to the product
$\nu^{\text{E}}_{k}\Phi_{k}^{\text{E}}$ in the constraint. A well known
convexification remedy for (44) is to consider optimizing over the marginal
and joint distribution instead of the marginal and the policy [22, 2]. Using
our identified structure from Section 4.2 we construct the following joint
distribution (written in matrix form):

$\displaystyle
J_{k}=\text{diag}(\nu^{\text{E}}_{k})\Phi^{\text{E}}_{k}\in\mathbb{R}^{|{\sf
X}|\times 2|{\sf X}|}.$ (48)

By construction, we have that
$\nu^{\text{E}}_{k+1}=\mathbb{1}^{T}J_{k}G^{\text{E}}_{k}$ and
$(\nu^{\text{E}}_{k})^{T}=J_{k}\mathbb{1}$ since
$\mathbb{1}^{T}\text{diag}(\nu^{\text{E}}_{k})=\nu^{\text{E}}_{k}$ and
$\mathbb{1}=\Phi^{\text{E}}_{k}\mathbb{1}$. It is straightforward to convert
the constraint set $\Phi^{\text{E}}_{k}\in\varPhi$ to the new decision
variables. For the equality constraints in $\varPhi$ if we have that
$\phi^{\text{E}}_{\text{\footnotesize{off}}}(u\ |\ j,l)=\kappa$, then in the
decision variables $J_{k}$ and $\nu_{k}^{\text{E}}$ we will have a linear
constraint of the form

$\displaystyle{\sf P}\left(m_{k+1}=u,\ I_{k}=j,\ {\sf{L}}_{k}=l,\
m_{k}=\text{off}\right)$
$\displaystyle=\kappa\nu_{\text{\footnotesize{off}}}[\lambda^{j},l,k],$ (49)

where the LHS of the above is some element in the matrix $J_{k}$. In addition
to the above equality constraints, requiring both $J_{k}$ and
$\nu_{k}^{\text{E}}$ to be within $[0,1]$ and the constraint
$(\nu^{\text{E}}_{k})^{T}=J_{k}\mathbb{1}$ will allow one to reconstruct a
policy $\Phi_{k}^{\text{E}}\in\varPhi$ from $J_{k}$ and $\nu_{k}^{\text{E}}$
(described shortly in Lemma 5.2). We denote the transcription of
$\Phi^{\text{E}}_{k}\in\varPhi$ to the new variables as
$(J_{k},\nu_{k}^{\text{E}})\in\bar{\varPhi}$. Optimizing over $J_{k}$ and
$\nu_{k}^{\text{E}}$ yields the convex program:

$\displaystyle\begin{split}\eta^{*}=&\min_{\nu^{\text{E}}_{k},J_{k}}\
\eta(\hat{\nu})=\sum_{k\in{\sf
T}}\Big{(}r^{\text{BA}}_{k}-\gamma^{\text{E}}_{k}\Big{)}^{2}\\\
\text{s.t.}\quad&\nu^{\text{E}}_{k+1}=\mathbb{1}^{T}J_{k}G^{\text{E}}_{k},\quad\nu^{\text{E}}_{{\sf
T}(0)}=\hat{\nu},\quad\gamma^{\text{E}}_{k}=\nu_{k}^{\text{E}}C^{\text{E}},\\\
&\nu^{\text{E}}_{k},J_{k}\in[0,1],\ (\nu^{\text{E}}_{k})^{T}=J_{k}\mathbb{1},\
(J_{k},\nu_{k}^{\text{E}})\in\bar{\varPhi}.\end{split}$ (50)

Once the convex problem is solved, the grid support control policies need to
be recovered from it by using the relation (48). If the matrix
$\text{diag}(\nu_{k}^{\text{E}})$ is invertible, then the policy can be
obtained trivially from inversion of $\text{diag}(\nu_{k}^{\text{E}})$. If
$\text{diag}(\nu_{k}^{\text{E}})$ is not invertible, then slight care is
required when reconstructing a policy from the solution of (50). We describe
this in the following Lemma.

######  Lemma 5.2.

Suppose for all $k\in{\sf T}$ that $\nu_{k}^{\text{E}}$ and $J_{k}$ satisfy
the constraints in problem (50). Then, there exists matrices
$H_{k}=H_{k}(\nu_{k}^{\text{E}})$ and $W_{k}=W_{k}(\nu_{k}^{\text{E}})$ so
that for all $k\in{\sf T}$ the quantity $\Phi_{k}^{\text{E}}=H_{k}J_{k}+W_{k}$
satisfies (48) and $\Phi_{k}^{\text{E}}\in\varPhi$.

######  Proof 5.3.

See Appendix A.3.

Exact construction of $H_{k}$ and $W_{k}$ is given in the proof of Lemma 5.2.
_Hence, the proof of Lemma 5.2 provides an algorithm for computing grid
support control policies that are feasible for the problem (44) from the
solutions of the convex problem (50)._ Further the two problems have a certain
equivalence described here in the following Theorem.

######  Theorem 1.

Denote $\eta^{*}_{\text{CVX}}$ the optimal cost for (50) and
$\eta^{*}_{\text{NCVX}}$ the optimal cost for (44) we have that
$\eta^{*}_{\text{CVX}}=\eta^{*}_{\text{NCVX}}$.

######  Proof 5.4.

See Appendix A.4.

This result, for a similar problem setup, is also reported in [2]. While we
have no guarantee on the difference of the argument minimizers (and hence the
policies obtained from both), Theorem (1) says that the policies will produce
the same tracking performance. Further, from Lemma 5.2, the policies produced
from either problem are guaranteed to ensure TCL QoS.

####  5.3.2 Computational considerations

The dimension of the program (50) can be quite large, so that even though it
is convex obtaining a solution requires care. We discuss now some practical
considerations that we found necessary to consider when solving the problem
(50).

Due to the structure of $\Phi^{\text{E}}_{k}$, we do not need to declare every
element in the matrix $J_{k}$ as a decision variable since many of these
elements will be zero. For instance, we see that
$\text{diag}(\nu^{\text{E}}_{k})\Phi^{\text{E}}_{k}$ is a block matrix, where
further each matrix block is diagonal. We express this as:
$\text{diag}(\nu^{\text{E}}_{k})\Phi^{\text{E}}_{k}=$

$\displaystyle\begin{bmatrix}B_{\text{\footnotesize{off}},\text{\footnotesize{off}}}[k]&B_{\text{\footnotesize{off}},\text{\footnotesize{on}}}[k]&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&B_{\text{\footnotesize{on}},\text{\footnotesize{off}}}[k]&B_{\text{\footnotesize{on}},\text{\footnotesize{on}}}[k]\end{bmatrix}$
$\displaystyle\mathrel{\ensurestackMath{\stackunder[1pt]{=}{\scriptstyle\triangledown}}}\
\text{sparse}(J_{k})$

where, e.g.,
$B_{\text{\footnotesize{off}},\text{\footnotesize{off}}}[k]=\text{diag}(\nu^{\text{E}}_{\text{\footnotesize{off}}}[k])(I-\Phi^{\text{E}}_{\text{\footnotesize{off}}}[k])$.
The other diagonal matrices appearing in (5.3.2) can be inferred by carrying
out the matrix multiplication.

If $J_{k}$ was declared directly as a decision variable the problem (50) would
have $(8N^{2}+2N(\tau+1))T_{\text{plan}}$ primal variables, whereas the
problem with $\text{sparse}(J_{k})$ as a decision variable only has
$2NT_{\text{plan}}(\tau+3)$ primal variables. As an example, consider $N=12$,
$T_{\text{plan}}=360$, and $\tau=5$, which are values used in numerical
results reported later. The problem (50) without the structure exploited has
$\approx 0.5$ million decision variables, but only $\approx 75,000$ when the
structure is exploited.

We also have found it helpful to include constraints of the form,

$\displaystyle\phi^{\text{GS}}_{\text{\footnotesize{off}}}(\text{on}\ |\
j-1)\nu_{\text{\footnotesize{off}}}[\lambda^{j-1},0,k]\leq\phi^{\text{GS}}_{\text{\footnotesize{off}}}(\text{on}\
|\ j)\nu_{\text{\footnotesize{off}}}[\lambda^{j},0,k],$ (51)
$\displaystyle\phi^{\text{GS}}_{\text{\footnotesize{on}}}(\text{off}\ |\
j+1)\nu_{\text{\footnotesize{on}}}[\lambda^{j+1},0,k]\leq\phi^{\text{GS}}_{\text{\footnotesize{on}}}(\text{off}\
|\ j)\nu_{\text{\footnotesize{on}}}[\lambda^{j},0,k],$ (52)

so to suggest that the switching on (resp., switching off) probability
increases as temperature increases (resp., decreases). Adding the constraints
(51)-(52) to the problem (50) is straightforward as both the LHS and RHS of
the inequalities are elements in the matrix $J_{k}$.

Matlab implementation of (50) and the algorithm to extract the policies from
$J_{k}$ (described in the proof of Lemma 5.2) is available at [8].

####  5.3.3 Communication burden

Once solved, the policies obtained from (50) need to be sent to each
individual TCL. Many of the policy state values are constrained to either zero
or one, which could be pre-programmed into each TCL. At each time index, $q-2$
(for the on to off policy) plus $N-m-1$ (for the off to on policy) numbers are
not constrained and need to be sent from the BA to each TCL. Recall that the
numbers $m$ and $q$ are temperature bin indices (see Figure 2) and $N$ is the
number of temperature bins. For illustrative purposes, consider the values
used in numerical experiments reported in the sequel: $N=12$ with $q=10$ and
$m=2$ and a time discretization $\Delta t=1$ minute. Since $N=q+m$, then the
BA has to broadcast $2(q-1)=18$ numbers every 1 minute to the TCLs. Each TCL
receives the same 18 numbers.

Communication from TCLs to the BA - about their temperature and on/off state -
is needed at the beginning of every planning period so that the BA can
determine the initial condition $\hat{\nu}$ in (50). The frequency of this
feedback is a design choice. In our numerical simulations reported later, a
planning horizon of 6 hours was used, and this feedback was necessary only
once in six hours.More frequent loop closure may be needed for higher
robustness to uncertainty in weather prediction etc., a topic outside the
scope of this paper.

##  6 Numerical experiments

Simulation involving coordination of ${\sf{N_{tcl}}}=20,000$ TCLs through our
proposed framework is presented here. Recall the two parts of the coordination
architecture shown in Figure 1: (i) planning and (ii) real time control.
Planning refers to the solution of the problem (50) at the BA to compute the
following two things for the planning period ${\sf T}$:

  1. 1.

$r_{k}$: the reference power consumption of the TCL collection, given the
problem data $r^{\text{BA}}_{k}$.

  2. 2.

$\phi^{\text{GS},*}_{\text{\footnotesize{off}}}[k]$ and
$\phi^{\text{GS},*}_{\text{\footnotesize{on}}}[k]$: grid support control
policies for each TCL.

This computation is performed at ${\sf T}(0)$. Real time control is then the
implementation of the grid support policies by each TCL to make on/off
decisions in real time. We imagine the BA broadcasts the policies
$\phi^{\text{GS},*}_{\text{\footnotesize{off}}}[k]$ and
$\phi^{\text{GS},*}_{\text{\footnotesize{on}}}[k]$ at each $k$, though it can
also broadcast all the policies, for all $k\in{\sf T}$, at ${\sf T}(0)$ and
not broadcast again until the beginning of the next planning horizon.

The goal of the numerical simulations of real time control is to show the
following.

  1. 1.

When each TCL uses the policies
$\phi^{\text{GS},*}_{\text{\footnotesize{off}}}[k]$ and
$\phi^{\text{GS},*}_{\text{\footnotesize{on}}}[k]$ to decide on/off actuation,
the collection’s power demand indeed tracks $r_{k}$.

  2. 2.

Every TCL’s QoS constraints - both temperature and cycling - are satisfied at
all times.

Temperature of each TCL is computed in these simulations with the ODE model
(1).

Table 1: Simulation Parameters Par. | Unit | value | Par. | Unit | value
---|---|---|---|---|---
${\sf{N_{tcl}}}$ | N/A | 2$\times 10^{4}$ | $\eta$ | $\frac{\text{kW-e}}{\text{kW-th.}}$ | $2.5$
$C$ | kWh$/^{\circ}$C | 1 | $P_{0}$ | kW | 5.5
$\lambda^{\text{min}}$ |  ∘C | 20 | $\lambda^{\text{max}}$ |  ∘C | $22$
$(\Delta t)\tau$ | Mins. | 5 | $P_{\text{\footnotesize{agg}}}$ | MW | 110
$R$ |  ∘C$/$kW | 2 | $\Delta t$ | Mins. | 1
$q$ | N/A | 10 | $m$ | N/A | 2
$N$ | N/A | 12 | $T_{\text{plan}}$ | N/A | 360

###  6.1 Planning

The demand needed for demand-supply imbalance at the BA, $r^{\text{BA}}_{k}$,
is chosen arbitrarily, and shown in Figure 7 (top). It is infeasible for the
collection: sometimes negative and sometimes far higher than the maximum power
demand of the collection. This is done to simulate a realistic scenario in
which many sources of demand and generation, not just TCLs, are managed by the
BA.

The baseline demand trajectory is defined by the equation (4), which is
approximately the power consumption for this collection of air conditioners
under thermostat control. The ambient air temperature is time varying and is
obtained from wunderground.com for a typical summer day in Gainesville,
Florida, USA. The other parameters that affect the Markov model are shown in
Table 1.

Planning computations are done with Matlab and CVX [14] using a desktop Linux
machine, with $N=12$, and for a six hour planning horizon with 1 minute
discretization ($T_{\text{plan}}=360$). The problem (50) takes about a minute
to solve. The quantity $r^{\text{BA}}_{k}$, the baseline power $\bar{P}_{k}$,
and the reference signal $r_{k}$, obtained from solving (50), are shown in
Figure 7 (top). The optimal reference for the collection, $r_{k}$, is as close
to $r^{\text{BA}}_{k}$ as the dynamics of TCLs allows without violating their
QoS constraints; recall Remark 5.1. Figure 7 (bottom) shows the two grid
support control polices for one time instant.

###  6.2 Real time control

The power consumption of the collection making on/off decisions according to
the obtained policies is shown in Figure 8 (top). The figure shows that the
TCLs are able to collectively track the reference signal $r_{k}$. We emphasize
that the computational effort at each TCL is negligible. Recall Remark 4.5:
once a TCL receives a grid support policy ($\approx$ 18 floating point
numbers, see Section 5.3.3) it only has to measure its current state
(temperature and on/off mode) and generate a uniformly distributed random
number in $[0,1]$ to implement the policy.

Verification of the grid support policies in ensuring QoS is shown in Figure
8. The bottom plots shows a histogram of the times between switches for 300
randomly chosen TCLs. The middle plot shows a histogram of temperature from
200 randomly chosen TCLs’ temperature trajectories. The histograms show that
the policies designed with (50) indeed satisfy the QoS constraints, which is
specified by the vertical lines in the figures. Some TCLs do escape the
temperature deadband by a little bit, which is expected and occurs also in
thermostatic control: the sensor must _first_ register a value outside the
deadband in order decide to switch the on/off state.

Figure 7: (Top): The quantity $r_{k}$ obtained from solving (50), the dashed
horizontal lines represent all of the TCLs on (top line) and off (bottom
line). (Bottom): Grid support control policies, obtained from solving (50), at
one time instance.

Figure 8: (Top): Reference tracking results for the TCLs under the influence
of the grid support control policies obtained by solving (50). (Middle):
Histogram of the 200 TCL’s temperature trajectories over the entire simulation
horizon. (Bottom): Histogram of the time between switches over 3000 TCLs with
the vertical line representing the minimum allowable time between switches.

##  7 Conclusion

In this work we present a unified framework for the distributed control of
TCLs. The framework enables: (i) reference planning for a collection of TCLs
and (ii) design of a randomized control policy for the individual TCLs, so
that both the BA’s requirement and consumers’ QoS are satisfied. The resulting
framework is (i) scalable to an arbitrary number of loads and is implemented
through _local_ feedback and minimal communication, (ii) able to guarantee
both temperature and cycling constraints maintenance in each TCL, and (iii)
based on convex optimization. Matlab/cvx implementation is publicly available
[8].

There are several avenues for future work. The optimal control problem is
solved in an open-loop fashion here. Feedback from TCLs is used only to
compute an initial condition that is needed as problem data for the off- line
planning problem. It is straightforward to close the loop between the TCL
collection and the BA with greater frequency for robustness to uncertainty in
weather forecast and TCL parameters. It will be of interest to identify
scenarios where closing loop, say, by using Model Predictive Control, is (i)
necessary, and (ii) at what frequency should information be communicated from
the TCLs to the BA. Another avenue is to investigate how the problem (50)
could be solved at each TCL, intermittently, instead of at the BA. Since the
computational power of the processor at each TCL is lower than that of the
processor at the BA, online distributed algorithms for convex optimization
could play a role. The Fokker-Planck equations from [21] we used here are
convenient for modeling TCL populations with a small deegree of heterogeneity.
Distributed computation of optimal policies locally at each TCLs may help
extend the method to a highly heterogeneous population of TCLs.

## References

  * [1] Saeid Bashash and Hosam K Fathy.  Modeling and control of aggregate air conditioning loads for robust renewable power management.  IEEE Transactions on Control Systems Technology, 21(4):1318–1327, 2013.
  * [2] Emilio Benenati, Marcello Colombino, and Emiliano Dall’Anese.  A tractable formulation for multi-period linearized optimal power flow in presence of thermostatically controlled loads.  In 2019 IEEE 58th Conference on Decision and Control (CDC), pages 4189–4194. IEEE, 2019.
  * [3] Ana Bušić and Sean Meyn.  Distributed randomized control for demand dispatch.  In IEEE conference on decision and control, pages 6964–6971, 2016\.
  * [4] D.S. Callaway and I.A. Hiskens.  Achieving controllability of electric loads.  Proceedings of the IEEE, 99(1):184–199, 2011.
  * [5] Yue Chen, Ana Bušić, and Sean P. Meyn.  State estimation for the individual and the population in mean field control with application to demand dispatch.  IEEE Transactions on Automatic Control, 62(3):1138–1149, 2017.
  * [6] Yue Chen, Md Umar Hashmi, Joel Mathias, Ana Bušić, and Sean Meyn.  Distributed control design for balancing the grid using flexible loads.  In IMA Volume on the Control of Energy Markets and Grids, pages 1–26. 2017.
  * [7] Chi-Min Chu and Tai-Lang Jong.  A novel direct air-conditioning load control method.  IEEE transactions on power systems, 23(3):1356–1363, 2008.
  * [8] Austin Coffman.  Distributed control of TCLs through convex optimization.  https://gitlab.com/austinrc925/distributed-control-of-tcls-through-convex-optimization/-/tree/main, 2021\.
  * [9] Austin Coffman, Ana Bušić, and Prabir Barooah.  Virtual energy storage from TCLs using QoS preserving local randomized control.  In 5th ACM International Conference on Systems for Built Environments (BuildSys), page 10, November 2018.
  * [10] Austin Coffman, Neil Cammardella, Prabir Barooah, and Sean P. Meyn.  Aggregate capacity of TCLs with cycling constraints.  ArXiV.org, October 2020.  arXiv:1909.11497.
  * [11] Austin R. Coffman, Ana Bušić, and Prabir Barooah.  Aggregate capacity for TCLs providing virtual energy storage with cycling constraints.  In IEEE Conference on Decision and Control, December 2019.
  * [12] Austin R. Coffman, Ana Bušić, and Prabir Barooah.  Control oriented modeling of TCLs.  In American Control Conference, May 2021.  accepted, extended version in arXiv:2009.12960.
  * [13] Austin R. Coffman, Zhong Guo, and Prabir Barooah.  Characterizing capacity of flexible loads for providing grid support.  IEEE Transactions on Power Systems, 36:2428 – 2437, May 2021\.
  * [14] M. Grant and S. Boyd.  CVX: Matlab software for disciplined convex programming, version 1.21.  http://cvxr.com/cvx, February 2011.
  * [15] H. Hao, B. M. Sanandaji, K. Poolla, and T. L. Vincent.  Aggregate flexibility of thermostatically controlled loads.  IEEE Transactions on Power Systems, 30(1):189–198, Jan 2015.
  * [16] Emre C Kara, Mario Bergés, and Gabriela Hug.  Impact of disturbances on modeling of thermostatically controlled loads for demand response.  IEEE Transactions on Smart Grid, 6(5):2560–2568, 2015.
  * [17] Adil Khurram, Roland Malhamé, Luis Duffaut Espinosa, and Mads Almassalkhi.  Identification of hot water end-use process of electric water heaters from energy measurements.  Electric Power Systems Research, 189:106625, 2020.
  * [18] Zachary E. Lee, Qingxuan Sun, Zhao Ma, Jiangfeng Wang, Jason S. MacDonald, and K. Max Zhang.  Providing Grid Services With Heat Pumps: A Review.  ASME Journal of Engineering for Sustainable Buildings and Cities, 1(1), 01 2020.  011007\.
  * [19] M. Liu and Y. Shi.  Model predictive control of aggregated heterogeneous second-order thermostatically controlled loads for ancillary services.  IEEE Transactions on Power Systems, 31(3):1963–1971, May 2016.
  * [20] M. Liu, Y. Shi, and X. Liu.  Distributed MPC of aggregated heterogeneous thermostatically controlled loads in smart grid.  IEEE Transactions on Industrial Electronics, 63(2):1120–1129, 2016\.
  * [21] Roland Malhame and Chee-Yee Chong.  Electric load model synthesis by diffusion approximation of a high-order hybrid-state stochastic system.  IEEE Transactions on Automatic Control, 30(9):854–860, 1985.
  * [22] Alan S. Manne.  Linear programming and sequential decisions.  Management Science, 6(3):259–267, 1960.
  * [23] J. L. Mathieu, S. Koch, and D. S. Callaway.  State estimation and control of electric loads to manage real-time energy imbalance.  IEEE Transactions on Power Systems, 28:430–440, 2013.
  * [24] Md Salman Nazir and Ian Hiskens.  Analysis of synchronization in load ensembles.  Electric Power Systems Research, 190:106779, 2021.
  * [25] Dario Paccagnan, Maryam Kamgarpour, and John Lygeros.  On the range of feasible power trajectories for a population of thermostatically controlled loads.  In 2015 54th IEEE Conference on Decision and Control (CDC), pages 5883–5888, 2015.
  * [26] L. C. Totu, R. Wisniewski, and J. Leth.  Demand response of a TCL population using switching-rate actuation.  IEEE Transactions on Control Systems Technology, 25(5):1537–1551, 2017.
  * [27] Henk Kaarle Versteeg and Weeratunge Malalasekera.  An introduction to computational fluid dynamics: the finite volume method.  Pearson education, 2007.
  * [28] 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.

##  Appendix A Proofs

###  A.1 Proof of Lemma 1

See Appendix B before reading this proof. Property (ii) is a consequence of
the upwind difference scheme used. We see that for the internal CVs we have

off CVs:
$\displaystyle\quad-\Big{(}F^{i,+}_{\text{\footnotesize{off}}}+D\Big{)}$
(A.61) on CVs:
$\displaystyle\quad\Big{(}F^{i,-}_{\text{\footnotesize{on}}}-D\Big{)}.$ (A.62)

From Assumption A.2 we have that $F^{i,-}_{\text{\footnotesize{on}}}\leq 0$
and $F^{i,-}_{\text{\footnotesize{off}}}\geq 0$ so that both of the above
terms are negative. The upwind scheme is what ensured appropriate sign was
added to the terms $F^{i,-}_{\text{\footnotesize{on}}}$ and
$F^{i,-}_{\text{\footnotesize{off}}}$ so that the above coefficients are
negative. Similar arguments can be applied for the off diagonal terms of the
internal CVs and the boundary CVs.

To show property (i) we consider solely an internal CV for the off state as
the arguments for all other CVs are identical in structure. Note that showing
$A(t)\mathbb{1}=\mathbf{0}$ is equivalent to
$\mathbb{1}^{T}\mathcal{A}(t)=\mathbf{0}^{T}$. Hence we need to show, for an
arbitrary $i$ that all coefficients acting on
$\nu_{\text{\footnotesize{off}}}(\lambda^{i},t)$ sum to $0$. We collect the
coefficients corresponding to
$\nu_{\text{\footnotesize{off}}}(\lambda^{i},t)$:

$\displaystyle\text{From CV($i$)}:$
$\displaystyle\quad-F^{i,+}_{\text{\footnotesize{off}}}(t)-D.$
$\displaystyle\text{From CV($i-1$)}:$ $\displaystyle\quad\frac{D}{2}.\quad$
$\displaystyle\text{From CV($i+1$)}:$
$\displaystyle\quad\frac{D}{2}+F^{i+1,-}_{\text{\footnotesize{off}}}(t).\quad$

We then require the sum of these coefficients to be zero for all $t$ and any
index $i$ for the internal off CVs, adding yields

$\displaystyle
F^{i+1,-}_{\text{\footnotesize{off}}}(t)-F^{i,+}_{\text{\footnotesize{off}}}(t)=\frac{f_{\text{\footnotesize{off}}}(\lambda^{i+1,-},t)-f_{\text{\footnotesize{off}}}(\lambda^{i,+},t)}{\Delta\lambda}=0$

since by construction $\lambda^{i+1,-}=\lambda^{i,+}$ for the off CV’s. This
procedure can be repeated for $\nu_{\text{\footnotesize{off}}}(\lambda^{i},t)$
with $i\in\\{1,m,N\\}$, i.e., the boundary CVs in the off state and all of the
on CVs in a similar fashion.

###  A.2 Proof of Lemma 3

If $\alpha=(\Delta t)^{-1}$, the diagonal elements of $A_{k}$ with $\alpha$ in
them will go to zero and the non diagonal elements will go to 1. These non-
diagonal elements with value $1$ are the red dots in Figure 3 and encapsulate
the thermostat control law. Thus the construction of
$\Phi^{\text{\footnotesize{TS}}}$ with the canonical basis vectors. Now,
multiplying out the matrix we have,

$\displaystyle\Phi^{\text{\footnotesize{TS}}}G_{k}=\begin{bmatrix}\big{(}I-\Phi^{\text{\footnotesize{TS}}}_{\text{\footnotesize{off}}}\big{)}P_{k}^{\text{\footnotesize{off}}}&\Phi^{\text{\footnotesize{TS}}}_{\text{\footnotesize{off}}}S_{k}^{\text{\footnotesize{off}}}\\\
\Phi^{\text{\footnotesize{TS}}}_{\text{\footnotesize{on}}}S_{k}^{\text{\footnotesize{on}}}&\big{(}I-\Phi^{\text{\footnotesize{TS}}}_{\text{\footnotesize{on}}}\big{)}P_{k}^{\text{\footnotesize{on}}}\end{bmatrix}$
(A.63)

where
$\big{(}I-\Phi^{\text{\footnotesize{TS}}}_{\text{\footnotesize{off}}}\big{)}P_{k}^{\text{\footnotesize{off}}}$
(respectively,
$\big{(}I-\Phi^{\text{\footnotesize{TS}}}_{\text{\footnotesize{on}}}\big{)}P_{k}^{\text{\footnotesize{on}}}$)
is the matrix $P_{k}^{\text{\footnotesize{off}}}$ (respectively,
$P_{k}^{\text{\footnotesize{on}}}$) but with the last (respectively, first)
row zeroed out. The exact opposite statement is true for
$\Phi^{\text{\footnotesize{TS}}}_{\text{\footnotesize{off}}}P_{k}^{\text{\footnotesize{on}}}$
and
$\Phi^{\text{\footnotesize{TS}}}_{\text{\footnotesize{on}}}P_{k}^{\text{\footnotesize{off}}}$.
Hence, by definition of the matrices in $G_{k}$ we have
$P_{k}=\Phi^{\text{\footnotesize{TS}}}G_{k}$ where each non-zero element holds
the interpretation (23).

###  A.3 Proof of Lemma 4

We define the following transformation for $l\in\\{0,\dots,\tau\\}$ and
$j\in\\{1,\dots,N\\}$ as

$\displaystyle T(j,l)=lN+j$ (A.64)

that maps the integers $j$ and $l$ that label the state values to the absolute
index of either of the vectors $\nu^{\text{E}}_{\text{\footnotesize{off}}}$
and $\nu^{\text{E}}_{\text{\footnotesize{on}}}$. Now consider the following
two sets

$\displaystyle\mathcal{W}_{\text{\footnotesize{off}}}$
$\displaystyle\triangleq\Big{\\{}(u,j,l)\in{\sf X}\ \Big{|}\
\phi^{\text{E}}_{\text{\footnotesize{off}}}(u\ |\ j,\
l)=\beta_{\text{\footnotesize{off}}}(u,j,l)\Big{\\}}$ (A.65)
$\displaystyle\mathcal{W}_{\text{\footnotesize{on}}}$
$\displaystyle\triangleq\Big{\\{}(u,j,l)\in{\sf X}\ \Big{|}\
\phi^{\text{E}}_{\text{\footnotesize{on}}}(u\ |\ j,\
l)=\beta_{\text{\footnotesize{on}}}(u,j,l)\Big{\\}}.$ (A.66)

The values $\beta_{\text{\footnotesize{off}}}$ and
$\beta_{\text{\footnotesize{on}}}$ are chosen to ensure the structural
requirements in (31)-(32) and (35). For example, for $l=1$ and $u=\text{on}$
we have that
$\beta_{\text{\footnotesize{off}}}(\text{\footnotesize{on}},\cdot,1)=\phi^{\text{TS}}_{\text{\footnotesize{off}}}(\text{\footnotesize{on}}|\cdot)$
(and hence
$\beta_{\text{\footnotesize{off}}}(\text{\footnotesize{off}},\cdot,1)=1-\phi^{\text{TS}}_{\text{\footnotesize{off}}}(\text{\footnotesize{on}}|\cdot)$)
so to enforce the structural requirement in (35). Define for each $k\in{\sf
T}$ , $u,v\in\\{\text{\footnotesize{on}},\text{\footnotesize{off}}\\}$,
$j\in\\{1,\dots,N\\}$, and $l\in\\{0,\dots,\tau\\}$ the following vectors

$\displaystyle h^{u}_{k}[T(j,l)]$
$\displaystyle\triangleq\begin{cases}(\nu_{u}^{\text{E}}[\lambda^{j},l,k])^{-1}&\text{if}\
\nu_{u}^{\text{E}}[\lambda^{j},l,k]>0.\\\ 0&\text{otherwise}.\end{cases}$
$\displaystyle w^{u,v}_{k}[T(j,l)]$
$\displaystyle\triangleq\begin{cases}\beta_{v}(u,j,l)&\text{if}\
(u,j,l)\in\mathcal{W}_{v}\ \text{and}\\\
&\nu_{v}^{\text{E}}[\lambda^{j},l,k]=0.\\\ 0.5&\text{if}\
(u,j,l)\notin\mathcal{W}_{v}\ \text{and}\\\
&\nu_{v}^{\text{E}}[\lambda^{j},l,k]=0.\\\ 0&\text{otherwise}.\end{cases}$

where $T(\cdot,\cdot)$ is defined in (A.64),
$\mathcal{W}_{\text{\footnotesize{off}}}$ in (A.65), and
$\mathcal{W}_{\text{\footnotesize{on}}}$ in (A.66). Let
$W^{u,v}_{k}=\text{diag}(w^{u,v}_{k})$ and
$H^{u,v}_{k}=\text{diag}(h^{u,v}_{k})$, and construct the following matrices

$\displaystyle H_{k}$
$\displaystyle=\begin{bmatrix}H_{k}^{\text{\footnotesize{off}}}&\mathbf{0}\\\
\mathbf{0}&H_{k}^{\text{\footnotesize{on}}}\end{bmatrix},\quad\text{and}$
(A.67) $\displaystyle W_{k}$
$\displaystyle=\begin{bmatrix}W^{\text{\footnotesize{off}},\text{\footnotesize{off}}}_{k}&W^{\text{\footnotesize{off}},\text{\footnotesize{on}}}_{k}&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&W^{\text{\footnotesize{on}},\text{\footnotesize{off}}}_{k}&W^{\text{\footnotesize{on}},\text{\footnotesize{on}}}_{k}\end{bmatrix}.$
(A.68)

We first show that $\Phi_{k}^{\text{E}}=H_{k}J_{k}+W_{k}$ satisfies (48). Note
that $\text{diag}(\nu_{k}^{\text{E}})W_{k}=\mathbf{0}$ since by construction
if the $i^{th}$ row of $W_{k}$ has a non zero entry then the $i^{th}$ diagonal
entry of $\text{diag}(\nu_{k}^{\text{E}})$ is zero. The product
$\text{diag}(\nu_{k}^{\text{E}})H_{k}$ is a diagonal matrix with with entries
of either zero or one. The zero entries also correspond to the zero entries of
$\nu_{k}^{\text{E}}$. In this case, the respective entry in $J_{k}$ is also
zero so that $\text{diag}(\nu_{k}^{\text{E}})H_{k}J_{k}=J_{k}$ as desired.

We now show that $\phi_{k}^{\text{E}}=H_{k}J_{k}+W_{k}\in\varPhi$. First
consider an arbitrary state indexed by $(\text{\footnotesize{off}},j,l)$ at
time $k$, if the corresponding value in
$\nu_{\text{\footnotesize{off}}}[\lambda^{j},l,k]>0$ then the two policy
values are defined as

$\displaystyle\frac{{\sf P}\left(m_{k+1}=\text{on},\ I_{k}=j,\
{\sf{L}}_{k}=l,\
m_{k}=\text{off}\right)}{\nu_{\text{\footnotesize{off}}}[\lambda^{j},l,k]}$
(A.69) $\displaystyle\frac{{\sf P}\left(m_{k+1}=\text{off},\ I_{k}=j,\
{\sf{L}}_{k}=l,\
m_{k}=\text{off}\right)}{\nu_{\text{\footnotesize{off}}}[\lambda^{j},l,k]}.$
(A.70)

If either of the above values are fixed in the constraint set $\Phi$, then the
constraint (49) will ensure this. Further, since we have that
$\nu_{k}^{\text{E}},J_{k}\in[0,1]$ and that
$(\nu_{k}^{\text{E}})^{T}=J_{k}\mathbb{1}$ this ensures that the above policy
values are within $[0,1]$ and sum to 1. The above argument is valid for any
pair of state values such that the corresponding value of $\nu_{k}^{\text{E}}$
is non-zero. If the corresponding value of
$\nu_{\text{\footnotesize{off}}}[\lambda^{j},l,k]=0$ and the policy
(conditioned on this state value) has a constraint, the first if case in the
definition of $w^{u,v}_{k}$ ensures this constraint. Further, the constraint
values must also be chosen to ensure the respective policy values are in
$[0,1]$ and sum to one. Lastly, if
$\nu_{\text{\footnotesize{off}}}[\lambda^{j},l,k]=0$ and there is no
constraint for the policy conditioned on this state value the second if case
in the definition of $w^{u,v}_{k}$ ensures the policy value sums to 1 and the
respective elements are in $[0,1]$. Thus $\Phi_{k}^{\text{E}}\in\varPhi$ for
all $k\in{\sf T}$.

###  A.4 Proof of Theorem 1

The proof structure is similar to the one in [2]. The idea is to exploit the
fact that: (i) $\nu_{k}^{\text{E}}$ is a decision variable for both
optimization problems (50) and (44) and (ii) the objective function is the
same for both problems and solely a function of the marginal
$\nu_{k}^{\text{E}}$. We rewrite these problem compactly below,

$\displaystyle\eta^{*}_{\text{CVX}}$
$\displaystyle=\min_{(\nu^{\text{E}},J)\in X}\eta(\nu^{\text{E}}),$ (A.71)
$\displaystyle\eta^{*}_{\text{NCVX}}$
$\displaystyle=\min_{(\nu^{\text{E}},\Phi^{\text{E}})\in
Y}\eta(\nu^{\text{E}}),$ (A.72)

where the sets $X$ and $Y$ collect all of the relevant constraints for the
problems. The variables $\nu^{\text{E}}$, $\Phi^{\text{E}}$, and $J$ are
concatenated over the considered finite time horizon and hence are not sub-
scripted by $k$. We proceed by showing that
$\eta^{*}_{\text{CVX}}\leq\eta^{*}_{\text{NCVX}}$ and
$\eta^{*}_{\text{NCVX}}\leq\eta^{*}_{\text{CVX}}$ to give the desired result.

####  A.4.1 $\eta^{*}_{\text{CVX}}\leq\eta^{*}_{\text{NCVX}}$

Pick any argument minimizer that achieves value $\eta^{*}_{\text{NCVX}}$ and
denote the pair as
$(\nu^{\text{E}}_{\text{NCVX}},\Phi^{\text{E}}_{\text{NCVX}})$. Trivially
construct $J$ through the relation (48) so that this constructed $J$ and
$\nu^{\text{E}}_{\text{NCVX}}$ (that is optimal for (44)) are also feasible
for (50), i.e., $(\nu^{\text{E}}_{\text{NCVX}},J)\in X$. This is since
$\mathbb{1}^{T}\text{diag}(\nu^{\text{E}}_{k})=\nu^{\text{E}}_{k}$ and
$\mathbb{1}=\Phi^{\text{E}}_{k}\mathbb{1}$. Hence we have that

$\displaystyle\eta^{*}_{\text{CVX}}$
$\displaystyle=\min_{(\nu^{\text{E}},J)\in
X}\eta(\nu^{\text{E}})\leq\eta(\nu^{\text{E}}_{\text{NCVX}})=\eta^{*}_{\text{NCVX}}$
(A.73)

since by definition $\eta^{*}_{\text{CVX}}$ is the minimum value over the set
of feasible solutions.

####  A.4.2 $\eta^{*}_{\text{NCVX}}\leq\eta^{*}_{\text{CVX}}$

We take a pair $(\nu^{\text{E}}_{\text{CVX}},J_{\text{CVX}})$ that achieve
optimal cost $\eta^{*}_{\text{CVX}}$ and construct a feasible solution for
(44), denoted $(\eta^{\text{E}}_{\text{NCVX}},\Phi^{\text{E}}_{\text{NCVX}})$,
as follows (for each $k$)

$\displaystyle\Phi^{\text{E}}_{k,\text{NCVX}}$
$\displaystyle=H_{k}J_{k}+W_{k},\quad\text{and}$ (A.74)
$\displaystyle\nu^{\text{E}}_{k,\text{NCVX}}$
$\displaystyle=\nu^{\text{E}}_{k,\text{CVX}}.$ (A.75)

Where $H_{k}$ and $W_{k}$ are defined in Lemma 5.2. This constructed solution
is then feasible for (44) as the constraint
$\Phi^{\text{E}}_{\text{NCVX}}\in\varPhi$ is part of the result in Lemma 5.2
and

$\displaystyle\nu^{\text{E}}_{k,\text{NCVX}}\Phi^{\text{E}}_{k,\text{NCVX}}G_{k}^{\text{E}}$
$\displaystyle=\nu^{\text{E}}_{k,\text{NCVX}}\big{(}H_{k}J_{k}+W_{k}\big{)}G_{k}^{\text{E}}$
(A.76)
$\displaystyle=\mathbb{1}^{T}J_{k}G_{k}^{\text{E}}=\nu^{\text{E}}_{k+1,\text{NCVX}}.$
(A.77)

The fact that
$\nu^{\text{E}}_{k,\text{NCVX}}\big{(}H_{k}J_{k}+W_{k}\big{)}=\mathbb{1}^{T}J_{k}$
is since $\nu^{\text{E}}_{k,\text{NCVX}}W_{k}=\mathbf{0}$ and
$\nu^{\text{E}}_{k,\text{NCVX}}H_{k}J_{k}=\mathbb{1}^{T}J_{k}$. The matrix
$W_{k}$ only has non zero entries for row indices where the index of the row
vector $\nu^{\text{E}}_{k,\text{NCVX}}$ is zero so that the resulting product
is the zero vector. The product $\nu^{\text{E}}_{k,\text{NCVX}}H_{k}$ is a
vector of ones and zeros, specifically, if the $i^{th}$ element of this vector
is zero then the entire $i^{th}$ column of the matrix $J_{k}$ will be the zero
vector. Thus the equivalence between
$\nu^{\text{E}}_{k,\text{NCVX}}\big{(}H_{k}J_{k}+W_{k}\big{)}$ and
$\mathbb{1}^{T}J_{k}$. Since the constructed solution is feasible we have that

$\displaystyle\eta^{*}_{\text{NCVX}}$
$\displaystyle=\min_{(\nu^{\text{E}},\Phi^{\text{E}})\in
Y}\eta(\nu^{\text{E}})\leq\eta(\nu^{\text{E}}_{\text{CVX}})=\eta^{*}_{\text{CVX}}$
(A.78)

since by definition $\eta^{*}_{\text{NCVX}}$ is the minimum value over the set
of feasible solutions.

##  Appendix B PDE discretization

We denote the $i^{th}$ CV as CV($i$) and further adopt the following
notational simplifications,

$\displaystyle\mu_{\text{\footnotesize{off}}}(\lambda^{i},t)\triangleq\mu_{\text{\footnotesize{off}}}(\lambda^{i}_{\text{\footnotesize{off}}},t),\quad\text{and}\quad\mu_{\text{\footnotesize{on}}}(\lambda^{i},t)\triangleq\mu_{\text{\footnotesize{on}}}(\lambda^{i}_{\text{\footnotesize{on}}},t).$

Highlighted red in Figure 2 are the two control volumes to assist in enforcing
boundary conditions that coincide with the thermostat policy (6). This is
discussed further in Appendix B.2 when the boundary conditions CVs are
discretized.

###  B.1 Internal CV’s

Consider the RHS of the pde (9) integrated over CV($i$):

$\displaystyle\int_{\text{CV(i)}}\bigg{(}\frac{\sigma^{2}}{2}\frac{\partial^{2}}{\partial\lambda^{2}}\big{(}\mu_{\text{\footnotesize{on}}}(\lambda,t)\big{)}-\frac{\partial}{\partial\lambda}\big{(}f_{\text{\footnotesize{on}}}(\lambda,t)\mu_{\text{\footnotesize{on}}}(\lambda,t)\big{)}\bigg{)}d\lambda$
$\displaystyle=\bigg{(}\frac{\sigma^{2}}{2}\frac{\partial}{\partial\lambda}\mu_{\text{\footnotesize{on}}}(\lambda,t)-f_{\text{\footnotesize{on}}}(\lambda,t)\mu_{\text{\footnotesize{on}}}(\lambda,t)\bigg{)}\bigg{|}_{\lambda^{i,-}}^{\lambda^{i,+}},$
(B.79)

where equality is by the divergence theorem [27]. Note, the points
$\lambda^{i,-}$ and $\lambda^{i,+}$ are not control volume variables, but
rather the boundaries of a single control volume. Hence, quantities in (B.79)
need to be approximated in terms of the nodal points of the neighboring
control volumes. The approximations for the partial derivative are,

$\displaystyle\frac{\partial}{\partial\lambda}\mu_{\text{\footnotesize{on}}}(\lambda^{i,+},t)$
$\displaystyle\approx\frac{\mu_{\text{\footnotesize{on}}}(\lambda^{i+1},t)-\mu_{\text{\footnotesize{on}}}(\lambda^{i},t)}{\Delta\lambda},\quad\text{and}$
(B.80)
$\displaystyle\frac{\partial}{\partial\lambda}\mu_{\text{\footnotesize{on}}}(\lambda^{i,-},t)$
$\displaystyle\approx\frac{\mu_{\text{\footnotesize{on}}}(\lambda^{i},t)-\mu_{\text{\footnotesize{on}}}(\lambda^{i-1},t)}{\Delta\lambda}.$
(B.81)

For the integrated convective term, we use the so-called upwind scheme [27].
This scheme elects the FVM equivalent of a forward or backward difference
based on the sign of the convective velocity $f_{\text{on}}(\lambda,t)$. By
assumption A.2, $f_{\text{\footnotesize{on}}}(\lambda,t)\leq 0$ and the upwind
scheme prescribes:

$\displaystyle
f_{\text{\footnotesize{on}}}(\lambda^{i,-},t)\mu_{\text{\footnotesize{on}}}(\lambda^{i,-},t)$
$\displaystyle=f_{\text{\footnotesize{on}}}(\lambda^{i,-},t)\mu_{\text{\footnotesize{on}}}(\lambda^{i},t),\quad\text{and}$
$\displaystyle
f_{\text{\footnotesize{on}}}(\lambda^{i,+},t)\mu_{\text{\footnotesize{on}}}(\lambda^{i,+},t)$
$\displaystyle=f_{\text{\footnotesize{on}}}(\lambda^{i,+},t)\mu_{\text{\footnotesize{on}}}(\lambda^{i+1},t).$
(B.82)

When the TCL is off $f_{\text{\footnotesize{off}}}(\lambda,t)\geq 0$ (also by
Assumption A.2) the upwind scheme prescribes:

$\displaystyle
f_{\text{\footnotesize{off}}}(\lambda^{i,-},t)\mu_{\text{\footnotesize{off}}}(\lambda^{i,-},t)$
$\displaystyle=f_{\text{\footnotesize{off}}}(\lambda^{i,-},t)\mu_{\text{\footnotesize{off}}}(\lambda^{i-1},t),\
\text{and}$ $\displaystyle
f_{\text{\footnotesize{off}}}(\lambda^{i,+},t)\mu_{\text{\footnotesize{off}}}(\lambda^{i,+},t)$
$\displaystyle=f_{\text{\footnotesize{off}}}(\lambda^{i,+},t)\mu_{\text{\footnotesize{off}}}(\lambda^{i},t).$
(B.83)

Now returning to the discretization of the PDE (9) over an arbitrary internal
CV. We approximate the LHS of (9) integrated over the control volume as,

$\displaystyle\int_{\text{CV(i)}}\frac{\partial}{\partial
t}\mu_{\text{\footnotesize{on}}}(\lambda,t)d\lambda\approx\frac{d}{dt}\mu_{\text{\footnotesize{on}}}(\lambda^{i},t)\Delta\lambda=\frac{d}{dt}\nu_{\text{\footnotesize{on}}}(\lambda^{i},t),$

where we have defined

$\displaystyle\nu_{\text{\footnotesize{on}}}(\lambda^{i},t)\triangleq\mu_{\text{\footnotesize{on}}}(\lambda^{i},t)\Delta\lambda.$
(B.84)

Now, denote the following

$\displaystyle
D\triangleq\frac{\sigma^{2}}{(\Delta\lambda)^{2}},\quad\text{and}\quad
F_{\text{\footnotesize{on}}}^{i}(t)\triangleq\frac{f_{\text{\footnotesize{on}}}(\lambda^{i},t)}{\Delta\lambda},$
(B.85)

where the quantities $F_{\text{\footnotesize{off}}}^{i}(t)$,
$F_{\text{\footnotesize{on}}}^{i,+}(t)/F_{\text{\footnotesize{off}}}^{i,+}(t)$,
and
$F_{\text{\footnotesize{on}}}^{i,-}(t)/F_{\text{\footnotesize{off}}}^{i,-}(t)$
are defined similarly to $F_{\text{\footnotesize{on}}}^{i}(t)$, e.g.,
$F_{\text{\footnotesize{off}}}^{i,+}(t)\triangleq
f_{\text{\footnotesize{off}}}(\lambda^{i,+},t)/\Delta\lambda$. Now equating
the approximation of the RHS (9) with the approximation of the LHS of (9) we
have,

$\displaystyle\frac{d}{dt}\nu_{\text{\footnotesize{on}}}(\lambda^{i},t)$
$\displaystyle=\Big{(}F_{\text{\footnotesize{on}}}^{i,-}(t)-D\Big{)}\nu_{\text{\footnotesize{on}}}(\lambda^{i},t)+\frac{D}{2}\nu_{\text{\footnotesize{on}}}(\lambda^{i-1},t)$
$\displaystyle+\Big{(}\frac{D}{2}-F_{\text{\footnotesize{on}}}^{i,+}(t)\Big{)}\nu_{\text{\footnotesize{on}}}(\lambda^{i+1},t).$
(B.86)

The spatial discretization for the PDE (10) is similar and yields,

$\displaystyle\frac{d}{dt}\nu_{\text{\footnotesize{off}}}(\lambda^{i},t)$
$\displaystyle=\frac{D}{2}\nu_{\text{\footnotesize{off}}}(\lambda^{i+1},t)-\Big{(}F_{\text{\footnotesize{off}}}^{i,+}(t)+D\Big{)}\nu_{\text{\footnotesize{off}}}(\lambda^{i},t)$
$\displaystyle+\Big{(}\frac{D}{2}+F_{\text{\footnotesize{off}}}^{i,-}(t)\Big{)}\nu_{\text{\footnotesize{off}}}(\lambda^{i-1},t),$
(B.87)

where
$\nu_{\text{\footnotesize{off}}}(\lambda^{i},t)\triangleq\mu_{\text{\footnotesize{off}}}(\lambda^{i},t)\Delta\lambda.$

###  B.2 Boundary CV’s

The boundary CVs are the CVs associated with the nodal values:
$\lambda_{\text{\footnotesize{on}}}^{1}$,
$\lambda_{\text{\footnotesize{on}}}^{q}$,
$\lambda_{\text{\footnotesize{on}}}^{N}$,
$\lambda_{\text{\footnotesize{off}}}^{1}$,
$\lambda_{\text{\footnotesize{off}}}^{m}$, and
$\lambda_{\text{\footnotesize{off}}}^{N}$. The superscript, for example the
integer $q$ in $\lambda_{\text{\footnotesize{on}}}^{q}$ represents the CV
index. All boundary CVs can be seen in Figure 2. Discretization of the
boundary CVs requires care for atleast two reasons. First, this is typically
where one introduces the BCs of the PDE into the numerical approximation.
Secondly, on finite domains the endpoints present challenges, for example,
there is no variable $\mu_{\text{\footnotesize{on}}}(\lambda^{N+1},t)$ for
computation of the derivative values for node
$\lambda^{N}_{\text{\footnotesize{on}}}$.

The BC’s for the coupled PDEs (9)-(10) are [21]:

Absorbing Boundaries:
$\displaystyle\qquad\qquad\qquad\mu_{\text{\footnotesize{on}}}(\lambda^{\text{min}},t)=\mu_{\text{\footnotesize{off}}}(\lambda^{\text{max}},t)=0.$
(B.88) Conditions at Infinity:
$\displaystyle\qquad\qquad\qquad\mu_{\text{\footnotesize{on}}}(+\infty,t)=\mu_{\text{\footnotesize{off}}}(-\infty,t)=0.$
(B.89) Conservation of Probability:
$\displaystyle\frac{\partial}{\partial\lambda}\bigg{[}\mu_{\text{\footnotesize{on}}}(\lambda^{q,-},t)-\mu_{\text{\footnotesize{on}}}(\lambda^{q-1,+},t)-\mu_{\text{\footnotesize{off}}}(\lambda^{N-1,+},t)\bigg{]}=0.$
(B.90)
$\displaystyle\frac{\partial}{\partial\lambda}\bigg{[}\mu_{\text{\footnotesize{off}}}(\lambda^{m,+},t)-\mu_{\text{\footnotesize{on}}}(\lambda^{2,-},t)-\mu_{\text{\footnotesize{off}}}(\lambda^{m+1,-},t)\bigg{]}=0.$
(B.91) Continuity:
$\displaystyle\qquad\qquad\qquad\mu_{\text{\footnotesize{on}}}(\lambda^{q,-},t)=\mu_{\text{\footnotesize{on}}}(\lambda^{q-1,+},t).$
(B.92)
$\displaystyle\qquad\qquad\qquad\mu_{\text{\footnotesize{off}}}(\lambda^{m,+},t)=\mu_{\text{\footnotesize{off}}}(\lambda^{m+1,-},t).$
(B.93)

As we will see, implementation of some of the above conditions will require a
bit of care. However, some are quite trivial to enforce. For example, by
default, the continuity conditions (B.92) and (B.93) are satisfied due to our
choice of CV structure, since, for example, for any $i$ we have
$\lambda_{\text{\footnotesize{off}}}^{i,-}=\lambda_{\text{\footnotesize{off}}}^{i-1,+}$
and
$\lambda_{\text{\footnotesize{off}}}^{i,+}=\lambda_{\text{\footnotesize{off}}}^{i+1,-}$.

Now focusing on the conditions at infinity BC (B.89), we enforce instead the
following conditions:

$\displaystyle\frac{\partial}{\partial\lambda}\mu_{\text{\footnotesize{off}}}(\lambda^{1,-},t)=0,\quad\text{and}\quad\frac{\partial}{\partial\lambda}\mu_{\text{\footnotesize{on}}}(\lambda^{N,+},t)=0.$
(B.94)

Our computational domain cannot extend to infinity, where the BC (B.89) is
required to hold, but the temperature values
$\lambda_{\text{\footnotesize{off}}}^{1}$ and
$\lambda_{\text{\footnotesize{on}}}^{N}$ are quite far away from the deadband
and so the density here will be near zero.

Now, consider the spatial discretization of the CVs associated with the BC at
infinity. First considering the CV associated with the temperature
$\lambda^{1}_{\text{\footnotesize{off}}}$, we have that the differential
equation is

$\displaystyle\frac{d}{dt}\nu_{\text{\footnotesize{off}}}(\lambda^{1},t)$
$\displaystyle=\Big{(}-F_{\text{\footnotesize{off}}}^{1,+}(t)-\frac{D}{2}\Big{)}\nu_{\text{\footnotesize{off}}}(\lambda^{1},t)$
(B.95)
$\displaystyle+\Big{(}\frac{D}{2}+F_{\text{\footnotesize{off}}}^{2,-}(t)\Big{)}\nu_{\text{\footnotesize{off}}}(\lambda^{2},t).$

Considering the CV associated with the temperature $\lambda^{N}_{\text{on}}$,
we have

$\displaystyle\frac{d}{dt}\nu_{\text{\footnotesize{on}}}(\lambda^{N},t)$
$\displaystyle=\Big{(}F_{\text{\footnotesize{on}}}^{N,+}(t)-\frac{D}{2}\Big{)}\nu_{\text{\footnotesize{on}}}(\lambda^{N},t)$
(B.96)
$\displaystyle+\Big{(}\frac{D}{2}-F_{\text{\footnotesize{on}}}^{N,+}(t)\Big{)}\nu_{\text{\footnotesize{on}}}(\lambda^{N-1},t).$

In the above we make the assumption that
$\nu_{\text{\footnotesize{off}}}(\lambda^{1,-}-\Delta\lambda,t)=0$ and
$\nu_{\text{\footnotesize{on}}}(\lambda^{N,+}+\Delta\lambda,t)=0$.

Now focus on the absorbing boundary (B.88) and conservation of probability
(B.90)-(B.91) boundary conditions. These BCs have the following meaning. The
condition (B.88) clamps the density at the end of the deadband to zero. BC
(B.90) reads: the net-flux across the temperature value
$\lambda^{q}_{\text{\footnotesize{on}}}$ is equal to the flux of density going
from off to on. In order to enforce both (B.90) and (B.91) we will model the
flux of density due to the thermostat control policy as a source/sink. Before
doing this, we mention some issues with enforcing the BC (B.88).

A TCL’s temperature trajectory will not satisfy the BC (B.88) since to switch
its mode the TCL’s temperature sensor will have to register a value outside
the deadband. Therfore, we introduce two additional CV’s associated with the
temperatures $\lambda^{1}_{\text{\footnotesize{on}}}$ and
$\lambda^{N}_{\text{\footnotesize{off}}}$, which are shown in red in Figure 2.
We then transfer the BC (B.88) to one on the added CVs, which becomes:

$\displaystyle\mu_{\text{\footnotesize{on}}}(\lambda^{1,-},t)=\mu_{\text{\footnotesize{off}}}(\lambda^{N,+},t)$
$\displaystyle=0.$ (B.97)

As mentioned, to enforce the conservation of probability BC we use a
source/sink type argument, which we also enforce on the added CVs. To see what
we mean by source/sink argument, consider the following: some rate of density
is transferred out of the CV $\lambda^{N}_{\text{\footnotesize{off}}}$ and
into the CV $\lambda^{q}_{\text{\footnotesize{on}}}$ (as depicted in Figure 2)
due to thermostatic control. We model the sink as simply
$-\nu_{\text{\footnotesize{off}}}(\lambda^{N},t)$. The rate of the sink is
then given as $-\gamma\nu_{\text{\footnotesize{off}}}(\lambda^{N},t)$, where
$\gamma>0$ is a modeling choice and a constant of appropriate units that
describes the discharge rate. We shortly give insight on how to select a value
for $\gamma$. Now discretizing the CV corresponding to the nodal value
$\lambda^{N}_{\text{\footnotesize{off}}}$ subject to the BC (B.97) and the
sink $-\nu_{\text{\footnotesize{off}}}(\lambda^{N},t)$ we obtain,

$\displaystyle\frac{d}{dt}\nu_{\text{\footnotesize{off}}}(\lambda^{N},t)$
$\displaystyle=\Big{(}\frac{D}{2}+F_{\text{\footnotesize{off}}}^{N,-}(t)\Big{)}\nu_{\text{\footnotesize{off}}}(\lambda^{N-1},t)$
(B.98) $\displaystyle-\alpha\nu_{\text{\footnotesize{off}}}(\lambda^{N},t),$

where $\alpha\triangleq\big{(}\gamma+D\big{)}$. In obtaining the above, we
have made the reasonable assumption that
$\nu_{\text{\footnotesize{off}}}(\lambda^{N,+}+\Delta\lambda,t)=0$. The
quantity $\alpha\nu_{\text{\footnotesize{off}}}(\lambda^{N},t)$ represents the
rate of change of density from the CV
$\lambda^{N}_{\text{\footnotesize{off}}}$ to the CV
$\lambda^{q}_{\text{\footnotesize{on}}}$, as depicted in Figure 2.
Consequently, to conserve probability, we must add this quantity as a source
to the ode for the CV $\lambda^{q}_{\text{\footnotesize{on}}}$, i.e.,

$\displaystyle\frac{d}{dt}\nu_{\text{\footnotesize{on}}}(\lambda^{q},t)$
$\displaystyle=\dots+\alpha\nu_{\text{\footnotesize{off}}}(\lambda^{N},t).$
(B.99)

The dots in equation (B.99) represent the portion of the dynamics for the
standard internal CV (i.e., the RHS of (B.86)) for the temperature node
$\lambda^{q}_{\text{\footnotesize{on}}}$. A similar argument is used for the
BC (B.91) with the CV’s $\lambda^{1}_{\text{\footnotesize{on}}}$ and
$\lambda^{m}_{\text{\footnotesize{off}}}$, and the corresponding differential
equations are,

$\displaystyle\frac{d}{dt}\nu_{\text{\footnotesize{on}}}(\lambda^{1},t)$
$\displaystyle=\Big{(}\frac{D}{2}-F_{\text{\footnotesize{on}}}^{1,+}(t)\Big{)}\nu_{\text{\footnotesize{on}}}(\lambda^{2},t)-\alpha\nu_{\text{\footnotesize{on}}}(\lambda^{1},t),$
(B.100)
$\displaystyle\frac{d}{dt}\nu_{\text{\footnotesize{off}}}(\lambda^{m},t)$
$\displaystyle=\dots+\alpha\nu_{\text{\footnotesize{on}}}(\lambda^{1},t).$
(B.101)

To better understand the role of $\gamma$ consider the following example.
Electing $\gamma$ in the above so that $\alpha=(\Delta t)^{-1}$, where $\Delta
t$ is a time increment, has the following interpretation: all mass starting in
state $\nu_{\text{\footnotesize{off}}}(\lambda^{N},\cdot)$ at time $t$ is
transferred out by time $t+\Delta t$ into the state
$\nu_{\text{\footnotesize{on}}}(\lambda^{q},\cdot)$.

####  B.2.1 Additional conditions

Two additional conditions are enforced, namely that once mass is transferred
to the nodes $\lambda^{N}_{\text{\footnotesize{off}}}$ or
$\lambda^{1}_{\text{\footnotesize{on}}}$ it cannot “travel backwards.” For
example, mass is transferred from $\lambda^{N}_{\text{\footnotesize{off}}}$
entirely to the corresponding on temperature bin and no mass is transferred
backwards to $\lambda^{N-1}_{\text{\footnotesize{off}}}$. This corresponds to
setting: (i) the coefficient on
$\nu_{\text{\footnotesize{off}}}(\lambda^{N},t)$ in the ode for
$\nu_{\text{\footnotesize{off}}}(\lambda^{N-1},t)$ to zero and (ii) the
coefficient on $\nu_{\text{\footnotesize{on}}}(\lambda^{1},t)$ in the ode for
$\nu_{\text{\footnotesize{on}}}(\lambda^{2},t)$ to zero.

###  B.3 Overall system

Now, combining the odes–(B.86) and (B.87) for all of the internal CVs and
(B.95), (B.96), (B.98), (B.99), (B.100), (B.101) for the BC CVs–we obtain the
linear time varying system,

$\displaystyle\frac{d}{dt}\nu(t)=\nu(t)A(t).$ (B.102)