1

A Unified Stochastic Hybrid System Approach to
Aggregate Modeling of Responsive Loads

arXiv:1503.06911v3 [cs.SY] 18 Aug 2017

Lin Zhao and Wei Zhang

Abstract—Aggregate load modeling is of fundamental importance for systematic analysis and design of various demand
response strategies. Instead of keeping track of the trajectories
of individual loads, the aggregate modeling problem focuses
on characterizing the density evolution of the load population.
Most existing models are only applicable to Thermostatically
Controlled Loads (TCL) with first-order linear dynamics. This
paper develops a unified aggregate modeling approach that can
be used for general TCLs as well as deferrable loads. We propose
a deterministic hybrid system model to describe individual load
dynamics under demand response rules, and develop a general
stochastic hybrid system (SHS) model to capture the population
dynamics. We also derive a set of partial differential equations
(PDE) that governs the probability density evolution of the SHS.
Our results cannot be obtained using the exiting SHS tools in the
literature as the proposed SHS model involves both random and
deterministic switchings with general switching surfaces in multidimensional domains. The derived PDE model includes many
existing aggregate load modeling results as special cases and can
be used in many other realistic modeling scenarios that have not
been studied in the literature.

I. I NTRODUCTION
Many electric loads, such as residential HVACs (Heating,
Ventilation, and Air Conditioning) and PEVs (Plug-in Electric
Vehicles), can be equipped with smart local controllers that
can modify their power consumptions in response to some
coordination signals (e.g. price or power system frequency).
Coordination of a large population of such responsive loads is
an important paradigm of demand response programs and can
provide various services to the power grid, such as power capping, energy arbitrage, and frequency regulation [1], [2], [3],
[4], [5], [6], [7]. Aggregate load modeling is concerned with
characterizing the aggregate dynamics of responsive loads. The
problem is of fundamental importance for systematic analysis
and design of demand response strategies. It is challenging
due to the various kinds of uncertainties about individual load
models and the complex interplay between the load dynamics
and their local response/control rules.
One of the earliest results on aggregate load modeling was
given by Malhamé and Chong in [8]. They considered a population of thermostatically controlled loads (TCLs), where each
TCL was modeled by a two-mode stochastic hybrid system
(SHS) with continuous dynamics in each mode governed by a
scalar linear stochastic differential equation (SDE). Through a
This work was supported in part by the National Science Foundation under
grant ECCS-1309569 and grant CNS-1552838.
L. Zhao and W. Zhang (Corresponding Author) are with the Dept. of
Electrical and Computer Engineering, The Ohio State University, Columbus,
OH, 43210. Email: {zhao.833, zhang.491}@osu.edu

complicated probabilistic argument, a set of coupled FokkerPlanck equations with boundary conditions was derived, which
characterizes the time evolution of the hybrid-state probability
density function (p.d.f.).
The PDE (Partial Differential Equation) model derived
in [8] plays a fundamental role in the field of aggregate load
modeling. Many recent studies in this area can be viewed as
modifications, discretizations, approximations, or applications
of the PDE model in [8]. For example, the author in [1] derived
the stationary solution of the coupled Fokker-Planck equations.
The solution was then utilized to perform a set point control of
TCLs for demand response. Following this work, the authors in
[9] derived a scalar transport equation for deterministic hybrid
TCL models using the control-volume method from continuum
mechanics. The control design was based on an ordinary
differential equation (ODE) model, which is obtained from
the space-discretization of the PDE. Similar works on deriving
PDE models for TCLs were also reported in [10], [11], [12],
where diffusion terms were used to account for parameter
heterogeneity. In particular, a linear integral output feedback
control algorithm was designed directly based on the PDE
model in [12]. Moreover, Markov chain models were proposed
in [2], [3], [5] to capture the evolution of the temperature
distribution of a TCL population. These models can be viewed
as some approximations of the space-discretization of the
underlying PDE model.
The aforementioned works are only applicable to TCLs
with linear first-order temperature dynamics. Their extensions
to more complex load dynamics are not well understood.
The authors in [4] studied aggregation of second-order TCL
models that involve coupled air and mass temperature dynamics. It was shown that the consideration of the secondorder effect can improve the modeling performance. However,
the aggregate model in [4] was developed mainly based on
heuristic arguments. In fact, when considering higher-order
load dynamics, the boundary conditions of the coupled FokkerPlanck equations can be challenging to obtain, especially when
there are diffusion terms in load models. In addition, many
demand response applications involve deferrable loads, such
as PEVs, washers, dryers, among others. These loads are dramatically different from TCLs. They are similar to computer
jobs, characterized by job size and deadline. When imposing
hard constraints on the deadlines, deferrable loads also exhibit
dynamic behaviors. Aggregate modeling of dynamic deferrable
loads has not been adequately studied in the literature.
This paper studies aggregate modeling of responsive loads.
Different from most existing works that focus on first-order
TCLs models, we develop a unified framework that can

2

be used to obtain aggregate models for general TCLs as
well as deferrable loads. In particular, we propose a general
(deterministic) hybrid system model to capture individual load
dynamics. The proposed hybrid system model is convenient
to describe the complex multi-modal dynamics of responsive
loads induced by demand response strategies. We then consider a large population of such responsive loads described by
hybrid systems. To account for various kinds of uncertainties
at the population level, we develop a SHS model to capture
the aggregate population dynamics. Each mode of the SHS
model is governed by a multi-dimensional nonlinear SDE, and
mode transitions can be triggered by both deterministic and
random switchings. The main technical result of this paper is
the derivation of a set of coupled PDEs (later referred to as
forward equations) and their boundary conditions that govern
the evolution of the hybrid-state p.d.f. of the proposed SHS.
The boundary conditions are particularly challenging to derive
as the proposed SHS involves deterministic mode transitions
characterized by a fairly general class of switching surfaces.
The proposed SHS and the derived forward equations constitute a general aggregate modeling framework for responsive
loads. The main contribution of our result is two-fold. First,
it provides a unified way to obtain aggregate models for a
variety of responsive loads under different demand response
strategies. In particular, it contains the famous result of Malhamé and Chong [8] as a special case. It can also be used
to directly obtain the aggregate PDE models for second-order
TCLs ([4], [13]), which has not been formally derived in the
literature. In addition, it can incorporate new scenarios such as
nonlinear load dynamics, random mode switching (to account
for unmodeled uncertainties on mode transitions), aggregation
with deferrable loads, among others. These features can enable
a larger range of applications for aggregate load modeling.
The second contribution of this paper is on SHS. Our result,
though derived in the context of aggregate load modeling, can
also be used to obtain forward equations of a general SHS
with both deterministic and random switchings. Although SHS
has been studied extensively in the literature [14], [15], [16],
[17], [18], [19], its probability density evolution cannot be
explicitly characterized using the existing results, especially
when there are deterministic switchings with general switching
surfaces in multi-dimensional spaces. It is well known that
for a standard diffusion process, the adjoint of its (strong)
generator determines the form of the forward equation and
the domain of the generator (indirectly) affects the boundary
conditions. Unfortunately, the (strong) generator, including its
domain, of a general SHS is not available in the literature.
The extended generator of a generalized SHS derived in [14],
[19] cannot be directly used to obtain the forward equations
either (See Section V for technical details). Therefore, our
work also contains nontrivial and important extensions of the
existing theoretical results in the field of SHS.
Finally, it is worth mentioning the difference of our result
with respect to two closely related works ([20] and [21]) which
also derived the forward equations for certain classes of SHS
processes. The SHS model considered in [20] is a special
case of Piecewise-Deterministic Markov Process (PDMP) with
a unbounded continuous state space. It does not contain

Brownian motion terms and does not involve deterministic
switchings. As a result, the author does not need to derive
boundary conditions, which are the key challenge of our work.
Reference [21] focuses on a measure-valued formulation of the
forward equation based on the Levy’s identity [22], [23]. Their
derivation does not provide explicit characterization of the
generator boundary conditions, which, however, is the key to
correctly derive the PDE boundary conditions for the forward
equations.
The rest of the paper is organized as follows. Section II
provides an illustrating example for the aggregate load modeling problem. A unified hybrid system model is proposed for
individual responsive loads in Section III. Section IV proposes
a SHS for modeling the load population considering various
uncertainties. Section V develops the theory of the SHS and
derives the corresponding forward equation and boundary
conditions. Section VI discusses the applications of the main
theorems through several examples. The paper is concluded in
Section VII.
Notation: As usual, (Ω, F, P ) denotes the underlying probability space of a stochastic process, where ω ∈ Ω represents a
sample path. E denotes the expectation operation with respect
to an appropriate probability measure. For a set U ⊂ Rn ,
U ◦ , ∂U , Ū , and 1U represent the interior, closure, boundary,
and indicator function of U , respectively. For a matrix M ,
T r(M ), M T , and Mi are the trace, transpose, and ith row of
M , respectively. For vectors x, y ∈ Rn , we denote by x · y
the Euclidean inner product. For a function f : Rn 7→ R,
∂f
interchangeably for the first
we will use Dxi f , fxi , and ∂x
i
order partial derivatives with respect to xi . In addition, we use
∇f and ∇2 f to denote the gradient and Hessian matrix of f ,
respectively. The divergence of a vector field g : Rn 7→ Rn is
denoted by ∇ · g; Similarly, the divergence of a matrix-valued
function M : Rn 7→ Rm×n is denoted by ∇ · M , which is
a vector-valued function whose ith element is ∇ · Mi . Let
C k be the class of real-valued kth-continuously differentiable
functions in some open subset U ∈ Rn , and Cbk ⊂ C k denote
those functions in C k with bounded partial derivatives of order
up to k. For f ∈ C k , we assume as a convention that the value
of f and its kth order partial derivatives on the boundary of
U are defined by their continuous extensions.
II. M OTIVATING E XAMPLE
Consider a population of residential HVAC systems. Each
HVAC can be described by the second-order ETP (Equivalent
Thermal Parameter) model [24], The ETP model is a hybrid
system, which has two modes, representing the “ON/OFF”
power states of the device, and a two-dimensional continuous
state vector, representing the air and mass temperatures of
the house. Each mode is governed by a linear differential
equation and the mode transition is triggered when the air
temperature hits the boundary of the temperature dead band.
Although each HVAC has simple hybrid dynamics, a large
number of these loads may exhibit rather complex aggregate
dynamics under demand management strategies. Fig. 1 shows
the aggregate power response of 2000 HVACs under the socalled thermostat setback program [25], [4]. The parameters

