# Power Tracking Control of Heterogeneous Populations of TCLs with Partially
Measured States

Zhenhe Zhang1<EMAIL_ADDRESS>Jun Zheng1,2<EMAIL_ADDRESS>Guchuan Zhu1<EMAIL_ADDRESS>1Department of Electrical Engineering,
Polytechnique Montréal, P.O. Box 6079, Station Centre-Ville, Montreal, QC,
Canada H3T 1J4
2School of Mathematics, Southwest Jiaotong University, Chengdu 611756,
Sichuan, China

###### Abstract

This paper presents a new aggregate power tracking control scheme for
populations of thermostatically controlled loads (TCLs). The control design is
performed in the framework of partial differential equations (PDEs) based on a
late-lumping procedure without truncating the infinite-dimensional model
describing the dynamics of the TCL population. An input-output linearization
control scheme, which is independent of system parameters and uses only
partial state measurement, is derived, and a sliding model-like control is
applied to achieve finite-time input-to-state stability for tracking error
dynamics. Such a control strategy can ensure robust performance in the
presence of modeling uncertainties, while considerably reducing the
communication burden in large scale distributed systems similar to that
considered in the present work. A rigorous analysis of the closed-loop
stability of the underlying PDE system was conducted, which guaranteed the
validity of the developed control scheme. Simulation studies were performed
while considering two TCL populations with a significant difference in their
size, and the results show that the developed control scheme performs well in
both cases, thereby confirming the effectiveness of the proposed solution.

###### keywords:

Aggregate power tracking control, finite-time input-to-state stability, input-
output linearization, partial differential equations, thermostatically
controlled loads.

, ,

##  1 Introduction

In the context of today¡¯s smart grids, it is widely recognized that demand
response (DR) programs have great potential in dealing with ongoing demands,
while enhancing the energy efficiency and resilience of the power grid [4, 28,
7, 40]. As a promising demand-response enabled resource, thermostatically
controlled loads (TCLs), such as air conditioners (ACs), space heating
devices, refrigerators, and water heaters, are attracting increasing
attention. Although a single TCL unit has very limited power regulation
capability, ensembles of a large number of TCLs, when managed in an orderly
and controllable manner, can have a significant impact on the entire power
grid [26, 38, 12]. It has been shown that a large TCL population can be
managed to support demand response tasks, including peak load shaving and load
following [34, 9, 25], and to provide ancillary services, such as primary or
secondary frequency controls [36, 21, 20, 33].

The present work focuses on load tracking control, which allows the aggregate
power of a TCL population to follow a desired consumption profile. The control
design is based on a model of the dynamics of the TCL population described by
partial differential equations (PDEs). Specifically, we consider a set of TCLs
in which the dynamics of every individual device are modeled by a lumped
stochastic hybrid system (SHS) operated through thermostat-based deadband
control. The aggregate dynamics of such a TCL population can be modeled by two
coupled Fokker-Planck equations (see, e.g, [19, 2, 37]) describing the
evolution of the probability distribution of TCLs in the ON and OFF states
over the temperature. Note that the same form of PDE-based models can also be
derived by assuming that the dynamics of individual TCLs are described by
deterministic systems while considering population heterogeneity [1, 5, 23].

Another widely adopted method to build the aggregate dynamical model of TCL
populations is to divide a fixed range of temperatures into several segments,
called state-bins, each of which is associated with the number of TCLs with
their temperature fitting in this bin. The dynamics of state-bin transactions
can be described by a Markov chain (see, e.g., [35, 27, 14, 18, 29]) or state
queue (see, e.g., [17, 32]), which leads to finite-dimensional state-space
models. It is worth noting that discretizing a PDE with respect to (w.r.t.)
the space variable (temperature) also leads to a finite-dimensional state-
space model. However, as the considered Fokker-Planck equation is a semi-
linear time-varying PDE, its discretization results in a finite-dimensional
nonlinear time-varying system. Consequently, a model described by the linear
time invariant (LTI) system, which is the most used state-bin model in the
existing literature, may be equivalent to that derived from PDEs only locally
around particular equilibrium points and operational conditions (e.g.,
temperature set-point, ambient temperature, deadband), even with a variety of
extensions. Therefore, the PDE provides a more generic framework for modeling
the aggregate dynamics of TCL populations, which allows handling nonlinearity,
time-varying operational conditions, and parametric uncertainties with often
very simple control algorithms. However, the PDE control system design
procedure generally involves more complex mathematical analysis and is more
challenging.

The main objective of TCL population control is to manipulate the total power
consumption of the entire population, which can be achieved by changing the
temperature set-point, moving the deadband, or interfering with the
probability distributions of the TCLs via forced switches (see, e.g., [36, 1,
2, 18, 30, 39]). Because a TCL population usually contains a large number of
units that may spread over a large geographical area, only decentralized or
distributed schemes are applicable control strategies. In fact, a remarkable
amount of work on the control of TCL populations has been reported in the
literature, and the majority of the proposed solutions are based on lumped
models by applying optimization theory and optimal control techniques, in
particular model predictive control (see, e.g, [36, 20, 35, 27, 14, 18, 29,
17, 32, 1, 30]). It should be noted that, owing the nature of the considered
problem, control schemes requiring the state measurement of the entire
population in real-time are practically infeasible (see, e.g., [31] and the
references therein). This problem can be addressed using state observers [20,
22]. Nevertheless, it is still very challenging to assess the performance of
model-based state estimation algorithms because it depends heavily on the
accuracy of the system parameters.

The load tracking control algorithm developed in the present work is a
decentralized scheme in which the rates for set-point temperature adjustment
generated by a central unit are broadcast to the TCLs over the population.
Emphasis is placed on solving issues arising in practical applications,
particularly communication restrictions and modeling uncertainties for large
scale TCL populations. The control system design is carried out in the
framework of PDE-based modeling and control techniques. It should be noted
that the two basic paradigms in PDE control system design and implementation,
namely early-lumping and late-lumping procedures, have all been applied to the
control of the coupled Fokker-Planck equations associated with TCL
populations. The early-lumping method discretizes the underlying PDEs to
obtain a lumped model, and then applies the techniques for finite-dimensional
control system design [1, 2, 5, 23, 30]. In contrast, with the late-lumping
method, the controller is designed using the PDE model and then discretized
for implementation [6, 39]. A significant advantage of the late-lumping method
is that it can preserve the essential properties of the PDE model and no
approximation is required in the control design. However, some issues remain
open. More specifically, the schemes developed in [39] and [6] are based on
input-output linearization by state feedback control, which may incur a
communication burden. In addition, these control schemes require an accurate
knowledge of the system parameters, for example, the diffusion coefficient in
Fokker-Planck equations, which are not easy to determine from both theoretical
and practical viewpoints considering the nature of the problem under
investigation. Finally, although taking a weighted power load as the system
output proposed in [39] can avoid the controllability issue introduced by the
use of the total power load of the in-band TCLs as the system output in [6],
such a choice lacks physical interpretation and is unsuitable for practical
operation.

In this paper, we developed a new control algorithm based on the input-output
linearization technique, which results in a system composed of finite-
dimensional input-output dynamics and infinite-dimensional internal dynamics.
The control design amounts then to finding a robust closed-loop control law
that stabilizes the finite-dimensional input-output dynamics while
guaranteeing the stability of the infinite-dimensional internal dynamics.
Specifically:

  * •

A new system output for power tracking control is proposed that can guarantee
the controllability of the input-output dynamics.

  * •

A linearization control law, which is independent of system parameters, e.g.,
the diffusion coefficient, while requiring only knowledge of the states of
TCLs near the deadband boundaries, is derived.

  * •

To tackle modeling uncertainties while making the control scheme
computationally tractable, a sliding model-like tracking control scheme that
can achieve finite-time input-to-state stability (FTISS)[16, 8], is designed.

  * •

The non-negativeness of the solution to the Fokker-Planck equations under the
developed control law and other properties required to ensure closed-loop
stability are rigorously validated.

The main contribution of the present work lies in the simplicity, scalability,
and applicability of the control strategy developed under a generic framework.
In addition, it is worth noting that as the developed control algorithm
requires only measuring the state of the TCLs on the end-points of the
deadband, TCLs need to notify their state only when switching occurs. Because
the cyclic rate of the TCLs is much slower than the controller sampling rate,
the communication burden can be significantly reduced. Obviously, it is very
difficult for state feedback control schemes based on lumped aggregate models
to achieve such features, which is critical for practical implementations.

The remainder of this paper is organized as follows. Section 2 introduces the
notations used in the study and preliminaries on FTISS. Section 3 presents the
first-order equivalent thermal parameter (ETP) model for a single TCL unit and
the coupled Fokker-Planck model for the aggregate dynamics of the TCL
population. Section 4 presents the power tracking control design and closed-
loop stability analysis. The experimental validation of the developed control
strategy and the simulation results are reported in Section 5, followed by
concluding remarks in Section 6. Finally, the proof of one of the main
theoretical result is presented in the appendix.

##  2 Notations and preliminaries

###  2.1 Notations

Let $\mathbb{R}:=(-\infty,+\infty),$ $\mathbb{R}_{\geq 0}:=[0,+\infty)$,
$\mathbb{R}_{>0}:=(0,+\infty)$, and $\mathbb{R}_{\leq 0}:=(-\infty,0]$. Denote
by $\partial_{s}f$ the derivative of the function $f$ w.r.t. argument $s$.
Note that, for notation simplicity, we may omit the arguments of functions if
there is no ambiguity.

By convention, we denote by $|\cdot|$ the module of a function. For positive
integers $m,n$ and a given (open or closed) domain
$\Omega\subset\mathbb{R}^{n}$, let
$L^{\infty}(\Omega;\mathbb{R}^{m}):=\\{\phi:\Omega\rightarrow\mathbb{R}^{m}|~{}\phi$
is measurable in $\Omega$ and satisfies $\text{ess
sup}_{s\in\Omega}|\phi(s)|<+\infty\\}$. For $\phi\in
L^{\infty}(\Omega;\mathbb{R}^{m})$, the norm of $\phi$ is defined by
$\|\phi\|_{L^{\infty}(\Omega)}:=\text{ess sup}_{s\in\Omega}|\phi(s)|$. Let
$L^{\infty}_{loc}(\Omega;\mathbb{R}^{m}):=\\{\phi:\Omega\rightarrow\mathbb{R}^{m}|~{}\phi\in
L^{\infty}(\Omega^{\prime};\mathbb{R}^{m})$ for any
$\Omega^{\prime}\subsetneqq\Omega\\}$

For given (open or closed) domains
$\Omega_{1},\Omega_{2}\subset\mathbb{R}^{n}$ and
$\Omega_{3}\subset\mathbb{R}$, let
$C\left(\Omega_{1};\Omega_{3}\right):=C^{0}\left(\Omega_{1};\Omega_{3}\right):=\\{\phi:\Omega_{1}\rightarrow\Omega_{3}|~{}\phi$
is continuous w.r.t. its all augments in $\Omega_{1}\\}$. For positive
integers $i,j$, let
$C^{i}\left(\Omega_{1};\Omega_{3}\right):=\\{\phi:\Omega_{1}\rightarrow\Omega_{3}|~{}\phi$
has continuous derivatives up to order $i$ w.r.t. its all augments in
$\Omega_{1}\\}$, and
$C^{i,j}\left(\Omega_{1}\times\Omega_{2};\Omega_{3}\right):=\\{\phi:\Omega_{1}\times\Omega_{2}\rightarrow\Omega_{3}|~{}\phi$
has continuous derivatives up to order $i$ w.r.t. its augments in $\Omega_{1}$
and up to order $j$ w.r.t. its augments in $\Omega_{2}\\}$. In particular, if
$\Omega_{3}=\mathbb{R}$, we denote
$C\left(\Omega_{1}\right):=C^{0}\left(\Omega_{1};\mathbb{R}\right)$ and
$C^{i}\left(\Omega_{1}\right):=C^{i}\left(\Omega_{1};\mathbb{R}\right)$ for
$i>0$.

As in [16] and [11], we define the following sets of comparison functions. Let
$\mathcal{K}:=\\{\vartheta:\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}_{\geq 0}|\
\vartheta(0)=0,\vartheta$ is continuous, strictly increasing$\\}$;
$\mathcal{L}:=\\{\vartheta:\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}_{\geq 0}|\
\vartheta$ is continuous, strictly decreasing,
$\lim_{s\rightarrow+\infty}\vartheta(s)=0\\}$;
$\mathcal{KL}:=\\{\beta:\mathbb{R}_{\geq 0}\times\mathbb{R}_{\geq
0}\rightarrow\mathbb{R}_{\geq 0}|\ {\beta(\cdot,t)}\in\mathcal{K},\forall
t\in\mathbb{R}_{\geq 0}$, and $\beta(s,\cdot)\in\mathcal{L},\forall
s\in{\mathbb{R}_{+}}\\}$; $\mathcal{K}_{\infty}:=\\{\vartheta:\mathbb{R}_{\geq
0}\rightarrow\mathbb{R}_{\geq
0}|~{}\vartheta\in\mathcal{K}~{}\text{and}~{}\lim_{s\rightarrow+\infty}\vartheta(s)=+\infty\\}$;
$\mathcal{G}\mathcal{K}\mathcal{L}:=\\{\beta:\mathbb{R}_{\geq
0}\times\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}_{\geq
0}|~{}\beta(\cdot,0)\in\mathcal{K}$, and for each fixed $s\in\mathbb{R}_{>0}$
there exists $\widetilde{T}(s)\in\mathbb{R}_{\geq 0}$ such that $\beta(s,t)=0$
for all $t\geq\widetilde{T}(s)\\}$.

###  2.2  Finite-time input-to-state stability of finite dimensional systems

Consider the following nonlinear system

$\displaystyle\dot{z}(t)=$ $\displaystyle f(z(t),d(t)),\ \ \forall
t\in\mathbb{R}_{\geq 0},$ (1a) $\displaystyle z(0)=$ $\displaystyle z_{0},$
(1b)

where $z:=[z_{1},z_{2},...,z_{n}]^{T}\in\mathbb{R}^{n}$ is the state,
$z_{0}\in\mathbb{R}^{n}$ is the initial datum,
$d\in\mathcal{D}:=L^{\infty}_{\text{loc}}(\mathbb{R}_{\geq 0};\mathbb{R}^{m})$
is the input (disturbance) to the system,
$f:\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$ is a nonlinear
function that is continuous w.r.t. $(z,d)$, ensures the forward existence of
the system solutions, at least locally, and satisfies $f(0,0)=0$, and $m\geq
1$ and $n\geq 1$ are integers.

######  Definition 2.1.

System (1) is said to be finite-time input-to-state stable (FTISS) if there
exist functions $\vartheta\in\mathcal{K}$ and $\beta\in\mathcal{GKL}$ such
that for any $x_{0}\in\mathbb{R}^{n}$ and $d\in\mathcal{D}$ its trajectory
satisfies

$\displaystyle|z(t)|\leq\beta(|z_{0}|,t)+\vartheta(\|d\|_{{L^{\infty}(0,t)}}),\
\ \forall t\in\mathbb{R}_{\geq 0}.$ (2)

######  Remark 2.1.

