NREL is a national laboratory of the U.S. Department of Energy

Office of Energy Efficiency & Renewable Energy

Operated by the Alliance for Sustainable Energy, LLC

This report is available at no cost from the National Renewable Energy

Laboratory (NREL) at www.nrel.gov/publications.

Contract No. DE-AC36-08GO28308

Conference Paper

NREL/CP-5D00-76334

July 2020

Data-Driven Linear Parameter-Varying 

Modeling and Control of Flexible Loads 

for Grid Services

Preprint

Yue Chen and Andrey Bernstein

National Renewable Energy Laboratory

Presented at the 2020 American Control Conference (ACC) 

July 1–3, 2020 
NREL is a national laboratory of the U.S. Department of Energy

Office of Energy Efficiency & Renewable Energy

Operated by the Alliance for Sustainable Energy, LLC

This report is available at no cost from the National Renewable Energy

Laboratory (NREL) at www.nrel.gov/publications.

Contract No. DE-AC36-08GO28308

National Renewable Energy Laboratory

15013 Denver West Parkway

Golden, CO 80401

303-275-3000 • www.nrel.gov

Conference Paper

NREL/CP-5D00-76334

July 2020

Data-Driven Linear Parameter-Varying 

Modeling and Control of Flexible Loads 

for Grid Services

Preprint

Yue Chen and Andrey Bernstein 

National Renewable Energy Laboratory 

Suggested Citation 

Chen, Yue, and Andrey Bernstein. 2020. Data-Driven Linear Parameter-Varying Modeling 

and Control of Flexible Loads for Grid Services: Preprint. Golden, CO: National 

Renewable Energy Laboratory. NREL/CP-5D00-76334. 

https://www.nrel.gov/docs/fy20osti/76334.pdf. 
NOTICE

This work was authored in part by the National Renewable Energy Laboratory, operated by Alliance for Sustainable 

Energy, LLC, for the U.S. Department of Energy (DOE) under Contract No. DE-AC36-08GO28308. This 

work was supported by the Laboratory Directed Research and Development (LDRD) Program at NREL. The 

views expressed herein do not necessarily represent the views of the DOE or the U.S. Government. 

This report is available at no cost from the National Renewable 

Energy Laboratory (NREL) at www.nrel.gov/publications. 

U.S. Department of Energy (DOE) reports produced after 1991 

and a growing number of pre-1991 documents are available 

free via www.OSTI.gov. 

Cover Photos by Dennis Schroeder: (clockwise, left to right) NREL 51934, NREL 45897, NREL 42160, NREL 45891, NREL 48097, 

NREL 46526.

NREL prints on paper that contains recycled content.
Data-Driven Linear Parameter-Varying Modeling and Control

of Flexible Loads for Grid Services

Yue Chen, Andrey Bernstein

Abstract— Flexible loads have great potential to improve the

electric grid’s flexibility a nd s tability. T o e ffectively control

large ensembles of heterogeneous loads, reliable models thereof

are required. This paper presents a data-driven modeling and

control approach to manage flexible l oads f or p roviding grid

services. We leverage a linear parameter-varying autoregressive

moving average (LPV-ARMA) model to describe the aggregate

load response, where the parameters in the model are used

to capture external environmental impacts (e.g., weather). A

gain-scheduling feedback controller is then developed to adapt

to environmental variations. This data-driven approach can

be easily applied to different types of loads in various environmental conditions. In addition to the ensemble controller,

distributed load controllers are designed to deliver grid services,

while maintaining the quality of service of inherent load tasks.

We demonstrate the work on the IEEE 37-node distribution

system for real-time power regulation services through control

of thermostatically controlled loads.

I. INTRODUCTION

Modern power systems require increased flexibility to

adapt to the widespread penetration of volatile renewable

energy resources. It has been widely acknowledged that the

flexibility i n electric l oads c an b e h arnessed t o p rovide grid

services. Compared to expensive and limited-lifetime batteries, demand-side loads are more economic and sustainable

[1]; however, the performance of demand response is largely

restricted by the nature of working patterns of each type

of load, and current demand response programs are mainly

applied to event-based load shedding services [2].

Extensive work in the literature explores the potential

