1

State Estimation for the Individual and the
Population in Mean Field Control
with Application to Demand Dispatch

arXiv:1504.00088v3 [cs.SY] 30 May 2016

Yue Chen, Ana Bušić, and Sean Meyn

Abstract—This paper concerns state estimation problems in a
mean field control setting. In a finite population model, the goal
is to estimate the joint distribution of the population state and
the state of a typical individual. The observation equations are
a noisy measurement of the population.
The general results are applied to demand dispatch for
regulation of the power grid, based on randomized local control
algorithms. In prior work by the authors it is shown that local
control can be designed so that the aggregate of loads behaves
as a controllable resource, with accuracy matching or exceeding
traditional sources of frequency regulation. The operational cost
is nearly zero in many cases.
The information exchange between grid and load is minimal,
but it is assumed in the overall control architecture that the
aggregate power consumption of loads is available to the grid
operator. It is shown that the Kalman filter can be constructed
to reduce these communication requirements, and to provide the
grid operator with accurate estimates of the mean and variance
of quality of service (QoS) for an individual load.

I. I NTRODUCTION
Mean field models are a valuable tool for design and
performance approximation for certain classes of interacting
systems [1]–[3]. The infinite-population mean-field equations
provide tremendous insight, but ultimately we must translate
this insight to address a finite-population reality. In this paper
we propose algorithms based on the Kalman filter to obtain
estimates of first and second order statistics of the population
and a typical individual, based on noisy observations of the
population.
While the potential applications are far broader than power
systems, for ease of exposition it is convenient to restrict
attention to one application.
Renewable energy sources such as wind and solar power
have a high degree of unpredictability and time variation,
which complicates balancing supply and demand. One possible
way to address this challenge is to harness the inherent
flexibility in demand of many types of loads.
Demand Response is traditionally meant as a reduction in
load in response to some grid-level event. It is in use today
for peak-shaving (smoothing demand), and for contingency
reserves (load-shedding following generation loss).
It is argued in [4]–[6] that the value of demand-side flexibility is far greater than this. Loads can supply a range of grid
Research supported by NSF grants CPS-0931416 and CPS-1259040, and
the French National Research Agency grant ANR-12-MONU-0019
Y.C. and S.M. are with the Department of Electrical and Computer Engg.
at the University of Florida, Gainesville. A.B. is with Inria and the Computer
Science Dept. of École Normale Supérieure, Paris, France.

services, such as the balancing reserves required at BPA, or
the Reg-D/A regulation reserves used at PJM [5]. These grid
services can be obtained without impacting quality of service
(QoS) for consumers [3], [7]. This is only possible through
design. The term Demand Dispatch, introduced in [8], is used
to emphasize the difference between the goals of our own work
and traditional demand response.
The application in this paper concerns a large collection of
loads whose power consumption is not continuously variable.
Examples include thermostatically controlled loads (TCLs), as
considered in [2], [9], and irrigation or pool-pumps [3], [6]. In
these papers and [10], [11] it is argued that randomization at
the load is valuable to avoid synchronization, and to simplify
control at the grid level.
In much of this prior work, a mean-field model is obtained
for control design at the grid level – this is a deterministic
model of the aggregate of loads, obtained as a law of large
numbers limit as the population of loads tends to infinity. The
control solution adopted in [2] is based on state-feedback for
a linear state space model with partial observations. Although
the mean-field model is bi-linear, it is represented as a linear
model by treating the product of inputs and states as a new
input; this is why state estimation is needed for implementation
of the algorithm. In [6] a randomized policy is designed for
each load so that the mean-field model is an input-output
system that is easily controlled without the use of state
estimates.
For simplicity, in this paper attention is restricted to the
setting of [6], in which each load evolves as a controlled
Markov chain. The transition probability is determined by its
own state, and a scalar signal ζ broadcast from a balancing
authority (BA). The extension to vector inputs, as in [2],
requires only changes in notation.
The common dynamics are defined by a controlled transition
matrix {Pζ : ζ ∈ R}. For the ith load, there is a state process
X i whose transition probability is defined by,
i
Pζ (x− , x+ ) = P{Xt+1
= x+ | Xti = x− , ζt = ζ}

(1)

where x− and x+ are possible state-values. In the case of a
water heater, the state x ∈ X might represent temperature of
the water, and whether the unit is operating or not.
If there are N loads operating independently, conditional on
the common signal ζ, then the empirical distribution (i.e., the

2

histogram of state values) is defined as the average,
µN
t (x) =

1
N

N
X

I{Xti = x} ,

i=1

x∈X

Viewed as a row vector, the following recursion is central to
the analysis in [3], [6]:
N
T
µN
t+1 = µt Pζt + Wt+1 .

An observation model is also linear in the state,
X
Yt =
µN
t (x)U (x) + Vt

(2)

(3)

x

where U : X → R. In applications to demand dispatch, U (x)
represents
P power consumption of a load when its state is x,
so that x µN
t (x)U (x) is the average power consumption at
time t.
It is established in [3], [6] that W := {Wt : t ≥ 1} is a
d-dimensional martingale-difference sequence (and hence uncorrelated). The i.i.d. sampling model assumed in the present
paper implies that V = {Vt : t ≥ 1} is a martingale-difference
sequence that is also uncorrelated with W .
Prop. 2.1 contains full details of this system description.
The Kalman filter is developed in two settings, each with
the same observation process: The first is constructed to obtain
estimates of µN
t . The second filter obtains estimates of the joint
statistics of a larger state that includes both µN
t and the state
of a typical individual. The main conclusions are summarized
here:
(i) A measurement architecture is proposed in which each
load broadcasts its state only occasionally – say, once per
day. The observation equations in the aggregate model
then include white noise, whose conditional variance is
computed. The state equations for the population/individual
dynamics evolve as a linear stochastic system with white
noise disturbance, whose conditional covariance matrix is
also computable.
(ii) In the examples considered, the observability Grammian
is not full rank, and an approximate time-invariant model is
also unobservable. Moreover, Prop. 4.1 demonstrates that
a symmetry property (that holds in the example of [6]),
implies that the model cannot be observable.
However, in numerical experiments it is found that the
Kalman filter remains valuable for reducing the impact
of measurement noise, and for estimating the distribution
of quality of service. In particular, estimates of certain
first and second order statistics of an individual load are
remarkably accurate, even though the measurements are a
noisy sequence of samples from the population.
(iii) In the face of un-modeled dynamics such as load heterogeneity, or additional “opt-out” control used to enforce
QoS bounds, the Kalman filter combined with PI control
continues to perform nearly perfectly, even with 0.1%
sampling of loads.
There are in fact two general formulations of the Kalman
filter. In the first, most typical setting, the sequence of Kalman
gains is deterministic; obtained through a Riccati equation (a
recursive equation driven by the covariance matrices for the

state and observation noise). This is known to be L2 -optimal
over all estimators that are linear functions of the observations.
A second formulation of the Kalman filter uses conditional
covariance matrices to define the Kalman gain — see (22,23)
and surrounding discussion. If the state/observation noise
is conditionally Gaussian, then this Kalman filter coincides
with the nonlinear filter, which is L2 -optimal over all causal
estimators [12]. Because the Riccati equation is a nonlinear
function of the covariance matrices, this version of the Kalman
filter may be a nonlinear function of the observations.
The second is attractive because it is easy to compute
formulae for the conditional covariance matrices, while the
unconditional covariance matrices only admit approximations.
Moreover, when considering the dynamics of the aggregate, a
Gaussian approximation of the noise is justified by the Central
Limit Theorem (CLT).
Related research: In addition to the references cited
above, there are many papers on demand dispatch based on
centralized control, or relying on real-time prices to solve the
control problem of interest. Much of the latter is closer to
demand response, and has little intersection with the research
summarized here.
An application of this formulation of the Kalman filter
was considered previously in [13] for a single Markov chain
without control, with measurements subject to Gaussian error.
There are several recent papers with similar goals in the literature on mean-field models. Most closely related is [14], [15]
which concerns partially observed LQG mean field games,
with several classes of players. The state estimation problem is
from the point of view of the individual – each “minor agent”
obtains noisy and partial observations, and wishes to estimate
the “major state” as well as the aggregate. The solution is
obtained through the construction of a Kalman filter. This prior
work is also motivated by application to power systems.
The remainder of the paper consists of four sections
organized as follows. The following section describes the
stochastic model on which the estimation algorithms are based.
Filtering equations are derived in Section III for a collection of
controlled Markov models with partial observations. The algorithms have been tested in various different settings – results
for demand dispatch using residential pools and also TCLs
are summarized in Section IV. Conclusions and directions for
future research are contained in Section V.
Further results may be found in the dissertation [16].
II. M EAN F IELD M ODEL
It is assumed throughout the paper that a family of Markov
transition matrices {Pζ : ζ ∈ R} is given that is continuous
in the parameter ζ. The finite state space is denoted X =
{x1 , . . . , xd }, so that each Pζ is a d × d matrix.
The mean-field model is defined as the approximation of
(2) obtained as N → ∞. This is the deterministic recursion,
µt+1 = µt Pζt ,