Note that FTISS is defined in a similar way to the definition of input-to-
state stability (ISS) in [11, Chapter 4] via the norm of $d$ over the interval
$(0,t)$ rather than $(0,+\infty)$. Thus, the FTISS presented here is a refined
notion of the one introduced in [16, 8], where the second term in the right-
hand side of (2) is under the form $\vartheta(\|d\|_{L^{\infty}(0,+\infty)})$,
which describes the influence of the global bounds of $d$ instead of the
bounds of $d$ over the finite time interval $(0,t)$.

######  Definition 2.2.

A continuously differentiable function
$V:\mathbb{R}^{n}\rightarrow\mathbb{R}_{\geq 0}$ is said to be an FTISS
Lyapunov function for system (1) if there exist functions
$\mu_{1},\mu_{2}\in\mathcal{K}_{\infty}$, $\chi\in\mathcal{K}$ and constants
$c>0$ and $\theta\in(0,1)$ such that for all $x\in\mathbb{R}^{n}$ and all
$d\in\mathcal{D}$ it holds that

$\displaystyle{\mu_{1}({|x|})\leq V(x)\leq\mu_{2}({|x|})},$
$\displaystyle|z|\geq\chi(|d|)\Rightarrow DV(z)\cdot f(z,d)\leq-
cV^{\theta}(z),$

where $DV(z):=\begin{bmatrix}\frac{\partial V}{\partial
z_{1}},\ldots,\frac{\partial V}{\partial z_{n}}\end{bmatrix}$ .

The following Lyapunov-like lemma gives a sufficient condition for the FTISS.

######  Lemma 2.1.

System (1) is FTISS if it admits a finite-time ISS Lyapunov function.

Proof. Setting $\mathcal{V}:=\\{z|V(z)\leq\mu_{2}(\chi(|d|))\\}$ in the proof
of [8, Theorem 1(a)], the lemma statement follows immediately.
$\hfill\blacksquare$

##  3 Mathematical model and problem specification

###  3.1 Dynamics of individual TCLs

In the present work, we focus on modeling the population of residential air
conditioners (ACs). While, its extension to other cooling and heating devices
is straightforward. We consider the case where all ACs are operated by
thermostats hence, every AC switches between the ON and OFF states whenever it
reaches the prescribed lower or upper temperature bounds. For simplicity, we
ignore the solar irradiation and internal heat gains and assume that the ACs
operate at a fixed frequency. Then, the dynamics of the indoor temperature,
denoted by $x$, for a representative load can be modeled by the following SHS
(see, e.g., [19, 2, 30]):

$\displaystyle\operatorname{d}\\!{x}(t)=\dfrac{1}{CR}\left(x_{a}(t)-x(t)-s(t)RP\right)\operatorname{d}\\!{t}+\sigma\operatorname{d}\\!{w}(t),$
(3)

where $x_{a}(t)$ is the ambient temperature, $R$, $C$, and $P$ are the thermal
resistance, capacitance, and power, respectively, and $s(t)$ is the switching
signal. In (3), $w(t)$ is a standard Wiener process, which, along with the
parameter $\sigma$, represents modeling uncertainties, such as unaccounted
heat loss or heat gain, parameter variations, and disturbances.

For a thermostat-controlled AC, the switching signal $s(t)$ takes a binary
value from $\\{0,1\\}$, representing the OFF and ON states. We consider a
hybrid control scheme, as shown in Fig. 1, in which the device always switches
at the endpoints of the deadband. In addition, forced switches at any moment,
denoted by $r(t)$, may also occur to alert the probability distributions of
the TCL population. Let $r(t)$ take a binary value from $\\{0,1\\}$, with 1
representing the occurrence of switching and 0 otherwise. Letting
$\underline{x}$ and $\overline{x}$ be the prescribed lower and upper
temperature bounds, respectively, the deadband control for an AC can then be
expressed as

$\displaystyle s(t)=\begin{cases}1,&\hbox{if~{}}x\geq\overline{x};\\\
0,&\hbox{if~{}}x\leq\underline{x};\\\ (s(t^{-})\wedge r(t))+(s(t^{-})\vee
r(t)),&\hbox{otherwise};\end{cases}$

where “$+$” is the one-bit binary addition with overflow. In addition, the
notations $(\cdot)^{-}$ and $(\cdot)^{+}$ denote the left and right limits of
the scalar variable, respectively. Note that different actions, such as random
switches to avoid power demand oscillations due to synchronization within a
TCL population, mechanisms for blocking the switches to protect the ACs, etc.,
can be integrated in the design of forced switching schemes.

Figure 1: Hybrid thermostat-based deadband control scheme.

###  3.2 Dynamics of aggregate TCL population

As mentioned previously, the dynamics of an aggregate TCL population can be
characterized by the evolution of the distributions of the TCLs over
temperature. When the number of TCLs in the population tends to be infinite,
this population can be modeled as a continuum whose temperature distribution
is governed by the coupled Fokker-Planck equations [19, 23, 5, 2].
Specifically, we denote by $f_{1}(x,t)$ and $f_{0}(x,t)$ the probability
density functions (PDFs) of the TCLs in the ON and OFF states at temperature
$x$ and time $t$, respectively. As illustrated in Fig. 2, we assume that all
the loads are confined in a fixed temperature range $(x_{L},x_{H})$ along all
possible operations, where $x_{L}$ and $x_{H}$ are constants, which is a
reasonable assumption for practical application. Moreover, owing to the nature
of thermostat-based control, there must be that $f_{1}(x,t)=0$ for all
$x\leq\underline{x}$ and $t\in\mathbb{R}_{>0}$, and that $f_{0}(x,t)=0$ for
all $x\geq\overline{x}$ and $t\in\mathbb{R}_{>0}$. Therefore, we can divide
the range $(x_{L},x_{H})$ into three segments:

$I_{a}:=(x_{L},\underline{x}),I_{b}:=(\underline{x},\overline{x}),I_{c}:=(\overline{x},x_{H}),$

which will be used in the upcoming study.

Figure 2: Illustration of probability density functions of a TCL population at
a given time.

Suppose that the dynamics of each load in the TCL population are described by
(3). Let further

$\displaystyle\alpha_{0}(x,t):=$
$\displaystyle\frac{1}{CR}\left(x_{a}(t)-x\right),$
$\displaystyle\alpha_{1}(x,t):=$
$\displaystyle\frac{1}{CR}\left({x_{a}(t)}-x-RP\right).$

The evolutions of $f_{0}(x,t)$ and $f_{1}(x,t)$ are governed by the following
coupled Fokker-Planck equations [19, 2, 30]:

$\displaystyle\partial_{t}f_{0}=$
$\displaystyle{\partial_{x}}\\!\left(\dfrac{\sigma^{2}}{2}{\partial_{x}f_{0}}-(\alpha_{0}-u)f_{0}\right)\
\text{in~{}}I_{a}\times\mathbb{R}_{>0},$ (4a)
$\displaystyle\partial_{t}f_{0}=$
$\displaystyle{\partial_{x}}\\!\left(\dfrac{\sigma^{2}}{2}{\partial_{x}f_{0}}-(\alpha_{0}-u)f_{0}\right)-g(f_{0},f_{1})\
\text{in~{}}I_{b}\times\mathbb{R}_{>0},$ (4b)
$\displaystyle\partial_{t}f_{1}=$
$\displaystyle{\partial_{x}}\\!\left(\dfrac{\sigma^{2}}{2}{\partial_{x}f_{1}}-(\alpha_{1}-u)f_{1}\right)+g(f_{0},f_{1})\
\text{in~{}}I_{b}\times\mathbb{R}_{>0},$ (4c)
$\displaystyle\partial_{t}f_{1}=$
$\displaystyle{\partial_{x}}\\!\left(\dfrac{\sigma^{2}}{2}{\partial_{x}f_{1}}-(\alpha_{1}-u)f_{1}\right)\
\text{in~{}}I_{c}\times\mathbb{R}_{>0},$ (4d)

where $g(f_{0},f_{1})$ represents the net probability flux due to the switches
occurring over segment $I_{b}$, that is, the so-called forced switches. Hence,
the signs of $g(f_{0},f_{1})$ in (4b) and (4c) should be opposite to each
other, which implies a mass conservation property as claimed in Theorem 4.3 in
Section 4.3. Note that (4b) and (4c) have a general form compared to that
given in [30] (see (19a) and (19b) of that paper), where an explicitly linear
function $g(f_{0},f_{1})$ was used to model a switching rate control scheme.

Following [30], we introduce the notation of probability flows
$\mathcal{F}_{i}$. When there is no additional flux from the forced switches,
i.e., $g=0$, $\mathcal{F}_{i}$ is the integral of the probability fluxes
$\partial_{t}f_{i}$ over the temperature ($x$-) coordinate:

$\mathcal{F}_{i}(x,t):=\dfrac{\sigma^{2}}{2}\partial_{x}f_{i}(x,t)-(\alpha_{i}(x,t)-u(t))f_{i}(x,t),i=0,1.$

The boundary conditions can then be written as

$\displaystyle\mathcal{F}_{0}(x_{L}^{+},t)=$ $\displaystyle 0,\ \ \forall
t\in\mathbb{R}_{>0},$ (5a)
$\displaystyle\mathcal{F}_{0}(\underline{x}^{-},t)=$
$\displaystyle\mathcal{F}_{0}(\underline{x}^{+},t)+\mathcal{F}_{1}(\underline{x}^{+},t),\
\ \forall t\in\mathbb{R}_{>0},$ (5b) $\displaystyle
f_{0}(\underline{x}^{-},t)=$ $\displaystyle f_{0}(\underline{x}^{+},t),\ \
\forall t\in\mathbb{R}_{>0},$ (5c) $\displaystyle f_{0}(\overline{x}^{-},t)=$
$\displaystyle 0,\ \ \forall t\in\mathbb{R}_{>0},$ (5d) $\displaystyle
f_{1}(\underline{x}^{+},t)=$ $\displaystyle 0,\ \ \forall
t\in\mathbb{R}_{>0},$ (5e) $\displaystyle f_{1}(\overline{x}^{-},t)=$
$\displaystyle f_{1}(\overline{x}^{+},t),\ \ \forall t\in\mathbb{R}_{>0},$
(5f) $\displaystyle\mathcal{F}_{1}(\overline{x}^{+},t)=$
$\displaystyle\mathcal{F}_{0}(\overline{x}^{-},t)+\mathcal{F}_{1}(\overline{x}^{-},t),\
\ \forall t\in\mathbb{R}_{>0},$ (5g)
$\displaystyle\mathcal{F}_{1}(x_{H}^{-},t)=$ $\displaystyle 0,\ \ \forall
t\in\mathbb{R}_{>0},$ (5h)
$\displaystyle\mathcal{F}_{0}(\underline{x}^{-},t)>$
$\displaystyle\mathcal{F}_{0}(\underline{x}^{+},t),\ \ \forall
t\in\mathbb{R}_{>0},$ (5i) $\displaystyle\mathcal{F}_{1}(\overline{x}^{+},t)<$
$\displaystyle\mathcal{F}_{1}(\overline{x}^{-},t),\ \ \forall
t\in\mathbb{R}_{>0}.$ (5j)

The initial data of $f_{0}$ and $f_{1}$ defined over
$\overline{I}_{a0}:=[x_{L},\underline{x}(0)],\overline{I}_{b0}:=[\underline{x}(0),\overline{x}(0)]$,
and $\overline{I}_{c0}:=[\overline{x}(0),x_{H}]$ are given by

$\displaystyle f_{0}(0,x)=$ $\displaystyle f_{0}^{a0}(x),\ \ \forall
x\in\overline{I}_{a0},$ (6a) $\displaystyle f_{0}(0,x)=$ $\displaystyle
f_{0}^{b0}(x),\ \ \forall x\in\overline{I}_{b0},$ (6b) $\displaystyle
f_{1}(0,x)=$ $\displaystyle f_{1}^{b0}(x),\ \ \forall x\in\overline{I}_{b0},$
(6c) $\displaystyle f_{1}(0,x)=$ $\displaystyle f_{1}^{c0}(x),\ \ \forall
x\in\overline{I}_{c0}.$ (6d)

The total power demand of the TCL population at time $t\in\mathbb{R}_{\geq 0}$
is given by

$\displaystyle
y_{\text{total}}(t):=\frac{P}{\eta}\int_{\underline{x}(t)}^{x_{H}}f_{1}(x,t)\operatorname{d}\\!{x},$
(7)

where $\eta$ is the load efficiency coefficient.

######  Remark 3.1.

We provide remarks on the boundary conditions presented in (5).

  1. (i)

For continuous functions $\alpha_{0},\alpha_{1}$, and $u$, the boundary
conditions in (5) are equivalent to:

$\displaystyle\dfrac{\sigma^{2}}{2}{\partial_{x}f_{0}}(x_{L}^{+},t)=$
$\displaystyle(\alpha_{0}(x_{L}^{+},t)-u(t))f_{0}(x_{L}^{+},t),$ (8a)
$\displaystyle{\partial_{x}f_{0}}(\underline{x}^{-},t)=$
$\displaystyle{\partial_{x}f_{0}}(\underline{x}^{+},t)+{\partial_{x}f_{1}}(\underline{x}^{+},t),$
(8b) $\displaystyle f_{0}(\underline{x}^{-},t)=$ $\displaystyle
f_{0}(\underline{x}^{+},t),$ (8c) $\displaystyle f_{0}(\overline{x},t)=$
$\displaystyle 0,$ (8d) $\displaystyle f_{1}(\underline{x},t)=$ $\displaystyle
0,$ (8e) $\displaystyle f_{1}(\overline{x}^{-},t)=$ $\displaystyle
f_{1}(\overline{x}^{+},t),$ (8f)
$\displaystyle{\partial_{x}f_{1}}(\overline{x}^{+},t)=$
$\displaystyle{\partial_{x}f_{0}}(\overline{x}^{-},t)+{\partial_{x}f_{1}}(\overline{x}^{-},t),$
(8g) $\displaystyle\dfrac{\sigma^{2}}{2}{\partial_{x}f_{1}}(x_{H}^{-},t)=$
$\displaystyle(\alpha_{1}(x_{H}^{-},t)-u(t))f_{1}(x_{H}^{-},t),$ (8h)
$\displaystyle\partial_{x}f_{1}(\underline{x}^{+},t)>$ $\displaystyle 0,$ (8i)
$\displaystyle\partial_{x}f_{0}(\overline{x}^{-},t)<$ $\displaystyle 0.$ (8j)

  2. (ii)

It is worth noting that this set of boundary conditions ((5) or (8)), with
possible variations, is commonly used in the literature [19, 2, 30], which
captures the basic properties of the considered problem, for example,
impenetrable wall reflections ((8a) and (8h)), absorbing actions due to
thermostat switching ((8d)) and (8e)), and probability conservation at the
boundaries of the deadband ((8b) and (8g)). Note that because of the absorbing
property and the continuity of the PDFs on the boundaries of the deadband, the
conditions (8b) and (8g) remain the same as those originally derived in [19],
even though the considered problem in the present work contains control
actions.

###  3.3 Problem statement and basic assumptions

In this work, we study the dynamics described by the PDE model (4) under the
boundary and initial conditions (5) and (6). Based on (7), a new output
function will be defined and specified in Section 4. With these dynamics, a
continuous time controller that considers the convergence time and robustness
is designed to stabilize the tracking process.