3

10

Total Power (MW)

8

6

4

2

0
2

3

4

5

6
7
Time (Hour)

8

9

10

Fig. 1. GridLAB-D simulation of thermostat setback program of 2000 HVACs

of the ETP model used in the simulation are generated using
GridLAB-D [4]. At the beginning of the simulation, the initial
temperature setpoints of the HVACs are uniformly distributed
in [70◦ F, 78◦ F]. During the setback event starting at time t = 3
hour, all the setpoints are instructed to increase by 1◦ F, which
reduces the steady state power by about 10%. However, after
the setback control is released at time t = 6, a large rebound
is observed which may potentially damage the grid. Such
phenomenon is common in demand response problems. In fact,
the setpoint change described above can also be thought of as
being triggered by a price change for price-responsive loads
or a frequency drop for frequency-responsive loads. Therefore,
developing an aggregate model that can accurately capture the
collective dynamics of a population of responsive loads is of
fundamental importance for demand response. This will be the
main focus on this paper.

III. A U NIFIED H YBRID S YSTEM M ODEL FOR
R ESPONSIVE L OADS
Responsive loads typically have multiple discrete operation
modes. Transitions among these modes are governed by certain
switching logic rules that depend on the evolution of the
local continuous state variables (e.g. temperature of a TCL
or time to completion of a PEV load) as well as the external
control signals (e.g. price or direct load control). The switching
logic rules are designed to meet end-user performances or
energy consumption constraints, such as keep the temperature
within a desired range or finish charging before a given
deadline. To capture the interactions between the discrete and
the continuous states, we propose the following hybrid system
as a unified model for responsive loads:


(continuous dynamics)
ẋ(t) = f (q(t), x(t); θ),
−
−
−
q(t) = φ(q(t ), x(t ), υ(t ); θ), (mode transition)


y(t) = h(q(t); θ),
(output)
(1)
where the continuous state x(t) ∈ Xq given the discrete mode
q(t) ∈ Q , and υ(t) ∈ R is the external control input to
the system which is assumed to only affect the discrete mode
transitions. We assume that Xq◦ is an open subset of Rn

with boundary ∂Xq , Q is a finite subset of Z, and υ ∈ V
which is the space of piecewise constant functions on t ≥ 0.
We define the hybrid state space X := ∪q {q} × Xq . The
boundary of X is defined as ∂X = ∪q {q} × ∂Xq and
the closure of X is X̄ := X ∪ ∂X. For each q ∈ Q,
f (q, ·; θ) is a vector field of the continuous dynamics in mode
q. The mode transition is governed by a transition function
φ, where q(t− ) := lims↑t q(s) and x(t− ) and υ(t− ) are
defined similarly. The output function h(q(t); θ) represents
the power consumption of the load, which typically does not
depend on the continuous state. The hybrid system model is
parameterized by θ ∈ Θ ⊆ Rnθ .
The mode transition function φ is determined by the local
control logics of the responsive load. We assume that the
switching logic is characterized by switching surfaces. The
external control input υ(t) may directly modify the switching
surfaces. We assume that the change of υ(t) is sufficiently slow
as compared to the dynamics of the hybrid system. To avoid
further complicating the discussion, we will not explicitly
model the impact of the dynamics of υ(t). Nevertheless, it
will be shown through the numerical simulation. Therefore, we
will fix υ(t) from now on. A mode transition occurs when the
continuous state x(t) hits the corresponding switching surface
from within Xq . Let ν(q, x) denote the outer unit normal
vector on ∂Xq . We define the outflow switching surface Gq
as
Gq (θ) = {x ∈ ∂Xq : f (q, x; θ) · ν(q, x) > 0} .
(2)
For the responsive loads, we assume that the continuous state
remains the same after the mode switching. In addition, we
assume without loss of generality that for each mode q, there
are at most one pre-jump mode, denoted by q − , and at most
one post-jump mode, denoted by q + . The extension to the case
of multiple pre-jump and post-jump modes is straightforward.
Then we will also call Gq the inflow switching surface to the
mode q + . Moreover, we denote by G = ∪q {q}×Gq the overall
outflow switching surface and S = ∪q {q + } × Gq the inflow
switching surface of the hybrid system. Therefore φ : G 7→ Q
can be simply defined as q(t) = q + when x(t− ) ∈ Gq and
q(t− ) = q, otherwise q(t) = q(t− ).
The hybrid system model described above can be used
to capture the dynamics of various responsive loads under
different demand response strategies. In the following, we will
give two representative examples to illustrate its applications.
Example 1. HVACs
The second-order HVAC model in the thermostat setback
program described in [24], [4] is a special case of the proposed
hybrid sytsem (1). It has two discrete modes (Q = {0, 1}),
representing the “OFF” and “ON” power states of the HVAC,
respectively. Each mode has a linear continuous dynamics:
ẋ(t) = Ax(t) + Bq(t) ,

(3)

where x(t) = [x1 (t), x2 (t)]T are the air and mass temperatures and q(t) ∈ Q is the discrete mode. The matrix A and
Bq are determined by the ETP parameters such as heat flux,
thermal mass of the air, thermal mass of the inner solid mass,

4

Maximum

‫ݔ‬ଶ Deferrable Time

Mode 3

waiting

Mode 2

ON

completed
OFF

Fig. 2. Switching surfaces of the HVAC example

charging

0

‫ݔ‬ଵ

Fig. 3. Hybrid state trajectory of a PEV charging
load with no external control input

and so on. Given a temperature setpoint uset and a deadband
size δ, the continuous state spaces are given by
(

X1 = x ∈ R2 : x1 > uset − δ , ON mode,

X0 = x ∈ R2 : x1 < uset + δ , OFF mode,
and the corresponding switching surfaces are given by G1 =
∂X1 for the ON mode, and G0 = ∂X0 for the OFF mode
(see Fig. 2). The overall system is thus parameterized by θ =
(A, B1 , B2 , uset , δ). The system will switch between “ON”
and “OFF” upon reaching G. Thus the mode transition function
can be written as,


x1 (t− ) ≥ uset + δ,
1,
(4)
q(t) = 0,
x1 (t− ) ≤ uset − δ,

 −
q(t ), otherwise.
For the thermostat setback program, the external control
υ(t) will be the setpoint change command. It modifies the
switching surfaces and thus affects the mode transition. For
price responsive loads, υ(t) can also model the price signals
which trigger the setpoint change. Other TCLs under different
demand response strategies can also be modeled similarly
using the proposed hybrid system framework.
Example 2. Plug-in Electric Vehicle (PEV) Charging
The timing dynamics of the PEV charging tasks can also
be accurately represented by the hybrid system model (1).
For example, the timing dynamics of a PEV charging job
can be described as a hybrid system with three discrete
modes {0, 1, 2}, representing that the load is waiting to be
processed, being actively charging, and has been completed,
respectively. The continuous state is two-dimensional with
x1 (t) representing the remaining time to finish the load if
it is running, and x2 (t) representing the maximum time the
load can be further deferred. The vector field in each mode is
given by f (0, x) = [0, −1]T , f (1, x) = f (2, x) = [−1, 0]T .
Here, we assume that a negative value of x1 (t) indicates that
the job has been completed before.
 The 2continuous state spaces are
 given by X0 =
x ∈ R : x1 > 0, x2 > 0 , X1 = x ∈ R2 : x1 > 0 , and
X2 = x ∈ R2 : x1 ≤ 0 . Transition from “waiting” to
“charging” will be forced to take place when a PEV approaches G0 = {x ∈ R2 : x1 > 0, x2 = 0} from X0 in order
to meet the deadline (see Fig. 3). A PEV finishes charging
when it approaches G1 = {x ∈ R2 : x1 = 0} from X1 . The

Mode 1

Remaining
Charging Time
Fig. 4.
process

Sample trajectories of a typical SHS

external control υ(t) may directly trigger mode transitions at
an earlier stage. For example, a frequency responsive PEV
may stop charging if the frequency deviation (which is the
external input) exceeds a certain threshold. Other deferrable
loads under different local response rules can be modeled in
a similar way.
IV. AGGREGATE L OAD M ODELING U SING S TOCHASTIC
H YBRID S YSTEMS
Consider a large number N of responsive loads, each of
which is modeled as a hybrid system of the form (1) but may
have different parameters θi , 1 ≤ i ≤ N . The population
dynamics are described by
 i
i
i
i

ẋ (t) = f (q (t), x (t); θ ),

q i (t) = φ(q i (t− ), xi (t− ), υ(t− ); θi ),
(5)
y i (t) = h(q i (t); θi ),


PN
 Aggr
y
(t) = i=1 y i (t),
where y Aggr (t) represents the aggregate power output. The
load population is called homogeneous if all the parameters
θi are the same, otherwise it is called heterogeneous.
Due to issues like privacy, limited communication resources,
computational complexities, and infrastructure costs, it is
often difficult to keep track of the hybrid state trajectory
(xi (t), q i (t)) for all the individual loads. Alternatively, it is
natural to take a probabilistic (or distribution) viewpoint for
the overall population because many information regarding
individual loads is uncertain to the aggregator. In the rest of
this section, we will model the population uncertainties and
propose an SHS model for the load population.
A. Population Uncertainties
1) Modeling Errors: The continuous dynamics of responsive loads may depend on many factors that cannot be accurately modeled. One way to account for these unmodeled
dynamics is to add a noise process to the nominal model.
Therefore we assume that the ith load is modeled by
dxi (t) = f (ξ i (t); θi )dt + σ(ξ i (t); θi )dW i (t),

(6)

where ξ i (t) = (q i (t), xi (t)), W i (t) is the standard mdimensional Wiener process, and σ is the dispersion matrix
with appropriate dimensions. In addition, we define the diffusion matrix Σ := σσ T .

5

