1

Architecture and Algorithms for Privacy Preserving
Thermal Inertial Load Management by
A Load Serving Entity

arXiv:1606.09564v3 [cs.SY] 29 Nov 2016

Abhishek Halder, Member, IEEE, Xinbo Geng, Student Member, IEEE, P.R. Kumar, Fellow, IEEE,
and Le Xie, Senior Member, IEEE

Abstract—Motivated by the growing importance of demand
response in modern power system’s operations, we propose an
architecture and supporting algorithms for privacy preserving
thermal inertial load management as a service provided by the
load serving entity (LSE). We focus on an LSE managing a
population of its customers’ air conditioners, and propose a
contractual model where the LSE guarantees quality of service
to each customer in terms of keeping their indoor temperature
trajectories within respective bands around the desired individual
comfort temperatures. We show how the LSE can price the
contracts differentiated by the flexibility embodied by the width
of the specified bands. We address architectural questions of
(i) how the LSE can strategize its energy procurement based
on price and ambient temperature forecasts, (ii) how an LSE
can close the real time control loop at the aggregate level while
providing individual comfort guarantees to loads, without ever
measuring the states of an air conditioner for privacy reasons.
Control algorithms to enable our proposed architecture are given,
and their efficacy is demonstrated on real data.
Index Terms—Load serving entity, demand response, thermal
inertial load, privacy, contracts.

I. I NTRODUCTION

I

N traditional power systems, the role of feedback control
has primarily been concentrated on the generation side
where the demand side variability is dynamically matched
via operating reserve. As more renewables are being integrated into the power grid to reduce carbon footprint, a
new operational paradigm is emerging. Unlike fossil fuels,
generation from renewables such as solar and wind, cannot
be directly controlled. Thus, large scale renewable integration
requires moving from control of the supply to control of the
demand. This necessitates a collection of mechanisms, termed
as demand response, for manipulating aggregate demand to
offset the supply side variability. Thermal inertial loads such
as air conditioners (ACs) are one of the primary candidates
[1] for control on the demand side. They contribute to 48%
of U.S. annual residential consumption [2], and are potential
enablers for demand response due to their ability to store
energy and defer consumption without perceptible change in
end-user functionality.
The authors are with the Department of Electrical and Computer Engineering, Texas A&M University, College station, TX, 77843 USA, e-mail:
{ahalder,gengxbtamu,prk,le.xie}@tamu.edu. This material is
based upon work partially supported by NSF under Contract Nos. ECCS1546682, CPS-1239116, Science & Technology Center Grant CCF-0939370,
DGE-1303378, and ECCS-1150944.

In this paper, we consider an aggregator or load serving
entity (LSE) managing a finite population of residential ACs
through direct load control [1], [3]. We address three distinct
problems that are important for designing a direct load control
architecture, viz. (1) how to control a population of thermal
inertial loads to guarantee quality of service such as comfort
bounds for homes, (2) how to achieve such control while
preserving the privacy of the loads’ states from the LSE,
and (3) what to control the aggregate consumption to, so
that the cost of serving the collection of loads is minimized.
Our strategy effectively produces demand response from the
collection of homes according to how their temperatures are
distributed within their comfort zones at any given time,
without attempting to squeeze it from each and every home.
Holistically maintaining individual home privacy while still
generating demand response that respects home comfort constraints does not seem to have been addressed in the literature.
Previous works on thermal inertial load management [4]–
[11] have focused on obtaining an aggregate model for a
population of loads, and then used the resulting model, or an
approximation of it, to track a given reference power trajectory.
The purpose of this paper is to present a comprehensive
architecture and algorithms (Section III) so that the LSE can
manage a population of thermal inertial loads while addressing
the three questions mentioned in the preceding paragraph.
We address the first two problems by proposing a modelfree setpoint velocity control algorithm (Section III-1). We
formulate and solve (Section III-5) the third as an open-loop
optimal control problem. We also quantify (Section III-3) the
statistical performance of the LSE’s strategy under random
perturbations that entail different real-time energy usage than
planned in the day-ahead market. We show (Section IV) how
the LSE can price the contracts differentiated by the flexibility
provided by the customers with respect to their comfort
requirements, and illustrate (Section V) all our algorithms on
day-ahead price forecast data from the Electric Reliability
Council of Texas (ERCOT), and ambient temperature data
from a weather station in Houston, Texas. Some preliminary
results using these ideas applied to ERCOT real-time energy
market price forecast data, appeared in [12].
Notations: The notation x ∼ ρ means that the random
vector x follows joint density function ρ(x). We use the
shorthand Lap(κ) to denote the zero mean Laplace density
function with scale parameter κ > 0, Gamma(a, b) to denote
the Gamma density with parameters a, b > 0, Exp(λ) to denote

2

exponential density with rate λ > 0, and N (µ, Σ) to denote
the Gaussian density with mean µ and covariance Σ. We use
N (µ, σ, a, b) to denote a truncated univariate Gaussian density
with mean µ and standard deviation σ, supported on [a, b].
The symbols ∨ and ∧ denote the maximum and minimum,
respectively, and [r]+ := 0∨r, for any r ∈ R. The notation 1S
denotes the indicator function of a set S, while the abbreviation
i.i.d. stands for “independent and identically distributed”.

II. M ODEL

Fig. 2. A schematic of the proposed architecture for thermal inertial load
management by the LSE.

Fig. 1. Thermal inertial dynamics of an AC, not under direct load control,
is shown with fixed setpoint s0 , and a comfort temperature interval [L0 , U0 ]
with range 2∆, where U0 = s0 + ∆ and L0 = s0 − ∆. The ambient
temperature trajectory (green) is θa (t). The indoor temperature trajectory
θ(t) ∈ [L0 , U0 ] consists of alternating OFF (blue, up-going) and ON (red,
down-going) segments, as shown, where the boundaries L0 and U0 act as
reflecting barriers for the ON and OFF segments, respectively.

A. Dynamics of Thermal Inertial Load
In this paper, we consider the ambient temperature is high
enough that ACs are cooling homes, a situation of the sort
where demand response is appropriate. We denote the indoor
temperature at time t by θ(t), and the ambient temperature by
θa (t). At time t = 0, an occupant privately sets a temperature
s0 , called setpoint, with a willingness to tolerate at most ±∆
temperature deviation from s0 , thus defining a temperature
comfort range [L0 , U0 ] := [s0 − ∆, s0 + ∆]. The flexibility
of an AC is measured by 2∆; larger the flexibility, greater
the potential for demand response. The goal of the proposed
architecture and algorithms is to enable the LSE exploit this
flexibility without the LSE being aware of the load temperature
θ(t) or the setpoint s0 , or even ∆.
The dynamics of an AC not under direct load control, is
shown in Fig. 1. When the AC is OFF, the indoor temperature
θ(t) rises exponentially toward θa (t) until it hits the upper
comfort limit U0 , at which point the AC turns ON. Once ON,
θ(t) decreases exponentially until it hits the lower comfort
limit L0 , at which point the AC turns OFF. While this
qualitative behavior shown in Fig. 1 is valid for any ambient
temperature trajectory satisfying θa (t) > U0 , quantitatively
speaking, temporal variations in θa (t) produce time-varying
rates of heating (when AC is OFF) and cooling (when AC
is ON). Let us denote the mode of an AC at time t by