In the sequel, we assume that $x_{a}\in{C(\mathbb{R}_{\geq 0})}$
$\underline{x},\overline{x}\in C^{1}(\mathbb{R}_{\geq 0};\mathbb{R}_{>0})$,
and denote

$\displaystyle S_{ab}:=$
$\displaystyle\left({C^{2,1}}(I_{a}\times\mathbb{R}_{>0})\cap
C(\overline{I}_{a}\times\mathbb{R}_{\geq
0})\right)\cup\left({C^{2,1}}(I_{b}\times\mathbb{R}_{>0})\cap
C(\overline{I}_{b}\times\mathbb{R}_{\geq 0})\right),$ $\displaystyle S_{bc}:=$
$\displaystyle\left({C^{2,1}}(I_{b}\times\mathbb{R}_{>0})\cap
C(\overline{I}_{b}\times\mathbb{R}_{\geq
0})\right)\cup\left({C^{2,1}}(I_{c}\times\mathbb{R}_{>0})\cap
C(\overline{I}_{c}\times\mathbb{R}_{\geq 0})\right).$

Based on the physical properties of the problem, we impose the following
structural conditions and basic assumptions on the solution and control for
the system:
$\bullet$ The function of net probability flux $g$ belongs to
$C^{1}(\mathbb{R}^{2};\mathbb{R})$ and satisfies

  1. (G1)

$g(0,\tau)\leq 0$ for all $\tau\in\mathbb{R}$;

  2. (G2)

$g(s,0)\geq 0$ for all $s\in\mathbb{R}$;

  3. (G3)

$|g_{s}(s,\tau)|+|g_{\tau}(s,\tau)|\leq 1$ for all
$(s,\tau)\in\mathbb{R}^{2}$.

$\bullet$ The pair of solution $(f_{0},f_{1})$ and the control $u$ satisfy

  1. (U)

$u\in C(\mathbb{R}_{\geq 0};\mathbb{R})$ such that
$\dot{\underline{x}}=\dot{\overline{x}}=u$ in $\mathbb{R}_{\geq 0}$;

  2. (F1)

$f_{0}^{a0}\in C(\overline{I}_{a0};\mathbb{R}_{\geq 0})$, $f_{0}^{b0}\in
C(\overline{I}_{b0};\mathbb{R}_{\geq 0})$, $f_{1}^{b0}\in
C(\overline{I}_{b0};\mathbb{R}_{\geq 0})$, $f_{1}^{c0}\in
C(\overline{I}_{c0};\mathbb{R}_{\geq 0})$;

  3. (F2)

$f_{0}\in S_{ab}$ and has derivatives $\partial_{x}f_{0}(x_{L}^{+},t)$,
$\partial_{x}f_{0}(\underline{x}^{\pm},t)$ and
$\partial_{x}f_{0}(\overline{x}^{-},t)$ for any fixed $t\in\mathbb{R}_{>0}$;

  4. (F3)

$f_{1}\in S_{bc}$ and has derivatives $\partial_{x}f_{1}(x_{H}^{-},t)$,
$\partial_{x}f_{1}(\overline{x}^{\pm},t)$ and
$\partial_{x}f_{1}(\underline{x}^{+},t)$ for any fixed $t\in\mathbb{R}_{>0}$.

######  Remark 3.2.

It should be mentioned that for $f_{0}=0$ (or $f_{1}=0$), condition (G1) (or
(G2)) guarantees $-g(f_{0},f_{1})\geq 0$ (or $g(f_{0},f_{1})\geq 0$) in (4b)
(or (4c)) . This indicates that forced switching, which generates additional
fluxes, is only possible from the $f_{1}$ system into the $f_{0}$ system when
$f_{0}$ is zero.

Condition (G3) indicates that the change in the probability density of the
additional flux cannot be too fast for practical applications. This is in
accordance with the suggestion in [30].

Condition (F1) indicates that the initial data are assumed to be nonnegative
and continuous over the given domains. Conditions (F2) and (F3) describe the
regularity of the solutions at the endpoints of the given domains at any time
$t$.

##  4 Control design and stability analysis

In this section, we design a feedback control to ensure that the output of the
system (4)-(6) tracks a reference power curve, and assess the stability of the
error dynamics in the framework of FTISS theory. Moreover, we study the mass
conservation and non-negativeness properties of the solutions to the
considered system, which allows further clarification of the physical meanings
of the mathematical model.

###  4.1 Control design

The control objective is to drive the power consumption of the population to
track the desired regulation signal. To this end, we choose an output of the
power tracking control scheme as

$\displaystyle y(t):=$ $\displaystyle
y_{\text{total}}(t)+\frac{P}{\eta}\int_{\overline{x}(t)}^{x_{H}}f_{1}(x,t)\operatorname{d}\\!{x}-\frac{P}{\eta}\int_{x_{L}}^{\underline{x}(t)}f_{0}(x,t)\operatorname{d}\\!{x},\
\ t\in\mathbb{R}_{\geq 0}.$ (9)

It is worth noting that, as the probability flows of $f_{0}$ and $f_{1}$
always move towards the deadband, $y(t)$ defined in (9) converges to the
aggregated power demand $y_{\text{total}}(t)$ in the steady state. The
motivation to add two extra terms to $y_{\text{total}}(t)$ is to ensure the
controllability of the input-output dynamics.

The regulation of power consumption of the TCL population is achieved by
moving the mass of the temperature distribution, and the control signal is
chosen to be the set-point temperature variation rate $\dot{x}_{sp}$, which
may induce a change in the probability flux [2, 37]. As we consider a control
scheme with a fixed deadband width, denoted by $\delta_{0}$, we have
$\overline{x}=x_{sp}-\frac{\delta_{0}}{2},\;\underline{x}=x_{sp}+\frac{\delta_{0}}{2}$.
Thus, the actual control signal is given by
$u(t):=\dot{x}_{sp}=\dot{\underline{x}}=\dot{\overline{x}}$.

Let $y_{d}:{\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}}$ be the desired power
profile, which is sufficiently smooth, and define the power tracking error as

$e(t):=y(t)-y_{d}(t).$

In what follows, we introduce a nonlinear control law and derive the
corresponding tracking error dynamics.

######  Theorem 4.1.

Consider the system given in (4) and (9) under the boundary conditions in (5)
(or equivalently (8)). Let the control input be defined as

$\displaystyle{u(t):=}\dfrac{k|e(t)|^{\gamma}\mathrm{sgn}(e(t))+\Phi(t)}{2\left(f_{1}(\overline{x},t)+f_{0}(\underline{x},t)\right)},$
(10)

where $k\in\mathbb{R}_{>0}$ and $\gamma\in(0,1)$ are constants,
$\mathrm{sgn}(e)$ is the sign function defined by

$\displaystyle\mathrm{sgn}(e)$ $\displaystyle:=\begin{cases}-1,&e<0,\\\
0,&e=0,\\\ 1,&e>0,\end{cases}$

and

$\displaystyle\Phi(t):=$ $\displaystyle-\frac{\eta}{P}\dot{y}_{d}(t).$ (11)

Then, the power tracking error dynamics are given by

$\displaystyle\dot{e}(t)=-\frac{P}{\eta}k|e(t)|^{\gamma}\mathrm{sgn}(e(t))+\Gamma(t),$
(12)

where

$\displaystyle\Gamma(t):=$
$\displaystyle\dfrac{P}{\eta}\left(\alpha_{1}(\overline{x},t)f_{1}(\overline{x},t)+\alpha_{0}(\underline{x},t)f_{0}(\underline{x},t)\right)-\dfrac{\sigma^{2}P}{2\eta}\left({\partial_{x}f_{1}(\underline{x}^{+},t)}+{\partial_{x}f_{1}(\overline{x}^{+},t)}\right)$
$\displaystyle-\dfrac{\sigma^{2}P}{2\eta}\left({\partial_{x}f_{0}(\underline{x}^{-},t)}+{\partial_{x}f_{0}(\overline{x}^{-},t)}\right)+\frac{P}{\eta}\int_{\underline{x}(t)}^{\overline{x}(t)}g(f_{0},f_{1})\operatorname{d}\\!{x}.$
(13)

######  Remark 4.1.

$\Gamma(t)$ defined in (4.1) captures the terms depending on the diffusion
coefficient or requiring instantaneous state measurements and will be treated
as a disturbance thereafter. Moreover, the control law given in (10) involves
only the measurement of the states (probability distributions $f_{0}$ and
$f_{1}$) on the end-points of the deadband ($\underline{x}$ and
$\overline{x}$), which results in a control scheme with significantly reduced
communication burden compared to control schemes that require full-state
measurements.

Proof of Theorem 4.1. Note that

$\displaystyle\dot{e}(t)=$ $\displaystyle\dot{y}(t)-\dot{y}_{d}(t)$
$\displaystyle=$
$\displaystyle\frac{\text{d}}{\text{d}t}\left(\frac{P}{\eta}\int_{\underline{x}(t)}^{x_{H}}f_{1}(x,t)\operatorname{d}\\!{x}+\frac{P}{\eta}\int_{\overline{x}(t)}^{x_{H}}f_{1}(x,t)\operatorname{d}\\!{x}-\frac{P}{\eta}\int_{x_{L}}^{\underline{x}(t)}f_{0}(x,t)\operatorname{d}\\!{x}\right)-\dot{y}_{d}(t)$
$\displaystyle=$
$\displaystyle\frac{P}{\eta}\frac{\text{d}}{\text{d}t}\int_{\underline{x}(t)}^{x_{H}}f_{1}(x,t)\operatorname{d}\\!{x}+\frac{P}{\eta}\frac{\text{d}}{\text{d}t}\int_{\overline{x}(t)}^{x_{H}}f_{1}(x,t)\operatorname{d}\\!{x}-\frac{P}{\eta}\frac{\text{d}}{\text{d}t}\int_{x_{L}}^{\underline{x}(t)}f_{0}(x,t)\operatorname{d}\\!{x}-\dot{y}_{d}(t).$

Hence, we decompose the whole computation process into three steps.

_Step 1:_ Compute
$\frac{\text{d}}{\text{d}t}\int_{\overline{x}(t)}^{x_{H}}f_{1}(x,t)\operatorname{d}\\!{x}$.
It follows immediately from Leibniz’s integral rule and (4d) that

$\displaystyle\frac{\text{d}}{\text{d}t}\int_{\overline{x}(t)}^{x_{H}}f_{1}(x,t)\operatorname{d}\\!{x}$
$\displaystyle=$ $\displaystyle
0-\dot{\overline{x}}(t)f_{1}(\overline{x},t)+\int_{\overline{x}(t)}^{x_{H}}\partial_{t}f_{1}(x,t)\operatorname{d}\\!{x}$
$\displaystyle=$
$\displaystyle-u(t)f_{1}(\overline{x},t)+\int_{\overline{x}(t)}^{x_{H}}{\partial_{x}}\\!\left(\dfrac{\sigma^{2}}{2}{\partial_{x}f_{1}(x,t)}-(\alpha_{1}(x,t)-u(t))f_{1}(x,t)\right)\operatorname{d}\\!{x}$
$\displaystyle=$
$\displaystyle-u(t)f_{1}(\overline{x},t)+\left(\dfrac{\sigma^{2}}{2}{\partial_{x}f_{1}(x_{H}^{-},t)}-(\alpha_{1}(x_{H}^{-},t)-u(t))f_{1}(x_{H}^{-},t)\right)-\left(\dfrac{\sigma^{2}}{2}{\partial_{x}f_{1}(\overline{x}^{+},t)}-(\alpha_{1}(\overline{x},t)-u(t))f_{1}(\overline{x},t)\right).$

Using boundary condition (8h), it follows that

$\displaystyle\frac{\text{d}}{\text{d}t}\int_{\overline{x}(t)}^{x_{H}}f_{1}(x,t)\operatorname{d}\\!{x}=$
$\displaystyle-2u(t)f_{1}(\overline{x},t)-\dfrac{\sigma^{2}}{2}{\partial_{x}f_{1}(\overline{x}^{+},t)}+\alpha_{1}(\overline{x},t)f_{1}(\overline{x},t).$
(14)

_Step 2:_ Compute
$\frac{\text{d}}{\text{d}t}\int_{\underline{x}(t)}^{x_{H}}f_{1}(x,t)\operatorname{d}\\!{x}$.
Since

$\displaystyle\frac{\text{d}}{\text{d}t}\int_{\underline{x}(t)}^{x_{H}}f_{1}(x,t)\operatorname{d}\\!{x}=$
$\displaystyle\frac{\text{d}}{\text{d}t}\int_{\underline{x}(t)}^{\overline{x}(t)}f_{1}(x,t)\operatorname{d}\\!{x}+\frac{\text{d}}{\text{d}t}\int_{\overline{x}(t)}^{x_{H}}f_{1}(x,t)\operatorname{d}\\!{x},$

and
$\frac{\text{d}}{\text{d}t}\int_{\overline{x}}^{x_{H}}f_{1}(x,t)\operatorname{d}\\!{x}$
is given by (14), we only need to compute
$\frac{\text{d}}{\text{d}t}\int_{\underline{x}(t)}^{\overline{x}(t)}f_{1}(x,t)\operatorname{d}\\!{x}$.
It follows from (4c) and (8e) that

$\displaystyle\frac{\text{d}}{\text{d}t}\int_{\underline{x}(t)}^{\overline{x}(t)}f_{1}(x,t)\operatorname{d}\\!{x}=$
$\displaystyle
u(t)f_{1}(\overline{x},t)+\int_{\underline{x}(t)}^{\overline{x}(t)}g(f_{0},f_{1})\operatorname{d}\\!{x}+\int_{\underline{x}(t)}^{\overline{x}(t)}\\!{\partial_{x}}\\!\left(\dfrac{\sigma^{2}}{2}{\partial_{x}f_{1}(x,t)}-(\alpha_{1}(x,t)-u(t))f_{1}(x,t)\right)\operatorname{d}\\!{x}$
$\displaystyle=$ $\displaystyle
u(t)f_{1}(\overline{x},t)+\int_{\underline{x}(t)}^{\overline{x}(t)}g(f_{0},f_{1})\operatorname{d}\\!{x}+\left(\dfrac{\sigma^{2}}{2}{\partial_{x}f_{1}(\overline{x}^{-},t)}-(\alpha_{1}(\overline{x}^{-},t)-u(t))f_{1}(\overline{x}^{-},t)\right)$
$\displaystyle-\left(\dfrac{\sigma^{2}}{2}{\partial_{x}f_{1}(\underline{x}^{+},t)}-(\alpha_{1}(\underline{x}^{+},t)-u(t))f_{1}(\underline{x}^{+},t)\right)$
$\displaystyle=$ $\displaystyle
u(t)f_{1}(\overline{x},t)+\dfrac{\sigma^{2}}{2}{\partial_{x}f_{1}(\overline{x}^{-},t)}-(\alpha_{1}(\overline{x},t)-u(t))f_{1}(\overline{x},t)-\dfrac{\sigma^{2}}{2}{\partial_{x}f_{1}(\underline{x}^{+},t)}$
$\displaystyle+\int_{\underline{x}(t)}^{\overline{x}(t)}g(f_{0},f_{1})\operatorname{d}\\!{x}$
$\displaystyle=$ $\displaystyle
2u(t)f_{1}(\overline{x},t)+\dfrac{\sigma^{2}}{2}{\partial_{x}f_{1}(\overline{x}^{-},t)}-\alpha_{1}(\overline{x},t)f_{1}(\overline{x},t)-\dfrac{\sigma^{2}}{2}{\partial_{x}f_{1}(\underline{x}^{+},t)}+\int_{\underline{x}(t)}^{\overline{x}(t)}g(f_{0},f_{1})\operatorname{d}\\!{x}.$
(15)