2) Model Parameter Uncertainties: In contrast to a homogeneous load population, a heterogeneous one with welldiversified load parameters will result in a natural damping
and a more stable aggregate response [1], [2], [9]. From the
aggregator’s perspective, the real model parameters θi in (5)
are not known precisely, but their distribution may be acquired
or estimated. Without loss of generality, we assume θ ∼ pθ .
3) Random Switching: Users’ interference can cause the
spontaneous mode switchings. For example, a user may decide
to turn on/off the HVACs anytime without abiding the local
control rule. In the PEV charging example, the user may
start charging the PEVs before the deadline and drive it
away without being fully charged. This phenomenon can
affect the aggregate power output significantly. It can be
modeled by the so-called random jump mechanism [20], [16],
[15], which is determined by a transition intensity function
λ(q, x) : Q × Rn 7→ R+ , where R+ denotes the non-negative
real numbers. In particular, the probability of a random jump
happening in an infinitesimal time interval (t, t + dt] is given
by λ(q(t), x(t))dt .
B. SHS Model for Load Population
Although each individual load is modeled as a deterministic hybrid system (1), the overall population dynamics are
stochastic due to the aforementioned uncertainties intrinsic at
the aggregator level. In this subsection, we will introduce an
SHS model to describe the aggregate population dynamics.
We assume in (6) that all the initial conditions ξ i (0)0 s
are i.i.d. and {W i (t), 1 ≤ i ≤ N } are independent Wiener
processes which are also independent of {ξ i (0), 1 ≤ i ≤ N }.
Under this assumption, the loads are modeled by independent
stochastic processes. Thus if N is sufficiently large, the
population dynamics can be approximated by the following
SHS process,

dx(t) = f (ξ(t), θ(t))dt + σ(ξ(t), θ(t))dW (t),



dθ(t) = 0,
(7)

q(t) = φ(q(Tk−1 ), x(Tk− ), υ(Tk− ), θ(Tk− )),



∀t ∈ [Tk , Tk+1 )
where θ(0) ∼ pθ , ξ(0) = (q(0), x(0)) ∼ p0 , and we assume
without loss of generality that pθ and p0 are independent.
The special features of the proposed SHS model are worth
a detailed explanation. The first one is on incorporating θ
as part of the continuous state with dθ(t) = 0. In this way,
the parameter heterogeneity can be addressed equivalently by
specifying the initial distribution pθ . Furthermore, the SHS
incorporates both deterministic and random mode switchings
through the transition function φ, where Tk , ∀k ∈ N are
random variables denoting the kth jump instant of the SHS
process. Specifically, φ : X̄ × V × Θ 7→ Q can be simply
defined as q(t) = q + , ∀t ∈ [Tk , Tk+1 ), if q(Tk−1 ) = q, where
Tk depends on x, υ, and θ.
The jump instants Tk , ∀k ∈ N are a special class of random
variables called the stopping time. The following definitions
for Tk extends that of the piecewise deterministic Markov
process (PDMP) [14] to the case of the SHS process. We
first define the stopping time t∗ which triggers a deterministic

jump, and then define Tk by incorporating the random jump.
Let ω k := (q k (ω), xk (ω)) be the trajectory of the SHS
process between [Tk (ω), Tk+1 (ω)), where xk is the continuous
component in ω k . Then the stopping time t∗ (ω k ) which
triggers a deterministic jump in mode q k can be defined as

t∗ (ω k ) = inf t > 0; xk (t, ω) ∈ Gqk .
Now we incorporate the random jump. Let S k be the dwell
time of the process in mode q k , which is defined by
S k (ω k+1 ; ω k ) = inf{t : β(t; ω k ) ≤ u(ω k+1 )},

(8)

where β is the survival function of S k such that β(t) =
P (S k > t), and u is a uniformly distributed random variable
between [0, 1]. The survival function is given by [14],

 ˆ t
k
k
k
λ(q , x (s, ω))ds , (9)
β(t; ω ) = I{t<t∗ (ωk )} exp −
0

where λ is the transition intensity function. Thus the mode
transition of the SHS process occurs at the following stopping
times
Tk+1 = Tk + S k , T0 = 0.
(10)
Without loss of generality, we assume that λ(q, x) is piecewise
continuous, ∀q ∈ Q. This is to ensure that S k > 0 almost
surely (a.s.). Furthermore, it can be seen that the longer the
process has been staying in a mode, the more likely the mode
will change.
A sample path of (7) can be viewed as the trajectory of
a randomly selected load in the population. Fig. 4 illustrates
the sample path of a typical SHS process with three mode.
The boundary curves of the shaded area represent the outflow
switching surfaces, and the dashed curves represent the inflow
switching surfaces. We see that a random jump happens in
mode 1 before the load reaches the outflow switching surface,
and in the other two modes, deterministic jump happens upon
reaching the outflow switching surfaces. Suppose p(q, x, t) is
the p.d.f. of the hybrid state of the SHS process, and Q0 ⊂ Q
are the modes under which the loads are consuming energy.
Then the aggregate power response can be calculated by
X ˆ ˆ
y Aggr (t) = N
h(q, θ)p(q, x, t)pθ dxdθ.
q∈Q0

Θ

Xq

V. PDE C HARACTERIZATION OF AGGREGATE DYNAMICS
As discussed in the previous section, the load population
can be modeled as an SHS of the form (7). For systematic
analysis, it is often desired or even necessary to characterize
the evolution of the hybrid-state probability density function
of the SHS, which is the main focus of this section. In fact,
most of the existing aggregate load models in the literature
are essentially characterizing the density evolutions of some
simplified versions of our general SHS model (7).
Our main results will be derived for the case of homogeneous parameters first. The extension to the case of
heterogeneous parameters will be discussed at the end of this
section. We therefore drop the dependence of θ for now and
focus on the following SHS,

6



dx(t) = f (q(t), x(t))dt + σ(q(t), x(t))dW (t),
q(t) = φ(q(Tk−1 ), x(Tk− ), υ(Tk− )),


∀t ∈ [Tk , Tk+1 )

(11)

where note that the definition of Tk depends on the transition
intensity function λ and the outflow switching surface G.
For a standard diffusion process, its density evolution is
characterized by the Fokker-Planck equation, which can be
derived using the infinitesimal generator and the associated
Dynkin’s formula [26]. A similar idea can be used to derive
a set of coupled PDEs characterizing the density evolution
for a switching diffusion process [8]. However, characterizing
the density evolution for the general SHS model (11) is much
more challenging. The main challenge lies in the deterministic
switchings that are forced to occur when the continuous
state hits certain switching surfaces. These switchings, along
with the possible random switchings and diffusion noises,
significantly complicate the boundary conditions of the resulting density PDEs. In the rest of the section, we will first
review some important concepts and the key steps in deriving
forward equations of a general Markov process. Then we will
generalize the existing results and derive the density evolution
PDEs for the homogeneous SHS (11). The last subsection will
introduce a method to deal with the general heterogeneous
case.
A. Preliminaries on Markov Processes
We first review some classical concepts and results that
are useful in deriving the Kolmogorov forward equations of a
general Markov process. The discussion in this subsection will
not only set up stages for our main results, but also allow us
to see the subtle technical challenges in deriving the forward
equations for the general SHS model.
Let ξt , t ≥ 0 be a (time homogeneous) Markov process
defined on X, whose transition probability measure is denoted
by % : R+ × X × B(X) 7→ [0, 1]. Let Pξ0 be the Wiener
probability measure of ξt such that Pξ0 (ξt ∈ A) = %(t, ξ0 , A),
∀A ∈ B(X). Given µ the initial distribution of ξ0 , the
´ abstract
probability measure P can be related to Pξ0 by P = X Pξ0 dµ
. If µ and P admits a probability density function p0 and
p respectively, then we have the relation dµ = p0 (ξ)dξ and
dP = p(ξ)dξ.
For a real-valued bounded Borel measurable function ψ
on X, we can define the semi-group Zt ψ := Eξ0 ψ(ξt ),
t ≥ 0, where the expectation Eξ0 is with respect to Pξ0 . The
(infinitesimal or strong) generator of the process ξt is defined
as follows:
Definition 1. (Strong Generator) Let D(L) denote the set of
bounded Borel measurable functions ψ :X 7→ R with the
property that ∀ψ ∈ D(L), the limit
Lψ := lim
t↓0

Zt ψ − ψ
,
t

(12)

exists in the supremum norm k·k := supξ∈X |·|. Then we call
(L, D(L)) the strong generator of the process ξt .
Furthermore, from the Markov process theory [14], we
know:

Lemma 1. Each ψ ∈ D(L) is associated with a martingale
defined as,
ˆ t
ψ
Ct = ψ(ξt ) − ψ(ξ0 ) −
Lψ(ξs )ds,
(13)
0

which satisfies the Dynkin’s formula,
ˆ t

Eξ0 ψ(ξt ) = ψ(ξ0 ) + Eξ0
Lψ(ξs )ds .

(14)

0

The Dynkin’s formula defines the time evolution of the
expectation of a function ψ of the stochastic processes. It
plays a major role in deriving the PDE model. Let Cc∞ (X)
denote the space of real-valued smooth functions on X with
compact support, h, i denote the inner product of L2 (X) of
square integrable functions on X. For a (linear) differential
operator V , we will say V ∗ is the L2 formal adjoint of L if
V ψ1 , ψ2 = ψ1 , V ∗ ψ2 ,

(15)