σ(t) = 1(0), signifying if it is ON (OFF). Its hysteretic
switching trajectory is


if θ(t) = U0 ,
1
(1)
σ(t) := 0
if θ(t) = L0 ,


−
σ (t )
otherwise.
The indoor temperature dynamics of the home is governed by
θ̇(t) = −α (θ(t) − θa (t)) − βP σ(t),

(2)

where the parameters α, β > 0 denote the heating time constant and thermal conductivity, respectively. Here, P denotes
the amount of thermal power drawn by the AC when it is
ON (σ(t) = 1). The electrical power drawn, Pe , is related
to the thermal power P , via the formula Pe = Pη , where the
parameter η > 0 denotes the efficiency of the AC.
B. Direct Load Control Strategy
In this paper, we propose an architecture where the setpoint
of an individual AC will be changed over time by the LSE
through a velocity control strategy. The rate of change of
setpoint s(t) will constitute a “control variable” that can be
used by the LSE. This velocity command, suitably modulated
so that the comfort zone is never violated, can be followed
by an AC locally without the LSE being aware of the load’s
chosen values of L0 , U0 or its current temperature θ(t). Yet,
by squeezing the band close to the upper or lower limits,
it can cause a certain fraction of loads to turn ON/OFF in
response, thereby collectively invoking demand response. ACs
that are ON with temperature trajectory headed downward
can be intercepted by a rising lower boundary (Lt ) to turn
OFF. Similarly, ACs that are OFF with temperature trajectory
headed upward can be intercepted by a falling upper boundary
(Ut ) and to force them turn ON. Thus, the demand response
can induce less electricity consumption by issuing positive
velocity commands, and it can also induce more cooling which
allows storage of energy in the homes. This will be detailed
in Section III-1.

3

Fig. 3. The two layer hierarchical strategy of the LSE.

If a home agrees to allow the LSE to manage the energy
consumption of its AC, the LSE is obligated to respect privacy
in the sense that it may influence the consumption but should
do so without ever measuring the state {s(t), θ(t), σ(t)} of
any individual AC. Furthermore, the LSE must guarantee that
the indoor temperature θ(t) lies within the home’s comfort
range [L0 , U0 ], at all times. These privacy and comfort range
guarantees will constitute a contract between the LSE and
an individual home. Larger the value of ∆, greater is the
flexibility of the load that an LSE can potentially exploit in
shaping demand response. The value of the contract therefore
depends on the magnitude of ∆. A home may choose its
tolerance ∆, and is charged lesser for electricity the larger
∆ it can tolerate. In Section IV, we will show how the LSE
can determine the cost of providing electricity to a customer
as a function of its flexibility ∆, which it can then use to
determine how to price contracts.

in the day-ahead market so that it can support the aggregate
load of the homes under contract. We term this operational
planning. For t ∈ [0, T ], the LSE takes the measured real time
total power consumption Ptotal (t) as feedback to compute a
control signal v(t) that it broadcasts to all homes, requesting
the setpoint of the ith AC at time t, denoted as si (t), be moved
according to the velocity control policy
ṡi (t) = ∆i v(t),

i = 1, . . . , N.

(3)

Homes with different contractual ∆i correspond to different
comfort tolerances, and experience different setpoint velocities. The operational planning, followed by setpoint velocity
control, comprise a two layer hierarchical strategy (Fig. 3) for
the LSE, which we expound next.
6

Ptotal (t)
Power (kW)

C. Contracts

4
ref
Ptotal
(t)

2

D. Planning in the Day-Ahead Market
0

III. M EASUREMENT AND C ONTROL A RCHITECTURE AND
S TRATEGY FOR THE LSE
We propose a novel architecture (Fig. 2) and operational
planning and control strategies for the LSE to manage a
population of N ACs over a finite time interval [0, T ], subject
to the contractual obligations, privacy constraints, and load
dynamics constraints, and to do so at low cost to the LSE.
Our strategy consists of a planning phase where the LSE
contracts to purchase the optimal amount of energy at low cost

0

2000

4000

6000

8000

10000

12000

14000

10000

12000

14000

t (seconds)
30

Temperature (◦ C)

Much of the electricity purchase takes pace in the dayahead market. In Section III-5, we will show how the LSE
can operationally plan the amount of energy to purchase in
each of the 24 one-hour periods of the day-ahead market, so
that it can support the loads at least cost. The goal of such
planning is to ensure that the LSE purchases more energy at
times when the day-ahead prices are low, and less energy in
periods when day-ahead prices are high. The LSE must do so
taking into consideration the flexibilities of the collection of
homes and the need to keep their temperatures within their
comfort zones.

U0

28

Ut

θ(t)

26
24

Lt

s(t)

22

L0

20
0

2000

4000

6000

8000

t (seconds)

Fig. 4. Top: Using setpoint velocity control, the LSE makes the real time
consumption Ptotal (t) (solid line), shown here for a single AC with θa (t) =
32◦ C and t ∈ [0, 4 hours], track the planned optimal reference consumption
ref (t) (dashed line). Bottom: corresponding nonlinear deformation of the
Ptotal
hysteretic band [L0 , U0 ](thick black lines) 7→ [Lt , Ut ] ([blue line, red line]),
shown as the gray area, in real time under the setpoint velocity control. Also
shown are the θ(t) (solid green) and s(t) (dashed brown) trajectories.

1) Privacy Preserving Load Control: We suppose that the
day-ahead operational planning, to be detailed in Section III-5,
has led to an optimal reference aggregate power consumption
ref
trajectory Ptotal
(t) for 0 ≤ t ≤ T . We now focus on the

4

design of a real time privacy preserving setpoint velocity
controller for the LSE so as to make the actual total power
consumption Ptotal (t) track the planned optimal target conref
(t). This corresponds to the “second layer” in
sumption Ptotal
Fig. 3. Assuming that the LSE can measure or estimate the
aggregate output Ptotal (t) in real time, define the error signal
ref
(t). Then, a simple PID controller
e(t) := Ptotal (t) − Ptotal
can generate a common
velocity command for all homes
Rt
v(t) := kp e(t)+ki 0 e(ς)dς +kd de
dt , with kp , ki , and kd being
the proportional, integral, and derivative gains, respectively.
We choose PID controller since it is easy to implement.
Alternatively, one may use a model-based feedback controller
for which the “number density” [10] of ON/OFF ACs is
approximated in real time. As a model-free implementation,
the PID controller obviates this difficulty.
The LSE broadcasts the common velocity command v(t) to
all homes. Upon receiving the LSE’s broadcast, an individual
home’s AC moves its setpoint according to (3). At any given
time, a home with larger comfort range moves its setpoint
faster than the one with smaller range, and hence provides
more flexibility to the LSE for reference tracking. To respect
the contractual comfort guarantee Li0 ≤ θi (t) ≤ Ui0 for each
i = 1, . . . , N , the setpoint controller (3), at any time t, utilizes
the time-varying band [Lit , Uit ], where
Lit
Uit