Combining (14) and (15) we obtain by (8h)

$\displaystyle\frac{\text{d}}{\text{d}t}\int_{\underline{x}(t)}^{x_{H}}f_{1}(x,t)\operatorname{d}\\!{x}=$
$\displaystyle-\dfrac{\sigma^{2}}{2}{\partial_{x}f_{1}(\underline{x}^{+},t)}-\dfrac{\sigma^{2}}{2}{\partial_{x}f_{0}(\overline{x}^{-},t)}+\int_{\underline{x}(t)}^{\overline{x}(t)}g(f_{0},f_{1})\operatorname{d}\\!{x}.$
(16)

_Step 3:_ Compute
$\frac{\text{d}}{\text{d}t}\int_{x_{L}}^{\underline{x}(t)}f_{0}(x,t)\operatorname{d}\\!{x}$.
According to (4a) and (8a), we have

$\displaystyle\frac{\text{d}}{\text{d}t}\int_{x_{L}}^{\underline{x}(t)}f_{0}(x,t)\operatorname{d}\\!{x}=$
$\displaystyle
u(t)f_{0}(\underline{x},t)+\int_{x_{L}}^{\underline{x}(t)}\partial_{t}f_{0}(x,t)\operatorname{d}\\!{x}$
$\displaystyle=$ $\displaystyle
u(t)f_{0}(\underline{x},t)+\int_{x_{L}}^{\underline{x}(t)}{\partial_{x}}\\!\left(\dfrac{\sigma^{2}}{2}{\partial_{x}f_{0}}-(\alpha_{0}-u)f_{0}\right)\operatorname{d}\\!{x}$
$\displaystyle=$ $\displaystyle
u(t)f_{0}(\underline{x},t)+\left(\dfrac{\sigma^{2}}{2}{\partial_{x}f_{0}(\underline{x}^{-},t)}-(\alpha_{0}(\underline{x},t)-u(t))f_{0}(\underline{x},t)\right)$
$\displaystyle-\left(\dfrac{\sigma^{2}}{2}{\partial_{x}f_{0}(x_{L}^{+},t)}-(\alpha_{0}(x_{L}^{+},t)-u(t))f_{0}(x_{L}^{+},t)\right)$
$\displaystyle=$ $\displaystyle
2u(t)f_{0}(\underline{x},t)+\dfrac{\sigma^{2}}{2}{\partial_{x}f_{0}(\underline{x}^{-},t)}-\alpha_{0}(\underline{x},t)f_{0}(\underline{x},t).$
(17)

Finally, by combining (14), (16), and (17), we obtain:

$\displaystyle\dot{e}(t)=$
$\displaystyle\frac{P}{\eta}\left(-\dfrac{\sigma^{2}}{2}{\partial_{x}f_{1}(\underline{x}^{+},t)}-\dfrac{\sigma^{2}}{2}{\partial_{x}f_{0}(\overline{x}^{-},t)}\right)+\frac{P}{\eta}\left(-2u(t)f_{1}(\overline{x},t)-\dfrac{\sigma^{2}}{2}{\partial_{x}f_{1}(\overline{x}^{+},t)}+\alpha_{1}(\overline{x},t)f_{1}(\overline{x},t)\right)$
$\displaystyle-\frac{P}{\eta}\left(2u(t)f_{0}(\underline{x},t)+\dfrac{\sigma^{2}}{2}{\partial_{x}f_{0}(\underline{x}^{-},t)}-\alpha_{0}(\underline{x},t)f_{0}(\underline{x},t)\right)-\dot{y_{d}}(t)+\frac{P}{\eta}\int_{\underline{x}(t)}^{\overline{x}(t)}g(f_{0},f_{1})\operatorname{d}\\!{x}$
$\displaystyle=$
$\displaystyle-\frac{2P}{\eta}u(t)\left(f_{1}(\overline{x},t)+f_{0}(\underline{x},t)\right)-\dfrac{\sigma^{2}P}{2\eta}{\partial_{x}f_{1}(\underline{x}^{+},t)}-\dfrac{\sigma^{2}P}{2\eta}{\partial_{x}f_{1}(\overline{x}^{+},t)}-\dfrac{\sigma^{2}P}{2\eta}{\partial_{x}f_{0}(\underline{x}^{-},t)}$
$\displaystyle-\dfrac{\sigma^{2}P}{2\eta}{\partial_{x}f_{0}(\overline{x}^{-},t)}+\frac{P}{\eta}\int_{\underline{x}(t)}^{\overline{x}(t)}g(f_{0},f_{1})\operatorname{d}\\!{x}+\frac{P}{\eta}\alpha_{1}(\overline{x},t)f_{1}(\overline{x},t)+\frac{P}{\eta}\alpha_{0}(\underline{x},t)f_{0}(\underline{x},t)-\dot{y}_{d}(t).$

The error dynamics can then be expressed as

$\dot{e}(t)=-\frac{2P}{\eta}u(t)\left(f_{1}(\overline{x},t)+f_{0}(\underline{x},t)\right)+\frac{P}{\eta}\Phi(t)+\Gamma(t).$

Let

$\displaystyle
u(t):=\dfrac{v(t)+\Phi(t)}{2\left(f_{1}(\overline{x},t)+f_{0}(\underline{x},t)\right)},$

where $v(t)$ is an auxiliary control input, then

$\displaystyle\dot{e}(t)=-\frac{P}{\eta}v(t)+\Gamma(t).$ (18)

Considering an auxiliary control of the form:

$\displaystyle v(t):=k|e(t)|^{\gamma}\mathrm{sgn}(e(t)),$ (19)

the tracking error dynamics in the closed loop are then given by (12).
$\hfill\blacksquare$

######  Remark 4.2.

Note that for the given initial data (see (F1)), it can be shown that the term
$f_{1}(\overline{x},t)+f_{0}(\underline{x},t)$ is strictly positive (see
Theorem 4.4 (iii) in Section 4.2). Therefore, the control signal $u$, given in
(10) is well-defined. In addition, $u$ is continuous due to the fact that
$\gamma\in(0,1)$ and the assumptions on the continuity of $\dot{y}_{d}(t)$ and
$f_{1}(\overline{x},t)+f_{0}(\underline{x},t)$ (see (F2) and (F3)). It is also
worth noting that, as $f_{1}(\overline{x},t)$ and $f_{0}(\underline{x},t)$
describe the probability density of TCLs in the ON and OFF states at the
prescribed upper and lower temperature boundaries $\overline{x}$ and
$\underline{x}$, respectively, it is impossible in practice that
$f_{1}(\overline{x},t)+f_{0}(\underline{x},t)\rightarrow 0$ as $t\to+\infty$.

###  4.2 Finite-time input-to-state stability of the tracking error dynamics

In this section, we assess the robust stability of the tracking error dynamics
in the sense of FTISS, with $\Gamma$ as the input (disturbance). One of the
main properties of the closed-loop system is stated below.

######  Theorem 4.2.

The power tracking error dynamics (12) under the control law given in (10) are
FTISS w.r.t. $\Gamma(t)$ for any $\gamma\in(0,1)$.

Proof. Consider a Lyapunov candidate of the form $V(e)=\frac{1}{2}e^{2}$. The
time derivative of $V$ along the trajectory of the tracking error dynamics
(12) is given by:

$\displaystyle\dot{V}$
$\displaystyle=e\dot{e}=e\left(-\frac{P}{\eta}k|e|^{\gamma}\mathrm{sgn}(e)+\Gamma\right)=-\frac{P}{\eta}k|e|^{1+\gamma}+e\Gamma=-\frac{P}{\eta}k\left(\sqrt{2V}\right)^{1+\gamma}+e\Gamma=-\frac{P}{\eta}k\left(2V\right)^{\frac{1+\gamma}{2}}+e\Gamma,$

which implies that

$DV(e)\cdot
f(e,\Gamma)\leq-\frac{P}{\eta}k(2V)^{\frac{1+\gamma}{2}}+|e||\Gamma|$ (20)

with
$f(e,\Gamma):=-\frac{P}{\eta}k|e(t)|^{\gamma}\mathrm{sgn}(e(t))+\Gamma(t)$.

Let $C_{0}\in(0,k)$ be a constant. Then, for any
$|e|\geq\left(\frac{\eta}{PC_{0}}|\Gamma|\right)^{\frac{1}{\gamma}}$, i.e.,
$|\Gamma|\leq\frac{P}{\eta}C_{0}|e|^{\gamma}$, we deduce by (20) that

$\displaystyle DV(e)\cdot
f(e,\Gamma)\leq-\frac{P}{\eta}k(2V)^{\frac{1+\gamma}{2}}+\frac{P}{\eta}C_{0}|e|^{1+\gamma}=-\frac{P}{\eta}k(2V)^{\frac{1+\gamma}{2}}+\frac{P}{\eta}C_{0}(2V)^{\frac{1+\gamma}{2}}=-\frac{P}{\eta}(k-C_{0})2^{\frac{1+\gamma}{2}}V^{\frac{1+\gamma}{2}}.$

Note that $\frac{P}{\eta}(k-C_{0})2^{\frac{1+\gamma}{2}}>0$,
$\frac{1+\gamma}{2}\in(\frac{1}{2},1)$, and that
$\chi(s):=(\frac{\eta}{PC_{0}}s)^{\frac{1}{\gamma}}$ is a
$\mathcal{K}$-function w.r.t. $s\in\mathbb{R}_{\geq 0}$. The FTISS of system
(12) is then guaranteed by Lemma 2.1. $\hfill\blacksquare$

###  4.3 Properties of the governing PDEs

In practice, we can assume that the number of TCLs in a population remains
unchanged within a specific DR control period. Therefore, the mass
conservation property of the solutions to the system (4)-(6) should be
verified under the imposed boundary conditions, thereby conforming the
compliance of the mathematical model with the imposed condition. Moreover,
non-negativeness of the solutions is also required.

######  Theorem 4.3 (Mass conservation property).

The solution to the initial-boundary value problem (IBVP) (4)-(6) is
conservative in the sense that

$\displaystyle\int_{x_{L}}^{\overline{x}(t)}f_{0}(x,t)\operatorname{d}\\!{x}+\int_{\underline{x}(t)}^{x_{H}}f_{1}(x,t)\operatorname{d}\\!{x}=1\
\ \forall t\in\mathbb{R}_{\geq 0},$ (21)

provided that

$\displaystyle\int_{x_{L}}^{\underline{x}(0)}f_{0}^{a0}(x)\operatorname{d}\\!{x}+\int_{\underline{x}(0)}^{\overline{x}(0)}f_{0}^{b0}(x)\operatorname{d}\\!{x}+\int_{\underline{x}(0)}^{\overline{x}(0)}f_{1}^{b0}(x)\operatorname{d}\\!{x}+\int_{\overline{x}(0)}^{x_{H}}f_{1}^{c0}(x)\operatorname{d}\\!{x}=1.$
(22)

Proof. Using (4a), (4b), (8a), (8b), and (8d), and noting (U) and (F2), we
have

$\displaystyle\frac{\text{d}}{\operatorname{d}\\!{t}}\left(\int_{x_{L}}^{\overline{x}(t)}f_{0}(x,t)\operatorname{d}\\!{x}\right)$
$\displaystyle=$
$\displaystyle\frac{\operatorname{d}\\!{}}{\operatorname{d}\\!{t}}\left(\int_{x_{L}}^{\underline{x}(t)}f_{0}(x,t)\operatorname{d}\\!{x}+\int_{\underline{x}(t)}^{\overline{x}(t)}f_{0}(x,t)\operatorname{d}\\!{x}\right)$
$\displaystyle=$
$\displaystyle\int_{x_{L}}^{\underline{x}(t)}\partial_{t}f_{0}(x,t)\operatorname{d}\\!{x}+f_{0}(\underline{x}(t),t)\dot{\underline{x}}(t)+\int_{\underline{x}(t)}^{\overline{x}(t)}\partial_{t}f_{0}(x,t)\operatorname{d}\\!{x}+f_{0}(\overline{x}(t),t)\dot{\overline{x}}(t)-f_{0}(\underline{x}(t),t)\dot{\underline{x}}(t)$
$\displaystyle=$
$\displaystyle\int_{x_{L}}^{\underline{x}(t)}\partial_{x}\\!\left(\frac{\sigma^{2}}{2}\partial_{x}f_{0}(x,t)-(\alpha_{0}(x,t)-u(t))f_{0}(x,t)\right)\operatorname{d}\\!{x}+\int_{\underline{x}(t)}^{\overline{x}(t)}\partial_{x}\\!\left(\frac{\sigma^{2}}{2}\partial_{x}f_{0}(x,t)-(\alpha_{0}(x,t)-u(t))f_{0}(x,t)\right)\operatorname{d}\\!{x}$
$\displaystyle-\int_{\underline{x}(t)}^{\overline{x}(t)}g(f_{0},f_{1})\operatorname{d}\\!{x}$
$\displaystyle=$
$\displaystyle\left(\frac{\sigma^{2}}{2}\partial_{x}f_{0}(x,t)-(\alpha_{0}(x,t)-u(t))f_{0}(x,t)\right)\bigg{|}_{x_{L}^{+}}^{\underline{x}^{-}(t)}\
+\left(\frac{\sigma^{2}}{2}\partial_{x}f_{0}(x,t)-(\alpha_{0}(x,t)-u(t))f_{0}(x,t)\right)\bigg{|}_{\underline{x}^{+}(t)}^{\overline{x}^{-}(t)}$
$\displaystyle-\int_{\underline{x}(t)}^{\overline{x}(t)}g(f_{0},f_{1})\operatorname{d}\\!{x}$
$\displaystyle=$
$\displaystyle\frac{\sigma^{2}}{2}\partial_{x}f_{0}(\underline{x}^{-}(t),t)-(\alpha_{0}(\underline{x}(t))-u(t))f_{0}(\underline{x}(t),t)-0+\frac{\sigma^{2}}{2}\partial_{x}f_{0}(\overline{x}^{-}(t),t)-(\alpha_{0}(\overline{x}(t))-u(t))f_{0}(\overline{x}(t),t)$
$\displaystyle-\left(\frac{\sigma^{2}}{2}\partial_{x}f_{0}(\underline{x}^{+}(t),t)-(\alpha_{0}(\underline{x}(t))-u(t))f_{0}(\underline{x}(t),t)\right)-\int_{\underline{x}(t)}^{\overline{x}(t)}g(f_{0},f_{1})\operatorname{d}\\!{x}$
$\displaystyle=$
$\displaystyle\frac{\sigma^{2}}{2}(\partial_{x}f_{0}(\underline{x}^{-}(t),t)-\partial_{x}f_{0}(\underline{x}^{+}(t),t))+\frac{\sigma^{2}}{2}\partial_{x}f_{0}(\overline{x}^{-}(t),t)-\int_{\underline{x}(t)}^{\overline{x}(t)}g(f_{0},f_{1})\operatorname{d}\\!{x}$
$\displaystyle=$
$\displaystyle\frac{\sigma^{2}}{2}\partial_{x}f_{1}(\underline{x}^{+},t)+\frac{\sigma^{2}}{2}\partial_{x}f_{0}(\overline{x}^{-}(t),t)-\int_{\underline{x}(t)}^{\overline{x}(t)}g(f_{0},f_{1})\operatorname{d}\\!{x}.$
(23)

