Estimation and Control of
Quality of Service in Demand Dispatch
Yue Chen, Ana Bušić, and Sean Meyn∗†

arXiv:1609.00051v1 [] 31 Aug 2016

March 13, 2018

Abstract
It is now well known that flexibility of energy consumption can be harnessed for the purposes
of grid-level ancillary services. In particular, through distributed control of a collection of loads, a
balancing authority regulation signal can be tracked accurately, while ensuring that the quality of
service (QoS) for each load is acceptable on average. In this paper it is argued that a histogram of
QoS is approximately Gaussian, and consequently each load will eventually receive poor service.
Statistical techniques are developed to estimate the mean and variance of QoS as a function
of the power spectral density of the regulation signal. It is also shown that additional local
control can eliminate risk: The histogram of QoS is truncated through this local control, so that
strict bounds on service quality are guaranteed. While there is a tradeoff between the grid-level
tracking performance (capacity and accuracy) and the bounds imposed on QoS, it is found that
the loss of capacity is minor in typical cases.
Index Terms—Demand dispatch, demand response, ancillary services, mean field control.

1

Introduction

The power grid requires regulation to ensure that supply matches demand. Regulation is required by
each balancing authority (BA) on multiple time-scales, corresponding to the time-scales of volatility
of both supply and demand for power. Resources that supply these regulation services are collectively
known as ancillary services. FERC orders 755 and 745 are examples of recent federal policy intended
to provide incentives for the provision of these services.
A number of papers have explored the potential for extracting ancillary service through the
inherent flexibility of loads. Examples of loads with sufficient flexibility to provide service to the
grid are aluminum manufacturing, plug-in electric vehicles, heating and ventilation (HVAC), and
water pumping for irrigation [1–5]. Even with direct load control, there may be delay and dynamics,
so harnessing ancillary services from flexible loads amounts to a control problem: The BA wishes
to design some signal to be broadcast to loads, so that deviation in power consumption tracks a
reference signal. It has been argued that a randomized control architecture at each load can simplify
this control problem [4, 6, 7].
Fig. 1 shows a schematic of the Demand Dispatch control architecture adopted in [4], in which
each load operates according to a randomized policy based on its internal state, and a common control signal ζ. Theoretical results and examples in this prior work demonstrate that local randomized
policies can be designed to simplify control of the aggregate to provide reliable ancillary services.
∗

This research is supported by the NSF grants CPS-0931416 and CPS-1259040, the French National Research
Agency grant ANR-12-MONU-0019, and US-Israel BSF Grant 2011506.
†
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.

1

Load 1

r

BA

Gc

ζ

Load 2

+

...

Reference (MW)

Power
Consumption (MW)

Load N

Figure 1: Control architecture for Demand Dispatch
Absent in prior work is any detailed analysis of risk for an individual load (with the exception of
the preliminary work [8] on which the present paper is built). In the setting of [4] it can be argued
that the quality of service (QoS) for each load is acceptable only on average. Strict bounds on QoS
are addressed in [9] for a deterministic model. The service curves considered there are of similar
flavor to the QoS metrics used in the present work.
The main contributions of this paper are summarized as follows: QoS metrics are proposed in the
first part of the paper, along with techniques to approximate their mean and variance. Estimates
of second order statistics are a function of the power spectral density (PSD) of the grid regulation
signal. These theoretical results are developed in Section 3.
A simple approach is proposed to restrict QoS to pre-specified bounds: A load will opt-out
of service to the grid temporarily, whenever its QoS is about to exit pre-specified bounds. This
essentially eliminates risk since the histogram of QoS is restricted to these bounds.

Fraction of loads

×10 -3 4

Histogram QoS
Gaussian approximation

3
2
1
0

-250

0

QoS without opt out

250

-250

0

QoS with opt out

250

Figure 2: Discounted QoS (9) with and without local opt-out control.
Numerical results surveyed in Section 4 confirm that the histogram of QoS is truncated through
this local control, so that strict bounds on service quality are guaranteed. As long as the QoS bounds
are not overly restrictive, the impact on the grid-level performance is insignificant.
Fig. 2 shows histograms of QoS based on simulation experiments described in Section 4. The
plot on the left hand side shows that a Gaussian approximation is a good fit with empirical results
when there is no local opt-out control. The figure on the right shows how the histogram is truncated
when opt-out control is in place.
The methodology used in this work employs a linear representation of the controlled Markov
chain that represents the state of an individual load [10]. Similar representations are used in [11,12]
2