for all smooth functions ψ1 , ψ2 ∈ Cc∞ (X). Note that V ∗
can be calculated by shifting the differential operator from ψ1
to ψ2 using integration by parts or divergence theorem (see
Lemma 5 in the Appendix D).
Lemma 2. Let (L, D(L)) be the strong generator of the
Markov process ξt on X and p(ξ, t) be its probability density
function. Suppose L is a differential operator and X has
a piecewise smooth boundary. Then, in the sense of weak
derivatives, p(ξ, t) satisfies the PDE
∂p
= L∗ p, t ≥ 0,
∂t
with boundary conditions uniquely determined by D(L).
Proof: For smooth functions ψ ∈ D(L), using the
divergence theorem, we can write formally
∂p
∂p
− L∗ p, ψ =
, ψ − p, Lψ + Γ(ψ, p),
∂t
∂t
where Γ is the surface integral on ∂X induced by the differential operator L. Then the PDE holds in the sense of weak
derivatives [27], if
(
∂p
∂t , ψ = p, Lψ ,
(16)
Γ(ψ, p) = 0,
for all the smooth functions ψ ∈ D(L). Let Λ(L) ⊃ D(L)
denote the set of the boundary conditions satisfied by ψ ∈
D(L), then the boundary conditions Λ(L∗ ) is defined as the
minimal set of (homogeneous) conditions such that Γ(ψ, p) =
0 for all ψ ∈ Λ(L) and p ∈ Λ(L∗ ) (see [28, page 103]). Note
that Λ(L∗ ) is uniquely determined by Λ(L) by this definition.
Furthermore, we know from Lemma 1 that the Dynkin’s
formula holds for all ψ ∈ D(L). Taking the time derivative of
both sides of (14), we have ∂p
∂t , ψ − p, Lψ = 0. Thus it
completes the proof.
We see from Lemma 2 that the PDE is determined by the
formal adjoint operator of L, and the associated boundary
conditions are uniquely determined by the boundary conditions
in D(L). Lemma 2 has been used to obtain the PDEs for
diffusion processes [26] on subsets of Rn with different kinds

7

of boundaries. Similar idea has also been used to derive the
PDE for the jump-diffusion processes [29].
For the SHS process defined on a hybrid state space X,
characterizing the strong generator (L, D(L)) can be very
difficult. However, for the purpose of deriving the PDE, we
can instead characterize an extended operator of the strong
generator. Before proceeding, we will introduce the following
assumptions for the later discussion.
B. Standing Assumptions
In this section, we will introduce some regularity conditions.
The first two assumptions are imposed to ensure that the
SHS processes behave nicely, while Assumptions 3 and 4 are
introduced to facilitate the derivation of the PDE.
Denote by ξt = ξ(t) = (q(t), x(t)) the hybrid state of the
SHS. Let fi be the ith entry of f , and σij be the (i, j)-entry
of σ. For a function ψ defined on the hybrid state space X,
we denote ψ q = ψ Xq the restriction of ψ to Xq .
Assumption 1. For the SDE in (11), we assume that fiq ∈ Cb1 ,
q
σij
∈ Cb1 ∩ C 2 , ∀q ∈ Q. In addition, the initial state ξ0 is
2
independent of {W (t), t ≥ 0} and satisfies E(kξ0 k ) < ∞.
Remark 1. This is a standard assumption to guarantee that
the initial value problem of (11) has a unique continuous
solution
(i.e., the Itô diffusion), and the solution satisfies
´t
2
E( 0 kx(s)k ds) < ∞, ∀t > 0 [30, Theorem 5.2.1]. The
smoothness of f and σ is imposed for obtaining a PDE model
in the classical sense.
Assumption 2. There is no Zeno execution for the hybrid
system of (11).
Remark 2. Let Nt (ω) := sup {k ∈ N : Tk (ω) ≤ t},
which is the number of jumps happened before time t. Then
Assumption 2 implies Eξ0 (Nt ) < ∞, ∀ξ0 ∈ X. This is a
common assumption for the SHS.
Regarding the structure of the hybrid state space X of the
SHS process, we have,
Assumption 3. ∀q ∈ Q, X̄q is a connected and oriented C 2
manifold with corners.
The above assumption requires that the state space is locally
C 2 -diffeomorphic to [0, ∞)n . It unifies the description of the
state spaces of both the TCLs (e.g. HVACs) and defferable
loads (e.g. PEV) and facilitates our proof later. More details
on manifold with corners can be found in [31].
Recall that Gq is the outflow switching surface in mode q
and G := ∪q {q} × Gq . For the SHS process, we can define Gq
explicitly as
n
Gq := x ∈ ∂Xq : f (q, x) · ν(q, x) > 0
o
or ν T (q, x)Σ(q, x) 6= 0 , (17)
where ν(q, x) is the outer unit normal vector on ∂X. This
extends the definition in (2) for the deterministic system. It is
possible that the outflow switching surface Gq may not contain
all of ∂Xq , i.e., not the entire boundary can be reached from

some interior point of Xq . For example, in the deterministic
case (i.e., σ ≡ 0), Gq is the forward reachable boundary
defined in (2). Moreover, we assume that:
Assumption 4. The outflow switching surface Gq is an open
subset of ∂Xq .
Remark 3. Assumptions 3 and 4 imply that Gq is also oriented [31, Chapter 15]. This is to guarantee that the integrals
over manifolds are defined consistently. In particular, this is
a standing assumption of the divergence theorem that will be
used in the proof of our Theorem 1.
The following defines a metric on the hybrid state space
X [14, page 58]. A distance function between the hybrid states
0
0
can be defined by ρ(ξ, ξ ) = 1 if q 6= q 0 and ρ(ξ, ξ ) = χ(x −
2
−1
0
0
x ) if q = q , where χ(x) = π tan (kxk), ∀x ∈ Rn . Now
X can be endowed with the Borel σ−algebra B(X) generated
by its metric topology, where B(X) = σ {∪q {q} × Bq } and
Bq is the σ−algebra on Xq .
For clarity, we introduce the following assumptions on the
switching surfaces and define the partitions of the hybrid state
space. Since the SHS process is defined on X, we assume
that the post-jump position (φ(q, x), x) ∈ X, for a jump from
(q, x) ∈ X̄. Recall that S := ∪q {q + } × Gq is the union of
all the inflow switching surfaces. Therefore, we have S ⊂ X.
In addition, it is assumed that S ∩ G = ∅. However, note
that S ∩ ∂X may not be empty. Moreover, define the surface
E := ∪q {q + } × ∂Xq . Clearly, by Assumption 4, S is an open
subset of E. We assume that E forms a partition of X̄, that
is, X̄ = ∪q∈Q,i∈Nq {q} × X̄qi , such that the each partition
{q} × Xqi satisfies {q} × (X̄qi ∩ X̄qj ) ⊂ E, ∀q ∈ Q, i, j ∈ Nq ,
i 6= j, where Nq is the index set of the partitions of {q} × Xq .
C. Dynkin’s Formula for the SHS Processes
As discussed before, the PDE is completely determined
by the strong generator (L, D(L)). However, generally it can
be very difficult to directly characterize D(L). Note that in
order to apply the idea of Lemma 2, we only need the
Dynkin’s formula hold, and it might be easier to characterize
the generator boundary condition by characterizing a super set
of D(L). To this end, we introduce the notion of the extended
generator (c.f. [14, Definition 14.15]).
Definition 2. (Extended Generator) Let D(L̂) denote the set
of Borel measurable functions ψ : X 7→ R with the following
property: there exists a measurable function h : X 7→ R, such
that the function t 7→ h(ξt ) is integrable Pξ0 a.s. for each
ξ0 ∈ X and the process
ˆ t
Ctψ = ψ(ξt ) − ψ(ξ0 ) −
h(ξs )ds,
0

is a local martingale. Then we write h = L̂ψ and call
(L̂, D(L̂)) the extended generator of the process ξt .
Since a martingale is also a local martingale, we see that
D(L) ⊂ D(L̂), and L̂ψ = Lψ for ψ ∈ D(L), and therefore the
name “extended generator”. More importantly, it can be easily
verified that set Dm := {ψ ∈ D(L̂) : Ctψ is a martingale} is
precisely the largest class of functions for which the Dynkin’s

8

formula (14) holds. Clearly, we have D(L) ⊂ Dm ⊂ D(L̂).
Since our purpose is to derive the PDE using the Dynkin’s
formula, it is our major focus to characterize Dm now.
Before proceeding, we introduce the following notations [14]. For a measurable function
g : X × R+ 7→ R,
´
we say g(ξ, s) ∈ L1 (τ ) if E( |g| dτ ) < ∞, where τ is a
stochastic process that can be identified by a random measure
and thus defines a Stieljes integral in the expression. Typically,
we will denote τ as the counting process defined by
X
τ (t, A) :=
1{Tk ≤t} 1{ξTk ∈A} ,
(18)
which counts the number of jumps of the GSHS process ξt .
Then we say g ∈ Lloc
1 (τ ) if there is a sequence of stopping
times γn with γn ↑ ∞ a.s. such that g1{s≤γn } ∈ L1 (τ ) for
∀n ∈ N+ . For ψ : X 7→ R, let Bψ : X × R+ × Ω 7→ R be
defined as Bψ(ξ, s, ω) := ψ(ξ)−ψ(ξs− (ω)). The function Bψ
is introduced to evaluate the increments of ψ due to the jumps.
Finally, we assume as a convention that for a continuous
function on X, its value on the boundary G is defined by
its continuous extension.
We first cite a characterization result of the extended generator of a class of general SHS (GSHS) processes proposed
in [15].
Lemma 3. (GSHS extended generator [15, Theorem 4.11])
The extended generator (L̂, D(L̂)) of a GSHS process ξt
satisfies ∀ψ ∈ D(L̂),
1
(L̂ψ)(ξ) = ∇ψ(ξ) · f (ξ) + T r(Σ(ξ)∇2 ψ(ξ))
2
ˆ
+ λ(ξ)
(ψ(ζ) − ψ(ξ)) R(ξ, dζ), (19)
X

with D(L̂) including at least those functions ψ : X 7→ R such
that:
1. ψ q ∈ C 2 , ∀q ∈ Q.
´
2. (boundary condition) ψ(ξ) = X ψ(ζ)R(ξ, dζ), ξ ∈ G.
3. Bψ ∈ Lloc
1 (τ ).
where R : X̄ × B(X) 7→ [0, 1] is a transition measure on the
post-jump positions given the pre-jump position.
Remark 4. An exact characterization of the extended generator
has been obtained for the PDMP [14, Theorem (26.14)]. It
hinges on the local martingale representation theorem for the
general jump process. However, such a theorem is unavailable
for the GSHS process, and only the sufficiency of these
conditions can be claimed. In Lemma 3, the key information
of the boundary condition is obtained through characterizing
the extended generator. It is worth mentioning that this is the
same boundary condition satisfied by the extended generator
of the PDMP.
Next, we will characterize a subset of D(L̂) on which Ctψ
is actually a martingale (rather than a local martingale) with
respect to the natural filtration Ft generated by the GSHS
process.
Lemma 4. (Martingale characterization)Let (L̂, D(L̂)) be the
extended generator of a GSHS process ξt and ψ ∈ D(L̂). If
ψ q is bounded ψ q ∈ Cb1 for all q ∈ Q, then Ctψ is an Ft martingale.

Proof: See Appendix A.
Based on Lemmas 3 and 4, we are now in the position to
establish the Dynkin’s formula for the proposed SHS process
in (11). Note that the post-jump position of the proposed
SHS process is specified by the mode transition function φ
deterministically, that is, there is a one-to-one correspondence
between the pre-jump position and post-jump position. This
will give us a local boundary condition for the extended
generator.
Proposition 1. (SHS extended generator) Under the Assumptions 1 and 2, the extended generator (L̂, D(L̂)) of the SHS
process (11) satisfies ∀ψ ∈ D(L̂),
1
(L̂ψ)(q, x) = ∇ψ(q, x) · f (q, x) + T r(Σ(q, x)∇2 ψ(q, x))
2

+ λ(q, x) ψ(q + , x) − ψ(q, x) , (20)
with D(L̂) including at least those functions ψ : X 7→ R
such that: (1) ψ q ∈ C 2 , ∀q ∈ Q, (2) (boundary condition)
ψ(q + , x) = ψ(q, x), ∀(q, x) ∈ G, and (3) Bψ ∈ Lloc
1 (τ ).
Proof: See Appendix B.
For all ψ ∈ D(L̂) characterized in Proposition 1, Ctψ is a
local martingale on the SHS process (11). We define,
n
D = ψ : X 7→ R ψ is bounded, ψ q ∈ Cb1 ∩ C 2 ,
o
and satisfies ψ(q + , x) = ψ(q, x), ∀(q, x) ∈ G .
Then combining Lemma 4 and Proposition 1, we can obtain:
Proposition 2. (SHS Dynkin’s formula) For all ψ ∈ D, Ctψ is
a martingale and the Dynkin’s formula (14) holds.
Proof: See Appendix C.
Note that the characterized set D is a subset of Dm .
However, it contains enough functions for deriving the PDE
model.
D. PDE Model
In this subsection, we derive the PDE model based on the
Dynkin’s formula (14) proved in Proposition 2.
In view of the partition of the state space X, we will start by
assuming that the p(q, x, t) is piecewise C 2,1 in the interior
of Xq × R+ , ∀q ∈ Q, where C 2,1 means that p(q, x, t) is
C 2 with respect to x and C 1 with respect to t. Of course,
certain smoothness conditions on the initial distribution and
coefficients are required to make this happen, otherwise p is
only a weak solution in the sense of (16) (i.e., in the sense
of weak derivatives, see [27]). The detailed discussion will be
lengthy and out of the scope of this paper. Hence, we will
restrict ourselves to the smooth assumption.
Let us denote p(q, x, t) = pi (q, x, t) in Xqi , ∀q ∈ Q, i ∈ Nq .
Theorem 1. (PDE characterization) Under the Assumptions 1-4, for all q ∈ Q, i ∈ Nq , x ∈ Xqi , the hybrid-state
probability density function pi (q, x, t) satisfies the following
PDE:
T [pi ] = 0,
(21)

9

where
∂p(q, x, t)
T [p] =
+ ∇ · γ(q, x, t)
∂t
+ λ(q, x)p(q, x, t) − λ(q − , x)p(q − , x, t),
and γ is known as the probability flux given as
1
γ(q, x, t) = f (q, x)p(q, x, t) − ∇ · (p(q, x, t)Σ(q, x)).
2
Proof: See Appendix D.
We see that the above PDEs are coupled in the state spaces
Xqi due to the source terms λp caused by the random jumps.
Noting that without the random jump, i.e. λ ≡ 0, equation
(21) reduces to the well-known Fokker-Planck equations.
Moreover, they are also coupled on the boundaries G and
S due to the deterministic jumps. Let ν(q, x) denote the
−
outer unit normal vector on Gq , and define →
p (q + , x, t) :=
→
−
+
+
limε↑0 p(q , x + εν, t) and γ (q , x, t) := limε↑0 γ(q + , x +
εν, t). These are respectively the continuous extensions of the
probability flux p and density γ on S from the side of S that
has the same outer normal direction as G. Similarly, define
←
−(q + , x, t) := lim γ(q + , x + εν, t) and ←
−
γ
p (q + , x, t) :=
ε↓0
+
limε↓0 p(q , x + εν, t) the quantities continuously extended
from the opposite side of S. We have:
Theorem 2. (PDE boundary conditions) The PDE model in
(21) satisfies the following boundary conditions, for all q ∈ Q,
x ∈ Gq ,
p(q, x, t)ν T (q, x)Σ(q, x) = 0,
(22)
 T
→
−
←
−
+
+
+
p (q , x, t) − p (q , x, t) ν (q, x)Σ(q , x) = 0, (23)
→
−
−(q + , x, t) + γ(q, x, t) · ν(q, x) = 0. (24)
γ (q + , x, t) − ←
γ
Proof: See Appendix E.
A closer look at (22)-(24) will give us insight into the patterns of the boundary conditions. Physically, ν T (q, x)Σ(q, x)
represents the projection of the multi-dimensional Brownian
motion along the normal direction of the outflow switching
surface. In particular, if ν T (q, x)Σ(q, x) 6= 0, then
1) Condition (22) implies that the p.d.f. at any point of the
outflow switching surfaces must be 0;
2) Condition (23) implies that the p.d.f. is continuous
across the inflow switching surfaces;
3) Condition (24) implies that the sum of the probability
fluxes along the normal direction of the inflow switching
surface is zero.
Therefore, following the convention of [8], we will call (22)
the absorbing condition, (23) the continuity condition, and (24)
the probability conservation condition. In particular, note that
for the deterministic continuous systems where there is no
Brownian motion, (i.e., Σ = 0) the results in Theorems 1
and 2 still hold.
Remark 5. It is worth mentioning that the author in [21] derived a measure-valued formulation of the forward equation for
the GSHS based on the Levy’s identity [22], [23]. In contrast,
our method is based on a sufficient characterization of the SHS
martingales (or the Dynkin’s formula). This martingale characterization is equivalent to the Levy’s identity in evaluating
the expected increments caused by the jumps. Nevertheless,