:= Ui0 ∧ [Li0 ∨ (si (t) − ∆i )] ,

:= Li0 ∨ [Ui0 ∧ (si (t) + ∆i )] .

(4)
(5)

Thus, the width of the time-varying bands expand (contract)
when the aggregate power consumption needs to be increased
(reduced). In Fig. 4, we show the deformation of the timevarying band [Lit , Uit ] ⊆ [Li0 , Ui0 ] for N = 1. Furthermore,
we notice that only the computation of v(t) is centralized, but
the computation (3), (4) and (5) are all decentralized.
From an architectural standpoint, there are two ways to interpret the computation in equations (3)–(5). One may interpret
them as part of control, and in that case the proposed real-time
control algorithm becomes partially centralized (computation
of v(t)) and partially decentralized (computation of (3)–(5)).
This is unlike existing direct load control programs [13]–[16],
where the controls are fully centralized while the actuation
are decentralized. Alternatively, one can consider only the
centralized computation of v(t) as control, and regard the
decentralized computation of (3)–(5) as part of smart actuation. In this interpretation, the proposed method differs from
the existing direct load control programs in that traditional
thermostatic actuation simply realize the control commands,
whereas in our case, additional computation are performed at
the actuators to modify the received control commands before
realizing them.
2) Estimating the Total Power Consumption: We note first
that the setpoint control loop is closed only over the aggregate
total power Ptotal (t) and not the individual homes’ power
draws, thus requiring no intrusive sensing of every home. The
LSE only needs to measure or estimate Ptotal (t) in real time,
for the AC population it manages.
Using the technique of differential privacy [17], we now
show how this can be estimated even if the population of

Fig. 5. A schematic for the privacy preserving sensing architecture for the
LSE to measure the real-time aggregate power Ptotal (t) for controlling the
setpoints of N ACs it manages. Each AC, at any given time, consumes either
zero or Pe power. For small  > 0 and fixed p ∈ (0, 1), each home adds local
1
i.i.d. noise from Gamma( pN
, pP ), to the individual consumption at time t,
e
and then transmits the noisy data with probability p. The sum of the reported
1
consumptions, is first scaled by p , then perturbed further by subrtactive noise
ν ∼ Exp( P ), and the resulting value Petotal (t) is sensed by the LSE. This
e
guarantees (Appendix A) that at the LSE level, the individual consumption
data remains -differentially private, while the data transmitted by each AC
remains secure against eavesdropping.

customers for the LSE constitutes only a subset of the total
population whose consumption is being measured at a distribution substation. The ith customer’s home, with fixed probability
p, reports a noise corrupted version Pbi (t) := Pi (t) + ni of
its power consumption Pi (t) to the LSE at time t, where
1
ni ∼ Gamma( pN
, pP e )’s are chosen independently and identically distributed by the homes. The probability p ∈ (0, 1)
b , where N
b is the number of customers’
satisfies pN = N
homes that report consumption to the LSE. From this, the LSE
generates an estimate of total consumption for its population
of customers Petotal (t) = Pbtotal (t) − ν, where Pbtotal (t) :=
PNb b

N
i=1 Pi (t), and ν ∼ Exp( Pe ). It can be shown (Appendix
b
N
P
b
N
Pi + n, where the additive zero
A) that Petotal (t) := N
b
N

i=1

mean noise random variable n ∼ Lap( Pe ), and the standard
differential privacy formulation for the sum query ensures that
an adversary would not be able to infer the N ×1 consumption
>
vector (P1 (t), . . . , PN (t)) ∈ {0, Pe }N , at any time t, if it
were to snoop Petotal (t). Also, even if Pbi (t) are eavesdropped,
the same -level of global differential privacy can be attained
by this local implementation of differential privacy (different
from the local differential privacy [18], [19]), as shown in Fig.
5. In Section V, we will computationally analyze the tradeoff
between privacy and tracking performance.
3) Performance of the LSE’s Strategy: Since computing
ref
Ptotal
(t) depends on both price and ambient temperature
forecasts π
b(t) and θba (t), the LSE’s plan, and consequently
the tracking performance during real time setpoint velocity
control, are affected by the uncertainties in these forecasts,
as well as the uncertainties in real time ambient temperature
θa (t). Let Ω be the set of all tuples of forecasted price, forecasted and real time ambient temperatures. Then each forecast
scenario ω ∈ Ω refers to a triplet of correlated trajectories
(b
π (ω, t), θba (ω, t), θa (ω, t)), and engenders an effective width
stochastic process for the ith AC’s band at time t, denoted

5

hereafter as weff (i, ω, t), given by
+

Z t
, (6)
v(ω, ς)dς
weff (i, ω, t) := Uit − Lit = ∆i 2 −
0

where the last equality results from substituting (4) and
(5) for Lit and Uit respectively, followed by some algebra
(Appendix B). To statistically quantify the performance of the
LSE’s control strategy, we observe that the limit of control
performance for the ith AC, is
ξT (ω) :=

1
1
lim
T ↓0 

Z
0

Normalized histogram

0.04

T

1{0≤weff (i,ω,t)<} dt,

(7)
0.01

which is the normalized local time [21] for the effective
width stochastic process. In words, the random variable ξT (ω)
measures the fraction of time ith AC remains inflexible over
the fixed horizon [0, T ], for a scenario ω ∈ Ω. Notice that
(7) does not depend on i, which is a consequence of (3), i.e.,
zero effective width epochs are synchronized among all ACs.
In Section V, we will compute ξT (ω) for different forecast
scenarios of (b
π (ω, t), θba (ω, t), θa (ω, t)), available from the
historical data.
4) Assumptions: We make the following assumptions.
ba (t), θa (t), and
• The ambient temperature trajectories θ
the price forecast π
b(t) are positive smooth (continuously
differentiable) functions of time t.
• ACs are necessary, i.e., for all t ∈ [0, T ], we have
θba (t), θa (t) > maxi=1,...,N Ui0 .
• Without loss of generality, the initial indoor temperatures
θi0 := θi (0) ∈ [Li0 , Ui0 ], for all i = 1, . . . , N .