to obtain nonlinear filters for state estimation.
We begin with a brief survey of a portion of results from [4], and a precise definition of QoS for
a load.

2
2.1

Randomized control and mean-field models
Randomized control

The system architecture considered in this paper is illustrated in Fig. 1, based on the following
components:
(i) There are N homogeneous loads that receive a common scalar command signal from the balancing
authority (BA), denoted ζ = {ζt } in the figure.

(ii) Each load evolves as a controlled Markov chain on the finite state space X = {x1 , . . . , xd }. Its
transition probability is determined by its own state, and the BA signal ζ. 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 given by,
i
P{Xt+1
= x0 | Xri , ζr : r ≤ t} = Pζ (x, x0 ),

(1)

for each x0 ∈ X, Xti = x ∈ X and ζt = ζ ∈ R.
(iii) The BA has measurements of the other two scalar signals shown in the figure: The aggregate
power consumption y and desired deviation power consumption r.
An approach to construction of {Pζ : ζ ∈ R} was proposed in [4] based on information-theoretic
arguments. The nominal behavior is defined as the dynamics with ζ ≡ 0.
In the present paper we do not require a specific construction of Pζ , but the following assumptions
are imposed in our main results.
A1: The transition matrix Pζ is twice continuously differentiable (C 2 ) in a neighborhood of ζ = 0,
and the second derivative is Lipschitz continuous. In addition, the nominal transition matrix P0
is irreducible and aperiodic.
The first and second order derivatives of the transition matrix at ζ = 0 are denoted,
E=

d
Pζ
dζ

ζ=0

E (2) =

,

d2
Pζ
dζ 2

ζ=0

.

(2)

A2: ζt = εζt1 , where 0 ≤ ε < 1 and ζ 1 = {ζt1 : t ∈ Z} is a real-valued stationary stochastic process
with zero mean. The following additional assumptions are imposed:
(i) It is bounded, |ζt1 | ≤ 1 for all t with probability one. Hence σζ2 = E[(ζt )2 ] ≤ ε2 .

(ii) Its auto-covariance satisfies, for each t,

|Σζ (t)| ≤ ε2 bρ|t| , with b < ∞, and |ρ| < 1.
It is assumed that the power consumption at time t from load i is equal to some function of the
state, denoted U (Xti ). The normalized power consumption is denoted,
ytN =

N
1 X
U (Xti ).
N
i=1

3

(3)

0
P Under Assumption A1, P0 has a unique pmf (probability mass function) π0 . The value y :=
x π0 (x)U (x) is interpreted as the average nominal power usage. On combining the ergodic theorem
for Markov chains with the Law of Large Numbers for i.i.d. sequences we can conclude that ytN ≈ y 0
when both N and t are large, and ζ ≡ 0.
It is assumed that the signal r is also normalized so that tracking amounts to choosing the signal
ζ so that ỹtN ≈ rt for all t, where ỹtN = ytN −y 0 is the deviation from nominal behavior. For example,
we might use error feedback of the form,

et = rt − ỹtN ,

ζt = Gc et ,

(4)

where Gc is the control transfer function [4, 13, 14]. Adopting terminology from control engineering,
r will be called the reference signal.

2.2

Mean-field model

The mean-field model is based on the empirical pmfs:
µN
t (x)

N
1 X
:=
I{Xti = x},
N
i=1

x ∈ X.

(5)

Each entry of the vector µN
t represents the fraction of loads in a particular state. Under very general
conditions on the input sequence ζ, it can be shown that the empirical distributions converge to a
solution to the nonlinear state space model equations,
µt+1 = µt Pζt ,

(6)

P
where µt is considered as a row vector. The output is denoted yt = x µt (x)U (x), which is the limit
of the {ytN } (defined in (3)) as N → ∞. See [4] for details.
The unique equilibrium with ζ ≡ 0 is µt ≡ π0 and yt ≡ y 0 . The analysis in [4] can be extended
to obtain a linearization about this equilibrium. This is described by the linear state space model,
Φt+1 = AΦt + Bζt
γt = CΦt

(7)

where A = P0T , C is a row vector of dimensionPd = |X| with Ci = U (xi ) for each i, and B is a
d-dimensional column vector with entries Bj = x π0 (x)E(x, xj ), where the matrix E is defined in
(2). In the state equations (7), the state Φt is d-dimensional, and Φt (i) is intended to approximate
µt (xi ) − π0 (xi ) for 1 ≤ i ≤ d. The output γt is an approximation of ỹt = yt − y 0 .