it extracts more information from the individual load model,
namely the boundary condition satisfied by the generator of the
SHS. The generator boundary condition directly determines
the PDE boundary conditions through the adjoint relation.
However, this information was not employed in [21], and
therefore the results there do not directly apply to the PDE
modeling of responsive load aggregation. In addition, it is
worth mentioning that the presence of the generator boundary
condition eliminates the need of proving the existence of the
mean jump intensity as required by [21].
E. Load Population with Heterogeneous Parameters
In this subsection, we consider the case of heterogeneous
parameters. As shown in (7), theoretically we can treat the parameter distribution in the same way as the initial distribution.
However, the increase in dimensions will cause major computational issues. Some of the existing methods rely on the noise
process to account for the parameter heterogeneity [10], [11].
Nonetheless, this is generally not an accurate characterization,
since the randomly distributed load parameters θ may not be
equivalent to an additive drift term to the vector field f that is
normally distributed. Therefore, the corresponding collective
behaviors can be quite different.
Since parameter heterogeneity cannot be captured well by
diffusion terms, we approximate the parameter heterogeneity
using several homogeneous population whose parameters are
obtained from clustering over the samples drawn from the
distribution pθ . Denote p(q, x, t) the p.d.f. of the heterogeneous population dynamics (7). Then the conditional p.d.f.
p(q, x, t|θi ) represents the p.d.f. of the homogeneous population with parameters θi , that is, the SHS model (7) with
fixed load parameters θi , assuming that it has the same initial
p.d.f p0 as the heterogeneous population. Intuitively, p(q, x, t)
can be estimated by p.d.f. of sufficiently many homogeneous
populations. Specifically, we have
Proposition 3. Let θi , i = 1, 2, . . . n, be the samples drawn
from the distribution pθ : Θ 7→ R+ , and assume that the set
Θ are compact. If p(q, x, t|θ) is continuously dependent on θ,
then ∀q, x, and t > 0,
n

1X
p(q, x, t|θi ) = p(q, x, t), almost surely
n→∞ n
i=1
lim

(25)

Proof: See Appendix F.
Remark 6. The continuous dependence of the p.d.f. on the
initial data may be most appropriately investigated from the
corresponding PDE. However, it is out of the scope of this
paper to go into that direction to determine the conditions
that guarantee the well-posedness of the PDE. Instead, for the
practical problem considered in this paper, we shall assume
this property carried over from the SHS modeling. Later on,
we will demonstrate the effectiveness of our method through
simulation on a realistic example.
We can further P
approximate the empirical distribution funcnc
wk 1{θ̄k ≤θ} using much fewer samples
tion Fn by F̃c := k=1
Pnc
k
θ̄ , where nc  n and k=1
wk = 1 such that fixing nc ,
the samples θ̄k and the index sets Ik minimize the sum of

10

Pnc P
k
i
the within-cluster distances (n) := k=1
i∈Ik θ̄ − θ ,
i
where θ , i ∈ {1, 2, . . . n} are the original samples. Hence,
Ik contains the indices of θi which are closest to θ̄k , and
wk = nnk , nk = |Ik | is the weight of the kth cluster. In
this paper, we adopt the k−means clustering algorithm which
efficiently calculates a suboptimal solution to the minimization
of (n). Furthermore, by assuming the continuous dependence
on the parameters, the approximation error can be bounded by
ˆ
ˆ
p(q, x, t|θ)dFn − p(q, x, t|θ)dF̃c
nc X
1X
p(q, x, t|ζ̄ k ) − p(q, x, t|ζ i )
n
k=1 i∈Ik
c
≤ (n),
n
where c > 0 is a constant depending on the set Θ, and
the specific SDEs. Since the approximation error decreases
linearly with (n), we may increase the number of the clusters
nc to reduce the error. However, it is a trade-off between the
accuracy and the computational complexity.

≤

VI. A PPLICATIONS
The proposed SHS approach provides a unified framework
to obtain aggregate models of responsive loads. It contains
many existing models as special cases. For example, the
pioneer work by [8] can be easily obtained from Theorems 1
and 2 (See the following example in Section VI-A). In this
section, we will discuss three aggregate modeling examples
and use the proposed SHS framework to derive the associated density evolution PDEs. These examples have not been
formally studied in the literature. Hence, this section not only
demonstrates how to use our SHS framework, but also contains
new contributions in the field of aggregate load modeling.
A. PDE Model for HVACs
Consider the HVAC aggregation problem similar to Example 1. To draw connections with other results in the literature,
we assume a homogeneous population where each HVAC is
modeled by a 2D SDE,