5) Operational Planning: We now address how the refref
erence trajectory Ptotal
(t) can be chosen so as to minimize
the cost to the LSE. We suppose that the LSE is exposed
to a price forecast π
b(t) over [0, T ] from a day-ahead market.
Suppose that the LSE serves N homes with comfort tolerances
{∆i }N
i=1 . Further, for simplicity of exposition, we assume that
each AC in the population, when ON, draws same thermal
power P . The homes may have different thermal coefficients
(α, β), and the LSE only knows a joint distribution ρα,β across
the population of its customers. To provision an optimal aggreref
gate consumption profile Ptotal
(t) over the time window [0, T ],
the LSE makes its decision based on a price forecast π
b(t)
available from the energy market, energy budget E available
from the load forecast, and predicted ambient temperature
trajectory θba (t) available from the weather forecast. This is
depicted as the “first layer” in Fig. 3.
For i = 1, . . . , N , let the indicator variable ui (t) ∈ {0, 1}
denote whether the ith AC is ON/OFF at time t. Let nON (t)
denote the number of ACs that are ON at time t. Since
the aggregate electrical power drawn byPthe AC population
N
at time t is Pe (t) = Pη nON (t) = Pη i=1 ui (t), and the
procurement
cost for total energy consumption over [0, T ] is
R T ref
P
(t)b
π
(t)dt,
the LSE’s objective is to
total
0
P
u1 (t),...,uN (t)∈{0,1}N η

Z

minimize

0

T

π
b(t)

N
X
i=1

ui (t) dt,

(8)

3.7

4.6

↵ (⇥10

3

seconds

1

)

7.4

9.2
(⇥10

3

C
)
kW seconds

Fig. 6. Normalized histograms for α ∼ N (µα , 0.1µα , 0.9µα , 1.1µα ) (left),
1
and for β ∼ N µβ , 0.1µβ , 0.9µβ , 1.1µβ (right), with µα = RC
h-1 ,
1 ◦
C/kWh, used to generate the heterogeneous population of thermal
µβ = C
inertial loads.

Fig. 7. Initial conditions for N = 500 homes, for four different contractual
comfort tolerance (∆) distributions, shown left to right. The red circles denote
ON (σ0 = 1), and the blue circles denote OFF (σ0 = 0) initial states.

subject to the constraints


θ̇i (t) = −αi θi (t) − θba (t) − βi P ui (t),
Z
N
P TX
ui (t) dt = E,
η 0 i=1
Li0 ≤ θi (t) ≤ Ui0 ,

(9)
(10)
(11)

where [Li0 , Ui0 ] := [si0 −∆i , si0 +∆i ], and i = 1, . . . , N . We
notice that the dynamics constraints (9) and the comfort range
constraints (11) are decoupled, while the cost function (8) and
the energy budget constraint (10) are coupled. The constraint
(10) provides a modeling flexibility for the LSE, in case a total
energy budget E > 0 is available from load forecast.
Equations (8)–(11) represent a continuous time deterministic optimal control problem. Denoting its solution as
{u∗i (t)}N
i=1 , expressed as an open-loop optimal control, the
ref
optimal reference consumption is given by Ptotal
(t) =
PN ∗
P
u
(t).
In
Section
V,
we
will
follow
a
“discretize-theni=1 i
η
optimize” approach to numerically solve (8)–(11).
6) Feasibility of the Planning Problem: Let τ := ηE
P ,
and notice that constraint (10) imposes a necessary feasibility

6

shown; left bottom: the corresponding
Fig. 8. Left top: Increasing day-ahead price forecast π
bDA (t) as a piecewise constant function for each hour, with τ := ηE
P
optimal indoor temperature θ∗ (t), and ON/OFF control u∗ (t) obtained by first discretizing (8)–(11) for a single AC, and then solving the resulting MILP
(using Gurobi) and its LP relaxation (using MATLAB). In this simulation, [L0 , U0 ] = [20◦ C, 30◦ C], the ambient forecast is constant, i.e., θba (t) ≡ 32◦ C,
ref (t) computed from the MILP and its LP
and the parameters α, β, P, η are from Table 1, p. 1392 in [6]. Right: The optimal consumption trajectory Ptotal
relaxation, for the planning problem shown in the left.