of flexible l oads t o e nhance t he fl exibility an d st ability of

power systems [3]–[10]. The main idea is to transform the

flexibility f rom l oads t o t he g rid b y s hifting l oad power

consumption. Electric loads, such as plug-in electric vehicles

(EVs), heating, ventilation, and air conditioning (HVAC)

systems, data centers, residential thermostatically controlled

loads (TCLs), pool pumps, and agricultural irrigation systems

have sufficient fl exibility to su pport an d pr ovide se rvices to

the grid.

Y. Chen and A. Bernstein are with the Power Systems Engineering

Center, National Renewable Energy Laboratory, Golden, CO 80401, USA

{yue.chen, andrey.bernstein}@nrel.gov

Because the contribution from a single load is negligible

to the grid (especially for low-power residential loads), it

is important to develop aggregate load models that will be

used in demand response control schemes. An ensemble of

loads can be modeled as a dynamic system whose dynamics

depend on the internal load operation physics as well as the

external signal that is used to engage them. Most existing

works consider or approximate the aggregate system by a

linear time-invariant (LTI) model, e.g., [4], [11], [12], where

it is assumed that the load working environment is time

invariant. LTI systems have attractive properties that make

the controller design easier. In many situations, however, the

LTI model is not sufficient to capture the load response.

As an example, the power consumption of a water heater

is largely related to the incoming water temperature.

This paper proposes linear parameter-varying (LPV) modeling and control of the aggregate loads, where the parameters capture the environment variables that affect load behavior. The model is identified using a data-driven approach so

that the general framework can be applied to different loads

in various environments. The LPV model presents a set of

linear approximations for the system with respect to varying

parameters. As a result, the standard analysis tools for the

LTI systems are still applicable to the LPV system. In this

regard, we design a gain-scheduling feedback controller to

deal with the variability of the environmental parameters.

In addition to the aggregate controller, local load controllers need careful design to ensure that loads complete

their routine operation tasks while providing grid services.

The decision-making methods for local controllers can be

classified into two categories: centralized methods and distributed methods. The centralized methods usually need

full load information to make decisions for each individual

load. In contrast, the distributed methods allow each local

controller to make its own decision to balance the grid

objective and its load quality of service (QoS). It relieves

the computation and communications requirements in the

centralized methods and protects the local privacy. This paper

adopts the distributed approach for local load control. In

addition, we employ a randomized policy in the local control

to desynchronize load responses.

The proposed hierarchical control architecture is illustrated

in Fig. 1. At the upper level, there is a central controller at

a utility that receives grid measurements (e.g., frequency,

voltage) and computes a control common for flexible loads.

At the local level, we propose a parallel distributed control

architecture where each load has its own controller. Each

local controller is designed to balance the utility’s grid

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 1
Load 

1

Local Stochastic Decision Making

Power Grid

...

Grid Services

Measurements

Load 

2

Load 

N

Grid-Level Controller 

Utility

Fig. 1. Schematic diagram of the control architecture.

service and the local objective that maintains the QoS to the

load owner. To desynchronize loads and unlock flexibility, a

stochastic method is introduced to randomize load behaviors

so that two identical loads could arrive at different decisions.

In summary, the contributions of this paper are: 1) we

develop an LPV model to capture the aggregate load response through a data-driven approach; 2) we design a gainscheduling feedback controller based on the developed LPV

model; and 3) we develop distributed local controllers to

respond to the grid-level controller and maintain its own load

QoS within a strict bound.

The rest of the paper is organized as follows. Section II

introduces the load QoS metric and proposes a distributed

randomized local control. Section III presents the data-driven

LPV modeling and control techniques. Section IV presents

simulation results on the IEEE 37-node distribution system

with 500 residential homes. Section V concludes the paper

and discusses future research.

II. INDIVIDUAL LOAD QOS AND CONTROL

Under nominal operation, loads are controlled to complete

their routine tasks, e.g., an air conditioner keeps the room

temperature within a comfortable range. We use the concept

of QoS to represent the load work performance on its routine

tasks. Without loss of generality, the QoS for an individual

load can be expressed as [5]:

Qt =

X

Tf

τ=0

β