dW1 (t)
dx(t) = Ax(t) + Bq(t) + σ
,
dW2 (t)
and the same local control rule (4), where Wi , i = 1, 2
are the standard 1D Brownian motions. For the convenience
of comparison, we do not consider the random jump in
this example. The hybrid state space X, outflow switching
surface G, and inflow switching surface S are the same as
in Example 1, which is illustrated in Fig. 2. We denote the
partitions
of X by Xqi = {q} × Γi (q), where Γ2 (q) :=

2
x ∈ R : uset − δ < x1 < uset + δ is the overlapped continuous state space, and Γ1 (q) represents the rest of the state
space given by Xq \Γ2 (q). Let pi (q, x, t) be the corresponding
p.d.f. defined on Γi (q). By Theorem 1, we can easily have for
all q = 0, 1, i = 1, 2, and x ∈ Γi (q),
Dt pi (q, x, t)+∇·(f (q, x)pi (q, x, t))−

σ2
∆pi (q, x, t) = 0,
2

where f (q, x) = Ax + Bq and ∆ is the Laplacian operator.
The PDE boundary conditions directly follow from (22)(24). We have the outer unit normal vector ν(0, x) = [1, 0]T
on the outflow switching surface G0 . Similarly, for the ON
mode ν(1, x) = [−1, 0]T on G1 . Therefore, the condition (22)
yields p2 (q, x, t) = 0, ∀q = 0, 1, x ∈ ∂Xq . Condition (23)
reduces to p1 (q, x, t) = p2 (q, x, t) on the inflow switching
surface S. Note that by the definition used in Theorem 2,
−
the probability flux →
γ (q + , x, t) is defined in Γ2 (q + ), while
γ(q, x, t) is the outflow probability flux from Γ2 (q), for
q = 0, 1. Furthermore, using the boundary conditions obtained
above and the continuity of the vector field f (q, x, t), we have
from (24) that for all q = 0, 1 and x ∈ Gq , the following
probability conservation condition holds:
Dx1 p2 (q + , x, t) − Dx1 p1 (q + , x, t) + Dx1 p2 (q, x, t) = 0.
The result in [8] is a simple special case of the above
obtained model (one immediately get the PDE model of [8] if
the x2 dynamic is 0). Many other PDE models for aggregating
HVACs [9], [10], [11], [12], [32], [33] can also be justified or
obtained conveniently from Theorems 1 and 2.

B. PDE Model for PEV Charging
We consider the aggregation of PEVs with the charging
dynamics as detailed in Example 2. Note that the heterogeneous parameters such as battery capacity, charging rate,
charging deadlines may all be normalized or translated into
the initial distribution of the timing dynamics. In addition, the
spontaneous mode switching behaviors will be modeled for
the realistic consideration. The PDE model can capture the
distribution evolution of the PEVs in different modes, which
will be useful to predict the aggregate power consumption.
By using Theorems 1 and 2, the major work to obtain the
PDE model is only left to identify the partitions of the state
spaces. We see that only X1 is partitioned into Γ11 := X1 \X0
and Γ12 := X0 . Let p1 (1, x, t) and p2 (1, x, t) be the p.d.f.
on Γ11 and Γ12 , respectively. Suppose the transition intensity
function λ(q, x) is given. Then we obtain directly from (21):
Dt p(0, x, t) − Dx2 p(0, x, t) + λ(0, x)p(0, x, t) = 0 in X0 ,
Dt p1 (1, x, t) − Dx1 p1 (1, x, t) + λ(1, x)p1 (1, x, t) = 0 in Γ11 ,
Dt p2 (1, x, t) − Dx1 p2 (1, x, t) + λ(1, x)p2 (1, x, t)
−λ(0, x)p(0, x, t) = 0 in Γ12 ,
Dt p(2, x, t) − Dx1 p2 (1, x, t) = 0 in X2 .
Since Σ is a zero matrix, there are only probability conservation conditions. On the outflow switching surface G0 ,
we have ν(0, x) = [0, −1]T . Let S1 = {x1 = 0, x2 < 0}
and S2 = {x1 = 0, x2 ≥ 0} be the switching surfaces
of Γ12 and Γ11 , respectively, on both of which we have
ν(1, x) = [−1, 0]T . Then we can obtain from (24) the
following PDE boundary conditions:
p(0, x, t) = 0, ∀x ∈ G0 ,
p(2, x, t) = ph (1, x, t), ∀x ∈ Sh , h = 1, 2.

11

C. Price Response of Aggregate HVACs

10
Monte Carlo
Weighted PDE
8
Total Power (MW)

In this example, we consider the price response of HVACs
with heterogeneous parameters θ = (A, B1 , B2 , uset , δ),
see (3) in Example 1 for the 2D hybrid system model of
the individual HVAC load. The price responsive HVACs can
adjust their power state according to a price signal. In addition
to the local control rules (4), the setpoint of the HVAC can
be modified by an external price signal υ(t) according to a
adjustable price response curve [34]. We will derive the PDE
model for a heterogeneous population of HVACs, and use it
to characterize the aggregate power response under the price
variation. Note that different from Section VI-A, we will not
consider the noise process in this example.
As discussed in Section V-E, we can cluster the heterogeneous parameters θ to obtain several homogeneous populations. Then the p.d.f. of the heterogeneous HVACs can be
approximated by the weighted sum of the solutions of the
PDE models with homogeneous parameters. Specifically, the
coefficients of each PDEs are evaluated at the cluster centers
of the sampled parameters. For example, the p.d.f. evolution
of the kth cluster is governed by:

Dt pki (q, x, t) + ∇ · f k (q, x)pki (q, x, t) = 0,
(26)

6

4

2

0
2

3

−f

k

(q, x)pk2 (q, x, t)

= 0,

for x = uset + δ if q = 0, and x = uset − δ if q = 1.
The above PDE and boundary conditions can be easily
obtained from the model derived in Section VI-A by setting
σ = 0. Once pk (q, x, t) is obtained for each cluster, the
aggregate power response can be calculated as
ˆ
nc
X
wk
y a (t) = N W
pk (1, x, t)dx,
k=1

X1

where N is the total number of the HVACs , W is the power
consumption of one HVAC, and wk is the weight of the kth
cluster.
Without loss of generality, we assume that the price changes
at t = 3 hour and t = 6 hour, which cause the setpoints
changes of all HVACs by 1◦ F and −1◦ F, respectively. For the
Monte Carlo simulation, the load parameters are generated
randomly using GridLAB-D [35], and then the ETP models
are derived using these parameters. The initial setpoints for the
2000 HVACs are generated uniformly between [70◦ F, 78◦ F],
and the initial temperatures for air and internal solid mass are
generated uniformly within the temperature deadband. For the
PDE simulation, the k-means clustering method is performed
over the system matrices (which are generated by the load
parameter and note that f is linear in x and q) to obtain 10
clusters. The initial distribution pξ0 is given arbitrarily since
the aggregate power output converges to the steady state very
fast before the setpoint changes. Then the PDEs in (26) are

5

6
7
Time (Hour)

8

9

10

Fig. 5. Comparison of aggregate response between Monte Carlo simulation
and weighted PDE using 10 clusters.

solved for each cluster using the Donor-Cell finite volume
scheme along with the dimensional splitting method for multidimensional problem [36]. We compare in Fig. 5 the aggregate
response of Monte Carlo simulation of 2000 HVACs and that
of the weighted PDEs of 10 clusters. Clearly, the weighted
PDE model captures the major oscillation cycles after the
price/setpoint changes very accurately.

in Γi (q), ∀q = 0, 1, i = 1, 2, (see Subsection VI-A for
the definition of Γi (q)), where f k (q, x) is the vector field at
the kth center, k = 1, 2, ..., nc . The boundary conditions are
obtained from Theorem 2 as
f k (q + , x)pk1 (q + , x, t) − f k (q + , x)pk2 (q + , x, t)

4

VII. C ONCLUSION
This paper developed a unified stochastic hybrid system
(SHS) framework for the aggregation of a large population
of responsive loads. General nonlinear stochastic differential
equations were used to describe the continuous state evolution
between discrete mode transitions, and the mode transition was
modeled by both deterministic and random jumps. The PDE
satisfied by the hybrid-state probability density function was
derived based on the adjoint relationship with the generator
of the SHS process. In particular, the PDE boundary conditions were uniquely determined from the generator boundary
condition. These results generalized many existing models and
can be directly applied to more general modeling scenarios.
Several examples were provided to illustrate the effectiveness
of the proposed modeling framework.
A PPENDIX
A. Proof of Lemma 4 (Martingale Characterization)
Following the same steps in [19, Proof of Theorem 2], we
can represent Ctψ as
ˆ t
ˆ
∇ψ(ξt ) · σ(ξt )dW (t) +
Ctψ =
Bψdq,
(27)
0

[0,t]×X

where q(t, A), ∀A ∈ B(X) is the stochastic process defined in
[14, page 67, Proposition 26.7]. We claim that both terms on
the right side of (27) are martingales with respect to the natural
filtration Ft generated by the GSHS process ξt . Note that all
the positive constants will be denoted by c in the following.

12

To show that the first term is a Ft -martingale, we will prove
that all the conditions of [30, page 33, Corollary 3.2.6] hold.
q
We know σij
∈ Cb1 by Assumption 1. Therefore, we can have
|σij (x)| ≤ c(1 + kxk),
for some constant c > 0. Since ψ q ∈ Cb1 , we have ∇ψ is
bounded. Note that the first term can be written as
ˆ tX
m X
n
(∇ψ)i σij dWj (t),
0 j=1 i=1

where (∇ψ)i and Wi (t) are the ith entry of ∇ψ and W (t),
respectively. Then for a fixed index j, ∀t ≥ 0, we have
"ˆ

n
X
(∇ψ)i σij

t

E
0

!2

#

t

ˆ

dt ≤ cE


(1 + kxk)2 dt < ∞,

0

i=1

(28)