condition
τ
ηE
0 ≤ τ :=
=
≤ 1.
(12)
NT
NPT
Given an ambient forecast θba (t) and parameters of the TCL
population, to respect constraint (11), τ is restricted to
0 ≤ τ ` ≤ τ ≤ τ u ≤ 1,

(13)

ηEu
`
where τ ` := Nτ`T = NηE
P T , and τ u = N P T . Here, E`
(resp. Eu ) is the aggregate energy consumed if the indoor temperatures for each home in the population were
to be restricted at their private upper (resp. lower) setpoint
boundaries, thus resulting in the lowest (resp. highest) total
energy consumption
while respecting (11). In other words,
PN R T
E` = Pη i=1 0 ui (t) dt, where the corresponding controls are ui (t) = βαi Pi (θba (t) − Ui0 ), and hence τ ` =
R
PN αi b
1
1 T b
b
i=1 βi (hθa i − Ui0 )), where hθa i := T 0 θa (t) dt.
NP (
PN
Similar calculation yields τ u = N1P ( i=1 αβii (hθba i − Li0 )).
Equation (13) gives the necessary and sufficient conditions
for feasibility of the planning problem (8)–(11).
that if
R T PN Notice
∗
the constraint (10) is inactive, then Pη 0
u
dt
=
E` .
i=1 i
Also, note that (13) is equivalent to the energy inequality
Emin ≤ E` ≤ E ≤ Eu ≤ Emax , where Emin := 0 and
Emax := N Pη T . For privacy preserving planning purpose, as
earlier, the LSE can use noisy estimates and random samples
for the initial conditions of the AC states (si0 , θi0 , σi0 )N
i=1 .
Similarly, the thermal coefficient parameters (αi , βi )N
i=1 can
be sampled from ρα,β .

customer of the given type. Assuming the LSE has an existing
customer
of N ACs, let the optimal value of (8) be
R Tpopulation
ref
J ∗ := 0 π
b(t)Ptotal
(t) dt. Recall that π
b(t) has unit $/MWh,
and for T = 24 hours, J ∗ (in $) denotes the daily energy
∗
procurement cost. For i = 1, . . . , N , let J−i
denote the optimal
th
value of (8) with the i AC removed from the population.
∗
∗
Clearly, J−i
> J ∗ , and the increase in per day cost J−i
− J∗
represents the marginal value ($/day) of an AC with comfort
∗
tolerance ∆i . Thus, the LSE can use the graph of J−i
− J ∗ as
function of ∆i , as a price chart for the thermal inertial load
management service. When an individual AC, as new customer
to the LSE, specifies its choice of comfort tolerance, this graph
will determine the price of the new contract between the LSE
and the AC. This will be illustrated further in Section V.
V. N UMERICAL R ESULTS
A. Parameters and Simulation Setup

IV. P RICING THE C ONTRACTS
To address the issue of how the LSE can estimate the cost
of a contract for an individual home with a given flexibility
∆, we employ a sensitivity analysis approach. The LSE can
simply calculate the cost (8) with and without an additional

Fig. 9. The day-ahead forecasted ambient temperature θba (t) (dashed blue),
real-time ambient temperature θa (t) (solid blue), the day-ahead forecasted
price of energy π
bDA (t) (dashed green), and real-time price of energy πRT (t)
(solid green) data used in Section V-B.

7

(kp , ki , kd )

(kp , ki , kd )

(a.1)

(a.2)

(b.1)

(b.2)

(a.3)

(a.4)

(b.3)

(b.4)

Fig. 10. (a) The setpoint tracking performance for the PID velocity controller implemented without (4) and (5); and (b) the same with (4) and (5), for different
combinations of PID gain tuple (kp , ki , kd ), shown as different colors in the legend box at the top. The four subfigures, denoted as (a.1)–(a.4) and (b.1)–(b.4),
refer to four different contractual ∆ distributions mentioned in Section V-A.

Following Table 1 (p. 1392) in [6], we fix the thermal
resistance R = 2◦ C/kW, the thermal capacitance C = 10
kWh/◦ C, the thermal power drawn by an ON AC as P =
14 kW, and the load efficiency η = 2.5. We choose the
distribution of thermal coefficient parameters as truncated
Gaussians (Fig. 6), given by α ∼ N (µα , 0.1µα , 0.9µα , 1.1µα )
1
and β ∼ N (µβ , 0.1µβ , 0.9µβ , 1.1µβ ), with µα = RC
h-1 ,
1 ◦
µβ = C C/kWh. We consider a population of N = 500
homes, a time horizon of T = 24 hours, and suppose that
the consumers have comfort tolerances (∆) between ∆min =
0.1◦ C and ∆max = 1.1◦ C. We investigate four cases: (i) the
entire population of customers has identical ∆ ≡ 1◦ C, (ii)
the population has i.i.d. uniform ∆ over [∆min , ∆max ], (iii)
the population has i.i.d. right triangular ∆ distribution over
[∆min , ∆max ] with peak at ∆max , (iv) same as (iii) but with
peak at ∆min . Intuitively, in case (iii) (case (iv)), there are
more homes with large (small) ∆, and the population offers
higher (lower) flexibility for the LSE’s control strategy.
We sample the initial conditions (si0 , θi0 ) from a bivariate
Gaussian N (µ0 , Σ0 ), with
µ0 = (20, 20)> and
 mean vector

1 0.5
covariance matrix Σ0 =
, subject to the constraint
0.5 3
si0 − ∆i ≤ θi0 ≤ si0 + ∆i , where the comfort tolerance
samples ∆i are drawn from distributions according to case (i)(iv), and i = 1, . . . , N . For each tuple (si0 , θi0 ) thus generated,
we assign σi0 = 1 or 0 with probability 0.5. The resulting
initial conditions are shown in Fig. 7.

step size, and use Euler discretization for the ODE (9), to
transcribe the optimal control problem (8)–(11) as a mixed
integer linear programming (MILP) problem in decision vector
{θ(1), u(1), . . . , θ(T ), u(T )}. In Fig. 8 left bottom and right,
we compare the MILP solution, computed using Gurobi [20],
with the solution of its linear programming (LP) relaxation,
computed using MATLAB. Intuitively, since π
bDA (t) is increasing, it is optimal to stay ON during [0, τ ], and OFF
∗
thereafter. We observe that (θMILP
(t), u∗MILP (t)) matches with
∗
∗
(θLP (t), uLP (t)) except at times t ∈ [thit , τ ], where thit is
the first hitting time to L0 . In particular, during [thit , τ ], the
∗
∗
solution θMILP
(t) (resp. θLP
(t)) chatters (resp. slides) along L0 .
In the chattering period, u∗MILP (t) switches between 0 and 1 at
every time step, and as shown in Fig. 8 right, u∗MILP (t) ∈ {0, 1}
averaged over the chattering period yields u∗LP (t) ∈ [0, 1], i.e.,
the proper fractional values of u∗LP (t) can be interpreted as the
fraction of time the AC remains ON during [thit , τ ].

B. Case Studies

2) Planning and Setpoint Control for Actual Data: We
now apply the LSE’s strategy proposed in this paper, to
the day-ahead price forecast π
bDA (t) data for Houston on
August 10, 2015, available [22] from ERCOT. Further, we
use Houston weather station data, for day-ahead forecasted
ambient temperature θba (t) on August 10, 2015, and for real
time ambient temperature θa (t) on August 11, 2015. With

1) Direct Numerical Solution of the Planning Problem: To
elucidate the numerical solution of (8)–(11), let us consider
N = 1 home with [L0 , U0 ] = [20◦ C, 30◦ C], θ0 = 25◦ C,
constant ambient temperature forecast θba (t) = 32◦ C, dayahead price forecast π
bDA (t) and τ as shown in Fig. 8,
left top. We discretize (8)–(11) with ∆t = 1 minute time

For N homes with 1 minute time step size, the planning
problem requires solving an MILP in 24×60×N ×2 variables,
and in our numerical experiments with Gurobi, solving the
MILP even for N = 2 was found to have more than 24
hours of CPU runtime. Armed with the physical meaning of
the LP relaxation discussed above, in this paper, we solve
(8)–(11) under control convexification (u1 (t), . . . , uN (t)) ∈
{0, 1}N 7→ (u1 (t), . . . , uN (t)) ∈ [0, 1]N . This LP relaxation
allows us to compute for large N (e.g. N = 500 in 10 minutes
CPU time of which only 36 seconds are needed to run the LP
solver), which is of interest from the LSE’s perspective.

TABLE I
E NERGIES ( IN MW H ) MENTIONED IN S ECTION III-6, ASSOCIATED WITH
THE OPTIMAL PLANNING PROBLEM IN S ECTION V-B2. T HE GRAY
COLUMNS SHOW THE REAL TIME CONTROLLED (Ec ) AND
UNCONTROLLED (Eunc ) ENERGY CONSUMPTIONS IN MW H .
∆ cases

Emin

E`

E

21.5088

Ec

Eunc

Eu

22.4681

25.6821

26.3328

Emax

1.0

22.1496
0.1

22.4485

25.6905

26.6778

1.1

0

0.1

1.1

0.1

1.1

22.4

67.2

21.9287

22.4690

25.7193

25.9663

22.5625

22.4299

25.0748

25.2126

Fig. 11. Tracking performance of the PID setpoint velocity controller with
local implementation of differential privacy, as described in Section III-2,
with  = 0.1, p = 0.9, and (kp , ki , kd ) = (10−4 , 10−6 , 10−4 ). The
ref (t) for ∆
black curve above is the same as that in Fig. 10(b.2), i.e., Ptotal
uniformly distributed in [0.1, 1.1]◦ C. Different colored Ptotal (t) correspond
to 10 different realizations of the noise random variable n.