τ

q(Xt−τ , ut−τ ) (1)

Xt+1 =f(Xt, ut) (2)

where Xt represents the load state; ut is the control signal

applied on the load; the function q extracts QoS features

from the load state and control; and the QoS is defined as

the discounted sum of the q values, in which the discount

factor β ∈ [0, 1] and Tf defines the moving window size.

With proper selection of Tf , β, and q(X, u), the QoS

function Qt can be designed to represent a broad range of

load QoS. For example, if β = 0 and/or Tf = 0 (we use

the convention 0

0 = 1), then Qt = q(Xt, ut) can be used

to define an instant load QoS, such as the temperature of

a TCL. When Tf > 0 and 0 < β ≤ 1, QoS value Qt

considers historical trajectory of relevant quantities, such as

energy consumption and cycling behavior of an on/off load.

Generally, the load QoS is determined by an adjustable

setpoint. For example, the cleanliness of a swimming pool

is controlled by the setpoint of daily operation time, and the

room temperature is maintained by the temperature setpoint

of an air conditioner. In this paper, we assume that the

load QoS has some level of flexibility, such that Qt ∈ Qt

is considered good QoS. Typically, Qt represents a box

constraint: Qt = [Q ̄

t, Qt

].

Using (1), the constraint Qt ∈ Qt can be transformed to

the control variable ut ∈ Ut, such that

Ut :=







u :

X

Tf

τ=1

β

τ

q(Xt−τ , ut−τ ) + q(f(Xt, u)) ∈ Qt







(3)

Thus, each load has its real-time control variable constraints,

according to its predefined QoS flexibility.

We propose a randomized control policy to avoid load synchronization. Assuming a load received a control command

u ̄t from the grid operator, it is randomized by a uniform

distribution:

ut = PUt {u ̄t(1 + αNt)} (4)

where PUtis the closest point projection operator that

ensures good load QoS; the random variable Nt follows

the uniform distribution on [−1, 1]; and α > 0 is a small

constant.

Example: Thermostatically Controlled Loads

This paper considers TCLs as the example to demonstrate

the proposed approach. TCLs such as air conditioners, electric water heaters, and refrigerators are the major power

consumers in residential households. The power flexibility

is extracted from the thermal inertia in TCLs that allows for

a slight shift in the power consumption without a noticeable

change in the temperature. Under proper control, TCLs are

great demand-side resources to improve the grid’s flexibility

and stability without significant impact on load owners’

comfort.

The dynamics for one TCL can be approximated by an

RC model [4], [13]:

Tt+1 = aTt + (1 − a)(T

o

t − mtRPtrans) + wt (5)

where Tt is the object temperature; a ∈ (0, 1) is the factor

of inertia, defined as a = e

−τ/(CR)

, where τ > 0 is the time

step, and C and R are the TCL’s thermal capacitance and

resistance, respectively; T

o

t

is the outdoor temperature; wt is

the modeling noise; and P

trans is the TCL’s energy transfer

rate, which is positive for cooling and negative for heating.

The internal control variable mt equals to 1 when the TCL

is ON and 0 when the TCL is OFF. It is defined as follows

for a cooling device:

mt =







0, Tt < min(Ωt)

1, Tt > max(Ωt)

mt−1, otherwise

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 2
where Ωt is the TCL’s temperature deadband, which is

defined by the temperature setpoint T

set

t

and the deadband

width δ as Ωt := [T

set

t − δ/2, Tset

t + δ/2].

Each TCL has its nominal temperature setpoint that

reflects the load owner’s preference. We define the QoS

function as the deviation from the preferred temperature

setpoint T

set

•

:

Qt := T

set

t − T

set

•

(6)

A large |Qt| indicates bad QoS.

In this particular TCL application, we define the control

variable the same as the deviation from the nominal QoS

setpoint:

T

set

t