where the last inequality is by Remark 1. Clearly, ∇ψ(ξt )σ(ξt )
is Ft -adapted. Together with the fact that ξt is càdlàg [19], it
can be concluded that ∇ψσ as a function of (t, ω) ∈ [0, ∞)×Ω
is B × F measurable [37, page 5, Proposition 1.13], where
B is the Borel σ-algebra on [0, ∞).
´ t Combining the above
results, it can be concluded that 0 ∇ψσdW (t) is an Ft martingale by the standard results of stochastic integration,
see for example [30, Corollary 3.2.6].
For the second term, the proof is similar to [14, page 272,
Theorem (A4.9)] for the single jump process and combines the
idea of [38, Lemma 3.1] to decompose a general jump process
into the summation of single jump processes. However, the
natural filtration and the probability measure P (ω) of the
GSHS are different from those of the general jump process
due to the additional Itô diffusions between jumps. To prove
the claim, we can express P (ω) explicitly as follows, for
t ∈ [Tk , Tk+1 ),
ˆ
ψ(ξt (ω))dP (ω) =
ˆΩ ˆ
ˆ
k
···
ψ(ξt (ω))dPξ(T
(ω k )µk (dYk ; ω k−1 ) · · ·
k)
Ω0

Ω1

Ωk

1
× dPξ(T
(ω 1 )µ1 (dY1 ; ω 0 )dPξ00 (ω 0 )µ0 (dξ0 ), (29)
1)

where ω is the concatenation of the trajectories of the Itô
diffusions ω k ∈ Ωk killed at time Tk+1 − Tk with ω k (0) =
k
ω(Tk ) = ξTk , ∀k ∈ N. The measure Pξ(T
is the Wiener
k)
probability corresponding to the kth Itô diffusion starting at
Yk := (ξTk , Tk ). The measure of the post-jump position and
the jump time of the kth jump conditioned on the previous Itô
diffusion ω k−1 is defined by ∀A ∈ B(X), ∀s > Tk−1 ,
µk (A×(s−Tk−1 , ∞); ω k−1 ) := R(ξs− , A)β(s−Tk−1 ; ω k−1 ),
where R is the transition measure defined in Lemma 3 and
β is the survival function defined in (9). Thus using (29),
several intermediate results in proving [14, Theorem (A4.9)]
and [38, Lemma 3.1] can be easily verified by direct calculations. Hence, in the following we will leave out the explicit
calculations, and instead, focus on the key arguments of the
proof.
Following [38, Lemma 3.1], we can decompose q(t, A)
into the sum of single jump processes conditioned on

their immediate precedent Itô diffusions, that is, q(t, A) =
P
k−1
(t, A; ω k−1 ). Then by direct calculation using
k≥1 q
(29) in a similar way to [14, Theorem (A4.4)], it can be
shown that each q k−1 is an Ft -martingales (which is often
called martingale-valued measures or simply martingale measures [39, page´ 105]). Since Bψ is bounded, it can be proved
similarly that [0,t]×X Bψdq k−1 is also an Ft -martingales. It
follows immediately that q(t∧Tk , A) is a uniformly
integrable
´
Ft -martingale for each k ∈ Z+ , and so is [0,t∧Tk ]×X Bψdq.
Now the´result in [40, Proposition 1.7] implies that in order
to show [0,t]×X Bψdq is a Ft -martingale, it suffices to prove
that Bψ1{s≤t} ∈ L1 (q) for all t ≥ 0. Through a similar
calculation to [14, Theorem (A4.5)], this condition can be
shown to be equivalent to Bψ1{s≤t} ∈ L1 (τ ) for all t ≥ 0,
where τ is a counting process defined earlier in (18). The latter
condition can be further evaluated as ∀t ≥ 0,

ˆ
Bψ1{s≤t} dτ =
Eξ0
!
Nt
X
Eξ0
ψ(ξTk ) − ψ(ξT − ) , (30)
k

k=1

which is finite since ψ is bounded and Eξ0 (Nt ) < ∞ by
Remark 2. Hence, Bψ1{s≤t} ∈ L1 (τ ) and we conclude that
the second term on the right side of (27) is a Ft -martingale.
This completes the proof of this lemma.

B. Proof of Proposition 1 (Extended Generator)
The formula in (20) is obtained by taking the transition
measure R in Lemma 3 as a Dirac measure (unit mass),
∀(q, x) ∈ X̄,
(
1, if ζ = (q + , x);
(31)
δ(q+ ,x) (ζ; (q, x)) =
0, otherwise.
Similarly, by using (31), the boundary condition in Lemma 3
reduces to ψ(q, x) = ψ(q + , x), ∀(q, x) ∈ G.

C. Proof of Proposition 2 (Dynkin’s Formula)
This directly follows from Corollary 1 and Lemma 4. The
boundedness of ψ guarantees that Bψ ∈ Lloc
1 (τ ) as can be
seen from (30) in the proof of Lemma 4. Therefore, the set D
satisfies all the conditions for ψ in Corollary 1 and Lemma 4.
Hence, Ctψ is a martingale. The Dynkin’s formula follows by
taking the expectation of Ctψ and noticing that it is 0.

D. Proof of Theorem 1 (PDE Characterization)
The following divergence theorem simplifies our main
proofs of obtaining the forward equation. It does not require
the vector field g to be compactly supported on M . Although
the original theorem given in [41] is stated for the manifold
with boundary, the same result holds for the manifold with
corners and can be proved in a similar way.
Lemma 5. (Divergence Theorem) Let M ⊂ Rn be an oriented
C 2 -manifold with corners and M is closed, and ν : ∂M 7→

13

Rn be the unit outward pointing normal to M . If g : M →
7 Rn
1
◦
is continuous on M and C in M , and
ˆ
ˆ
{|g| + |∇ · g|} dV +
|g · ν| dS < ∞,
(32)
M

then

one of the partition Xqi . Then all the surface integrals in (37)
vanishes, and we are left with
ˆ
T [p]ψdV = 0.
Xqi

∂M

ˆ

ˆ
g · νdS =

∂M

Since ψ is arbitrary, we must have T [p] = 0. Thus it
completes the proof.

∇ · gdV,
M

where S is the surface measure on ∂M and V is the volume
measure.
Proof of Theorem 1: By Proposition 2, ∀ψ ∈ D, the
Dynkin’s formula (14) holds. Evaluate the expectation using
p(q, x, t), and take the time derivative of both sides of (14),
yielding
ˆ
h
i
∂X
ψ(q, x)p(q, x, t)dV = E L̂ψ(q(t), x(t)) , (33)
∂t q Xq
By the bounded convergence theorem, we can exchange the
time derivative and the integral, and it follows that
Xˆ 
∂p(q, x, t)
− ∇ψ(q, x) · f (q, x)p(q, x, t)
ψ(q, x)
∂t
Xq
q

1
− T r Σ(q, x)∇2 ψ(q, x) p(q, x, t)
2


+
− λ(q, x) ψ(q , x) − ψ(q, x) p(q, x, t) dV = 0, (34)
By the chain rule of differentiation, the second term in (34)
can be written as
∇ψ · f p = ∇ · (f pψ) − ψ∇ · (f p) ,

(35)

and the third term can be written as

1
T r Σ∇2 ψ p
2

n

=
=

1X
pΣj (∇ψ)xj
2 j=1
n
1X

2 j=1

(pΣj ∇ψ)xj − (pΣj )xj ∇ψ

1
1
∇ · (pΣ∇ψ) − ∇ · (∇ · (pΣ)ψ)
2
2
1
(36)
+ ∇ · (∇ · (pΣ)) ψ.
2
where we dropped the dependence on the time and state where
no confusion arises.
Suppose that all the vector fields under the divergence
operation in (35) and (36) satisfy the absolutely integrable
condition (32), and then we can use Lemma 5 in (34) to shift
the differential operator from ψ to p. Since p is piecewise
C 2 in the interior of each Xq , we can apply the divergence
theorem to each partition of X. It then follows after grouping
the terms that


ˆ
Xˆ
1
T [p]ψdV −
γψ + pΣ∇ψ · ν(i)dS = 0,
2
i
∂Xqi
q,i Xq
(37)
where T [p] and γ is defined in (21) and ν(i) is the outer unit
normal vector on ∂Xqi . Now let ψ be compactly supported in
=

E. Proof of Theorem 2 (Boundary Conditions)
Note that in general the surface integral in (37) will be
integrated over ∂X once, and over E twice from both sides,
where recall that ∂X = ∪q {q} × ∂Xq and E = ∪q {q + } ×
∂Xq . We first claim that it is equivalent to evaluate the surface
integral in (37) over G and both sides of S only. Since ψ q ∈
C 2 , ∀q ∈ Q, ψ and ∇ψ are continuous over E\S. Hence,
the surface integrals on both sides of E\S cancel each other
out. By the definition of Gq in (17), we have ν · f ≤ 0 and
ν T Σ = 0 on ∂X\G. If ν · f = 0 and ν T Σ = 0, then clearly
the integral is 0 over∂X\G. If ν ·f < 0 and ν T Σ = 0, we have
two situations: 1) on (∂X\G)\S, we must have p(q, x, t) = 0
for all t > 0 (since all particles left the boundary immediately
after the starting time and no particle jumps to the boundary
from other modes), and hence, the surface integral is 0 on
(∂X\G)\S; 2) on (∂X\G) ∩ S: the integration over this set is
in fact included in the integration over E. Thus it proves the
claim.
Now using the generator boundary condition ψ(q, x) =
ψ(q + , x), ∀(q, x) ∈ Gq in Proposition 1, we can collect the
γψ terms which are integrated on G and on both sides of
S. Furthermore, since ∇ψ is continuous, the 12 pΣ∇ψ terms
integrated on both sides of S can be collected. The last group
is the 21 pΣ∇ψ integrated on G. Moreover, noting that the outer
normal directions of Gq as in G and in S are the same, we
can obtain from (37) that
Xˆ
q,i

T [p]ψ(q, x)dV −

X ˆ

Xqi

q

ˆ

ˆ
H∇ψ(q + , x)dS +

+

F ψ(q + , x)dS

Gq

Gq

!
G∇ψ(q, x)dS

= 0, (38)

Gq

where F , H, G denote the left sides of (24), (23) and (22),
respectively.
From Theorem 1 we know that the volume integral in (38) is
0. Since ψ(q + , x), ∇ψ(q, x), and ∇ψ(q + , x) can be arbitrary
continuous functions on Gq , then F , G, and H must vanish
on Gq , which are the boundary conditions in (22)-(24).

F. Proof of Proposition 3 (Clustering)
The density functions for the heterogeneous and the homogeneous populations can be related by
ˆ
p(q, x, t) =
p(q, x, t|θ)pθ dθ.
Θ