these price and ambient temperature data (Fig. 9), and using
the parameter values in Section V-A, we first solve the dayahead planning problem described in Section III-5 to design
the optimal aggregate consumption, and then track that reference consumption in real time via setpoint velocity control
described in Section III-1. As in Section V-B1, the planning
problem is solved with time step size equal to 1 minute,
and with τ = 31 , which is feasible from (13). The setpoint
velocity control is implemented with a smaller time step size
of 1 second, to emulate the actual dynamics possibly being
different from the model used in planning computation.
In Fig. 10, two sets of four subfigures, labeled as (a.1)–(a.4)
and (b.1)–(b.4), are shown. The subfigures (a.1)–(a.4) (resp.
(b.1)–(b.4)), show the tracking performance of the setpoint
velocity controller (3) in the absence of noise n, without (resp.
with) enforcing the comfort range constraints (11) in real time
via (4) and (5). The four subfigures in each set correspond
to the four contractual cases mentioned in Section V-A. The
ref
black curve in each subfigure denotes Ptotal
(t) obtained by
solving the discretized LP for the control convexified (see
Section V-B1) version of (8)–(11) with the corresponding
initial conditions and parameters. As shown in Table I, for

Ptotal (t) (kW)

8

t (hours)

t (hours)

(a)

(b)

Fig. 12. Tracking performance of the PID setpoint velocity controller with
ref (t) same as in
(kp , ki , kd ) = (10−4 , 10−6 , 10−4 ), and black curves Ptotal
Fig. 10(b.2) and 11, when individual consumption data at each time are
(a) corrupted by ±10% i.i.d. uniform noise about the true values, and (b)
both noise corrupted as in (a) as well as incomplete data reported due to
communication by (each time randomly selected) 400 out of N = 500 homes.
The blue curves are the controlled real-time aggregate consumption. There is
no significant difference in the tracking performance between those shown
above and the same in Fig. 10(b.2) and 11.

R T ref
(t) dt satisfies the energy
each contractual case, E := 0 Ptotal
inequality mentioned in Section III-6.R The gray columns in
T
Table I show the actual consumptions 0 Ptotal (t) dt with and
without PID control.
Different colored Ptotal (t) trajectories in Fig. 10 correspond
to different combinations of the PID gain tuple, shown in the
top legends. In both sets (a.1)–(a.4) and (b.1)–(b.4), the real
time aggregate consumption Ptotal (t) trajectories, in general,
ref
show more oscillation than the optimal Ptotal
(t), due to different time scales for operational planning (1 minute) and
real time tracking (1 second); the gain tuple (kp , ki , kd ) =
(10−4 , 10−6 , 10−4 ) (shown in maroon), in general, performs
better than others. Fig. 10 (a.1)–(a.4) show almost perfect real
time tracking for all four ∆ distributions, since the movement
of the comfort boundaries are unrestricted. In the presence
of (4) and (5), as shown in (b.1)–(b.4), the tracking accuracy
is limited compared to (a.1)–(a.4), and differs for different
contractual comfort (∆) distributions.
In subfigure (b.1), all ACs have identical ∆ ≡ 1◦ C, and
hence identical weff (t) for all t. Notice that about 16 hours
ref
onward, Ptotal
(t) is decreasing, and to follow it in real time, the
LSE moves all ACs’ setpoints up at identical rates, resulting
in all hitting their respective upper boundaries Ui0 simultaneously at around 20 hours. Between 20–23 hours, weff (t) = 0
ref
at Ui0 , and Ptotal (t) cannot track further decrease in Ptotal
(t),
as seen in subfigure (b.1). In subfigures (b.2)–(b.4), unlike
(b.1), the setpoint boundaries for different ∆i are moved at
different rates. In (b.3) (resp. (b.4)), more ACs have large (resp.
small) ∆i , thus allowing the LSE more (resp. less) flexibility in
real time tracking, and result in less (resp. more) overshoots
in Ptotal (t). Comparing Fig. 10(b.2) with Fig. 11, it is seen
that the privacy for individual consumptions can be preserved
via local noise injection described in Section III-2, without
affecting the tracking performance. Since the setpoint velocity
control is reactive to error in the aggregate consumption, hence
the tracking performance is robust against outliers, noise,
and/or incompleteness in individual consumption data. In other
words, the robustness observed in Fig. 11 is typical in the
sense that individual noise statistics or incomplete data does
not significantly affect the aggregate tracking performance, as
evidenced in Fig. 12.
To demonstrate the statistical performance of the LSE’s

9

5 PM

Fig. 13. We use historical data for 30 days in August 2015, to generate 30 scenarios of θba (ω, t) (top left), θa (ω, t) (top right), π
bDA (ω, t) (bottom left), for
ref (ω, t) sample paths are shown in the bottom right.
analyzing the performance of the LSE’s strategy. The corresponding Ptotal

strategy, as discussed in Section III-3, we take the historical
data for forecasted day-ahead price and ambient temperature
for 30 days in August 2015, and use these as the scenario
set Ω. One specific scenario ω ∈ Ω, refers to a trajecref
tory tuple (b
πDA (ω, t), θba (ω, t), θa (ω, t)) that yields Ptotal
(ω, t)
and ξT (ω) (see (7)). For ω ∈ Ω, Fig. 13 shows all scenarios (b
πDA (ω, t), θba (ω, t), θa (ω, t)), and the corresponding
ref
Ptotal
(ω, t) for ∆ uniformly distributed between [0.1, 1.1]◦ C
with respective initial conditions and parameters mentioned
in Section V-A, and τ = 0.45 feasible for all ω ∈ Ω. Fig.
14 shows that for most scenarios, the normalized local time
ξT (ω) is close to unity, meaning that under the proposed
strategy for the LSE, each AC would have zero effective
widths for more than 90% of times in a day, as expected in
hot days in Houston.
In Fig. 15, we plot results for the sensitivity analysis
described in Section IV, with the data shown in Fig. 9, for
∆ uniformly distributed in [0.1, 1.1]◦ C. The red line in Fig.
15 is the least squares linear fit for the scattered data, and can
be used to price the contract with a new home for its chosen
comfort tolerance ∆. As expected intuitively, a smaller ∆ leads

to higher per day service charge, and vice versa.
C. Computational Complexity
The computational complexity of the proposed method is
twofold: that of the optimal planning, and that of the setpoint
control. The first being an offline computation, scales as a function of both the number of homes N , and the discretizationT
of the time horizon [0, T ], where ∆t is the
length µ := ∆t
discretization step-size for the ODE (9). The setpoint control
being an online computation, scales as a function of N .
As described in Section V-B1, Euler discretization for ODE
(9) leads to an LP in decision vector {θ1 (1), u1 (1), θ2 (1),
u2 (1), . . . , θN (1), uN (1), . . . , θN (µ), uN (µ)} of length 2N µ
with N µ + 1 equality constraints and 2N µ inequality constraints. Here, θi (k) and ui (k) respectively denote the indoor
temperature and optimal control for the ith home at time t =
k∆t, where i = 1, . . . , N , and k = 1, . . . , µ. We solved this
LP via MATLAB linprog using Mehrotra’s implementation
of primal-dual interior point solver [23], that has worst-case
complexity polynomial in N µ. For our particular instance of
linear program and for ∆t = 1 minute (i.e., µ = 1440), the