2.3

QoS for an individual and the population

The QoS metrics considered in this paper are defined by a scalar-valued function ` : X → R, and a
stable transfer function HL . The QoS of the ith load at time t is then defined by Lit = HL Lit , where
Lit = `(Xti ). For example, the function ` may represent temperature, cycling, or power consumption
as a function of x ∈ X.
Two classes of transfer functions will be considered in numerical experiments:
(i) Summation over a finite time horizon Tf :
Lit

=

Tf
X

i
`(Xt−k
).

k=0

4

(8)

(ii) Discounted sum:
Lit

=

∞
X

i
β k `(Xt−k
),

(9)

k=0

where the discount factor satisfies β ∈ [0, 1).
In particular, setting Tf = 0 or β = 0 gives Lit = `(Xti ).
Unless elsewhere specified, the function ` is specialized to reflect the power consumption of a
load,
`(Xti ) = U (Xti ),
(10)
and its normalized form is also considered,
`(Xti ) = U (Xti ) − y 0

(11)

where y 0 is defined after (3).
The average QoS at time t is denoted,
N
1 X i
L̄t =
Lt ,
N
i=1

and the filtered signal is denoted Rt = HL rt . The following result follows from the definitions:
Proposition 2.1 Suppose there is perfect tracking: ỹtN = rt for all t. Then, under the definition of
`(Xti ) in (11),
L̄t = Rt .
(12)

In practice we can only expect the approximation ỹtN ≈ rt , which will imply the corresponding
approximation L̄t ≈ Rt . This leads to a useful heuristic notion of capacity of the aggregate of loads
as a function of time:
Battery analogy With ` defined in (11), the signal L̄t is similar to the SOC (state of charge) of
a battery. In particular, a large value of L̄t suggests loads have a large capacity for “discharge”.
Since Rt approximates L̄t by design, the former can be used as an indicator of SOC of loads.

2.4

Example: Intelligent pools

The paper [4] on which the present work is based considered an application of this system architecture
in which the loads are a collection of pools. The motivation for considering pools is the inherent
flexibility of pool cleaning, and because the total load in a region can be very large. The maximum
load is approximately 1GW in the state of Florida.
In this prior work, the state space was the finite set,
X = {(κ, j) : κ ∈ {⊕, }, j ∈ {1, 2, . . . , I}}.

(13)

The load samples the grid signal periodically (the sampling increments are assumed to be deterministic, or i.i.d. and distributed according to a geometric distribution). At the time of the kth sample,
if Xk = ( , j) then the load has remained off for the past j sampling times, and was turned off at
sampling time k − j. The interpretation of Xk = (⊕, j) is symmetrical, with “⊕” indicating that
the load is currently consuming power.
5

A technique to construct the state transition matrix Pζ was introduced in [4] using an optimalcontrol approach, and improved in [12, 15] to reduce the dimension of the state space, and extend
application beyond this class of loads. It has been demonstrated through extensive simulation
that tracking at the grid level is nearly perfect, provided that the reference signal is scaled and
filtered [4, 8, 12]. What was left out in prior work is any detailed consideration of the quality of
service to individual loads. The pool pump example is ideal for illustrating the possibility of risk
for an individual load, and illustrating how this risk can be reduced or even eliminated.
In a numerical experiment conducted with 104 pools, each pool pump was expected to clean the
pool 12 hour/day and assumed to consume 1 kW
Pwhen operating. Hence the function U that defines
(10) and (11) is the indicator function U (x) = j I{x = (⊕, j)}. Fig. 3 shows the histogram of the
QoS metric Lit based on (8) and (10), with Tf corresponding to 10 days.

Fraction of loads

5%

Histogram: σ 2 = 79.2

Approximation: σ 2 = 87.7

4%
3%
2%
1%

80

90

100

110

120

130

Operation hours

140

150

160

Figure 3: Histogram of the moving-window QoS metric (8).
The histogram appears Gaussian, with a mean of approximately 120 hours — this is consistent
with the time horizon used in this experiment, given the nominal 12 hour/day cleaning period. It
is evident that a fraction of pools are over-cleaned or under cleaned by 24 hours or more.
An analysis of QoS is presented in the next section based on a model of an individual load in
the mean field limit.

3

QoS analysis and opt-out control

In the mean-field limit, the aggregate dynamics are deterministic, following the discrete-time nonlinear control model (6). The behavior of each load remains probabilistic.

3.1

Mean field model for an individual load

N
The mean field model for one load is defined by replacing (µN
t , ζt ) with its mean-field limit (µt , ζt ).
The justification is that we have a very large number of loads, but our interest is in the statistics of
an individual.
The super-script i is dropped in our analysis of a single load. Hence X denotes the controlled
Markov chain whose transition probabilities are defined consistently with (1):

P{Xt+1 = x0 | Xr , ζr : r ≤ t} = Pζt (x, x0 ),
6

x = Xt

(14)

The construction of the mean field model (6) is based on lifting 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 Γt ∈ S
is the degenerate distribution whose mass is concentrated at x if Xt = x; that is, Γt = δx .
With this state description, the load evolves according to a random linear system, similar to (6):
Γt+1 = Γt Gt+1

(15)

in which ΓP
t is interpreted as a d-dimensional row vector. The d × d matrix Gt has entries 0 or 1
only, with x0 ∈X Gt (x, x0 ) = 1 for all x ∈ X. It is conditionally independent of {Γ0 , · · · , Γt }, given
ζt , with
E[Gt+1 |Γ0 , · · · , Γt , ζt ] = Pζt .
(16)
The random linear system (15) can be described as a linear system driven by “white noise”:
Γt+1 = Γt Pζt + ∆t+1

(17)

where, {∆t+1 = Γt (Gt+1 − Pζt ) : t ≥ 0} is a martingale difference sequence.
A Taylor-series approximation of Pζt leads to a useful approximation of (17). Recall that the
first and second derivatives E and E (2) were introduced in (2). The proof follows from the definitions.
Proposition 3.1 The nonlinear system (17) admits the LTI approximation,
Γt+1 = Γt P0 + Dt+1 + O(ε3 )

(18)

Dt+1 := BtT ζt + VtT ζt2 + ∆t+1 ,

(19)

in which,
with BtT = Γt E and VtT =

3.2

1
(2) .
2 Γt E

Steady-state QoS statistics

The main goal of this section is to estimate the second order statistics of {Lt } for arbitrary function
` and stable filter HL . These approximations are obtained for a stationary realization.
Theorem 3.2 Under assumptions A1 and A2, there exists a realization {ζt , Xt , Dt , Lt : −∞ < t <
∞} that is jointly stationary.