Let Fθ be the cumulative distribution
Pn function of θ such
that dFθ = pθ dθ. Let Fn (r) := n1 i=1 1{θi ≤r} be the corresponding empirical distribution function with independent

14

and identically distributed (θi ) ∼ pθ , i ≤ n. Then by the
Glivenko-Cantelli theorem [42, Theorem 2.4.7], as n goes to
infinity,
sup |Fn (r) − Fθ (r)| → 0,
r∈Θ

almost surely.
Since p(q, x, t|θ) is continuously dependent on θ, it is
bounded on the compact sets Θ. Then it follows from the
´dominated convergence theorem that, as n goes to infinity,
p(q, x, t|θ)dFn → p(q, x, t) a.s., which is exactly (25). 
R EFERENCES
[1] D. S. Callaway, “Tapping the energy storage potential in electric loads to
deliver load following and regulation with application to wind energy,”
Energy Conversion and Management, vol. 50, no. 5, pp. 1389 –1400,
May 2009.
[2] S. Koch, J. L. Mathieu, and D. S. Callaway, “Modeling and control of
aggregated heterogeneous thermostatically controlled loads for ancillary
services,” in 17th Power Systems Computation Conference, Stockholm,
Sweden, August 2011.
[3] J. L. Mathieu, S. Koch, and D. Callaway, “State estimation and control of
electric loads to manage real-time energy imbalance,” IEEE Transactions
on Power Systems, vol. 28, no. 1, pp. 430–440, Feb. 2013.
[4] W. Zhang, J. Lian, C.-Y. Chang, and K. Kalsi, “Aggregated modeling
and control of air conditioning loads for demand response,” IEEE
Transactions on Power Systems, vol. 28, no. 4, pp. 4655–4664, Nov.
2013.
[5] J. L. Mathieu, M. Kamgarpour, J. Lygeros, G. Andersson, and D. S.
Callaway, “Arbitraging intraday wholesale energy market prices with
aggregations of thermostatic loads,” IEEE Transactions on Power Systems, vol. 30, no. 2, pp. 763–772, March 2015.
[6] S. Li, W. Zhang, J. Lian, and K. Kalsi, “Market-based coordination of
thermostatically controlled loads-part i: A mechanism design formulation,” IEEE Transactions on Power Systems, vol. 31, no. 2, pp. 1170–
1178, 2016.
[7] ——, “Market-based coordination of thermostatically controlled loadspart ii: Unknown parameters and case studies,” IEEE Transactions on
Power Systems, vol. 31, no. 2, pp. 1179 – 1187, 2016.
[8] R. Malhamé and C.-Y. Chong, “Electric load model synthesis by
diffusion approximation of a high-order hybrid-state stochastic system,”
IEEE Transactions on Automatic Control, vol. 30, no. 9, pp. 854–860,
Sept. 1985.
[9] S. Bashash and H. K. Fathy, “Modeling and control of aggregate air
conditioning loads for robust renewable power management,” IEEE
Transactions on Control Systems Technology, vol. 21, no. 4, pp. 1318–
1327, July 2013.
[10] S. Moura, V. Ruiz, and J. Bendtsen, “Modeling heterogeneous populations of thermostatically controlled loads using diffusion-advection
PDEs,” in ASME Dynamic Systems and Control Conference, Stanford,
CA, 2013.
[11] S. Moura, J. Bendtsen, and V. Ruiz, “Parameter identification of aggregated thermostatically controlled loads for smart grids using PDE
techniques,” International Journal of Control, vol. 87, no. 7, pp. 1373–
1386, 2014.
[12] A. Ghaffari, S. Moura, and M. Krstić, “Analytic modeling and integral
control of heterogeneous thermostatically controlled load populations,”
in ASME Dynamic Systems and Control Conference, San Antonio, Texas,
2014.
[13] M. Liu and Y. Shi, “Model predictive control of aggregated heterogeneous second-order thermostatically controlled loads for ancillary
services,” IEEE Transactions on Power Systems, vol. 31, no. 3, pp.
1963–1971, 2016.
[14] M. H. A. Davis, Markov Models and Optimization, ser. Monographs on
Statistics and Applied Probability. Chapman & Hall, 1993, vol. 49.
[15] L. M. Bujorianu, Stochastic Reachability Analysis of Hybrid Systems.
Springer, 2012.
[16] J. Hu, J. Lygeros, and S. Sastry, “Towards a theory of stochastic hybrid
systems,” in Lecture Notes in Computer Science, N. Lynch and B. Krogh,
Eds. Springer Berlin Heidelberg, 2000, vol. 1790, pp. 160–173.
[17] H. A. Blom, “Stochastic hybrid processes with hybrid jumps,” in
Analysis and Design of Hybrid System. IFAC Press, 2003, pp. 319–324.
[18] G. G. Yin and C. Zhu, Hybrid Switching Diffusions: Properties and
Applications. Springer-Verlag New York, 2010.

[19] M. L. Bujorianu and J. Lygeros, “Toward a general theory of stochastic
hybrid systems,” in Lecture Notes in Control and Information Science,
H. Blom and J. Lygeros, Eds. Springer Berlin Heidelberg, 2006, vol.
337, pp. 3–30.
[20] J. P. Hespanha, “A model for stochastic hybrid systems with application
to communication networks,” Nonlinear Analysis: Theory, Methods &
Applications, vol. 62, no. 8, pp. 1353–1383, Sept. 2005.
[21] J. Bect, “A unifying formulation of the Fokker-Planck-Kolmogorov
equation for general stochastic hybrid systems,” Nonlinear Analysis:
Hybrid Systems, vol. 4, no. 2, pp. 357–370, May 2010.
[22] J. B. Walsh and M. Weil, “Représentation de temps terminaux et
applications aux fonctionnelles additives et aux systèmes de lévy,”
Annales scientifiques de l’École Normale Supérieure, vol. 5, no. 1, pp.
121–155, 1972.
[23] R. F. Bass, “Adding and subtracting jumps from markov processes,”
Transactions of the American Mathematical Society, vol. 255, pp. 363–
376, Nov. 1979.
[24] N. W. Wilson, B. S. Wagner, and W. G. Colborne, “Equivalent thermal
parameters for an occupied gas-heated house,” ASHRAE Transactions,
vol. 91, no. CONF-850606-, 1985.
[25] K. Kalsi, F. Chassin, and D. Chassin, “Aggregated modeling of thermostatic loads in demand response: A systems and control perspective,” in
50th IEEE Conference on Decision and Control and European Control
Conference (CDC-ECC), Dec 2011, pp. 15–20.
[26] Z. Schuss, Theory and Applications of Stochastic Processes An Analytical Approach, ser. Applied Mathematical Sciences. Springer, 2010,
vol. 170.
[27] L. C. Evans, Partial Differential Equations, 2nd ed.
American
Mathematical Society, 2010.
[28] P. DuChateau and D. W. Zachmann, Schaums outline of theory and
problems of partial differential equations. MacGraw-Hill, 1986.
[29] F. B. Hanson, Applied Stochastic Processes and Control for Jump
Diffusions: Modeling, Analysis, and Computation. Philadelphia, PA:
SIAM Books, 2004.
[30] B. Øksendal, Stochastic Differential Equations, 6th ed. Springer, 2003.
[31] J. M. Lee, Introduction to Smooth Manifolds. Springer, 2013.
[32] D. Paccagnan, M. Kamgarpour, and J. Lygeros, “On the range of feasible
power trajectories for a population of thermostatically controlled loads,”
in 54th IEEE Conference on Decision and Control (CDC), Dec. 2015,
pp. 5883–5888.
[33] L. C. Totu, R. Wisniewski, and J. Leth, “Demand response of TCL population using switching-rate actuation,” IEEE Transactions on Control
Systems Technology, vol. 25, no. 5, pp. 1537–1551, Sept. 2017.
[34] S. Li, W. Zhang, J. Lian, and K. Kalsi, “On market-based coordination
of thermostatically controlled loads with user preference,” in 53rd IEEE
Conference on Decision and Control (CDC), Dec. 2014, pp. 2474–2480.
[35] GridLAB-D residential module user’s guild. [Online]. Available: http://sourceforge.net/apps/mediawiki/gridlab-d/index.php?title=
Residential module user%27s guide
[36] R. J. LeVeque, Finite Volume Methods for Hyperbolic Problems, ser.
Cambridge Texts in Applied Mathematics. Cambridge University Press,
2002.
[37] I. Karatzas and S. E. Shreve, Brownian Motion and Stochastic Calculus,
2nd ed. Springer, 1991.
[38] R. J. Elliott, “Stochastic integrals for martingales of a jump process with
partially accessible jump times,” Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete, vol. 36, no. 3, pp. 213–226, 1976.
[39] D. Applebaum, Lévy Processes and Stochastic Calculus. Cambridge
University Press, 2009.
[40] D. Revuz and M. Yor, Continuous Martingales and Brownian Motion,
3rd ed., ser. Comprehensive Studies in Mathematics. Springer, 1999,
vol. 293.
[41] B. K. Driver. (2002, Mar.) Surfaces, surface integrals and integration by
parts. Lecture Notes on Partial Differential Equations. [Online]. Available: www.math.ucsd.edu/∼bdriver/231-02-03/Lecture Notes/pde8.pdf
[42] R. Durrett, Probability Theory and Examples, 4th ed.
Cambridge
University Press, 2010.

15

Lin Zhao received the B.E. and M.S. degrees in automatic control from the Harbin Institute of Technology, Harbin, China, in 2010 and 2012, respectively.
He is currently a Ph.D. student in the Department
of Electrical and Computer Engineering, The Ohio
State University, Columbus, OH, USA. His current
research interests include modeling and control of
large-scale complex systems with applications in
power systems.

Wei Zhang received a B.S. degree in automatic control from the University of Science and Technology
of China, Hefei, China, in 2003, and a M.S. degree in
statistics and a Ph.D. degree in electrical engineering
from Purdue University, West Lafayette, IN, USA,
in 2009. Between 2010 and 2011, he was a Postdoctoral Researcher with the Department of Electrical
Engineering and Computer Sciences, University of
California, Berkeley, CA, USA. He is currently an
Associate Professor in the Department of Electrical
and Computer Engineering, Ohio State University,
Columbus, OH, USA. His research focuses on control and game theory with
applications in power systems, robotics, and intelligent transportations.