Similarly, we infer from (4c), (4d), (8e), (8g), (8h), (U) and (F3) that

$\displaystyle\frac{\operatorname{d}\\!{}}{\operatorname{d}\\!{t}}\left(\int_{\underline{x}(t)}^{x_{H}}f_{1}(x,t)\operatorname{d}\\!{x}\right)=$
$\displaystyle\frac{\text{d}}{\text{d}t}\left(\int_{\underline{x}(t)}^{\overline{x}(t)}f_{1}(x,t)\operatorname{d}\\!{x}+\int_{\overline{x}(t)}^{x_{H}}f_{1}(x,t)\operatorname{d}\\!{x}\right)$
$\displaystyle=$
$\displaystyle-\frac{\sigma^{2}}{2}\partial_{x}f_{1}(\underline{x}^{+},t)-\frac{\sigma^{2}}{2}\partial_{x}f_{0}(\overline{x}^{-}(t),t)+\int_{\underline{x}(t)}^{\overline{x}(t)}g(f_{0},f_{1})\operatorname{d}\\!{x}.$
(24)

By (23) and (24), we obtain

$\displaystyle\frac{\text{d}}{\text{d}t}\left(\int_{x_{L}}^{\overline{x}(t)}f_{0}(x,t)\operatorname{d}\\!{x}+\int_{\underline{x}(t)}^{x_{H}}f_{1}(x,t)\operatorname{d}\\!{x}\right)=0\
\ \forall t\in\mathbb{R}_{\geq 0},$

which along with (22) implies (21). $\hfill\blacksquare$

######  Theorem 4.4 (Non-negativeness).

The following statements hold true for the solution to IBVP (4)-(6):

  1. (i)

$f_{0}(x,t)\geq 0$ for all $x\in[x_{L},\overline{x}(t)]$ and all
$t\in\mathbb{R}_{\geq 0}$;

  2. (ii)

$f_{1}(x,t)\geq 0$ for all $x\in[\underline{x}(t),x_{H}]$ and all
$t\in\mathbb{R}_{\geq 0}$;

  3. (iii)

$f_{0}(\underline{x}(t),t)+f_{1}(\overline{x}(t),t)>0$ for all
$t\in\mathbb{R}_{>0}$.

The proof of this theorem is provided in Appendix.

##  5 Experimental Validation

In this section, we present simulation results to demonstrate the
effectiveness of the proposed control scheme. Note that the control law given
in (10) is derived from the coupled Fokker-Planck equations, which assume a
population of an infinite number of TCLs. As the number of TCLs in a real-
world TCL population is always finite, and considering the fact that the
larger the population size, the more accurate the PDE model, we present a
comparative study of two heterogeneous populations with 1,000 and 100,000
TCLs.

###  5.1 Simulation setup

A numerical simulation is conducted to validate the proposed control scheme
and evaluate its performance. Table 1 lists the physical parameters of the AC
units utilized in the simulation, which are the same as those in [2]. The
thermal resistances and thermal capacitances are random variables following a
log-normal distribution with average mean values of $2$ ∘C/kW and $10$ kWh/∘C,
respectively. The level of heterogeneity is parameterized by the standard
deviation $\sigma_{p}$. In our experiment, the initial temperatures of the AC
units are uniformly distributed around the initial set-point $x_{sp}^{0}=20$
∘C over the deadband, and initially $40\%$ of the AC units are set randomly in
“ON”-state. This setting causes the population to begin running from an almost
steady state.

Table 1: Simulation parameter

Parameter | Description (Unit) | Value
---|---|---
$R$ | average thermal resistance (∘C/kW) | 2
$C$ | average thermal capacitance (kWh/∘C) | 10
$P$ | electric power (kW) | 14
$\eta$ | load efficiency | 2.5
$x_{sp}^{0}$ | initial temperature set-point (∘C) | 20
$\delta$ | temperature deadband width (∘C) | 0.5
$\sigma_{p}$ | standard deviation of lognormal distributions | 0.2
$p_{f}$ | forced switch probability per hour (%) | 3
$t_{ci}$ | control interval (second) | 30
$t_{\mathrm{lock}}$ | locked time of each TCL (minute) | 6

The disturbances brought into the system come mainly from the following three
sources. First, all AC units operate under the same varying outside
temperature, as depicted in Fig. 3, which rises from $30$∘C at 11:30 to $23$∘C
at 12:30 and then drops back from 14:30 to 15:30. Second, a forced random
switch mechanism is added to desynchronize AC operations. The number of forced
interrupts per hour can be adjusted through the hyper-parameter $p_{f}$.
Moreover, a safe border distance of $5\%$ of the deadband width is
incorporated to prevent forced switches from happening when an AC is around
$\overline{x}(t)$ and in “ON” state or around $\underline{x}(t)$ and in “OFF”
state. Finally, because frequent switching leads to reduced energy efficiency
and more rapid compressor wear out, a lockout time, $t_{\mathrm{lock}}$, is
included for each AC. Thus, an AC unit remains inactive to the control signals
when it is locked.

Figure 3: Ambient temperature.

The reference power is a predefined curve, as shown in Fig. 4. From 10:30 to
11:30, the normalized desired power is maintained constant at $0.4$. From
11:30 to 12:00, the reference power drops to $0.2$ and keeps constant for the
following two and a half hours. From 14:30, the desired power rises to $0.5$
in 30 minutes and remains constant until 16:30. During the rising and dropping
phases, the desired power is specified by a smooth polynomial with the
endpoint constraints given below:

$\displaystyle
y_{d}(t)=\left(y_{d}(t_{f})-y_{d}(t_{i})\right)\tau^{5}(t)\sum_{l=0}^{4}a_{l}\tau^{l}(t),t\in[t_{i},t_{f}],$
(25)
$\displaystyle{\dot{y}}_{d}(t_{i})={\dot{y}}_{d}(t_{f})={\ddot{y}}_{d}(t_{i})={\ddot{y}}_{d}(t_{i})={\overset{\ldots}{y}}_{d}(t_{i})={\overset{\ldots}{y}}_{d}(t_{f})=0,$
(26)

where $\tau(t)=(t-t_{i})/(t_{f}-t_{i})$. By a direct computation, the
coefficients can be determined as follows:

$\displaystyle a_{0}=126,\,a_{1}=420,\,a_{2}=540,\,a_{3}=315,\,\mathrm{and}\
a_{4}=70.$

Figure 4: Desired power profile.

In the simulation, the control signal is updated every 30 seconds ($t_{ci}$ in
Table LABEL:tab:table1). The control signal that every AC receives is the set-
point variation rate. Each AC computes then its set-point temperature offset
for the next control interval starting from $t_{k}$. To compute the
denominator of the controller given in (10), a mid-point rectangular method
with a temperature bin width $\delta_{x}$ is used to estimate
$f_{1}(\overline{x}(t_{k}),t_{k})$ and $f_{0}(\underline{x}(t_{k}),t_{k})$.
The percentage of ACs falling in the rectangular region is used as
$f_{1}(\overline{x}(t_{k}),t_{k})\times\delta_{x}$ or
$f_{0}(\underline{x}(t_{k}),t_{k})\times\delta_{x}$. In general, $\delta_{x}$
should not be too large because the underlying system has complex nonlinear
dynamics. On the other hand, considering the limited number of ACs involved in
the simulation, the bin width $\delta_{x}$ should not be too small, which may
introduce larger biases. In our implementation, histogram bin widths of
$0.008~{}^{\circ}$C, $0.004~{}^{\circ}$C, and $0.002~{}^{\circ}$C are used,
which are reasonable and provide reliable estimations of
$f_{1}(\overline{x}(t_{k}),t_{k})$ and $f_{0}(\underline{x}(t_{k}),t_{k})$.

Table 2: Tracking performance of 10 episodes for the population with 1,000
TCLs

Episode | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---
RMSE ($\%$) | 0.948 | 0.923 | 0.844 | 0.834 | 0.935
Episode | 6 | 7 | 8 | 9 | 10
RMSE ($\%$) | 0.880 | 0.923 | 0.890 | 0.925 | 0.861

###  5.2 Simulation results

First, we present the test results for the population with 1,000 TCLs. The
control cycle lasts for 6 hours, from 10:30 to 16:30. The test is performed
continuously for $10$ episodes, and the tracking performance is measured by
the root mean square error (RMSE), as reported in Table 2. In the test, the
controller parameters in (10) are set to be $k=8$ and $\gamma=0.5$,
respectively. The final result shows that the mean RMSE for this setting is
$0.896\%$, and the standard deviation (STD) of the dRMSEs is $0.040\%$.

Fig. 5 shows a sample of the control results corresponding to the episode with
an RMSE of $0.948\%$. It can be seen from Fig. 5a that the proposed control
strategy is effective. The temperature evolution of $10$ randomly selected ACs
in the population is presented in Fig. 5b. It can be observed that all of
them, unless forced switches occur, operate smoothly inside the deadband
between the turning on and turning off points. Fig. 5c shows the control
signal generated during this episode. During the first 30 minutes (from 10:00
to 10:30), the controller is inactive, and the system operates in an open-loop
mode.

Figure 5: Control performance for a population of 1,000 TCLs

(a) tracking performance

(b) temperature trajectories of $10$ ACs

(c) set-point variation rate

When the number of ACs increases, the model of the coupled Fokker-Planck
equations becomes more accurate. To evaluate the effectiveness of the proposed
control strategy, tracking control performance is examined for a population of
100,000 ACs. The RMSE values for 10 continuous tests are shown in Table 1,
which gives a mean RMSE of $0.497\%$ and an STD of $0.004\%$. In this test,
$k=15$ and $\gamma=0.5$ are used. Fig. 6 illustrates one of the control
samples corresponding to the episode with an RMSE of $0.505\%$. The normalized
power consumption is shown in Fig. 6a, and the temperature evolutions of $10$
ACs are shown in Fig. 6b. The control signal is shown in Fig. 6c.

Table 1: Tracking performance of 10 episodes for the population with 100,000 TCLs Episode | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---
RMSE ($\%$) | 0.505 | 0.500 | 0.496 | 0.491 | 0.490
Episode | 6 | 7 | 8 | 9 | 10
RMSE ($\%$) | 0.497 | 0.499 | 0.495 | 0.498 | 0.500

Figure 6: Control performance for a population of 100,000 TCLs

(a) tracking performance

(b) temperature trajectories of $10$ ACs

(c) set-point variation rate

The results of the comparative study show clearly that the tracking control
system performs better for the population of larger size with smaller RMSE,
smoother power trajectory, and less “noisy” control signals. This is
consistent with the nature of the PDE model on which the proposed control
scheme is based. Nevertheless, the performance is not significantly degraded
for a population with a significantly smaller size. This demonstrates the
robustness and potential applicability of the developed control strategy to
practical systems.

##  6 Conclusion

In this work, we have developed a strategy for the power tracking control of
heterogeneous TCL populations. The control scheme can ensure a robust
performance in the presence of modeling uncertainties in the sense of FTISS
and requires measuring the states of the system only on the end-points of the
deadband. The simulation results provided encouraging evidence that the
proposed control approach is highly effective. From a practical application
viewpoint, we can consider in our future work other types of devices, such as
battery charging systems, and other demand-response tasks, such as frequency
regulation or transaction control [13, 24, 10]. Control of populations of TCLs
described by the second-order equivalent thermal parameters model [15, 3] may
also be a research direction worthy of exploration.

## Appendix: Proof of Theorem 5

We first prove statement (i). Given any $T>0$, it suffices to show that
$f_{0}\geq 0$ over $[x_{L},\overline{x}(t)]\times[0,T]$ for all $t\in[0,T]$.

Indeed, the transformations of variable
$y:=\frac{x-x_{L}}{\overline{x}-x_{L}}:=\frac{x-x_{L}}{h}$ and
$f_{0}(x,t)=f_{0}(yh+x_{L},t):=\tilde{f}_{0}(y,t)$ yield

$\displaystyle\partial_{x}f_{0}$
$\displaystyle=\frac{1}{h}\partial_{y}\tilde{f}_{0},$
$\displaystyle\partial_{xx}f_{0}$
$\displaystyle=\frac{1}{h^{2}}\partial_{yy}\tilde{f}_{0},\partial_{t}f_{0}=\partial_{t}\tilde{f}_{0}+\partial_{y}\tilde{f}_{0}\frac{\partial
y}{\partial
t}=\partial_{t}\tilde{f}_{0}-(x-x_{L})\frac{\dot{\overline{x}}}{h^{2}}\partial_{y}\tilde{f}_{0}=\partial_{t}\tilde{f}_{0}-\frac{1}{h}yu\partial_{y}\tilde{f}_{0}.$

Note that

$\displaystyle x\in[x_{L},\underline{x}]\Leftrightarrow y\in[0,1],$
$\displaystyle 0<\delta_{0}\leq h(t)\leq x_{H}-x_{L},\forall t\in[0,T].$

The PDEs (4a) and (4b) are equivalent to

$\displaystyle\partial_{t}\tilde{f}_{0}-\frac{1}{h}\left(\frac{\sigma^{2}}{2h}\partial_{yy}\tilde{f}_{0}+\left((1+y)u-\tilde{\alpha}_{0}\right)\partial_{y}\tilde{f}_{0}-\tilde{\alpha}_{0y}\tilde{f}_{0}\right)=$
$\displaystyle 0,\ \ \forall y\in\left(0,z(t)\right),\forall t\in(0,T],$ (27a)
$\displaystyle\partial_{t}\tilde{f}_{0}-g(\tilde{f}_{0},\tilde{f}_{1})-\frac{1}{h}\left(\frac{\sigma^{2}}{2h}\partial_{yy}\tilde{f}_{0}+\left((1+y)u-\tilde{\alpha}_{0}\right)\partial_{y}\tilde{f}_{0}-\tilde{\alpha}_{0y}\tilde{f}_{0}\right)=$
$\displaystyle 0,\ \ \forall y\in\left(z(t),1\right),\forall t\in(0,T],$ (27b)

respectively, where $\tilde{\alpha}_{0}(y,t):=\alpha_{0}(yh(t)+x_{L},t)$,
$f_{1}(x,t)=f_{1}(yh+x_{L},t):=\tilde{f}_{1}(y,t)$, and
$z(t):=1-\frac{\delta_{0}}{h(t)}$.

Note that (8) is equivalent to (5), and (8a), (8b), and (8d) become

$\displaystyle\frac{\sigma^{2}}{2}\partial_{y}\tilde{f}_{0}(0^{+},t)-(\tilde{\alpha}_{0}(0^{+},t)-u(t))h(t)\tilde{f}_{0}(0^{+},t)=$
$\displaystyle 0,\ \forall t\in(0,T],$ (28a)
$\displaystyle\partial_{y}\tilde{f}_{0}\left(z^{-}(t),t\right)-\partial_{y}\tilde{f}_{0}\left(z^{+}(t),t\right)=$
$\displaystyle\sigma_{0}(t),\ \forall t\in(0,T],$ (28b)
$\displaystyle\tilde{f}_{0}(1^{-},t)=$ $\displaystyle 0,\ \forall t\in(0,T],$
(28c)