The proof of the existence of the stationary realization {ζt , Xt , Dt : −∞ < t < ∞} follows from
Proposition 2.3 of [10]. Stationarity of the joint process {ζt , Xt , Dt , Lt } follows from the assumption
that {Lt } is obtained from {`(Xt )} through a stable filter. Proposition 2.5 of [10] implies the
following approximation for the mean QoS.
Proposition 3.3 The mean QoS admits the approximation,
X
E[Lt ] = HL (1)E[Lt ], E[Lt ] =
`(x)π0 (x) + O(ε2 )
x

where HL (1) is the DC gain of the transfer function HL .



The PSD for the stationary realization of the stochastic process D can be approximated by a
d×d
second order Taylor series expansion: Applying Theorem 2.4 of [10], a function S(2)
D : [−π, π] → C
can be constructed such that
2
SD (θ) = S•D (θ) + ε2 S(2)
D (θ) + o(ε ),

θ ∈ [−π, π] ,

(20)

in which S•D is the PSD obtained with ζ ≡ 0. This is independent of θ since D is a white noise
sequence in this case. An approximation for the PSD of {Lt } follows from (20) and the following:
7

Dt

[Iz − P0 ]−1

Γt

Γt , 

Lt

HL (z)

Lt

Figure 4: Linear model of QoS evolution.
(i) The model (15) is approximated by the linear state space model,
b t+1 = Γ
b t P0 + Dt+1 .
Γ

(21)
(2)

(ii) Formulae for S•D and S(2)
D are given in [10]. The computation of SD in (20) requires the autocovariance of ζ. This must be estimated from data.
(iii) The approximation of the PSD SL for {Lt } is obtained using (20) and the linear model shown in
Fig. 4.
The variance is then obtained based on an inverse formula, or an approximation:
Theorem 3.4 The variance of QoS admits the following representation and approximations:
(i) The variance of QoS is the average of its PSD:
Z 2π
1
SL (θ) dθ.
VL =
2π 0

(22)

In the special case of the moving time-horizon (8), it admits the approximation
VL
≈ SL (0)
Tf

if Tf ≈ ∞,

(23)

where SL is the PSD of L.
(ii) A Taylor series of the right hand side of (22) implies the approximation,
(2)

VL = VL• + ε2 VL + o(ε2 ),

(24)

(2)

where the terms VL• and VL are obtained using the QoS variance formula in (22), and the
approximation of the PSD SL is obtained using (20).

Recall that the QoS histogram shown in Fig. 3 is based on the moving time-horizon (8), with
Tf corresponding to 10 days. The approximation of the QoS variance was obtained using (23) and
a Taylor series approximation of SL (0).

3.3

Opt-out control

An extra layer of control is required to truncate the two tails of the QoS histograms observed in
experiments.
The local control considered in this paper is a simple “opt-out” mechanism, based on pre-defined
upper and lower limits b+ and b− . A load ignores a command from grid operator at time t if it may
result in Lt+1 6∈ [b− , b+ ], and takes an alternative action so that Lt+1 ∈ [b− , b+ ] with probability
one. This ensures that the QoS metric of each load remains within the predefined interval for all
time.
In many cases, the poor QoS revealed by the two tails of the histogram represents only a small
portion of loads. Therefore, the impact from local opt-out control is insignificant at the grid level if
the QoS interval [b− , b+ ] is carefully chosen.
8

4

Numerics

Numerical experiments were conducted on the pool pump model to illustrate the main technical
conclusions. Simulation results are summarized below:
(i) Applications of Theorem 3.4 show that the approximations of QoS closely match observed QoS.
(ii) Opt-out control ensures that QoS lies within strict bounds, and tracking remains nearly perfect
in most cases. In some extreme cases, capacity is reduced with the introduction of opt-out
control. However, it is found through simulations that it is far less than might be predicted by
the approximation in (i).
(iii) The opt-out control can be applied to multiple QoS metrics. The capacity is further reduced with
the introduction of an additional QoS metric, but the reduction is found to be minor numerical
experiments.
(iv) Approaches are proposed and tested in Section 4.5 to condition the reference signal to reduce
potential capacity reduction caused by opt-out local control.

4.1

Simulation setup

The simulation used N = 104 homogeneous Markov models: Each pool pump is operated under a
12 hours/day cleaning cycle, and consumes 1 kW during operation.
Two QoS metrics are considered, differentiated by the function ` appearing in the definition
i
Lt = `(Xti ): In the first QoS function, the normalized power consumption (11) is considered so that
if Lit > 0 (Lit < 0) then the pool has been over-cleaned (under-cleaned).
i )=
The second QoS function is introduced to capture the on/off cycling of loads: `c (Xti , Xt+1
Xn
i
|I{Xt+1
= (⊕, j)} − I{Xti = ( , j)}|
j
(25)
o
i
i
+ |I{Xt+1 = ( , j)} − I{Xt = (⊕, j)}| .
These two QoS metrics can be applied to many other loads, such as air conditioners, refrigerators,
and water heaters [?, 15]. Consideration of the QoS metric (25) is postponed to Section 4.6.
The discounted sum (9) was used to define Lit in these experiments. The following interpretation
is used to obtain insight on the choice of the discount factor in this QoS metric.
Let ξ denote a random variable that is independent of X i . Its distribution is taken to be
geometric on Z+ , with parameter β. Its mean is thus
E[ξ] =

∞
X
k=0

P(ξ ≥ k) =

∞
X
k=0

βk =

1
.
1−β

By independence we also have,
E

ξ
hX
k=0

∞
i X
i
i
`(Xt−k
) =
E[`(Xt−k
)I(ξ ≥ k)]
k=0

=

∞
X
k=0

=

∞
X

i
E[`(Xt−k
)]P(ξ ≥ k)
i
E[`(Xt−k
)]β k = E[Lit ]

k=0

9

(26)

wt0

Gwr

rt•

GBP

rt

Gc
1 + G c Gp

ζt

Figure 5: The input ζ modeled as a stationary stochastic process
Hence the discounted sum (9) is similar to the moving window QoS metric (8). In the following
experiments we took E[ξ] = 10 days, which corresponds to 2880 time samples (given the 5 minute
sampling period). Solving (26) gives β = 1 − 1/2880.

4.2

Model for the stationary input

A linear model for the stationary input ζ was constructed based on the block diagram shown in
Fig. 5. As in the prior work [4,8], it is assumed that the signal r is obtained by filtering the regulation
signal r • ; the latter is modeled as filtered white noise with transfer function Gwr .
The BPA (Bonneville Power Authority [16]) balancing reserves, deployed in January 2015, were
taken for the regulation signal r • . The sampling time is 5 minutes.
Sr• (f )

0

Linear model
Matlab

PSD (dB)

-10
-20
-30
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Frequency (× 1/600 Hz)

1f

Figure 6: BPA signal: Power spectral density and its approximation.
It is assumed that r • evolves as the ARMA (autoregressive moving average) model,
•
•
rt• + a1 rt−1
+ a2 rt−2
= wt + b1 wt−1

(27)

2 . The extended least squares (ELS) algorithm was used
in which w is white noise with variance σw
2 based on the BPA balancing reserves.
to estimate the coefficients a1 , a2 , b1 , and the variance σw
T
2 = 4.36 × 10−3 . In the
The algorithm terminated at [a1 , a2 , b1 ] = [−1.16, 0.2301, −0.2489]T , and σw
z-domain, its transfer function is expressed

Gwr (z) =

1 − 0.2489z −1
.
1 − 1.6z −1 + 0.2301z −2

(28)

The dashed and solid lines in Fig. 6 represent the estimates of the spectrum given by |Sr• (ejθ )|2 =
2 |G (ejθ )|2 and Matlab’s psd command, respectively.
σw
wr
The transfer function GBP in Fig. 5 is a filter designed to smooth the balancing reserve signal.
A low pass filter was adopted with crossover frequency near the nominal period of a single load:
1/(24 hours) in this example. In most of the experiments reported here, GBP is the first-order
Butterworth low-pass filter with cut-off frequency fc = 1/(1000 minutes):
GBP (z) = 0.0155

1 + z −1
.
1 − 0.9691z −1

10

(29)

The reference signal must also be scaled so that the desired goal ỹtN ≈ rt is feasible for all t.
Denote by r 1 the signal obtained using the largest scaling, while also ensuring that this signal can
be tracked by the collection of pools. This was obtained through trial and error.
In the tracking plots, such as Fig. 11, the signals ỹtN and rt are re-scaled to their original units
in MWs.
The construction of a stationary model for the input process ζ was based on the linearized
mean-field model, and the scaled reference signal defined by rt = εrt1 , t ∈ Z. The linear state space
model (7) leads to a transfer function from ζ to γ that is denoted Gp (recall that γ is intended
to approximate ỹ). Based on this linear approximation, we obtain a linear mapping r → ζ via
the transfer function Gc /(1 + Gc Gp ), where Gc is the PI controller with proportional gain 50 and
integral gain 1.5.

4.3

Individual QoS

The QoS of pools without opt-out control is illustrated using the histogram on the left-hand-side
of Fig. 2 using the reference signal r 1 . The histogram is based on samples of the QoS function (9)
from each load over one month. The dashed lines represent Gaussian densities with zero mean, and
variances σL2 obtained using the approximation given in (22).

Variance of QoS

×10 3

7

Approximation
Experiment

6
5
4
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

ε2

Figure 7: Linear approximation for the variance of QoS.
The conclusions of Theorem 3.4 are illustrated in Fig. 7, where the QoS variances from simulation
and Gaussian approximations are presented for several values of ε. It is seen that the variance is
approximately linear in ε2 for small ε > 0, with slope as predicted in Theorem 3.4.
The role of bandpass filter A range of cut-off frequencies were considered in order to investigate
the impact of the bandpass filter that is used to obtain r. Fig. 8 shows a comparison of the variance
of the reference and the variance of QoS as a function of the cut-off frequency fc over a range of
frequencies; The linear growth in QoS variance compared to the slow growth of the variance of the
reference signal justifies a filter with fc < 10−2 .
The remaining numerical results that follow are based on the low pass filter (29) based on
fc = 1/(1000 minutes).

4.4

Grid level performance and opt-out control

We present in this section experimental results with opt-out control to ensure that QoS is subject
to strict constraints. Four QoS intervals were considered corresponding to error tolerances of, respectively, 5%, 10%, 15%, and 20%. Recall that the discount factor β was chosen by approximating

11

4

-2
4 ×10

3

3

2

2

1

1
1

2

3

4

5

6

7

8

9

Variance of reference

4

Variance of QoS

×10

10

Cut-off frequency fc of LP filter (cycles/minutes) ×10 -3

Figure 8: QoS variance increases with bandwidth of the reference signal
the moving window of 10 days, which corresponds to 2880 time samples. Based on this interpretation, these percentages are converted to intervals for local opt-out control (following the notation in
Section 3.3) as [b− , b+ ] = [−36, +36], [−72, +72], [−108, +108], and [−144, +144], respectively.
Fig. 9 shows four histograms of QoS with reference signal r 1 . These histograms are truncated
to the predefined QoS intervals as desired.
1.5%
Fraction of loads

QoS Constraint
5%

1.0%

10%
15%
20%

0.5%

-150

-100

-50

0

50

100

150

QoS

Figure 9: QoS histograms with opt-out control using different QoS bounds.
A normalized root mean square error (NRMSE) was adopted as the metric of grid level tracking
performance to study the impact of local opt-out control,
NRMSE =

RMS(e) − RMS(e• )
,
RMS(r)

(30)

where e• denotes the tracking error sequence obtained without opt-out control, and ε = 0 (so that
r ε ≡ 0). It is found in simulation RMS(e• ) ≈ 8.35 kW.
The grid level tracking performance with and without opt-out control is illustrated in Fig. 10.
The tracking performance with 10%, 15%, or 20% QoS interval remains nearly perfect. This is
surprising, given the improvement in QoS shown in Fig. 9. The explanation is that very few loads
opt out: For example, in the simulation with 10% QoS interval and reference signal r 1 , no more
than 1% of loads opt out at any time.
However, there are limitations on the capacity of ancillary service from a collection of loads, and
experiments reveal that opt-out control can reduce capacity. As seen from Fig. 10, the additional
12

Normalized RMS Error

0.3

QoS Constraint
5%

10%

0.2

15%
20%

0.1
0
0.1

Without opt-out

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

ε

Figure 10: Tracking NRMSEs (30), with and without local opt-out control.
opt-out control with 5% QoS interval degrades the grid level tracking performance, especially when
the reference signal is large.
We next apply Proposition 2.1 to better understand the grid-level impact of opt-out control,
and to design an algorithm to reshape the reference signal so that loads are less likely to opt out of
service.

4.5

Re-shaping the reference input

Recall the SOC heuristic introduced following Proposition 2.1. The proposition implies that L̄t ≈ Rt
for all t when there is near perfect tracking. Conversely, this result suggests that many loads will
opt-out, leading to poor tracking, if the SOC Rt is near the boundary of the QoS interval [b− , b+ ].

MW

5

rt

(Reference)

ỹtN (Power deviation)

0

-5

Rt

(Filtered ref signal)

L̄t

(Average QoS)

0

-7
350

5%

400

450

Hours

500

550

0

600

Opt out
tpercentage

7% Opt-out Interval

7

Figure 11: Poor tracking results when the SOC Rt reaches its boundaries.
To illustrate the application of these concepts, consider the case with 7% QoS interval and
reference signal r 1 . Results are provided in Fig. 11. Most of the time, the tracking results are nearly
perfect and hence the average QoS approximates the filtered reference signal, L̄t ≈ Rt . However, at
time ≈ 450 hr., Rt falls close to −50.4, which is the lower QoS bound for an individual load. During
this time period, many loads opted out, which resulted in degraded tracking.
In conclusion, to ensure good grid level tracking, the reference signal must respect any QoS
13

constraints. The grid operator should re-shape the reference signal to ensure that the SOC of loads
Rt does not approach its limits. A smooth transformation is required because the reference signal
for this example should not have significant high frequency content.
Here is one approach, based on two non-negative parameters: a threshold τ < 1 and a gain
parameter δ > 0. The following recursive algorithm is designed to increase (decrease) the reference
signal r when R reaches its lower (upper) threshold. The re-shaped reference signal is defined as
follows:


[rt − δ(Rt − τ b+ )]+ Rt > τ b+ and rt > 0
r̄t = [rt + δ(Rt − τ b− )]− Rt < τ b− and rt < 0
(31)


rt
otherwise
where, [b− , b+ ] defines the QoS interval, and [a]+ ([a]− ) denote the positive (negative) part of a real
number a.
When r 1 is transformed in this way using δ = 0.006 and τ = 0.65 it is found in numerical
experiments that this small change in the reference signal results in perfect tracking, and a much
smaller percentage of loads opting out – a complete simulation study will be included in [7].

4.6

Multiple QoS metrics

We next investigate the impact of opt-out control based on two QoS metrics: the cleaning QoS
function (11) and the cycling QoS function (25). The lower bound of cycling QoS was set to −∞.
The opt-out priority is assigned to the cleaning QoS metric when both QoS metrics reach their
bounds at the same time. Thus, the cycling QoS metric occasionally exceeds its predefined upper
bound.

Fraction of loads

0.3%

6%

without opt-out

5%

with opt-out

4%

0.2%

3%
2%

0.1%

1%
-300 -200 -100

0

100 200 300

15

Cleaning QoS

20

Cycling QoS

25

30

Figure 12: QoS bounds are maintained using local opt-out control.
Four error tolerances {5%, 10%, 15%, 20%} are considered for determining the opt-out interval
[b− , b+ ] for each of the two QoS metrics. Fig. 12 illustrates an example of QoS improvement based
on a 15% error tolerance on both QoS metrics, using the reference signal r 1 . Tracking in this case
was nearly perfect.
Tracking performance for a range of opt-out parameters is summarized in Fig. 13 using 16 colored
bars, distinguished by each pair of QoS constraints. Each bar represents tracking errors for different
reference signal scaling factors, 0.1 ≤ ε ≤ 1. The darkest color represents NRMSE (30) of 10% or
greater, and lighter colors represent smaller values of NRMSE (indicated on the color bar label).
Results in this figure show that opt-out control based on these two QoS metrics have little impact
on tracking error over a large range of opt-out intervals. For those cases that local opt-out largely
degrades grid-level tracking, we can either reduce the reference signal or relax QoS constraints to
maintain good tracking.
14

100%

8%

50%

7%
6%

0
20%

5%
4%

15%

3%

Cycling QoS 10%
5%

5%

10%

15%

20%

Tracking error

Capacity ε

≥ 10%
9%

2%
1%

Cleaning QoS

Figure 13: Tracking performance with two QoS constraints.

5

Conclusions

The main technical contribution of this paper is the approximation of QoS for an individual load. It
is remarkable that it is possible to obtain accurate estimates of first and second order statistics for
an individual load, taking into account second order statistics of exogenous inputs (in this case the
reference signal), along with correlation introduced by the Markovian model. It is also remarkable
that strict bounds on QoS can be guaranteed while retaining nearly perfect grid-level tracking.
Under certain conditions, the overall QoS of a collection of loads is predictable to the grid
operator. With this information, the grid operator can estimate the flexibility/capacity of loads and
modify the reference signal if necessary, to maintain high quality of both grid-level tracking and QoS
of loads.
Open problems remain in the area of estimation and control. It is hoped that estimates of the
statistics of QoS can be obtained in the presence of opt-out control. It is also likely that control
performance can be improved further with a more sophisticated approach to opt-out control.

References
[1] R. Malhame and C.-Y. Chong, “Electric load model synthesis by diffusion approximation of a
high-order hybrid-state stochastic system,” IEEE Trans. Automat. Control, vol. 30, no. 9, pp.
854 – 860, Sep 1985.
[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] A. Kizilkale and R. Malhame, “Mean field based control of power system dispersed energy
storage devices for peak load relief,” in 52nd IEEE Conference on Decision and Control (CDC),
2013, pp. 4971–4976.
[4] 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.
[5] 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.
15

[6] P. Barooah, A. Bušić, and S. Meyn, “Spectral decomposition of demand-side 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.
[7] 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.
[8] 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.
[9] J.-Y. Le Boudec and D.-C. Tomozei, “Demand response using service curves,” in International
Conference and Exhibition on Innovative Smart Grid Technologies (ISGT Europe), 2011, pp.
1–8.
[10] Y. Chen, A. Bušić, and S. Meyn, “Ergodic theory for controlled Markov chains with stationary
inputs,” ArXiv e-prints and submitted for publication, Annals of Applied Prob., April 2016.
[11] 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.
[12] Y. Chen, A. Bušić, and S. Meyn, “State estimation for the individual and the
population in mean field control with application to demand dispatch,” CoRR and to
appear, IEEE Transactions on Auto. Control, June 2017., 2015. [Online]. Available:
http://arxiv.org/abs/1504.00088v1
[13] R. Couillet, S. Perlaza, H. Tembine, and M. Debbah, “Electrical vehicles in the smart grid: A
mean field game analysis,” IEEE J. Selected Areas in Comm., vol. 30, no. 6, pp. 1086–1096,
2012.
[14] Z. Ma, D. Callaway, and I. Hiskens, “Decentralized charging control for large populations of
plug-in electric vehicles: Application of the Nash certainty equivalence principle,” in IEEE
International Conf. on Control Applications (CCA), 2010, pp. 191–195.
[15] A. Bušić and S. Meyn, “Distributed randomized control for demand dispatch,” March 2016,
submitted to IEEE Conference on Decision and Control, and arXiv:1603.05966v1.
[16] “BPA balancing authority,” Online, http://tinyurl.com/BPAgenload, http://tinyurl.com/
BPAbalancing.

16