(4)

with µ0 given, and where ζ is obtained via causal feedback.
This paper is concerned with the stochastic system (2), but
the steps used to justify the limit will lead to the second-order
statistics required to describe the Kalman filter.

3

A. Aggregate dynamics
The individual dynamics are described by the controlled
Markov model (1). We lift the state space from the d-element
set X = {x1 , · · · , xd }, to the d-dimensional simplex S. For
the ith load at time t, the element Γit ∈ S is the degenerate
distribution whose mass is concentrated at x if Xti = x; that
is, Γit = δx . These distributions evolve according to a random
linear system,
Γit+1 = Γit Git+1
(5)
in which Γit is interpreted as a d-dimensionalProw vector, Git is
a d × d matrix with entries 0 or 1 only, and l Git (xj , xl ) = 1
for all j.
The following assumptions are imposed throughout.
A1: The input is defined by causal output feedback: For a
continuous family of functions φt : Rt+1 → R,
ζt = φt (Y0 , . . . , Yt ) ,

Git = Ξ(ζt−1 , ξti ),

(6)
d

where
: t ≥ 1, i ≥ 1} are i.i.d. on [0, 1] with uniform
distribution.
A3: The initial conditions {Γi0 : 1 ≤ i ≤ N } are i.i.d., with
P{Γi0 (x) = 1} = µ0 (x),

k
Φt (k) = µN
t (x ),

1 ≤ k ≤ d, t ≥ 0.

Φit (k) = Γit (xk ),

1 ≤ k ≤ d, t ≥ 0.

Φt =

N
1 X i
Φ
N i=1 t

Proposition 2.1: The two state processes each evolve as
linear systems,
Φt+1

= At Φt + Wt+1

(12)

Φit+1

i
= At Φit + Wt+1

(13)

where At = PζTt , and for each i, t,
Wt+1
i
Wt+1

N
1 X k
Wt+1
N
k=1

T
= Γit (Git+1 − Pζt )

=

The first assumption is based on the control architecture
assumed in all prior work. Assumption A2 is standard in
this context [13], and can be assumed since any probability
mass function (pmf) is a function of a uniformly distributed
random variable.
Assumption A4 can be used to model random sampling with
or without replacement. Examples are given in Section IV.
This assumption leads to the measurement equation (3); additional additive noise can be included in the observation
equation, provided this measurement noise is independent of
the sampling process.
The filtration of observations is denoted,
t ≥ 0.

(7)

Under A1–A3, Git+1 is conditionally independent of
{Γi0 , · · · , Γit }, given Yt , with
E[Git+1 | Yt ] = Pζt .

(8)

(14)
(15)

The observation equation (3) can be written,
(16)

where Cj = U (xj ) for each j.
Moreover, each of the sequences W i , W , V is a martingale
difference sequence, and the sequence V is also uncorrelated
with each of the {W i } and also W .
Proof: The martingale difference property for each of the
{W i } (and hence also W ), follows immediately from (8); see
also [3]. The sequence V can be expressed,

i=1

Yt = σ{Yr , ζr : r ≤ t},

(11)

Yt = CΦt + Vt

Moreover, the distribution of γ is unchanged by permutations of its N components.

(10)

The two state processes are evidently related,

x∈X

A4: The measurements are obtained through sampling:
There is a bounded sequence {γt : t ≥ 1} of N dimensional vectors with non-negative entries, and independent of the {ξti }, such that E[γt (k)] = N −1 for each t
and k, and
N
X
Yt =
γt (i)U (Xti ).

(9)

In the second, the goal is to estimate the state for an individual
load. For the ith load, this is denoted

t ≥ 0.