where, for the given solution $f_{1}$,
$\sigma_{0}(t):=\frac{\sigma^{2}}{2}\partial_{x}f_{1}(\underline{x}^{+}(t),t)$
is a well-defined function w.r.t. $t$, and $\sigma_{0}(t)>0$ for all
$t\in[0,T]$ owing to (F3) and (8i).

The initial data of $\tilde{f}_{0}$ over the domain $\left[0,z(t)\right]$ and
$\left[z(t),1\right]$ are given by

$\displaystyle\tilde{f}_{0}^{a0}(y):=f_{0}^{a0}(yh(0)+x_{L})\geq 0,$

and

$\displaystyle\tilde{f}_{0}^{b0}(y):=f_{0}^{b0}(yh(0)+x_{L})\geq 0,$

respectively.

Let $\phi(y):=\e^{m(y-\frac{1}{2})^{2}}$ and $\tilde{f}_{0}:=\phi\e^{\gamma
t}\hat{f}_{0}$ with $m>0$ and $\gamma>0$ being constants that will be chosen
later. Then (27) and (28) lead to

$\displaystyle\partial_{t}\hat{f}_{0}-\frac{\sigma^{2}}{2h^{2}}\partial_{yy}\hat{f}_{0}+\mathcal{B}(y,t)\partial_{y}\hat{f}_{0}+\mathcal{C}(y,t)\hat{f}_{0}=$
$\displaystyle 0,\forall y\in\left(0,z(t)\right),\forall t\in(0,T],$ (29a)
$\displaystyle\partial_{t}\hat{f}_{0}-\frac{\sigma^{2}}{2h^{2}}\partial_{yy}\hat{f}_{0}+\mathcal{B}(y,t)\partial_{y}\hat{f}_{0}+\mathcal{C}(y,t)\hat{f}_{0}+\frac{e^{-\gamma
t}}{\phi(y)}g(\tilde{f}_{0},\tilde{f}_{1})=$ $\displaystyle 0,\forall
y\in\left(z(t),1\right),\forall t\in(0,T],$ (29b)
$\displaystyle\frac{\sigma^{2}}{2}\partial_{y}\hat{f}_{0}(0^{+},t)-k(t)\hat{f}_{0}(0^{+},t)=$
$\displaystyle 0,\forall t\in(0,T],$ (29c)
$\displaystyle\partial_{y}\hat{f}_{0}\left(z^{-}(t),t\right)-\partial_{y}\hat{f}_{0}\left(z^{+}(t),t\right)=$
$\displaystyle\hat{\sigma}_{0}(t),\forall t\in(0,T],$ (29d)
$\displaystyle\hat{f}_{0}(1^{-},t)=$ $\displaystyle 0,\forall t\in(0,T],$
(29e)

where

$\displaystyle\mathcal{B}(y,t):=$
$\displaystyle-\frac{1}{h}\left(\frac{\sigma^{2}}{2h}\frac{2\partial_{y}\phi}{\phi}+(1+y)u-\tilde{\alpha}_{0}\right),$
$\displaystyle\mathcal{C}(y,t):=$
$\displaystyle\frac{1}{h}\left(\gamma-\frac{\sigma^{2}}{2h}\frac{\partial_{yy}\phi}{\phi}-\frac{\partial_{y}\phi}{\phi}\left((1+y)u-\tilde{\alpha}_{0}\right)+\tilde{\alpha}_{0y}\right),$
$\displaystyle k(t):=$
$\displaystyle\frac{m\sigma^{2}}{2}+(\tilde{\alpha}_{0}(0^{+},t)-u(t))h(t),$
$\displaystyle\hat{\sigma}_{0}(t):=$ $\displaystyle\frac{\e^{-\gamma
t}}{\phi(1)}\sigma_{0}(t).$

The initial data for the $\hat{f}_{0}$-system over the domain
$\left[0,z(t)\right]$ and $\left[z(t),1\right]$ are given by

$\displaystyle\hat{f}_{0}^{a0}(y):=\frac{\tilde{f}_{0}^{a0}(y)}{\phi(y)}\geq
0\ \ \text{and}\ \
\hat{f}_{0}^{b0}(y):=\frac{\tilde{f}_{0}^{b0}(y)}{\phi(y)}\geq 0,$ (30)

respectively.

Note that $u,\tilde{\alpha}_{0}$, and $\tilde{\alpha}_{0y}$ are continuous in
$[0,1]\times[0,T]$. Letting first $m$ and then $\gamma$ be sufficiently large,
there must be positive constants $k_{0}$ and $c_{0}$ such that

$\displaystyle k(t)\geq$ $\displaystyle k_{0},\forall t\in(0,T],$ (31)
$\displaystyle\mathcal{C}(y,t)-1\geq$ $\displaystyle
c_{0},\forall(y,t)\in(0,1)\times(0,T].$ (32)

To prove the non-negativeness property of $f_{0}$, it suffices to show that
$\hat{f}_{0}\geq 0$ in $[0,1]\times[0,T]$. We now proceed with the proof by
contradiction. Assume that there exists a point
$(y_{0},t_{0})\in[0,1]\times[0,T]$ such that

$\displaystyle\hat{f}_{0}(y_{0},t_{0})=\min_{(y,t)\in[0,1]\times[0,T]}\hat{f}_{0}(y,t)<0.$

Considering (29e) and (30), we have $y_{0}\neq 1$ and $t_{0}\in(0,T]$.

_Case 1_ : $y_{0}\in\left(0,z(t_{0})\right)$. At point $(y_{0},t_{0})$, it
holds that

$\displaystyle\partial_{t}\hat{f}_{0}(y_{0},t_{0})\leq
0,\partial_{y}\hat{f}_{0}(y_{0},t_{0})=0,\partial_{yy}\hat{f}_{0}(y_{0},t_{0})\geq
0.$

Then (29a) and (32) imply that

$\displaystyle\begin{split}0>&\left(c_{0}+1\right)\hat{f}_{0}(y_{0},t_{0})\geq\partial_{t}\hat{f}_{0}(y_{0},t_{0})-\frac{\sigma^{2}}{2h^{2}(t_{0})}\partial_{yy}\hat{f}_{0}(y_{0},t_{0})+\mathcal{B}(y_{0},t_{0})\partial_{y}\hat{f}_{0}(y_{0},t_{0})+\mathcal{C}(y_{0},t_{0})\hat{f}_{0}(y_{0},t_{0})=0,\end{split}$

which leads to a contradiction.

_Case 2_ : $y_{0}\in\left(z(t_{0}),1\right)$. At the point $(y_{0},t_{0})$, it
also holds that

$\displaystyle\partial_{t}\hat{f}_{0}(y_{0},t_{0})\leq
0,\partial_{y}\hat{f}_{0}(y_{0},t_{0})=0,\partial_{yy}\hat{f}_{0}(y_{0},t_{0})\geq
0.$

In addition, using the Mean Value Theorem, (G1), and (G2), we obtain:

$\displaystyle\begin{split}g(\tilde{f}_{0}(y_{0},t_{0}),\tilde{f}_{1}(y_{0},t_{0}))=g(0,\tilde{f}_{1}(y_{0},t_{0}))+\tilde{f}_{0}(y_{0},t_{0})g_{s}(s,\tilde{f}_{1}(y_{0},t_{0}))|_{s=\xi}\leq|\tilde{f}_{0}(y_{0},t_{0})|,\end{split}$

where $\xi$ is between $0$ and $\tilde{f}_{0}(y_{0},t_{0})$.

It follows that

$\displaystyle\frac{\e^{-\gamma
t_{0}}}{\phi(y_{0})}g(\tilde{f}_{0}(y_{0},t_{0}),\tilde{f}_{1}(y_{0},t_{0}))\leq$
$\displaystyle|\tilde{f}_{0}(y_{0},t_{0})|\frac{\e^{-\gamma
t_{0}}}{\phi(y_{0})}=-\hat{f}_{0}(y_{0},t_{0}).$ (33)

From (29b), (32), and (33), we obtain:

$\displaystyle\begin{split}0>&c_{0}\hat{f}_{0}(y_{0},t_{0})\\\
\geq&\left(\mathcal{C}(y_{0},t_{0})-1\right)\hat{f}_{0}(y_{0},t_{0})\\\
\geq&\mathcal{C}(y_{0},t_{0})\hat{f}_{0}(y_{0},t_{0})+\frac{\e^{-\gamma
t_{0}}}{\phi(y_{0})}g(\tilde{f}_{0}(y_{0},t_{0}),\tilde{f}_{1}(y_{0},t_{0}))\\\
\geq&\partial_{t}\hat{f}_{0}(y_{0},t_{0})-\frac{\sigma^{2}}{2h^{2}(t_{0})}\partial_{yy}\hat{f}_{0}(y_{0},t_{0})+\mathcal{B}(y_{0},t_{0})\partial_{y}\hat{f}_{0}(y_{0},t_{0})+\mathcal{C}(y_{0},t_{0})\hat{f}_{0}(y_{0},t_{0})+\frac{\e^{-\gamma
t_{0}}}{\phi(y_{0})}g(\tilde{f}_{0}(y_{0},t_{0}),\tilde{f}_{1}(y_{0},t_{0}))\\\
=&0,\end{split}$

which leads to a contradiction.

_Case 3_ : $y_{0}=0$. It follows that
$\partial_{y}\hat{f}_{0}(0^{+},t_{0})\geq 0$, which, along with (29c) and
(31), yields

$\displaystyle\begin{split}0&<-k_{0}\hat{f}_{0}(0^{+},t_{0})\leq-k(t_{0})\hat{f}_{0}(0^{+},t_{0})\leq\frac{\sigma^{2}}{2}\partial_{t}\hat{f}_{0}(0^{+},t)-k(t_{0})\hat{f}_{0}(0^{+},t)=0.\end{split}$

We get a contradiction.

_Case 4_ : $y_{0}=1$. It follows that
$\partial_{y}\hat{f}_{0}(1^{+},t_{0})\leq 0$, which along with (29c) and (31)
yields

$\displaystyle\begin{split}0&<-k_{0}\hat{f}_{0}(0^{+},t_{0})\leq-k(t_{0})\hat{f}_{0}(0^{+},t_{0})\leq\frac{\sigma^{2}}{2}\partial_{y}\hat{f}_{0}(0^{+},t)-k(t_{0})\hat{f}_{0}(0^{+},t)=0.\end{split}$

We get a contradiction.

_Case 5_ : $y_{0}=z(t_{0})$. It follows that
$\partial_{y}\hat{f}_{0}(z^{-}(t_{0}),t_{0})\leq 0$ and
$\partial_{y}\hat{f}_{0}(z^{+}(t_{0}),t_{0})\geq 0$, which along with (29d)
and $\hat{\sigma}_{0}(t)>0$ yields

$\displaystyle
0\geq\partial_{y}\hat{f}_{0}(z^{-}(t_{0}),t_{0})-\partial_{y}\hat{f}_{0}(z^{+}(t_{0}),t_{0})=\hat{\sigma}_{0}(t_{0})>0,$

leading to a contradiction.

Because we always obtain a contradiction in each case, we have shown that
$\hat{f}_{0}\geq 0$ over the domain $[0,1]\times[0,T]$, which implies the non-
negativeness property of $f_{0}$ over the domain
$[x_{L},\overline{x}(t)]\times[0,T]$ for all $t\in[0,T]$ and all
$T\in\mathbb{R}_{>0}$.

Because the proof of statement (ii) can proceed in the same way as above, we
omit the details of the proof.

Finally, suppose that statement (iii) fails to be true; then, for any given
$T\in\mathbb{R}_{>0}$ there must be a $t_{0}\in(0,T]$ such that

$\displaystyle
f_{0}(\underline{x}(t_{0}),t_{0})+f_{1}(\overline{x}(t_{0}),t_{0})=0,$

which, along with the non-negativeness property of $f_{0}$ and $f_{1}$,
implies that $f_{0}$ and $f_{1}$ attain their minima at
$(\underline{x}(t_{0}),t_{0})$ and $(\overline{x}(t_{0}),t_{0})$,
respectively. Then, using the same argument as that in _Case 5_ , we obtain a
contradiction. Therefore, statement (iii) holds true. $\hfill\blacksquare$