10

runtime complexity scaling with N is shown in Fig. 16. It
is easy to verify that executing equations (3), (4) and (5) for
real-time setpoint control, requires O(N ) runtime complexity.
The LPs for computing optimal reference power consumption trajectories, were solved in MATLAB R2015a running
on Dell PC with 3.4 GHz Intel i7-2600 processor and 16
GB memory. The same platform was used for solving the
MILP mentioned in Section V-B1 using Gurobi. All other
simulations including setpoint velocity control were performed
in MATLAB R2013a running on a Macbook Pro with 2.6 GHz
Intel Core i5 processor and 8 GB memory.
VI. C ONCLUSION
We have presented an architecture and control algorithms
for an aggregator or load serving entity managing a population of thermal inertial loads. The proposed hierarchical
architecture solves the twin problems of minimum cost energy
procurement from the day-ahead market, and real-time setpoint
control to track the reference aggregate power trajectory that
corresponds to the optimal energy procurement. The control
algorithms respect privacy and contractual comfort at the
individual consumer level, but can be implemented at the
aggregate level. We demonstrate these algorithms using the
day-ahead price data from ERCOT, and using the forecasted
and real time ambient temperature data from a weather station
in Houston, Texas. Future work will investigate the potential of
such contracts as innovative business models for load serving
entities.
A PPENDIX A
D IFFERENTIALLY P RIVATE E STIMATE OF T OTAL P OWER
At time t, let us consider two consumption vectors
c(t), c0 (t) ∈ {0, Pe }N . We call two consumption vectors
neighboring if they differ by at most one element. Then,
for any neighboring consumption pair c(t), c0 (t), and for
sum query SUM : {0, Pe }N 7→ {0, Pe , 2Pe , . . . , N Pe }, we
have supc(t),c0 (t) |SUM(c(t)) − SUM(c0 (t))| ≤ Pe . Hence,
-differential
PNprivacy [17] at the LSE level is guaranteed if
Petotal (t) = i=1 Pi + n, where n ∼ Lap Pe .
To see that the strategy proposed in Section III-2 achieves
1
-differential privacy, notice that ni ∼ Gamma( pN
, pP e ) ⇒
1
1

b
e :=
p ni ∼ Gamma( pN , Pe ), and since N = pN , we have n
PNb 1

i=1 p ni ∼ Exp( Pe ). We thus arrive at

0.000

⇡

Fig. 14. The normalized local time random variable ξT (ω) (see (7)) is shown
as a function of the 30 scenarios described in Fig. 13.

Fig. 15. The scatterplot with blue circles denotes the increase in per day
energy procurement cost, if an AC with comfort tolerance ∆i , is removed
from the existing population of 500 ACs, as described in Section IV. The red
line is the least squares linear fit that the LSE can use for pricing the contract
of a new home.

and hence we can rewrite (4) and (5) as

b
N

1X
Pi (t) + (e
n − ν),
Petotal (t) = Pbtotal (t) − ν =
p i=1

(14)

A PPENDIX B
D ERIVING EQUATION (6)
For a, b, x ∈ R, recall that
(a − x) ∨ (b − x)
− (0 ∧ a)

=
=

=

(a ∧ b) + x,
(a ∨ b) − x,
+

[−a] ,

(16)

=

(15)

PNb
where the first term p1 i=1 Pi (t) provides an estimate for
PN
the true aggregate consumption i=1 Pi (t); the second term
n := n
e − ν ∼ Lap( Pe ) since n
e, ν are i.i.d. Exp( Pe ).

(a + x) ∧ (b + x)

Lit

(15)
(16)
(17)

=

Uit

(15)

=

(16)

=

(si0 + ∆i ) ∧ [(si0 ∨ si (ω, t)) − ∆i ]
{si0 ∧ [(si0 ∨ si (ω, t)) − 2∆i ]} + ∆i , (18)
(si0 − ∆i ) ∨ [(si0 ∧ si (ω, t)) + ∆i ]
{si0 ∨ [(si0 ∧ si (ω, t)) + 2∆i ]} − ∆i . (19)

This results weff (i, ω, t) := Uit − Lit = T1 − T2 − 2∆i ,
where T1 := {si0 ∨ [(si0 ∧ si (ω, t)) + 2∆i ]}, and T2 :=
{si0 ∧ [(si0 ∨ si (ω, t)) − 2∆i ]}.
Rt
If 0 v(ω, ς)dς > 0, then si (ω, t) > si0 , ⇒ T1 = si0 ∨
[si0 + 2∆
∧ [si (ω, t) − 2∆i ] =
h i ] = si0R+ 2∆i , and T2 = isi0 (15)
t
si0 ∧ si0 + ∆i 0 v(ω, ς)dς − 2∆i
= si0 + ∆i {0 ∧

11

[7] S. Bashash, and H.K. Fathy, “Modeling and Control Insights into
Demand-side Energy Management through Setpoint Control of Thermostatic Loads”. Proceedings of the 2011 American Control Conference, pp.
30
4546–4553, 2011.
[8] J.L. Mathieu, S. Koch, and D.S. Callaway, “State Estimation and Control of Electric Loads to Manage Real-time Energy Imbalance”. IEEE
25
Transactions on Power Systems, Vol. 28, No. 1, pp. 430–440, 2013.
[9] W. Zhang, J. Lian, C.-Y. Chang, and K. Kalsi, “Aggregate Modeling
20
and Control of Air Conditioning Loads for Demand Response”. IEEE
Transactions on Power Systems, Vol. 28, No. 4, pp. 4655–4664, 2013.
15
[10] A. Ghaffari, S. Moura, and M. Krstić, “Analytical Modeling and Integral
Control of Heterogeneous Thermostatically Controlled Load Populations”. ASME 2014 Dynamic Systems and Control Conference, Paper No.
10
DSCC2014-6022, pp. V002T22A002; 8 pages, 2014.
[11] S. Meyn, P. Barooah, A. Bušić, Y. Chen, and J. Ehren, “Ancillary Service
5
to the Grid using Intelligent Deferrable Loads”. IEEE Transactions on
Automatic Control, Vol. 60, No. 11, pp. 2847–2862, 2015.
[12] A. Halder, X. Geng, G. Sharma, L. Xie, and P.R. Kumar, “A Control
1
2
3
4
5
System Framework for Privacy Preserving Demand Response of Thermal
2
Inertial Loads”. Proceedings of the 2015 IEEE International Conference
N (×10 homes)
on Smart Grid Communications (SmartGridComm), pp. 181–186, 2015.
[13] Summer Discount Plan of Southern California Edison in California,
Fig. 16. CPU time (in seconds) needed to execute MATLAB linprog for
http://tinyurl.com/zly97mt, accessed on Sep. 21, 2016.
solving the planning problems for N = 100, 200, 300, 400 and 500 homes,
[14] Saver’s Switch Program of Xcel Energy in Colorado, Minnesota and
corresponding to case (iv) in Section V-B2. The CPU time data for other cases
Wisconsin, http://tinyurl.com/z4zayh7, accessed on Sep. 21,
are almost identical to the one shown, and omitted for brevity. The horizontal
2016.
axis is in logarithmic scale while the vertical axis is in linear scale.
[15] Energy Wise Rewards Program in Maryland and District of Columbia,
http://tinyurl.com/z5pov8s, accessed on Sep. 21, 2016.
[16] Power Partner Thermostats of Austin Energy in Texas, http://
Rt
tinyurl.com/zbvzhyb, accessed on Sep. 21, 2016.
[ 0 v(ω, ς)dς − 2]}. This yields
[17] C. Dwork and A. Roth, The Algorithmic Foundations of Differential