:= (

T

set

• − ut, for a cooling TCL

T

set

• + ut, for a heating TCL

(7)

The purpose of proposing two different interpretations for

the control signal is to ensure all TCLs respond in the

same direction towards a given control signal, i.e., a positive

(negative) u encourages power consumption (saving).

Based on (6) and (7), we obtain a simple expression for

QoS function (1):

Qt =

(

ut, for a cooling TCL

−ut, for a heating TCL

(8)

Therefore, the load QoS constraint can be easily transformed

to the control constraint to restrict the maximal temperature

setpoint deviation from the preferred value.

III. LPV SYSTEM MODELING AND CONTROL

In this section, we propose techniques to find an LPV

model for the grid-level load response, and we design the

controller for grid services. The grid operator needs only a

grid-level model that captures the aggregate load response,

without detailed individual dynamics. It largely reduces the

communications and computation requirements at the central

controller. We consider a data-driven approach to identify the

load response model so that it can be applied to various types

of loads.

To understand the load response over time, the inputoutput data are collected as a time series. The autoregressive

moving average (ARMA) model is a typical tool to model

time-series data. For a single-input-single-output system, the

ARMA model is given by

yt + a1yt−1 + a2yt−2 + · · · + amyt−m

= b1ut−1 + b2ut−2 + · · · + bnut−n

(9)

where ut and yt are system input and output at t, respectively; and {a1, . . . am, b1, . . . , bn} are the model coefficients. Here, the control signal ut is the command sent to

flexible loads, e.g., the change in temperature setpoint of

air conditioners; the output signal yt represents the power

system outputs, such as the feeder head power and node

voltages.

One benefit of using the ARMA model is that it can be

easily converted to frequency domain for analysis. By taking

the z-transform of (9), we obtain

H(z) = Y (z)

U(z)

=

b1z

−1 + b2z−2 + . . . bnz−n

1 + a1z−1 + a2z−2 + . . . amz−m

(10)

Because of the environmental variation, however, the

ARMA model cannot accurately represent the load response

in many situations. To address this issue, the ARMA model

is improved to include the impact from the external environment.

A. Linear Parameter-Varying Model

The following proposed method is motivated by [14].

In this paper, we consider multiple environment features

and propose coefficient-identification algorithms with a regularization term to prevent model overfitting. For practical

applications, we add the data normalization process so that

model coefficients are “fairly” treated.

Consider the following parameter-varying ARMA model

yt + a1(p)yt−1 + a2(p)yt−2 + · · · + am(p)yt−m

= b1(p)ut−1 + b2(p)ut−2 + · · · + bn(p)ut−n + ηt

(11)

where ηt is the modeling error; the vector p ∈ R

v

collects

v environmental features, and

ai(p) := ai,0 + ai,1p1 + · · · + ai,vpv, for i = 1, . . . , m

bi(p) := bi,0 + bi,1p1 + · · · + bi,vpv, for i = 1, . . . , n

The model (11) can be represented by a standard linear

regression form:

yˆt = θ

T φt−1 (12)

where yˆt is the estimate of system output; the matrix θ

collects all model coefficients to be identified; the vector φ

collects input-output data and environmental feature values.

The detailed representations of θ and φ can be found in

appendix.

Consider the loss function that penalizes the modeling

error and the L2 regularization term (standard ridge regression):

J(θ) = 1

N

X

N

t=1

1

2

(yt − yˆt)

2 +

λ

2

||θ||2

2

(13)

where λ is a weighting factor. It is not hard to find the

gradient

∇J(θ) = 1

N

X

N

t=1

(yt − θ

T φt−1)(−φt−1) + λθ.

Setting the gradient to zero, we obtain the least-square

solution:

θ =



1

N

φ ̄φ ̄T + λI−1

1

N

φ ̄y ̄ (14)

where φ ̄ = [φ0, φ1, . . . , φN−1] and y ̄ = [y1, y2, . . . , yN ]

T

.

Online Algorithm. For some applications, especially when

the system process is nonstationary, it is useful to update the

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 3
Fig. 2. Schematic of gain-scheduling feedback control.

model in an online fashion. Consider the loss function (13)

in one time step:

Jt =

1

2

(yt − yˆt)

2 +

γ

2

||θ||2

2

(15)

We obtain its gradient:

∇Jt = −(yt − yˆt)φt−1 + γθt−1

It motivates the following online stochastic gradient descent

algorithm:

et =yt − θ

T

t−1φt−1

θt =(1 − αγ)θt−1 + αetφt−1

This algorithm is also known as the on-line ridge regression.

Data Normalization. Because the environmental variables

are collected in different units with different ranges, data

normalization is required to ensure that all features are

weighted equally in the model identification. Typical normalization techniques include centering, min-max scaling,

standardization, and other transformations [15]. In many

applications, normalization also improves the running time

of data-driven optimization algorithms.

B. Gain-Scheduling Feedback Control

The LPV-ARMA model (11) can be viewed as a set

of LTI approximations with respect to the environmental

parameters. For fixed parameters, the model is simply LTI.

So the tools for LTI system analysis and controller design

are still applicable. Different from fixed controller gains for

an LTI system, the gains for an LPV system are adaptable

to environmental parameters p.

The structure for the gain-scheduling feedback control

is given in Fig. 2. Grid services are represented by the

reference signal, which is expected to be followed by the

distribution system. Based on the error between the reference

signal and the system outputs, the controller computes the

control signal for loads. We assume that the environmental

parameters (such as the weather data) are measurable so

that the controller can adjust the gains to adapt to the

environment.

In this paper, we apply classic frequency domain techniques (e.g., [16]) to analyze the system and design the controller. For a discrete-time system, the closed-loop transfer

function can be expressed in z-domain as

H(z, p) = Y (z, p)

R(z, p)

=

C(z, p)P(z, p)

1 + C(z, p)P(z, p)

(16)

where Y (z, p), R(z, p), C(z, p), P(z, p) are the z-domain

representations of the output, reference signal, controller, and

system plant, respectively, with respect to the environmental

parameter p.

It is suggested by (16) that Y (z, p)/R(z, p) ≈ 1 can be

achieved by choosing high gains in C(z, p); however, the

system would become unstable if the gains cross a certain

critical point. Note that the closed-loop system (16) has

infinite gain when C(z, p)P(z, p) = −1, which is equivalent

to |C(z, p)P(z, p)| = 1 and ∠C(z, p)P(z, p) = −180o.

Two commonly used stability criteria are the gain margin and

phase margin, which define the amount of gain and phase lag

that can be increased before instability results.

Therefore, the controller depends on the dynamics of the

grid-level loads response P(z, p), which is the z-transform

of the LPV model (11). The details of LPV modeling and

controller design are illustrated in Section IV.

IV. NUMERICAL DEMONSTRATION

In this section, we demonstrate the data-driven approach

for real-time power regulation through control of TCLs. The

numerical results described here are based on simulations

of 500 residential households on the balanced IEEE 37-node

distribution system [17]. It was observed from the simulation

results that a collection of residential TCLs can provide highquality grid services, with little impact on load QoS.

A. Data Acquisition

We simulated 500 residential households, where each

household has an air conditioner, a refrigerator, and a heat

pump. In total, 1, 500 TCLs were modeled by (5) with

the parameters provided in Table I, where U refers to the

uniform distribution. These 500 households were distributed

throughout the load nodes in the IEEE 37-node test feeder.

The aggregate load response was modeled using the LPV

model (11). Because the air conditioners and heat pumps do

not operate at the same time, we developed two different

LPV models for hot and cold weather, respectively. For

brevity, we show only the LPV modeling and controller

design process under hot weather, when air conditioners

and refrigerators are engaged. The modeling and control

techniques used for cold weather are the same.

To mimic a hot summer environment, we created a 24-hour

outdoor temperature profile, as shown on the top in Fig. 3.

The corresponding power at the feeder head is plotted on

the bottom in Fig. 3. It is found that the feeder head power

is closely related with the outdoor temperature because the

power is dominated by air conditioners and the contribution

from refrigerators is insignificant.

To collect input-output data for system identification, a

swept-sine signal (top in Fig. 4) was created to change

the TCL temperature setpoints. It allows us to observe the

system response to various frequencies of the input signal.

Under the same temperature in Fig. 3, the distribution system

feeder head power was plotted in the middle of Fig. 4. Its

deviation from the nominal power consumption in Fig. 3

was considered as the system output, which is shown on the

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 4
TABLE I

TCL PARAMETER VALUES.

Parameter Air Conditioners Heat Pumps Refrigerators

Dead-band width, δ(

◦C) ∼ U[0.5 1] ∼ U[0.5 1] ∼ U[0.5 1]

Temperature set point, T

set (◦C) ∼ U[23 25] ∼ U[19 21] ∼ U[1.7 3.3]

Thermal resistance, R(

◦C/kW) ∼ U[1.5 2.5] ∼ U[1.5 2.5] ∼ U[60 80]

Thermal capacitance, C(kWh/

◦C) ∼ U[1.5 2.5] ∼ U[1.5 2.5] ∼ U[.05 0.15]

Energy transfer rate, P

trans(kW) ∼ U[12 16] ∼ U[25 30] ∼ U[0.6 0.9]

Coefficient of performance 2.5 3.5 3

0 6 12 18 24

10

20

30

40

50

oC

Outdoor temperature

0 6 12 18 24

Time / hour

0

2

4

6

MW

Feeder head real power - nominal behavior (u=0)

Fig. 3. The outdoor temperature and the nominal feeder head power.

0 6 12 18 24

-0.4

-0.2

0

0.2

0.4

oC

Control input u

0 6 12 18 24

0

2

4

6

MW

Feeder head real power - with control u

0 6 12 18 24

Time / hour

-4

-2

0

2

4

MW

Feeder head real power deviation - with control u

Fig. 4. Feeder head power demand under temperature setpoint control, and

its deviation from nominal power in Fig. 3.

bottom in Fig. 4. The goal is to identify the LPV model using

input-output simulation data along with the temperature data.

B. System Identification

We used the following two polynomial features:

p1,t =

T

o

t − p ̄

σ1

, p2,t =

(T

o

t − p ̄)

2

σ2

where T

o

t

is the outdoor temperature; p ̄ = 20 oC; and σ1

and σ2 are standard deviations for the quantities in the

numerators, respectively.

0 6 12 18 24

-4

-2

0

2

4

MW

Power deviation - Modeling result

0 6 12 18 24

Time / hour

-4

-2

0

2

4

MW

Modeling error

RMSE = 0.02203

Fig. 5. Modeling result and error.

A 10-day simulation was performed to collect the inputoutput data. For brevity, we present only the first 24-hour

results in 10 days. The top subplot in Fig. 5 shows the LPV

modeling results (11), with modeling orders m = 31 and

n = 31. It closely matches the real system output that shown

on the bottom in Fig. 4, and the modeling error is given on

the bottom in Fig. 5.

C. Ensemble Controller Design

We designed the controller in the frequency domain. Fig.

7 shows the frequency response of the system with various

outdoor temperatures. When the temperature is 25 oC, because most air conditioners are idling, a slight change in

the temperature setpoint has little impact on the power consumption. So the response magnitude is very small when the

outdoor temperature is 25 oC. As the temperature increases,

air conditioners start to actively respond to the temperature

setpoint change, and they eventually saturate around 35 oC.

The resonant frequencies in these bode plots are in the

neighborhood of 5 × 10−2rad/s, which is consistent with

the nominal on-off cycle frequencies of the air conditioners

in simulation.

Based on the aggregate load response shown in Fig. 7 and

the proposed control structure shown in Fig. 2, we designed

the following gain-scheduling PI controller:

ut(T

o

t

) = Kp(T

o

t

)et + Ki(T

o

t

)

Xt

τ=0

eτ∆t

with proportional and integral gains

Kp(T

o

t

) = 70 − T

o

t

50

, Ki(T

o

t

) = 0.05Kp(T

o

t

)

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 5
-30

-20

-10

0

10

20

Magnitude (dB)

10-5 10-4 10-3 10-2 10-1 100

-270

-180

-90

0

90

Phase (deg)

25 oC 30 oC 35 oC 40 oC

Bode Diagram

Frequency (rad/s)

Fig. 7. Frequency response of the feeder head power, driven by the control

of TCL temperature setpoint.

where et is the error between the reference and the actual

system output; ∆t = 5 seconds is the simulation time step;

and T

o

t

is the outdoor temperature. The controller provides

gain margins ≥ 10.4 dB and phase margins ≥ 52.8

o

for

system stability.

The resulting closed-loop frequency response is presented

in Fig. 8. It indicates that these flexible loads can track

reference signals of periods as fast as a few minutes.

For the purpose of maintaining good load QoS, we restricted the control signal u within at most 1

oC deviation

from the nominal TCL temperature setpoints. In addition, we

restricted the integral error e

I

t within predetermined bounds

to prevent the integral windup situation in PI control.

D. Grid Service

To illustrate the control performance under a wide range

of outdoor temperatures, we created an artificial 48-hour

-40

-30

-20

-10

0

10

Magnitude (dB)

10-4 10-3 10-2 10-1 100

-270

-225

-180

-135

-90

-45

0

Phase (deg)

25 oC 30 oC 35 oC 40 oC

Bode Diagram

Frequency (rad/s)

Fig. 8. Frequency response of the closed-loop system, using the gainscheduling PI control.

temperature profile that spans [−20, 40] oC to engage all

considered TCLs. The temperature is shown on the top in

Fig. 6. The left-hand-side (LHS) in Fig. 6 (the first 24-hour

period) shows the nominal behavior of the TCLs without

control; the right-hand-side (RHS) in Fig. 6 (the second 24-

hour period) shows the control performance when the TCLs

were controlled to provide the grid service.

Based on the nominal power at the feeder head (on the

LHS of Subplot 2 in Fig. 6), we obtain the baseline of 24-

hour power, which is shown as the dotted line on RHS of the

Subplot 2. The goal is to control the TCLs so that the power

deviation from the baseline tracks some reference signal. We

consider the regulation signal from the Bonneville Power

Administration (BPA) [18] as the reference signal. Because

of the capacity limit, the BPA signal was scaled to a proper

magnitude that is trackable by the loads. It is shown as the

Fig. 6. Control performance of the proposed approach in providing power-balancing grid service through control of the TCL temperature setpoint.

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 6
red dashed line on the RHS of Subplot 3 in Fig. 6.

It is found that these TCLs can provide real-time regulation power in the range of [−1, 1] MW with power tracking

root-mean-square-error (RMSE) of 22 kW (see the RHS of

Subplot 3). Considering that the aggregate TCL power is

less than 4 MW, the regulation service extracted from these

TCLs is significant. In addition, note that the fluctuations that

appear in nominal power consumption are also suppressed

under the control.

While providing grid service, the impact to home comfort

is insignificant because the control signal is restricted to

[−1, 1] oC. It is verified from the comparison of two 24-hour

house temperatures on the bottom in Fig. 6. To conclude,

there is a trade-off between the grid service and load QoS:

if these TCLs allow more QoS flexibility, a greater regulation

service can be provided.

V. CONCLUSION AND FUTURE WORK

In this paper, we presented a data-driven LPV modeling

and control approach for flexible loads to provide grid

services. The approach enables the controller to tune its gains

to adapt to changes in the environment. It was demonstrated

on the IEEE 37-node distribution system for real-time power

regulation service through control of TCLs.

In the future, we plan to extend the application to a wider

range of loads that involve more environmental parameters.

Other grid services, e.g, voltage regulation, are also of

interest for future work. In addition, it is valuable to estimate

the grid service capacity from the loads. We aim to estimate

the capacity by a similar data-driven approach and use the

capacity information for power system planning.

ACKNOWLEDGEMENTS

We thank Emiliano Dall’Anese (CU Boulder) and Marcello Colombino (ETH) for valuable discussions.

REFERENCES

[1] N. J. Cammardella, R. W. Moye, Y. Chen, and S. P. Meyn, “An energy

storage cost comparison: Li-ion batteries vs distributed load control,”

in 2018 Clemson University Power Systems Conference (PSC), Sep.

2018, pp. 1–6.

[2] “Florida Power and Light: OnCall Program,” https://www.fpl.com/

save/programs/on-call.html.

[3] D. S. Callaway and I. A. Hiskens, “Achieving controllability of electric

loads,” Proceedings of the IEEE, vol. 99, no. 1, pp. 184–199, Jan 2011.

[4] J. L. Mathieu, S. Koch, and D. S. Callaway, “State estimation and

control of electric loads to manage real-time energy imbalance,” IEEE

Transactions on Power Systems, vol. 28, no. 1, pp. 430–440, Feb 2013.

[5] Y. Chen, A. Busic, and S. Meyn, “Estimation and control of quality

of service in demand dispatch,” IEEE Transactions on Smart Grid,

vol. 9, no. 5, pp. 5348–5356, Sep. 2018.

[6] Y. Chen, M. U. Hashmi, J. Mathias, A. Busiˇ c, and S. Meyn,  ́

Distributed Control Design for Balancing the Grid Using Flexible

Loads. New York, NY: Springer New York, 2018, pp. 383–411.

[Online]. Available: https://doi.org/10.1007/978-1-4939-7822-9 16

[7] H. Hao, Y. Lin, A. S. Kowli, P. Barooah, and S. Meyn, “Ancillary

service to the grid through control of fans in commercial building

hvac systems,” IEEE Transactions on Smart Grid, vol. 5, no. 4, pp.

2066–2074, July 2014.

[8] M. Almassalkhi, J. Frolik, and P. Hines, “Packetized energy management: Asynchronous and anonymous coordination of thermostatically

controlled loads,” in 2017 American Control Conference (ACC), May

2017, pp. 1431–1437.

[9] S. P. Meyn, P. Barooah, A. Busic, Y. Chen, and J. Ehren, “Ancillary

service to the grid using intelligent deferrable loads,” IEEE Transactions on Automatic Control, vol. 60, no. 11, pp. 2847–2862, Nov

2015.

[10] 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 2010 IEEE International

Conference on Control Applications, Sep. 2010, pp. 191–195.

[11] Y. Chen, A. Busic, 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.

[12] Y. Chen, A. Busic, and S. P. Meyn, “State estimation for the individual

and the population in mean field control with application to demand

dispatch,” IEEE Transactions on Automatic Control, vol. 62, no. 3,

pp. 1138–1149, March 2017.

[13] D. S. Callaway, “Tapping the energy storage potential in electric

loads to deliver load following and regulation, with application to

wind energy,” Energy Conversion and Management, vol. 50, no. 5,

pp. 1389 – 1400, 2009. [Online]. Available: http://www.sciencedirect.

com/science/article/pii/S0196890408004780

[14] B. Bamieh and L. Giarre, “Identification of linear parameter

varying models,” International Journal of Robust and Nonlinear

Control, vol. 12, no. 9, pp. 841–853, 2002. [Online]. Available:

http://dx.doi.org/10.1002/rnc.706

[15] S. Shalev-Shwartz and S. Ben-David, Feature Selection and Generation. Cambridge University Press, 2014, p. 309–322.

[16] G. F. Franklin, J. D. Powell, and A. Emami-Naeini, Feedback Control

of Dynamic Systems, 7th ed. Upper Saddle River, NJ, USA: Prentice

Hall Press, 2014.

[17] Distribution System Analysis Subcommittee, IEEE Power Engineering

Society, “IEEE 37 Node Test Feeder.” [Online]. Available: http:

//sites.ieee.org/pes-testfeeders/resources/

[18] “BPA balancing authority load and total wind, hydro, and thermal generation (website),” https://transmission.bpa.gov/business/operations/

Wind/baltwg.aspx.

APPENDIX I

LINEAR REPRESENTATION OF THE LPV MODEL

We represent the LPV model (11) in the following form

yˆt = trace(ΘT Ψt−1) (17)

with

Θ

T

:=











a1,0 . . . am,0 b1,0 . . . bn,0

.

.

.

.

.

.

.

.

.

.

.

. . . .

.

.

.

a1,v . . . am,v b1,v . . . bn,v











Ψt−1 :=





























−yt−1

.

.

.

−yt−m

ut−1

.

.

.

ut−n





























h

1 p1 . . . pv

i

.

We introduce a vectorization operator vec(·) that transforms a matrix to a column vector by stacking matrix

columns, e.g.,

vec









a1,1 a1,2

a2,1 a2,2







 =

h

a1,1 a2,1 a1,2 a2,2

iT

Letting θ = vec(Θ) and φt−1 = vec(Ψt−1), we obtain (12).

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 7