## References

  * [1] Saeid Bashash and Hosam K Fathy.  Modeling and control of aggregate air conditioning loads for robust renewable power management.  IEEE Trans. Control Syst. Technol., 21(4):1318–1327, Jul. 2013\.
  * [2] D. S. Callaway.  Tapping the energy storage potential in electric loads to deliver load following and regulation, with application to wind energy.  Energy Conv. Manag., 50(5):1389–1400, May 2009.
  * [3] Li-Min Cheng and Yu-Qing Bao.  A day-ahead scheduling of large-scale thermostatically controlled loads model considering second-order equivalent thermal parameters model.  IEEE Access, 8:102321–102334, May 2020.
  * [4] Clark W Gellings.  The Smart Grid: Enabling Energy Efficiency and Demand Response.  Taylor & Francis, London, England, Dec. 2020.
  * [5] A. Ghaffari, S. Moura, and M. Krstić.  Modeling, control, and stability analysis of heterogeneous thermostatically controlled load populations using partial differential equations.  J. Dyn. Syst. Meas. Control, 137:101009, 2015.
  * [6] M. Ghanavati and A. Chakravarthy.  Demand-side energy management by use of a design-then-approximate controller for aggregated thermostatic loads.  IEEE Trans. Control Syst. Technol., 26(4):1439–1448, Jul. 2018\.
  * [7] Haider Tarish Haider, Ong Hang See, and Wilfried Elmenreich.  A review of residential demand response of smart grid.  Renew. Sust. Energ. Rev., 59:166–178, June 2016.
  * [8] Yiguang Hong, Zhong-Ping Jiang, and Gang Feng.  Finite-time input-to-state stability and applications to finite-time control design.  SIAM J. Optim., 48(7):4395–4418, Jan. 2010.
  * [9] J. Hu, J. Cao, M. Z. Chen, J. Yu, Yang Yao, J., S., and T. Yong.  Load following of multiple heterogeneous TCL aggregators by centralized control.  IEEE Trans. Power Syst., 32(4):3157–3167, Jul. 2016.
  * [10] Sandro Iacovella, Frederik Ruelens, Pieter Vingerhoets, Bert Claessens, and Geert Deconinck.  Cluster control of heterogeneous thermostatically controlled loads using tracer devices.  IEEE Trans. Smart Grid, 8(2):528–536, Mar. 2017.
  * [11] H. K. Khalil.  Nonlinear Systems.  Prentice-Hall, Englewood Cliffs, NJ, 3rd edition, 2002.
  * [12] Xiangyu Kong, Bowei Sun, Jian Zhang, Shupeng Li, and Qun Yang.  Power retailer air-conditioning load aggregation operation control method and demand response.  IEEE Access, 8:112041–112056, May 2020.
  * [13] Jianming Lian, Huiying Ren, Yannan Sun, and Donald J Hammerstrom.  Performance evaluation for transactive energy systems using double-auction market.  IEEE Trans. Power Syst., 34(5):4128–4137, Sept. 2019.
  * [14] Mingxi Liu and Yang Shi.  Model predictive control of aggregated heterogeneous second-order thermostatically controlled loads for ancillary services.  IEEE Trans. Power Syst., 31(3):1963–1971, May 2016.
  * [15] Weiliang Liu, Changliang Liu, Yongjun Lin, Kang Bai, and Liangyu Ma.  Interval multi-objective optimal scheduling for redundant residential microgrid with vess.  IEEE Access, 7:87849–87865, May 2019.
  * [16] Francisco Lopez-Ramirez, Denis Efimov, Andrey Polyakov, and Wilfrid Perruquetti.  Finite-time and fixed-time input-to-state stability: Explicit and implicit approaches.  Syst. Control. Lett., 144:104775, 2020.
  * [17] N. Lu, D. P. Chassin, and S. E. Widergren.  Modeling uncertainties in aggregated thermostatically controlled loads using a state queuing model.  IEEE Trans. Power Syst., 20(2):725–733, May 2005.
  * [18] Nariman Mahdavi, Julio H Braslavsky, Maria M Seron, and Samuel R West.  Model predictive control of distributed air-conditioning loads to compensate fluctuations in solar power.  IEEE Trans. Smart Grid, 8(6):3055–3065, Nov. 2017.
  * [19] R Malhame and Chee-Yee Chong.  Electric load model synthesis by diffusion approximation of a high-order hybrid-state stochastic system.  IEEE Trans. Automat. Contr., 30(9):854–860, Sept. 1985.
  * [20] Johanna L Mathieu, Stephan Koch, and Duncan S Callaway.  State estimation and control of electric loads to manage real-time energy imbalance.  IEEE Trans. Power Syst., 28(1):430–440, Feb. 2012.
  * [21] W Mendieta and C A Cañizares.  Primary frequency control in isolated microgrids using thermostatically controllable loads.  IEEE Trans. Smart Grid, 12(1):93–105, Jan. 2021.
  * [22] S. Moura, J. Bendtsen, and V. Ruiz.  Observer design for boundary coupled PDEs: Application to thermostatically controlled loads in smart grids.  In 52nd IEEE Conference on Decision and Control, pages 6286–6291, Florence, Italy, 10–13 Dec. 2013.
  * [23] S. Moura, V. Ruiz, and J. Bendtsen.  Modeling heterogeneous populations of thermostatically controlled loads using diffusion-advection PDEs.  In ASME DSC 2013, page V002T23A001, Palo Alto, California, USA, 21–23 Oct. 2013.
  * [24] Md Salman Nazir and Ian A Hiskens.  A dynamical systems approach to modeling and analysis of transactive energy coordination.  IEEE Trans. Power Syst., 34(5):4060–4070, Sept. 2018.
  * [25] C Perfumo, E Kofman, J H Braslavsky, and J K Ward.  Load management: Model-based control of aggregate power for populations of thermostatically controlled loads.  Energy Convers. Manag., 55:36–48, Mar. 2012.
  * [26] N Ruiz, I Cobelo, and J Oyarzabal.  A direct load control model for virtual power plant management.  IEEE Trans. Power Syst., 24(2):959–966, May 2009.
  * [27] Borhan M Sanandaji, He Hao, and Kameshwar Poolla.  Fast regulation service provision via aggregation of thermostatically controlled loads.  In 2014 47th Hawaii International Conference on System Sciences, pages 2388–2397, Waikoloa, HI, USA, 6–9 Jan. 2014.
  * [28] Pierluigi Siano.  Demand response and smart grids: A survey.  Renew. Sust. Energ. Rev., 30:461–478, Feb. 2014.
  * [29] Meng Song, Ciwei Gao, Mohammad Shahidehpour, Zhiyi Li, Jianlin Yang, and Huaguang Yan.  State space modeling and control of aggregated TCLs for regulation services in power grids.  IEEE Trans. Smart Grid, 10(4):4095–4106, Jul. 2019.
  * [30] Luminita Cristiana Totu, Rafael Wisniewski, and John Leth.  Demand response of a TCL population using switching-rate actuation.  IEEE Trans. Control Syst. Technol., 25(5):1537–1551, Sept. 2017\.
  * [31] John S Vardakas, Nizar Zorba, and Christos V Verikoukis.  A survey on demand response programs in smart grids: Pricing methods and optimization algorithms.  IEEE Commun. Surv. Tutor., 17(1):152–178, Jul. 2015.
  * [32] Chon Hou Wai, Marc Beaudin, Hamidreza Zareipour, Antony Schellenberg, and Ning Lu.  Cooling devices in demand response: A comparison of control methods.  IEEE Trans. Smart Grid, 6(1):249–260, Jan. 2015.
  * [33] Ellen Webborn.  Natural heterogeneity prevents synchronization of fridges with deterministic frequency control.  IEEE Access, 7:130206–130214, June 2019.
  * [34] Helin Xu, Lin Cheng, Ning Qi, and Xuyan Zhou.  Peak shaving potential analysis of distributed load virtual power plants.  Energy Rep., 6:515–525, Dec. 2020.
  * [35] W Zhang, K Kalsi, J Fuller, M Elizondo, and D Chassin.  Aggregatemodel for heterogeneous thermostatically controlled loads with demand response.  In IEEE PES General Meeting, pages 1–8. San Diego, CA, 22–26 Jul. 2012.
  * [36] W. Zhang, J. Lian, C. Y. Chang, and K. Kalsi.  Aggregated modeling and control of air conditioning loads for demand response.  IEEE Trans. Power Syst., 28(4):4655–4664, Jul. 2013.
  * [37] L. Zhao and W. Zhang.  A unified stochastic hybrid system approach to aggregate modeling of responsive loads.  IEEE Trans. Autom. Control, 63(12):4250–4263, Dec. 2018.
  * [38] Lin Zhao, Wei Zhang, He Hao, and Karanjit Kalsi.  A geometric approach to aggregate flexibility modeling of thermostatically controlled loads.  IEEE Trans. Power Syst., 32(6):4721–4731, Nov. 2017.
  * [39] J. Zheng, G. Laparra, G. Zhu, and M. Li.  Aggregate power control of heterogeneous TCL populations governed by Fokker–Planck equations.  IEEE Trans. Control Syst. Technol., 28(5):1915–1927, Sept. 2020\.
  * [40] Xiaoming Zhou, Maosheng Sang, Minglei Bao, Sheng Wang, Wenqi Cui, Chengjin Ye, and Yi Ding.  Exploiting integrated demand response for operating reserve provision considering rebound effects.  IEEE Access, 10:15151–15162, Jan. 2022.


$N$-sociable Heronian triangles Nart Shalqini Franklin and Marshall College

###### Abstract

A _Heronian_ triangle is a triangle that has integer side lengths and integer
area. Praton and Shalqini [1] define amicable Heronian triangles to be two
Heronian triangles where the area of one equals the perimeter of the other,
and vice versa; analogous to the concept of amicable numbers. In this paper,
we generalize this notion to $n$_-sociable_ Heronian triangles and
characterize them for each $n$.

## Introduction

A _Heronian_ (or _Heron_) triangle is a triangle that has integer side lengths
and integer area. The name comes from the ancient Greek mathematician Heron of
Alexandria, who first provided the formula for the area of a triangle solely
in terms of its side lengths,

$\displaystyle A=\sqrt{s(s-a)(s-b)(s-c)},$ (1)

where $a,b,c$ are the side lengths of the triangle and $s$ is the semi-
perimeter i.e., $s=(a+b+c)/2$. We define a Heronian triangle to be _equable_
if its area equals its perimeter. It turns out that there are only five
equable triangles [2], namely the triangles with side lengths
$(5,12,13),(6,8,10),(6,25,29),(7,15,20)$, and $(9,10,17)$. The last equable
triangle will be of particular interest to us.

Equable triangles are sometimes called _superhero_ triangles as featured in
[3]. They are analogous to the concept of perfect numbers; numbers that are
equal to the sum of their proper divisors (aliquot sum), for example $6$ and
$28$. Praton and Shalqini [1] define an _amicable_ pair of triangles in
similar fashion to _amicable_ numbers, which are numbers whose aliquot sum
equal each other. The smallest pair of amicable numbers is $(220,284)$.
Although there are finitely many equable and amicable Heronian triangles [1],
the case for perfect and amicable integers remains an open question.

The next natural thing that comes after perfect and amicable numbers are
_sociable_ numbers. An $n$_-sociable_ cycle is a set of numbers whose aliquot
sums form a cyclic sequence. In other words, a set of integers
$a_{1},a_{2},\dots,a_{n}$ is an $n$_-sociable_ cycle if
$s(a_{1})=a_{2},s(a_{2})=a_{3},\dots,s(a_{n})=a_{1}$ where $s(a_{i})$
represents the aliquot sum of $a_{i}$. Note that the amicable numbers are a
special case that occurs when $n=2$. So far, research has largely focused on
perfect and amicable numbers, and not much is known about sociable numbers. In
this paper we define $n$_-sociable_ Heronian triangles analogously to
$n$_-sociable_ numbers. Surprisingly, the complete characterization of
$n$_-sociable_ Heronian triangles turns out be possible by only using
elementary techniques.

## Definition and proofs

###### Definition 1

Suppose we have $n$ Heronian triangles $H_{1},H_{2},\dots,H_{n}$ such that at
least one of them is not equable. We say that this set is an $n$_-sociable_
Heronian cycle if $A_{1}=p_{2},A_{2}=p_{3},\dots,A_{n}=p_{1}$ where $A_{i}$
and $p_{i}$ are the area and perimeter of $H_{i}$, respectively. Moreover, for
each $k$ such that $1\leq k\leq n$ the Heronian cycles
$H_{1},H_{2},\dots,H_{n}$ and $H_{k},H_{k+1},\dots,H_{n-k}$ are equivalent.

Figure 1: Left: the unique amicable pair. Right: a $3$_-sociable_ cycle.

The amicable pair defined in [1] is equivalent to a $2$_-sociable_ Heronian
cycle. Furthermore, [1] shows that there is a unique such pair of triangles;
namely, the triangles $H_{1}:(9,12,15)$ and $H_{2}:(3,25,26)$ where
$A_{1}=54=p_{2}$ and $A_{2}=36=p_{1}$. In this paper, we give a complete
characterization of $n$_-sociable_ Heronian cycles for all $n$. We show that
all $n$_-sociable_ cycles are only composed of a combination of our amicable
pairs with the equable triangle $(9,10,17)$. This in turn implies that there
are finitely many such cycles for each $n$. We begin by establishing a series
of lemmas and theorems necessary for our proof.

###### Lemma 1

Each $n$_-sociable_ Heronian cycle contains a triangle whose perimeter is
strictly greater than its area.

Proof: Let $H_{1},H_{2},\dots,H_{n}$ be a cycle of $n$_-sociable_ Heronian
triangles. By definition, there exists some $H_{i}$ such that $A_{i}\neq
p_{i}$. Without loss of generality, we can assume that $i=1$. If we have some
other triangle $H_{j}$ ($j\neq 1$) such that $p_{j}>A_{j}$ then we are done.
If not, then $A_{1}=p_{2}\leq A_{2}=p_{3}\leq A_{3}=\dots\leq A_{n}=p_{1}$.
This implies that $A_{1}<p_{1}$ since $H_{1}$ is not equable, proving our
claim. QED

###### Lemma 2

Let $H$ be a Heronian triangle in an $n$_-sociable_ Heronian cycle with
perimeter $p$ and area $A$. Then $A^{2}/p$ is an integer.

Proof: First we show that for any Heronian triangle, the semi-perimeter is an
integer. Let $a,b,c$ be the side lengths of $H$. Suppose $s=(a+b+c)/2$. We let
$x=s-a,y=s-b$ and $z=s-c$. We can further suppose $z\geq y\geq x$. Note that
$s=x+y+z$, and that $x+y=c$, $x+z=b$ and $y+z=a$. By Heron’s formula,
$\sqrt{sxyz}=A$ which implies $sxyz=A^{2}\in\mathbb{N}$. Note that $s$ must be
an integer. Otherwise $s$ would be a half-integer, meaning that $x,y,z$ are
all half integers. This implies that $sxyz$ equals an odd integer divided by
$16$; an impossibility given the area is an integer.

Assuming $H$ is in an $n$_-cycle_ , we know that $A=p^{\prime}$ where
$p^{\prime}$ is the perimeter of the proceeding Heronian triangle. Since
$p^{\prime}=2s^{\prime}$ and $s^{\prime}$ is an integer, then $A$ is even.
This implies at least one of $x,y,z$ is even. Therefore

$\displaystyle\frac{A^{2}}{p}=\frac{sxyz}{2s}=\frac{xyz}{2}\in\mathbb{N}.$ (2)

QED

###### Theorem 1

Let $H$ be a triangle like the one in Lemma 2. Then, $z$ divides
$2^{2^{n}}(x+y)^{2^{n}-1}$.

Proof: Suppose $H_{1},H_{2},\dots,H_{n}$ is an $n$_-sociable_ Heronian cycle.
Without loss of generality, assume that $H=H_{1}$. From Lemma 2 we know that
for each $i$, $A_{i}^{2}=c_{i}p_{i}$ for some $c_{i}\in\mathbb{N}$. Hence

$\displaystyle p_{1}^{2^{n}}$
$\displaystyle=\left(A^{2}_{n}\right)^{2^{n-1}}=\left(c_{n}p_{n}\right)^{2^{n-1}}=c^{2^{n}-1}_{n}\left(p^{2}_{n}\right)^{2^{n-2}}$
$\displaystyle=c_{n}^{2^{n}-1}\left(A^{2}_{n-1}\right)^{2^{n-2}}=c^{2^{n}-1}\left(c_{n-1}p_{n-1}\right)^{2^{n-2}}=\dots$
$\displaystyle=c^{2^{n-1}}_{n}c^{2^{n-2}}_{n-1}\cdots c_{2}^{2}A_{1}^{2}.$

Thus $p_{1}^{2^{n}}/A_{1}^{2}\in\mathbb{N}$, and we have

$\displaystyle\frac{p_{1}^{2^{n}}}{A_{1}^{2}}=\frac{2^{2^{n}}s^{2^{n}}}{sxyz}=2^{2^{n}}\frac{s^{2^{n}-1}}{xyz}\in\mathbb{N}$
$\displaystyle\Rightarrow 2^{2^{n}}\frac{s^{2^{n}-1}}{z}\in\mathbb{N}.$ (3)

If $c=x+y$, then by (3) the Binomial Theorem yields

$\displaystyle 2^{2^{n}}\frac{(c+z)^{2^{n}-1}}{z}\in\mathbb{N}$
$\displaystyle\Rightarrow
2^{2^{n}}\frac{\sum_{k=0}^{{2^{n}-1}}\binom{2^{n}-1}{k}z^{k}c^{2^{n}-1-k}}{z}\in\mathbb{N}$
$\displaystyle\Rightarrow
2^{2^{n}}\frac{\sum_{k=1}^{2^{n}-1}\binom{2^{n}-1}{k}z^{k}c^{2^{n}-1-k}}{z}+2^{2^{n}}\frac{c^{2^{n}-1}}{z}\in\mathbb{N}$
$\displaystyle\Rightarrow
2^{2^{n}}\sum_{k=1}^{2^{n}-1}\binom{2^{n}-1}{k}z^{k-1}c^{2^{n}-1-k}+2^{2^{n}}\frac{c^{2^{n}-1}}{z}\in\mathbb{N}.$