Z t

Privacy, Foundations and Trends in Theoretical Computer Science, Now
Publishers, Delft, The Netherlands, 2014.
weff (i, ω, t) = −∆i 0 ∧
v(ω, ς)dς − 2
[18] S.L. Warner, “Randomized Response: A Survey Technique for Elimi0
nating Evasive Answer Bias”. Journal of the American Statistical Asso
+
Z t
(17)
ciation, Vol. 60, No. 309, pp. 63–69, 1965.
= ∆i 2 −
(20) [19] J.C. Duchi, M.I. Jordan, and M.J. Wainwright, “Local Privacy and
v(ω, ς)dς .
0
Statistical Minimax Rates”. 54th IEEE Annual Symposium on Foundations
Rt
of Computer Science, pp. 429–438, 2013.
On the other hand, if 0 v(ω, ς)dς ≤ 0, then si (ω, t) ≤ [20] Gurobi Optimization, Inc., “Gurobi Optimizer Reference Manual”.
http://www.gurobi.com, 2015.
si0 ,R⇒ T1 = si0 ∨ (si (ω, t) + 2∆i ) = si0 ∨ (si0 +
t
[21] M. Marcus, and J. Rosen, Markov Processes, Gaussian Processes, and
∆i 0 v(ω, ς)dς + 2∆i ), and T2 = si0 ∧ (si0 − 2∆i ) =
Local Times. 1st ed., Cambridge University Press, 2006.
si0 − 2∆i , resulting
[22] Electric Reliability Council of Texas (ERCOT) day-ahead price data.
http://tinyurl.com/z3gmvt6, accessed on Feb. 23, 2016.
Z t
[23] S. Mehrotra, “On the Implementation of a Primal-Dual Interior Point
weff (i, ω, t) = [si0 ∨ (si0 + ∆i
v(ω, ς)dς + 2∆i )] − si0
Method”. SIAM Journal of Optimization, Vol. 2, No. 4, pp. 575–601,
0
1992.

CPU time (seconds)

35

(16)

=


+
Z t
∆i 2 +
v(ω, ς)dς .

(21)

0

Combining (20) and (21), we arrive at (6).
R EFERENCES
[1] D.S. Callaway, and I.A. Hiskens, “Achieving Controllability of Electric
Loads”. Proceedings of the IEEE, Vol. 99, No. 1, pp. 184–199, 2011.
[2] U.S. Energy Information Administration, “Residential Energy Consumption Survey 2009”. http://www.eia.gov/consumption/
residential/index.cfm, accessed on May 3, 2016.
[3] M.C. Caramanis, R.E. Bohn, and F.C. Schweppe, “Optimal Spot Pricing:
Practice and Theory”. IEEE Transactions on Power Apparatus and
Systems, Vol. PAS-101, No. 9, pp. 3234–3245, 1982.
[4] C.-Y. Chong, and A.S. Debs, “Statistical Synthesis of Power System
Functional Load Models”. Proceedings of the 18th IEEE Conference on
Decision and Control including the Symposium on Adaptive Processes,
No. 18, pp. 264–269, 1979.
[5] 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,
1985.
[6] 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, 2009.

Abhishek Halder (S’10-M’14) received the B.
Tech. and M. Tech. degrees in aerospace engineering
from IIT Kharagpur in 2008, and the Ph.D. degree
in aerospace engineering from Texas A&M, College
Station, TX, USA, in 2014. He is a Postdoctoral
Research Associate at the Department of Electrical
and Computer Engineering, Texas A&M University.
His research interest is in systems, control and
optimization with focus on large-scale cyberphysical
systems such as the smart grid and aerial swarm
robotics. He received the 2014 Outstanding Doctoral
Student Award from Texas A&M and the 2008 Best Thesis Award from IIT,
both in Aerospace Engineering.

Xinbo Geng (S’14) received the B.E. degree in electrical engineering in 2013 from Tsinghua University,
Beijing, China, and the Master of Science degree
from the Department of Electrical and Computer Engineering at Texas A&M University in 2015, where
he is currently working toward the Ph.D. degree. His
research interests include demand response, electricity markets and data analytics in power system.

12

P. R. Kumar (F’88) received the B. Tech. degree
from IIT Madras, Chennai, India, in 1973, and the
D.Sc. degree from Washington University, St. Louis,
MO, USA, in 1977. He was a faculty member
at UMBC (19771984) and University of Illinois,
Urbana-Champaign (19852011). He is currently at
Texas A&M University. His current research is focused on stochastic systems, energy systems, wireless networks, security, automated transportation,
and cyberphysical systems.
He is a member of the US National Academy
of Engineering and The World Academy of Sciences. He was awarded a
Doctor Honoris Causa by ETH, Zurich. He has received the IEEE Field
Award for Control Systems, the Donald P. Eckman Award of the AACC,
Fred W. Ellersick Prize of the IEEE Communications Society, the Outstanding
Contribution Award of ACM SIGMOBILE, the IEEE Infocom Achievement
Award, and the SIGMOBILE Test-of-Time Paper Award. He is a Fellow of
IEEE and an ACM Fellow. He was Leader of the Guest Chair Professor Group
on Wireless Communication and Networking at Tsinghua University, is a D.
J. Gandhi Distinguished Visiting Professor at IIT Bombay, and an Honorary
Professor at IIT Hyderabad. He was awarded the Distinguished Alumnus
Award from IIT Madras, the Alumni Achievement Award from Washington
Univ., and the Daniel Drucker Eminent Faculty Award from the College of
Engineering at the Univ. of Illinois.

Le Xie (M’10-SM’16) received the B.E. degree
in electrical engineering from Tsinghua University,
Beijing, China, in 2004, the S.M. degree in engineering sciences from Harvard University, Cambridge,
MA, USA, in June 2005, and the Ph.D. degree
from the Electric Energy Systems Group (EESG),
Department of Electrical and Computer Engineering,
Carnegie Mellon University, Pittsburgh, PA, USA,
in 2009. He is an Associate Professor in the Department of Electrical and Computer Engineering at
Texas A&M University, where he is affiliated with
the Electric Power and Power Electronics Group. His research interest includes
modeling and control of large scale complex systems, smart grid applications
in support of renewable energy integration, and electricity markets.