A2: For some function Ξ with domain R × [0, 1], and range
equal to the set of d × d matrices,
{ξti

Two filtering problems are considered in this paper. In the
first, the state is equal to the empirical distributions expressed
as column vectors:

Vt = −

X

k
k
µN
t (x )U (x )

k

+

N
X

γt (i)U (Xti )

i=1

This has mean zero, from the exchangeability assumption in
A4, and it is a martingale difference sequence because the
sequence γ is i.i.d. It is uncorrelated with {W i } and W under
the assumption that {γt } is independent of {ξti }.
A derivation of the conditional state covariances is given in
Section III-A.
For a linear-Gaussian model, the Kalman filter equations are
intended to approximate the conditional mean and covariance
of the state. In the first model (12) they are denoted,
b t = E[Φt | Yt ] ,
Φ

e tΦ
e T | Yt ]
Σt = E[Φ
t

(17)

e t = Φt − Φ
b t . For large N it can be argued via the
with Φ
CLT that (V , W ) is approximately conditionally Gaussian
given the observations. We might expect the Kalman filter to
approximate the optimal nonlinear filter in this case.
In Section III-C the state space is extended to obtain
estimates of QoS metrics for an individual load that may not
be a function of the respective state Φit .

4
10

B. Example: Intelligent pools
Here we focus on a single example in which each load is a
residential pool pump. Section IV contains extensions to other
loads.
In the original model of [6], the state space is taken to be
the finite set,
X = {(m, k) : m ∈ {⊕, }, k ∈ {1, . . . , I}}

(18)

where I > 1 is an integer. If Xti = (⊕, k), this means that the
pool pump is on at time t, and has remained on for the past
k time units. In this paper we take the same state space, but
with a new interpretation of each state.
1

...

2

On

...

I −1

I

2

1

...

I
Fig. 1.

...

I −1

State transition diagram for pool model.

The controlled transition matrix is of the form,
Pζ = (1 − δ)I + δ P̌ζ

(19)

in which P̌ζ is the transition matrix used in [6], and δ ∈ (0, 1).
At each time t, a weighted coin is flipped with probability of
heads equal to δ. If the outcome is a tail, then the state does
not change. Otherwise, a transition is made from the current
state x to a new state x+ with probability P̌ζt (x, x+ ).
A state transition diagram is shown in Fig. 1. The state
transition diagram for P̌ζ is identical, except that the self-loops
are absent.
Output deviation

Reference

3
0
-3

20

40

60

80

100
t/hour

120

140

160

Input ζt

0

-10
0

5

96 Eigenvalues of the
Observability Grammian

10 0

10 -5

λζ

λ0

10 -10
1

24

Fig. 2. The deviation in power consumption tracks well even with only 20
pools engaged.

The motivation comes from conflicting needs of the grid
and the load: a single load turns on or off only a few times
per day, yet the grid operator wishes to send a signal far more
frequently – In this example we assume every 5 minutes. If the
sampling increments for each load were taken to be 5 minutes,
then it would be necessary to take I very large in the approach
of [6].
In this paper P̌ζ is obtained using the optimal-control
approach of [6]; we take I = 48, and hence d = |X| = 96. It
is assumed that δ = 1/6, so that the pool state changes every
30 minutes on average. In [6] it is shown that the transition
matrix has desirable properties for control: the linearized

48

72

96

i

Fig. 3. Eigenvalues for the observability Grammian for the pool model in two
cases: The magnitude of eigenvalues decays rapidly for a typical sample-path
of ζ, and for the LTI model obtained with ζ ≡ 0.

dynamics are minimum phase, with positive DC gain. Hence,
for example, a persistent positive value of ζt will lead to an
increase in aggregate power consumption.
For sake of illustration, Fig. 2 shows tracking performance
of this scheme with only twenty pools. Each pool is assumed
to consume 1 kW when operating, and each has a 12 hour
cleaning cycle. The grid operator uses PI compensation to
determine ζ (see Section IV-B for details). With such a small
number of loads it is not surprising to see some evidence of
quantization. For 100 loads or more, and the reference scaled
proportionately, the tracking is nearly perfect.
Two QoS metrics have been considered for this model. First
is ‘chattering’ – a large number of switches from on to off. A
large value means poor QoS, but this is already addressed
through design of the controlled transition matrix [6]. The
design of Pζ also helps to enforce upper and lower bounds
on the duration of cleaning each day.
A second metric is total cleaning over a time horizon of one
week or more. This is the QoS metric considered in [3]. In
this paper we consider a discounted version: We assume that
P0 has a unique invariant pmf π0 . With
P ` : X → R a given
function with zero steady-state mean, x π0 (x)`(x) = 0, we
define for each i and t,
Lit =

kW

10

|λi |

t
X

β t−k `(Xki )

(20)

k=0

with β ∈ (0, 1] a constant. The function `(x) = I(m = ⊕) −
I(m = ) was used in [3] in the case of a 12 hour cleaning
cycle (recall the notation x = (m, k)).
Experimental results surveyed in Section IV demonstrate
that it is possible to estimate functionals of the state process
such as the QoS metric {Lit }, even if the observations are
subject to significant measurement noise.
In anticipation of the results to come we ask, what would
linear systems theory predict with respect to state estimation
performance? The observability Grammian associated with
(12) was computed for typical sample paths {At = PζTt :
1 ≤ t ≤ 2016}, where the value 2016 corresponds to one
week, and ζ was scaled to lie between ± 12 . Its rank was
found to be approximately 40, while the maximal rank is 96
(the dimension of the state). With ζt ≡ 0 the system is timeinvariant. In this case the rank of the observability Grammian
coincides with the rank of the observability matrix, which was
found to be 23.
However, these values were obtained using the “rank” command in Matlab, which relies on finite numerical precision.
A plot of the magnitude of the eigenvalues for the two

5

observability Grammians shown in Fig. 3 suggests that both
matrices are full rank.
Further analysis establishes that the LTI model obtained
with ζ ≡ 0 cannot be observable, due to a particular symmetry
found in this example. A general result given in Prop. 4.1
implies that λ0i = 0 for 50 ≤ i ≤ 96.
III. K ALMAN F ILTER E QUATIONS
The second order statistics for the disturbances appearing
in the linear model (12, 13, 16) are derived here. These
expressions are used to construct a Kalman filter that generates
approximations for the conditional mean and covariance (17).
Other statistics of interest are,
b it = E[Φit |Yt ] , Σit = E[Φ
e it (Φ
e it )T |Yt ]
Φ
i
bi
Φ
t+1|t = E[Φt+1 |Yt ] ,

e it+1 (Φ
e it+1 )T |Yt ]
Σit+1|t = E[Φ

b t+1|t = E[Φt+1 |Yt ] ,
Φ

e t+1 (Φ
e t+1 )T |Yt ]
Σt+1|t = E[Φ

b it .
e it = Φit − Φ
where again tildes represent deviations, such as Φ
Prop. 3.1 states that some statistics of the individual can be
expressed in terms of those of the population: It is not the
case that Σit = N Σt , since {Φit : 1 ≤ i ≤ N } are correlated.
Proposition 3.1: For each s, t, i, and any set S ⊂ Rd , the
conditional probability is independent of i:
P{Φis ∈ S | Yt } = P{Φ1s ∈ S | Yt }

(21)

and consequently, E[Φis | Yt ] = E[Φs | Yt ]. Moreover,
bi
b
Φ
t+1|t = At Φt , the state covariances for the individual are
Σit
Σit+1|t

b t) − Φ
b tΦ
b Tt
= diag (Φ
bT
b t+1|t ) − Φ
b t+1|t Φ
= diag (Φ
t+1|t ,

t+1|t

Proof: The proof of (21) follows from the symmetry and
independence conditions imposed in (A1–A4). The remaining
results follow from this, and the fact that Φit has binary entries
[in particular, Φit (Φit )T = diag (Φit )].
Recall from the introduction that two formulations of the
Kalman filter have been considered in this research. For a
conditionally Gaussian model, the Kalman filter equations
require the conditional covariances for the state noise,
i

i
i
T
ΣW
= E[Wt+1
(Wt+1
)T | Yt ] , ΣW
t
t = E[Wt+1 Wt+1 | Yt ]
(22)
and also the conditional covariance of the measurement noise,

(23)

Formulae for the state noise covariances can be obtained in full
generality. We require the distribution of the random vector γt
introduced in A4 to obtain a formula for ΣVt .
The Kalman filter that generates L2 -optimal estimates over
all linear functions of the observations uses instead the (unconditional) covariance matrices,
Wi

Σt

i
i
= E[Wt+1
(Wt+1
)T ],

W

T
Σt = E[Wt+1 Wt+1
],

A. State noise covariance
The following result provides formulae for the conditional
covariances for the state noise (22) as a function of the
b t.
conditional mean Φ
Proposition 3.2: Under Assumptions A1–A4,

1
T
b
b
diag
(A
Φ
)
−
A
diag
(
Φ
)A
(25)
ΣW
=
t
t
t
t
t
t
N
i
b i ) − At diag (Φ
b i )AT
ΣW
= diag (At Φ
(26)
t
t
t
t
The second covariance is independent of i, with common value
i
ΣW
= N ΣW
t
t .
Proof: Since {Wti : 1 ≤ i ≤ N } is uncorrelated,
ΣW
t =

N
1 X Wi
Σ
N 2 i=1 t

(27)

bi = Φ
b t , and from this or (11) it
Moreover, Prop. 3.1 gives Φ
t
is obvious that
N
X
bi
bt = 1
Φ
Φ
N i=1 t

and the cross covariances can be expressed,

 i
e t (Φ
e t )T | Yt = Σt
E Φ

 i
e
e t+1|t )T | Yt = Σt+1|t , 1 ≤ i ≤ N.
(Φ
E Φ

ΣVt = E[Vt2 | Yt−1 ]

V

and Σt = E[Vt2 ] (the notation used in standard texts is
W
V
Qt = Σt and Rt = Σt [12]). We show in Prop. 3.2 that
the two covariance matrices in (22) are linear functions of the
b t . Expressions for the two covariance
true conditional mean Φ
matrices in (24) follow from Prop. 3.2 and the smoothing
property of conditional expectation, provided we can compute
b t ] = E[Φt ]. The formula we obtain for ΣV in Prop. 3.3 is
E[Φ
t
a linear function of the conditional covariance matrices Σit+1|t
and Σt+1|t . It is unlikely we can obtain formula for the means
of these covariance matrices, and hence we do not expect to
V
obtain an exact formula for Σt .

(24)

Consequently, (25) follows from (26).
i
is similar to
The derivation of the formula (26) for ΣW
t
the Kalman filter construction in [13]. Given the larger sigmafield,
Yt+ = σ{Φir , Yr , ζr , Ar : r ≤ t, i ≤ N }
the smoothing property of conditional expectation implies,


i
i
i
ΣW
= E E[Wt+1
(Wt+1
)T | Yt+ ] | Yt
t
The inner conditional expectation is transformed using the
definition (15):
i
i
E[Wt+1
(Wt+1
)T | Yt+ ]

= E[Φit+1 (Φit+1 )T | Yt+ ]

− E[Φit+1 | Yt+ ]E[Φit+1 | Yt+ ]T

= diag (At Φit )

− At diag (Φit )ATt
where the final equation uses E[Φit+1 | Yt+ ] = At Φit , and the
fact that Φir has binary entries for each i and r.
Taking the conditional expectation given Yt gives (26).

6

B. Sampling and observation covariance
The observation model used in the numerical experiments
that follow is based on random sampling of loads: An integer
n < N is held fixed, and at each time instant t a distinct set
of n indices {k1 , . . . , kn } is chosen uniformly at random. The
observation is the average

Using familiar arguments, and applying (31),


e1
e2
(N − 1)E {C Φ
t+1|t }{C Φt+1|t } | Yt
=

k=2



e1
e
= N E {C Φ
t+1|t }{C Φt+1|t } | Yt


e1
}2 | Yt
− E {C Φ

n

Yt =

1X
U (Xtki ).
n i=1

(28)

Recalling the definition U (Xtki ) = CΦkt i , it follows from
Prop. 3.1 that this is an unbiased estimate of CΦt . An
expression for the conditional variance ΣVt defined in (23)
is obtained next.
Proposition 3.3: The conditional covariance is given by,

1 N −n  i
C Σt+1|t − Σt+1|t C T
ΣVt+1 =
n N −1
Proof: By definition,
n
2
i
h 1 X
V
i
CΦkt+1
− CΦt+1 | Yt
Σt+1 = E
n i=1
where the permutation is random, uniform, and independent
of Yt . It follows from symmetry of the model that we can
consider just the first n samples, replacing ki by i. We also
center the random variables to obtain
n
h 1 X
2
i
ei
e
| Yt
ΣVt+1 = E
CΦ
(29)
t+1|t − C Φt+1|t
n i=1
The quadratic is expanded as the sum of three terms,
ΣVt+1 = S1 − 2S2 + S3

(30)

in which the second two terms are transformed using the
following consequence of Prop. 3.1:


T
ei
e
E {C Φ
(31)
t+1|t }{C Φt+1|t } | Yt = CΣt+1|t C
The first term is the conditional variance,



e t+1|t 2 | Yt = CΣt+1|t C T
S1 := E C Φ
The second term is transformed using (31):
n
h 1 X


i
ei
e t+1|t | Yt
S2 := E
CΦ
CΦ
t+1|t
n i=1

(32)

(33)

The third term is another conditional variance,
n
2
i
1 hX e i
S3 := 2 E
C Φt+1|t | Yt
n
i=1
n

1 X  ei
ej
= 2
E {C Φt+1|t }{C Φ
} | Yt
t+1|t
n i,j=1

n
X


ei
ej
E {C Φ
t+1|t }{C Φt+1|t } | Yt

i,j=1



2
e1
= nE {C Φ
t+1|t } | Yt


e1
e2
+ (n2 − n)E {C Φ
t+1|t }{C Φt+1|t } | Yt

t+1|t

= N CΣt+1|t C T − CΣit+1|t C T

Putting these expressions together gives,
i
n2 − n h
1 
N Σt+1|t − Σit+1|t C T
S3 = 2 C nΣit+1|t +
n
N −1
The desired expression is obtained by substitution of these
terms in (30).
C. Estimation of the individual
A filter is introduced here to estimate the mean and variance
for an individual state Φi , and the mean and variance of the
QoS metric Lit for a typical load (recall (20)). Estimates of the
individual state and individual QoS can be used to estimate the
flexibility of loads, which may vary with time. For example,
if the grid operator believes that every water heater contains
water that is too cold, then it is unlikely that there is much
remaining flexibility to reduce power consumption from these
loads.
The Kalman filter must be modified to obtain estimates of
the QoS for individual loads. One approach is to restrict to
β = 1, and use the Markovian model (Xti , Lit ). This comes
with high complexity: In the case of the pool model, the first
component can take on d = 96 values, and the range of values
of Lit may be over 100. This means that the state space is
at least 104 , which is probably too large to be practical for
estimation.
Here we introduce a Kalman filter for the joint process Ψit =
(Φit ; Φt ; Lit ), which is of dimension 2d + 1. The construction
of a linear model for Ψi is based on (12, 13), and the onedimensional dynamics for QoS,
Lit+1 = βLit + C` Φit

= CΣt+1|t C T

Applying symmetry once more gives,

N
X


e1
ek
E {C Φ
t+1|t }{C Φt+1|t } | Yt

where C` Φit = `(Xti ). The measurements remain of the form
(3), which is why it is necessary to include Φt in the definition
of Ψit . The new A, C, state-noise covariance matrices are,

At
0
C`

0
At
0


0
0 ,
β


[0, C, 0] ,

i

ΣW
t
N −1 ΣW i
t
0

N −1 ΣW
t
ΣW
t
0

i


0
0
0

We thus have all of the system parameters needed to construct
the Kalman filter.
IV. E STIMATING Q O S AND C LOSING THE L OOP
We have conducted experiments in various settings, using
both formulations of the Kalman filter, differentiated by the
use of conditional or unconditional covariance matrices. For
estimating performance we obtained excellent results with the
Kalman filter whose gain is based on the conditional covariance matrices. The numerical results focus on the residential

7

Mean

50

Estimate

Empirical

0

Lt

−50
−100

Variance

x10 3 6

ΣL
t

4
2
0
0

100

200

300

400

500

600

700

t/hours

Fig. 4. Although the liner model is not observable, estimates of the mean
and variance of QoS for an individual are nearly perfect.

Kalman filter by doubling the dimension of the state space.
Instead, an approximate model was constructed in which the
parameters in the two state space models (differentiated by
cleaning cycle) were averaged; similar approximate models
are used in [2].
Estimation results for the case of two pool classes are shown
in Fig. 5. The QoS for the two classes of pools have different
means and variances. This difference cannot be captured using
this filter which is designed to estimate the overall mean and
variance. The filter generates estimates of the mean QoS that
lie between the two empirical means, and estimates of its
variance follow more closely the empirical variance of the
pools with 12hr/day cleaning cycle.

Mean

50

Empirical (12hr)

Empirical (8hr)

Estimate

0

Lt

−50
−100
x10 3 6

Variance

pool model, and closes with preliminary extension to TCL
models.
The Kalman filter that is optimal over all linear estimators
requires unconditional covariances. The state covariance matrix can be expressed as the mean of (25):
i
1 h
W
b t ) − At diag (Φ
b t )ATt
Σt = E diag (At Φ
N
b t ] = E[Φt ] does not have a closed form
The mean E[Φ
expression, so we make two approximations. First, we consider
the mean with ζ ≡ 0, and second we let t → ∞. The
covariance used in the filter is the resulting limit,
i
1h
W
(34)
Σ∞ =
diag (At π0 ) − At diag (π0 )ATt
N
where π0 is invariant for P0 . A similar approximation is used
V
for Σt .
The common features in all of the numerical experiments
that follow are listed here: The reference signal r was generated from the BPA balancing reserves deployed in the first
month of 2015 [17]. This signal was low pass filtered, and
then scaled to an appropriate magnitude to match the capacity
of the aggregate of loads, exactly as in [6].
The observation model was based on the sampling assumptions introduced in Section III-B, in which a fixed number of
loads n < N is sampled at each time. Recall from (28) that
the average of these n samples at time t is denoted Yt . Unless
stated otherwise, n is taken to be 0.1% of the total population.
Hence, for a population of N = 10, 000 pools, exactly 10 are
sampled at each time step.
It is assumed that each pool pump consumes 1kW while in
operation. The actual power consumption at time t is N CΦt .
b t are scaled by N
For this reason, plots of CΦt or Ybt = C Φ
to represent total power consumption or its estimate.
Finally, the discount factor β in the QoS metric (20) was
chosen as the value for which β k ≈ 1/2 when k corresponds
to one week. The value β = 0.9997 is obtained under the
assumption of five-minute sampling.

ΣL
t

4
2
0
0

100

200

300

400

500

600

700

t/hours

Fig. 5. Estimates of the mean and variance of QoS for an individual with a
heterogenous population. The un-modeled dynamics lead to some error.

A. Estimation of QoS
The Kalman filter described in Section III-C was used to
obtain estimates of the mean and variance of the QoS metric
Lit , with `(x) = Cx − y. These experiments were performed
with ζ equal to the exogenous “mean-field limit” used in the
linearized model of [3]:
ζt =

Gc
rt
1 + Gc Gp

where Gc is a PI compensator, and Gp is a transfer function
for a linearized model (see [3] for details).
Fig. 4 shows estimates of QoS based on a model with 10,000
pools, with 0.1% sampling rate. The accuracy of estimation is
remarkable, in spite of significant measurement noise.
Non-ideal settings were also considered. In one experiment
a population of 20,000 pools was divided into two classes
of equal size: 10,000 pools operated as previously with a
12hr/day cleaning cycle, and the other 10,000 operated with
a 8hr/day cleaning cycle. It is possible to construct the exact

B. Control of the heterogeneous population
Closed-loop experiments were conducted in which unmodeled dynamics are present due to both load heterogeneity
and additional local control. With large un-modeled dynamics
it is still possible to obtain estimates of the mean of QoS, but
the Kalman filter cannot be expected to provide estimates of
the variance.
Error feedback was used in all of the remaining experiments,
ζt = Gc et . A PI compensator was used to define Gc : A
proportional gain of 50, and integral gain of 1.5 worked well
in all examples. These values were based on a nominal LTI
model — this feedback design resulted in a phase margin of
about 70 degrees, and gain margin of about 10dB.
In the experimental results illustrated in Fig. 6, the error
signal was defined with respect to the raw measurements,

8

Cycle (hr/day)

20

18

16

14

12

10

8

6

4

Number (×104 )

1

1

1

3

5

10

5

3

1

Although the simulation uses these nine classes of pools,
for the purposes of filter design an approximate model was
constructed based on just two pool models, corresponding
to 8 and 12-hour cleaning cycles. Following [2], a convex
combination was chosen for the state matrix,
At =

αA8t

+ (1 −

α)A12
t

=

∗
ΣW
t

+

W
k t Σ∞ ,

Reference

50
0

3
0
-3

−100
0

Fig. 6.

20

40

60

80

t/hour

100

120

140

160

(36)

∗
where kt ≥ 0, ΣW
is a convex combination of the matrices
t
W
computed in Prop. 3.2, and Σ∞ is a convex combination of
the matrices (34). This structure ensures that ΣW
t is positive
semi-definite, and that ΣW
1
=
0
(consistent
with
the fact that
t
our state evolves in the simplex). The performance of the filter
was not very sensitive to the parameter kt , but the best value

Input ζt

−50

PI control with output-error feedback, with et = N −1 rt − Yt .

With the experimental setting fully described, we now
survey results from several control experiments.
We first show what happens if we forgo use of the Kalman
filter, and take et = N −1 rt − Yet , with Yet = Yt − y, and Yt
the noisy measurements obtained with 0.1% sampling. Results
from one experiment are shown in Fig. 6. The volatility of
the measurements resulted in similar volatility of input ζt .
This drove the pools to switch frequently, and also resulted
in degraded grid level tracking. This motivates the use of a
filter to smooth the measurements.
Results using the filtered error (35) are shown in Fig. 7.
This Kalman filter estimator based control gives nearly perfect
tracking performance, in spite of the significant modeling
error.
100

Output deviation

Reference

50
0

−50
3
0
-3

−100

in which 0 ≤ α ≤ 1 is independent of t, and the superscripts
8
and 12 refer to the respective pool classes.
The parameter α was chosen to be consistent with the hypothesis that the collection of pools consist of just two classes,
with 8 or 12 hr/day cycles. It is assumed that the nominal
average power consumption y is known to the grid operator –
this can be estimated easily from weekly measurements. The
2-class assumption would imply that y = αy 8 + (1 − α)y 12 ,
giving α = (y 12 − y)/(y 12 − y 8 ), with y 12 = 1/2 and
y 8 = 1/3.
The covariance matrix ΣW
t used in the Kalman filter was
modified to take into account the large un-modeled dynamics.
It was chosen as the sum of two terms,
ΣW
t

Output deviation

0

Fig. 7.

20

40

60

80

t/hour

100

120

140

160

Input ζt

TABLE I: D ISTRIBUTION OF POOLS

100

MW

1
bt − y
(35)
rt − Yet , Yet = C Φ
N
with y equal to the average nominal power consumption. The
b t is an
scaling of the reference by N −1 is required since C Φ
estimate of average power consumption.
The total number of pools was taken to be N = 300, 000.
Performance results using smaller values of N are discussed
at the close of Section IV-B, and also in Section IV-C.
The simulation model used in the remaining experiments
consisted of nine different classes of pools, distinguished
by cleaning cycle. The distribution of classes is shown in
Table I. Suppose that the grid operator had full information
regarding the number of pools in each class. The full state
space description would have dimension 9 × d, which is far
too complex. Moreover, in practice the grid operator will not
have complete system information.
Here we make some coarse approximations to simplify the
model. In addition to reducing complexity of the filter, our
goal is to investigate robustness of the estimation algorithms
in a closed-loop setting.
et =

was found to grow linearly with population size. The results
that follow used kt = N/300 (note that N cancels the factor
W
1/N appearing in the definition (34) of Σ∞ ).

MW

without use of the Kalman filter. In all other experiments, the
following definition was used:

PI control with smoothed error signal (35).

Tracking experiments were conducted in many other scenarios, distinguished by sampling rate and load population,
using the same distribution of load classes given in Table
1. Normalized tracking error was chosen as the performance
metric,
! 12
T
1X 2
yt − rt
n
, krk2 =
r
.
et =
krk2
T t=1 t
The full relationship between population size N , sampling
rate, and closed loop performance is illustrated in Fig. 8. The
vertical axis is the root mean square (RMS) of the tracking
error {ent }. For each of the 4 values of N considered, the
tracking error is less than 3.5% with 100% measurements. It is
found that the tracking error is best predicted by the number of
samples per time-step. Consequently, a smaller population of
loads requires a larger sampling rate to maintain good tracking
performance.

9

10 N = 30,000
5

Output deviation

50
MW

0

−5

Reference

0

opt out %

−50

0.5

−10
100

110

0
t/hour 130

120

0.5

−100
100

opt out %

MW

Closed-loop tracking

100 N = 300,000

110

120

0
t/hour 130

PI control and local opt-out control for the heterogeneous model with 0.1% sampling rate. Tracking performance improves as N is increased.

Fig. 9.

18

RMS Normalized Error (%)

is the d × d matrix whose kth row is equal to CAk−1 , and its
null space Σō is called the unobservable subspace. Suppose
that {Φt } and {Φ0t } are two solutions to the state space
equations with common input sequence {ut }, but different
initial conditions, Φ0 , Φ00 ∈ Rd . If Φ00 − Φ0 ∈ Σō , then for
each t ≥ 0,

Sampling Rate

16

0.1%
1.0%
10%
100%

14
12
10
8
6
4
2
0

3 × 10 3

3 × 10 4

3 × 10 5

3 × 10 6

N

Fig. 8. Tracking performance improves with increased sampling rate, or a
larger population of loads N .

C. Opt-out local control
In [3] an additional layer of control is introduced, based on
the QoS metric Lit introduced in (20): The ith load can opt-out
of service to the grid (ignore the signal ζt ) at any time t for
which Lit lies outside of pre-assigned bounds. The impact of
these additional un-modeled dynamics is investigated here.
At each load the QoS metric Lit is constrained to a common
interval [−50, +50]: if Lit ≥ 50, then the pool is turned off at
time t, regardless of the value of ζt . Similarly, if Lit ≤ −50
then it will be turned on.
In the construction of the Kalman filter, the state noise
covariance was taken as (36), but in this case the scaling factor
was chosen as a function of the QoS estimates, kt = k +b|Lbit |.
This is chosen to reflect the fact that un-modeled dynamics
increase with the percentage of pools that opt-out. The values
k = N/100 and b = 300 worked well, and the sensitivity to
these values was not high.
Closed loop performance remained nearly perfect, and the
percentage of opt out loads was observed to remain under
0.5% over the entire run. Fig. 9 shows results from two
experiments with N = 300, 000 and N = 30, 000, each with
sampling rate maintained at 0.1%.
It is seen again that reducing the number of loads results in
degradation of closed loop performance, since fewer loads are
sampled at each time-step; the two plots are nearly identical if
the sampling rate is increased to 1% for the smaller population.
D. Observability in Demand Dispatch systems
It is shown here that the lack of observability observed in
the pool model with 12 hr cleaning cycle can be attributed to
symmetry of the model.
A review of linear systems terminology is required. Consider the general LTI state-space model,
Φt+1 = AΦt + But ,

yt = CΦt

where A is a d×d matrix, B is a d-dimensional column vector,
and C is a d-dimensional row vector. The observability matrix

yt − yt0 = CAt (Φ0 − Φ00 ) = 0
That is, based on measurements {yt : t ≥ 0}, it is impossible
to know if the true state sequence is {Φt } or {Φ0t }.
For the time-varying system
Φt+1 = At Φt + But ,

yt = CΦt

the observability matrix is replaced by the observability Grammian, OG . The definition can be found in any introductory
text on state space control (or, see the encyclopedic text [12]).
A critical interpretation is this: If Φ0 , Φ00 are two initial
conditions, then with identical inputs
(Φ0 − Φ00 )T OG (Φ0 − Φ00 ) =

∞
X
t=0

(yt − yt0 )2

This is why the small eigenvalues observed in Fig. 3 are a
concern if we wish to reconstruct the entire state based on the
observations {yt }.
The proposition that follows concerns a class of LTI models
that have the following symmetry property:


Ag Ao
T
A=
,
C = [1 | 0T ]
(37)
Ao Ag
where A is a d × d matrix, and the decomposition is in terms
of two square matrices Ag , Ao . This requires that the integer
d is even.
An example is the pool pump model with 12 hour cleaning
cycle. Under certain conditions, a TCL model may be symmetric following a state transformation – an example is given
below. Regardless of the details of the model, symmetry rules
out observability.
Analysis of the symmetric model is based on the subspaces,
V = {v ∈ Rd : vi = vd/2+i , 1 ≤ i ≤ d/2}
T

V0 = {v ∈ V : 1 v = 0}

Proposition 4.1: Consider the symmetric LTI model (37).
If AT is a transition matrix, then the subspace V0 is contained
in the unobservable subspace. Consequently, the symmetric
model is never observable, and the rank of the observability
matrix is no greater than d/2 + 1.

10

T T

where w = (Ag + Ao )z. This shows that A : V →PV.
To complete
the proof, we next show that
wi = 0
P
whenever
zi = 0, so that A : V0 → V0 .
Since AT is a transition matrix, it follows that


1T = 1T A = 1T Ag + 1T Ao | 1T Ao + 1T Ag
T
That is, Q = (Ag + AP
o ) is a transition matrix: Q1 =
T
Consequently, if z 1 = i zi = 0, then

1.

1T w = 1T (Ag + Ao )z = z T Q1 = 0
Consider the following example in which the dimension of
the symmetric state space model is d = 8,


1−δ
0
0
0
 a21 1 − δ
0
0 

Ag = 
 0
a32 1 − δ
0 
0
0
a43 1 − δ


a15 a16 a17 δ
 0
0
0 0

and A0 = 
 0
0
0 0
0
0
0 0
Based on the proposition, the rank of the observability matrix
is no greater than 5.
The following values are consistent with the pool pump
model with 12 hour cycle considered previously: δ = 1/6,
a21 = 0.1654, a32 = 0.0840, a43 = 0.0026, a15 = 0.0013,
a16 = 0.0827, a17 = 0.1641. The 8×8 observability matrix O
was obtained using Matlab, and the null command results in
equality Σō = V0 . However, Matlab also returns these values
for the log of the modulus of the eigenvalues: {log |λi |} =
{1.4, 0.095, −1.8, −5.8, −19.7, −35.5, −36.5, −37.8}
Based on these values, it is not clear if the rank of O is equal
to 4, 5, or 8.
Temperature/ o C

Temp state Z

21
20.5
20
19.5

Power consumption / kW

15
10
5

0
0

0.2

0.4

0.6

0.8

1
1.2
Time/hour

1.4

1.6

1.8

2

Fig. 10. State transformation for the TCL model: the evolution of Z is
similar to the pool model.

A typical trajectory for an air-conditioning unit shown in
Fig. 10 reveals that the sample paths are very similar to the
pool model. If the Markovian dynamics of Z are symmetric,
then the LTI model is not observable.
E. Eigenvectors and reduced order observer
The structure of eigenvectors and reduced order observers
are described for both the pool model and a TCL model.
Observability of the pool model: We return to the homogeneous model with 12hr/day cleaning cycle. The observability Grammian for the time-varying linear system was found
to be highly ill-conditioned in all of our experiments – one
example is illustrated in Fig. 3. The LTI model obtained with
ζ ≡ 0 is not observable, by an application of Prop. 4.1.
v1

1
Real part

Av = (w | w )
T

Consider next a TCL that provides cooling (a refrigerator or
air-conditioner), with temperature θt constrained to the deadband [Θmin , Θmax ]. A new state process (Zt , mt ) is obtained
in which mt is the power state as before, and
(
θt
if mt = 0
Zt =
Θmin + Θmax − θt
if mt = 1

v2

v3

v4

0.5
0
−0.5
−1
1

Imaginary part

Proof: To prove the proposition we establish that
A : V0 → V0 , so that Ak v ∈ V0 for any k ≥ 1 and any
v ∈ V0 . By definition we also have Cv = 0 for v ∈ V0 .
It follows that CAk v = 0 for each k when v ∈ V0 , so that
V0 ⊆ Σō as claimed.
T
The subspace V is equal to vectors of the form vP
= (z T |
T
d/2
z ), with z ∈ R , and V0 is the subset satisfying
zi = 0.
Symmetry implies the identity

0.5
0
−0.5
−1

Fig. 11.

1

Pool is on

24 1

Pool is off

24 hrs

First four eigenvectors of A for the pool model.

Fig. 11 shows four eigenvectors of A (left eigenvectors of
P0 ), ordered with respect to the magnitude of the corresponding eigenvalues. To avoid redundancy we considered only
eigenvalues with non-negative imaginary part. The eigenvector
v 1 has real eigenvalue λ1 = 1: it is proportional to the invariant
pmf π0 , expressed as a column vector. The remaining three
eigenvalues are complex, as are the corresponding eigenvectors
{v 2 , v 3 , v 4 }.
A reduced-order observer was constructed using a state
transformation based on the LTI model. The construction can
be interpreted as a truncated representation of the state with
respect to the eigenvectors of A. In view of the sinusoidal form
shown in Fig. 11, this approximation is similar to a Fourier
approximation.
After several experiments it was found that a 7th order
approximation fit the dynamics well: the Bode plots nearly
agreed, and the reduced order observer tracked the true output
U (Xt ) nearly as well as the full order observer.
Experiments were conducted for the model described in
Section IV-A, with a 12hr/day cleaning cycle, and a sampling
rate of 0.1%. Fig. 12 shows a comparison of the state and its

11

Empirical dist. Φt
N = 10 5

N = 10 4

Conditional dist. est. Φt
N = 10 5 and 10 4

0

1

Pools is on

24 1

Pools is off

24

Fig. 12. The reduced-order observer yields good estimates of the state Φt
for typical state values.

estimate at a particular time in two simulations, distinguished
by load population.
The empirical distribution Φt is approximately the same in
each case, N = 104 and N = 105 ; it is slightly smoother for
the larger population. As seen in the figure, the state estimates
are indistinguishable in the two cases. The state Φt is far from
the steady-state (proportional to v 1 shown in Fig. 11), yet the
estimates are good in each experiment. The explanation for
success is most likely the similar continuity of the empirical
distribution and the first seven eigenvectors.
Observability of a TCL model: The pool example is
similar to the TCL model considered in [2] and several other
recent papers. The controlled Markov model for the ith load
again consists of two parts: a binary component indicating if
the unit is on or off, and a continuous component indicating
its temperature. The state of the ith load at time t is denoted
Xti = (mit , θti ).
The complete system description follows [2]: For the nominal dynamics considered in prior work, this is defined by a
dead-band [Θmin , Θmax ]. For cooling devices,

i
θt+1
< Θmin
 0,
i
i
1,
θt+1 > Θmax
mt+1 =
 i
otherwise
mt ,
The temperature is modeled as a linear system driven by white
noise:
i
θt+1
= ai θti + (1 − ai )(θt0 − mit Ri %etr ) + ηti ,

in which 0 < ai < 1.
The parameters summarized in Table II are taken from [2].
In this prior work a heterogeneous model was considered, in
which the values (ai , Ri , C i ) were sampled from a probability distribution. A homogeneous model is considered in the
experiments surveyed here, in which Ri = C i = 2 for each i.
TABLE II: TCL PARAMETERS

The state space for the TCL model is continuous. A finite
state-space approximation is obtained by binning the deadband interval into 40 bins of equal size. A Markov chain model
is desired with a total of 80 states, since the state captures
temperature as well as whether the load is on or off.
An empirical model was obtained via Monte-Carlo. Simulation experiments with 10,000 TCLs and 3,600 time steps
were run to obtain an approximation of the steady-state joint
distribution,
Π(xj , xk ) = P{Xt = xj , Xt+1 = xk },

1 ≤ j, k ≤ 80.

The approximate transition matrix is then obtained using
Bayes’ rule:
X
Π(xj , xk )
, with π(xi ) =
Π(xj , xk ).
P0 (xj , xk ) =
i
π(x )
k

The approximate Markovian dynamics are used in [2] to
define the heuristic mean-field model,
{Akj } = {P0 (xj , xk )}.

Φt+1 = AΦt ,

The observation matrix used in [2] is the same as in the pool
example, except for a scaling:
Ck = 105 × %etr × I{State k is ON}.
The observability Grammian for this LTI system is highly
ill-conditioned: Fig. 13 shows that the rapid decay of the
eigenvalues is similar to what was observed in Fig. 3 for the
pool model.
|λi |

105

80 Eigenvalues of the
Observability Grammian

100
10 -5

10 -10

0

10

20

30

40

50

60

70

80

i

Fig. 13. Magnitude of eigenvalues of observability Grammian for TCL model

The structure of eigenvectors is also similar. Fig. 14 shows
four eigenvectors of the matrix A, corresponding to the four
eigenvalues of maximum magnitude (ignoring duplication
from complex conjugate pairs).
v1

0.2
Real part

0.1

v2

v3

v4

0

−0.2

Value

τ : time step (secs)
[Θmin , Θmax ]: temperature deadband (o C)
θ0 : ambient temperature (o C)
R: thermal resistance (o C/kW)
C: thermal capacitance (kWh/o C)
η: modeling noise (o C)
%etr : energy transfer rate (kW)
ai = exp{−τ /(CR)}

2
[20, 21]
32
2
2
∼ N (0, 2.5 × 10−7 )
14

Imaginary part

0.2

Parameter

0

−0.2

20

Fig. 14.

TCL is on

21 20

TCL is off

21 oC

First four eigenvectors of A for a TCL model.

A finer analysis of the continuous time / continuous state
model for TCLs (as considered in [11], [18]) may give greater
insight on the eigenstructure observed for these models.

12

V. C ONCLUSIONS
The construction of a Kalman filter for the joint population/individual dynamics is possible, and performance is
remarkable in the test cases considered. In particular, it is
very surprising to obtain accurate tracking of both the mean
and variance of QoS for an individual, given extremely noisy
estimates of the population.
An open topic for future research is the state estimation
techniques that take into account opt-out, which is used to
ensure good QoS to loads [3]. It may be possible to obtain
a reduced order observer for the very complex model for
joint state-QoS dynamics. Alternatively, it may be possible
to exploit the structure of the distribution of QoS observed in
numerical experiments. When subject to opt-out control, the
histogram of QoS is similar to a conditional Gaussian that is
constrained to the interval specified by opt-out parameters. A
nonlinear filter may be constructed that takes this structure
into account.

R EFERENCES
[1] M. Huang, P. E. Caines, and R. P. Malhame, “Large-population costcoupled LQG problems with nonuniform agents: Individual-mass behavior and decentralized ε-Nash equilibria,” IEEE Trans. Automat. Control,
vol. 52, no. 9, pp. 1560–1571, 2007.
[2] J. Mathieu, S. Koch, and D. Callaway, “State estimation and control
of electric loads to manage real-time energy imbalance,” IEEE Trans.
Power Systems, vol. 28, no. 1, pp. 430–440, 2013.
[3] Y. Chen, A. Bušić, and S. Meyn, “Individual risk in mean field
control with application to automated demand response,” in 53rd IEEE
Conference on Decision and Control, Dec 2014, pp. 6425–6432.
[4] D. Callaway and I. Hiskens, “Achieving controllability of electric loads,”
Proceedings of the IEEE, vol. 99, no. 1, pp. 184 –199, jan. 2011.
[5] P. Barooah, A. Bušić, and S. Meyn, “Spectral decomposition of demandside flexibility for reliable ancillary services in a smart grid,” in
Proc. 48th Annual Hawaii International Conference on System Sciences
(HICSS), Kauai, Hawaii, 2015, pp. 2700–2709.
[6] S. Meyn, P. Barooah, A. Bušić, Y. Chen, and J. Ehren, “Ancillary service
to the grid using intelligent deferrable loads,” IEEE Trans. Automat.
Control, vol. 60, no. 11, pp. 2847–2862, Nov 2015.
[7] H. Hao, Y. Lin, A. Kowli, P. Barooah, and S. Meyn, “Ancillary service to
the grid through control of fans in commercial building HVAC systems,”
IEEE Trans. on Smart Grid, vol. 5, no. 4, pp. 2066–2074, July 2014.
[8] A. Brooks, E. Lu, D. Reicher, C. Spirakis, and B. Weihl, “Demand
dispatch,” IEEE Power and Energy Magazine, vol. 8, no. 3, pp. 20–29,
May 2010.
[9] A. Kizilkale and R. Malhame, “A class of collective target tracking
problems in energy systems: Cooperative versus non-cooperative mean
field control solutions,” in IEEE Conference on Decision and Control,
2014, pp. 3493–3498.
[10] K. Christakou, D.-C. Tomozei, J.-Y. Le Boudec, and M. Paolone,
“GECN: primary voltage control for active distribution networks via
real-time demand-response,” IEEE Trans. on Smart Grid, vol. 5, no. 2,
pp. 622–631, March 2014.
[11] L. C. Totu, “Large scale demand response of thermostatic loads,” Ph.D.
dissertation, Faculty of Engineering and Science, Aalborg University,
2015.
[12] P. E. Caines, Linear Stochastic Systems. New York: John Wiley &
Sons, 1988.
[13] N. V. Krylov, R. S. Lipster, and A. A. Novikov, “Kalman filter for
Markov processes,” in Statistics and Control of Stochastic Processes.
New York: Optimization Software, inc., 1984, pp. 197–213.

[14] P. Caines and A. Kizilkale, “Recursive estimation of common partially
observed disturbances in MFG systems with application to large scale
power markets,” in 52nd IEEE Conference on Decision and Control,
Dec 2013, pp. 2505–2512.
[15] N. Sen and P. Caines, “Mean field games with partially observed major
player and stochastic mean field,” in IEEE Conference on Decision and
Control, 2014, pp. 2709–2715.
[16] Y. Chen, “Markovian demand dispatch design for virtual energy storage
to support renewable energy integration,” Ph.D. dissertation, University
of Florida, Gainesville, FL, USA, July 2016.
[17] “BPA balancing authority,” Online, http://tinyurl.com/BPAgenload, http:
//tinyurl.com/BPAbalancing.
[18] R. Malhamé and C.-Y. Chong, “On the statistical properties of a
cyclic diffusion process arising in the modeling of thermostat-controlled
electric power system loads,” SIAM J. Appl. Math., vol. 48, no. 2, pp.
465–480, 1988.

Yue Chen (S’ 14) received the B.E. degree in
electronic and information engineering from Harbin
Engineering University, Harbin, China, in 2008, the
M.E. degree in electrical engineering from University of Detroit Mercy, Detroit, USA, in 2010, the
M.S. degree in Electrical and Computer Engineering
in 2012 from University of Florida, Gainesville,
USA. During the Fall of 2015 he was an intern at Los
Alamos National Laboratory conducting research on
control techniques for demand dispatch. In 2016 he
is expected receive the Ph.D. degree in electrical and
computer engineering at University of Florida.

Ana Bušić is a Research Scientist at Inria Paris
– Rocquencourt, and a member of the Computer
Science Department at Ecole Normale Supérieure,
Paris, France. She received a M.S. degree in Mathematics in 2003, and a Ph.D. degree in Computer Science in 2007, both from the University of Versailles.
She was a post-doctoral fellow at Inria Grenoble –
Rhône-Alpes and at University Paris Diderot-Paris
7. She joined Inria in 2009. She is a member of the
Laboratory of Information, Networking and Communication Sciences, a joint lab between AlcatelLucent, Inria, Institut Mines-Télécom, UPMC Sorbonne Universities, and
SystemX. Her research interests include stochastic modeling, simulation,
performance evaluation and optimization, with applications to communication
networks and energy systems.

Sean P. Meyn (S’85, M’87, SM’95, F’02) received
the B.A. degree in Mathematics Summa Cum Laude
from UCLA in 1982, and the PhD degree in Electrical Engineering from McGill University in 1987
(with Prof. P. Caines). He held a two year postdoctoral fellowship at the Australian National University
in Canberra, and was a professor at the University of
Illinois from 1989-2011. Since January, 2012 he has
been a professor and has held the Robert C. Pittman
Eminent Scholar Chair in Electrical and Computer
Engineering at the University of Florida. He is coauthor with Richard Tweedie of the monograph Markov Chains and Stochastic
Stability, Springer-Verlag, London, and received jointly with Tweedie the 1994
ORSA/TIMS Best Publication In Applied Probability Award. His research
interests include stochastic processes, optimization, complex networks, and
information theory, with applications to smarter cities and smarter grids.