The first summand is an integer, which implies that the second summand must
also be an integer. This is possible if and only if $z\ |\
2^{2^{n}}c^{2^{n}-1}$, completing our proof. QED

###### Lemma 3

Assume that $H$ is a Heronian triangle such that its perimeter $p$ is greater
than its area $A$. Then $x\leq 3$ and $y\leq 9$.

Proof: This follows from Lemma 3 in [1]. QED

Consequently, Theorem 1 implies that there are finitely many $n$-sociable
Heronian triangles whose perimeter is greater than their area, since $z$ has
to divide finitely many integers; namely the integers
$2^{2^{n}}(x+y)^{2^{n}-1}$ where $x\leq 3$ and $y\leq 9$. The next theorem
shows that there are at most four of them.

###### Theorem 2

Let $H$ be a Heronian triangle in an $n$_-sociable_ cycle such that its
perimeter is greater than its area. Then $H$ must be one of the following
triangles: $(3,4,5),(5,5,8),(3,25,26)$ and $(3,865,866)$

Proof: First, we should note that for $H$, we have

$\displaystyle p>A\Rightarrow 2s>A\Rightarrow 4s^{2}>A^{2}\Rightarrow
4s^{2}>sxyz\Rightarrow 4(x+y+z)>xyz.$ (4)

From Lemma 3 it follows that $x=1$, $x=2$, or $x=3$ so we treat each of these
cases separately. We assume, as previously, that $x\leq y\leq z$.

We begin with $x=3$. If $y\geq 4$, then by (4) we have $4(z+z+z)\geq
4(x+y+z)>xyz=3yz\geq 12z$. Hence $12z>12z$; a contradiction. So $y=3$, and
$z=3$ or $z=4$ since $4(3+3+z)>xyz=9z$. Neither of these values yield an
integer area.

Next, we consider $x=2$. For $y>5$, we have $xy$ is at least $12$, so
$12z>12z$ and a contradiction. If $y=5$, then $4(2+5+z)>10z$, so $28>6z$ and
$z<5$, a contradiction arises again. If $y=2$ then $sxyz=A^{2}$ which implies
$(4+z)4z=A^{2}$. Hence $16=(2z+4-A)(2z+4+A).$ After checking all possible
factorizations of $16$, we conclude that there are no positive integers
solutions for $z$. Supposing $y=3$ we get $z<10$. Then $z=4,5,8$ since $z$
divides $2^{2^{n}}(x+y)^{2^{n}-1}=2^{2^{n}}(5)^{2^{n}-1}$ by Theorem 1. None
of these values yield an integer area. Similarly if $y=4$ then $z<6$, which
means $z=4,5$. Neither of these values work either.

Suppose that $x=1$. By Lemma 3, $y\leq 9$. From (4), we get an upper bound on
$z$ such that $z<4(y+1)/(y-4)$ when $y>4$. If $y=9$, then $z<9=y$, a
contradiction. If $y=8$, then $z=8$ which yields an irrational area. If $y=7$,
then $7\leq z\leq 10$, none of which yield a rational area. If $y=6$ then
$6\leq z\leq 13$. The only values of $z$ that divide $2^{2^{n}}\cdot
7^{2^{n}-1}$ are $7$ and $8$, but both of them yield an irrational area.
Finally if $y=5$ then $5\leq z\leq 23$. The only values of $z$ that divide
$2^{2^{n+1}-1}3^{2^{n}-1}$ are $6,8,9,12,16$ and $18$. None of these values
work either.

We need to take more care for $y\leq 4$. Suppose that $y=4$. Again, we want
$4z(5+z)=A^{2}$ for $A,z\in\mathbb{N}$. As before, we get
$25=(2z+5-A)(2z+5+A)$. The only solution is $z=4$ and $A=12$, so $x=1,y=4,z=4$
is a candidate for $H$. If $y=1$ then $z(z+2)=A^{2}$ which is not possible
since $A^{2}=z(z+2)=(z+1)^{2}-1$ and no two perfect squares differ by one
except the trivial case $z=0$ and $A=0$.

Next, supposing $x=1$ and $y=3$ we get $3z(z+4)=A^{2}$. Also note that $z\ |\
2^{2^{n}}(1+3)^{2^{n}-1}\Rightarrow z\ |\ 2^{2^{n}}4^{2^{n}-1}$. Therefore,
$z=2^{a}$ for $a\in\mathbb{N}$. Clearly, $z\geq 3$, so $a\geq 2$. Hence
$3\cdot 2^{a}(2^{a}+4)=3\cdot 2^{a+2}(1+2^{a-2})=A^{2}$. Note that $a$ cannot
be odd; otherwise, the quantity $2^{a+2}(1+2^{a-2})$ would have an odd number
of twos and hence wouldn’t be a perfect square. If $a$ is even, then
$1+2^{a-2}=3c^{2}$ where $c$ is an integer. Thus $3c^{2}=1+2^{a-2}=1+4\cdot
2^{a-4}$ for $a\geq 4$ or $3c^{2}=2$ for $a=2$. Considering the equations
modulo $4$, we see each case is impossible. Hence $y=3$ also leads to a
contradiction.

Finally, suppose $x=1$ and $y=2$. Then $z$ has to divide
$2^{2^{n}}(1+2)^{2^{n}-1}=2^{2^{n}}3^{2^{n}-1}$, so $z=2^{a}3^{b}$ for some
$a,b\in\mathbb{N}$. Furthermore, $A^{2}=2z(z+3)=2\cdot
2^{a}3^{b}(2^{a}3^{b}+3)=2^{a+1}3^{b}(2^{a}3^{b}+3)$. If $b=0$ we get
$2^{a+1}(2^{a}+3)=A^{2}$, which implies that $a$ is odd and $2^{a}+3=c^{2}$.
We can further assume that $a\geq 2$ ($z=2$ does not yield a rational area).
The previous expression becomes $4\cdot 2^{a-2}+3=c^{2}$ which is impossible
modulo $4$, so we can suppose that $b\geq 1$ and get
$A^{2}=2^{a+1}3^{b+1}(2^{a}3^{b-1}+1)$. Now consider $b\geq 2$. Then,
$2^{a}3^{b-1}+1=c^{2}$ where $a,b$ are odd. Hence $2^{a}3^{b-1}=(c-1)(c+1)$
which implies $2^{p}3^{q}-2^{r}3^{s}=2$ where $p+r=a$ and $q+s=b-1$. So
$2^{p-1}3^{q}-2^{r-1}3^{s}=1$ where $p-1\geq 0$ and $r-1\geq 0$; otherwise we
have an even number equal an odd number. In turn, this means that $r-1=0$ and
$q=0$ ($2^{p-1}-3^{s}=1$) or $p-1=0$ and $s=0$ ($3^{q}-2^{r-1}=1$). Otherwise
the difference is either even or a multiple of three. Since $b\geq 2$, either
$s\geq 1$ or $q\geq 2$. The solution to such equations is well-known and first
proven by the medieval mathematician Gersonides [4] using arithmetic modulo
$8$. The only solution satisfying our conditions is $q=2$, $r-1=3$ or
$p-1=2,s=1$. The latter yields $a=4$ and $b=2$ which is not possible because
$a$ and $b$ are both odd. The former yields $a=5$ and $b=2$ so $z=864$.
Assuming $b=1$ we have $A^{2}=2^{a+1}3^{2}(2^{a}+1)$. This implies that either
$a=0$ which yields $z=3$ or that $a$ is odd and $2^{a}+1=c^{2}$ for some
$c\in\mathbb{N}$. This in turn implies that $a=3$ and $z=24$ using a similar
argument as before.

We have just shown that the only possible values of $x,y,z$ are $x=1,y=4,z=4$,
or $x=1,y=2,z=3$, or $x=1,y=2,z=24$, or $x=1,y=2,z=864$. These values yields
the following triangles: $(5,5,8),(3,4,5),(3,25,26)$ and $(3,865,866)$.

QED

From Lemma 1 we know that each $n$_-sociable_ Heronian cycle must contain a
triangle whose perimeter is greater than its area. Hence, theorem 2 provides
us with a possibility of characterizing all $n$_-sociable_ Heronian cycles.
The next step is to investigate each of these four triangles and see whether a
sociable cycle contains them, and if so, what kind of other triangles must it
contain. It turns out that three of the four cannot be part of this cycle,
except for the one that belongs to our amicable pair.

###### Theorem 3

There is a unique triangle $H$ that belongs to a $n$-_sociable_ cycle such
that its perimeter is greater than its area; namely, the triangle with side
lengths $3$, $25$ and $26$. Moreover, each $n$-_sociable_ cycle is composed
only of the amicable pairs and the triangle $(9,10,17)$.

Proof: Suppose $H$ belongs to an $n$-_sociable_ cycle such that its perimeter
is greater than its area, and assume it is the first member of our cycle i.e.,
$H=H_{1}$. Then there are finitely possibilities for $H_{1}$ as listed in
Theorem 2. Let $H_{1}$ be the triangle with side lengths $(5,5,8)$. This
triangle has area $12$ and perimeter $18$. The semi-perimeter of the successor
triangle $H_{2}$, if it exists, must be $6$. So for the successor triangle
$x=1,y=1,z=4$ or $x=1,y=2,z=3$, or $x=2,y=2,z=2$. The first and third do not
produce integer areas and the second triangle is our other candidate; namely,
the triangle with side lengths $(3,4,5)$. The successor of this triangle must
have perimeter $6$, and there is no Heronian triangle with such perimeter by a
similar argument as above. Note that these arguments are the same as the ones
used in [1]. The triangle with side lengths $(3,865,866)$ needs to be
considered with more care. Since this triangle has perimeter $1734$, then the
preceding triangle $H_{n}$ must have area $1734=2\cdot 3\cdot 17^{2}$. So, for
the preceding triangle, if it exists, $2^{2}\cdot 3^{2}\cdot 17^{4}=sxyz$. The
fact that $s=x+y+z$ helps cut down the cases considerably and we find that the
only possible values are $z=51,y=34$ and $x=17$. Hence $H_{n}$ has perimeter
$204$. So, $H_{n-1}$ must have area $204$, which by a similar argument leads
to the triangle $z=17,y=9$, and $x=8$ with perimeter $68$. This implies
$H_{n-2}$ has area $68$, but there is no such Heronian triangle. Since this
process stops, we cannot have a $n$-_sociable_ cycle that contains the
triangle $(3,865,866)$.

The only triangle left is the triangle with side lengths $(3,25,26)$. This
triangle has area $36$ and perimeter $54$. Similar computations show that
there is only one triangle with area $54$, namely, the triangle $(9,12,15)$.
This is the amicable pair found in [1]. Therefore the only triangle that
precedes $(3,25,26)$ is the triangle $(9,12,15)$ that has perimeter $36$. The
same kind of computations show that the only triangles with area $36$ are
$(3,25,26)$ and $(9,10,17)$. Thus, the only possibilities of the predecessor
of $(9,12,15)$ or $(9,10,17)$ are $(9,10,17)$ or the amicable pair. By Lemma
1, all $n$_-sociable_ Heronian cycles must contain the triangle $(3,25,26)$
since this is the only sociable triangle whose perimeter is greater than its
area. This concludes the fact that each $n$_-sociable_ cycle is only composed
of amicable pairs and the triangle $(9,10,17)$. QED

## Conclusion

Let $v$ represent the triangle $(3,25,26)$, $u$ represent the triangle
$(9,12,15)$, and $w$ represent the triangle $(9,10,17)$. Theorem 3 states that
any Heronian cycle is only composed of $u,v$ and $w$. All cycles must contain
$v$, which must be preceded by $u$. Hence, $u$ and $v$ come in pairs.
Additionally, a pair $u,v$ is proceeded or succeeded by the same pair or by
$w$. We immediately construct a Heronian cycle as follows:

$\displaystyle\underbrace{[u,v,w,w,\dots,w]}_{n},$ (5)

where the $i$th component denotes the $i$th triangle. Note that (5) represents
all cycles that have exactly one pair $u,v$ since Heronian cycles are
equivalent with respect to rotation; for example, $[u,v,w,w,w]$, $[w,w,u,v,w]$
and $[w,w,w,u,v]$ all represent the same cycle.

Figure 2: We can replace each _amicable_ pair with a pair of equable triangles
$(9,10,17)$. This is because the preceding triangle of the pair
$(9,12,15)\to(3,25,26)$ must have perimeter $36$ and similarly the successor
triangle must have area $36$

Next, we can replace a pair of equable triangles with our amicable pair (see
figure 2). Thus, we can form an additional $(n-2)/2$ (if $n$ is even) and
$(n-3)/2$ (if $n$ is odd) Heronian cycles by successive replacements in each
step,

$\displaystyle\underbrace{[u,v,u,v,\overbrace{w,\dots,w}^{n-4}]}_{n},\underbrace{[u,v,u,v,u,v,\overbrace{w\dots
w}^{n-6}]}_{n},\dots$

Starting with our cycle in (5) and applying this process we get $\lfloor
n/2\rfloor$ $n$_-sociable_ distinct Heronian cycles for each $n$. However,
this list does not contain all $n$_-cycles_ as for instance $[u,v,u,v,w,w]$ is
a different cycle from $[u,v,w,u,v,w]$. Nevertheless, we can use this process
to identify all $n$_-cycles_ for each $n$ by considering all the possibilities
that arise from having $k$ pairs $u,v$ and $n-2k$ $w$s for $1\leq k\leq\lfloor
n/2\rfloor$. For example, there are only two $5$-_sociable_ Heronian cycles:
$[u,v,u,v,w]$ and $[u,v,w,w,w]$, as shown in figure 3 .

Figure 3: All possible $5$-_sociable_ cycles.

One can claim that there essentially are no non-trivial $n$-_sociable_ cycles
except for $n=2$. All of the other cycles are formed from _amicable_ pairs and
the equable triangle with area and perimeter $36$. Surprisingly, the complete
characterization of $n$-_sociable_ cycles turns out to be provable by only
using classical number theory and elementary methods as opposed to some more
complex methods such as those involving algebraic number theory used in [5]
for other questions on Heronian triangles.

## References

  * [1] Praton and Shalqini, Amicable Heron Triangles, arXiv:2101.09768. Note: this paper has been accepted for publication and will appear in the second issue of Fibonacci Quarterly for 2021.
  * [2] L. Dickson, History of the Theory of Numbers, Volume II. Dover publications (2005).
  * [3] B. Haran, Superhero triangles,
https://www.youtube.com/watch?v=UIjeCKPHbso

  * [4] C. Merzbach and B. Boyer, A history of mathematics 3rd ed. John, Wiley, & Sons, Inc., (2011)
  * [5] Y. Hirakawa and H. Matsumura, A unique pair of triangles. J. Number Theory, 194(2019), 297-302

## About the author:

Nart Shalqini is an international student from Kosovo; a young and small
country in the Balkans. He was born in Istanbul in 1999 and it is his favorite
place to visit. Besides Mathematics, Nart enjoys playing board games and
reading mystery novels in his spare time.

### Primus Scriber

Franklin and Marshall College
Lancaster, Pennsylvania, 17603
<EMAIL_ADDRESS>