Distributed control design for balancing the grid using flexible - - PowerPoint PPT Presentation

Distributed control design for balancing the grid using flexible loads

YEQT XI: ”Winterschool on Energy Systems” Dec 11 - Dec 15 2017

Ana Buˇ si´ c Inria Paris D´ epartement d’Informatique de l’ENS Joint work with Prabir Barooah, Yue Chen, and Sean Meyn University of Florida

Motivation and challenges

Motivation: Climate change

Carbon dioxide concentrations are rising.

1769

Measurements of the CO2 concentration in the air from the year 1000AD to the present.

(Source: Fig. 1.4. from David MacKay, Sustainable Energy – without the hot air. 2009, and www.withouthotair.com)

Motivation and challenges

Motivation: Climate change

Carbon dioxide concentrations are rising.

1769

Measurements of the CO2 concentration in the air from the year 1000AD to the present.

(Source: Fig. 1.4. from David MacKay, Sustainable Energy – without the hot air. 2009, and www.withouthotair.com)

Industrial Revolution 1769 — the year in which James Watt patented his steam engine.

Motivation and challenges

Motivation: Climate change

Gigatons of CO2 per year into the atmosphere

1769

• Fig. 1.4

Motivation and challenges

Motivation: Climate change

Gigatons of CO2 per year into the atmosphere Up until 1769 there was a balance: Biosphere = 440 Oceans = 330 These 770 gigatons were absorbed by the planet

1769

• Fig. 1.4

Motivation and challenges

Motivation: Climate change

Gigatons of CO2 per year into the atmosphere Up until 1769 there was a balance: Biosphere = 440 Oceans = 330 These 770 gigatons were absorbed by the planet Current excess from burning fossil fuels: add 26

1769

• Fig. 1.4

Motivation and challenges

Motivation: Climate change

Gigatons of CO2 per year into the atmosphere Up until 1769 there was a balance: Biosphere = 440 Oceans = 330 These 770 gigatons were absorbed by the planet Current excess from burning fossil fuels: add 26 Small contribution in percentage terms, but the input-output balance is no longer maintained.

1769

• Fig. 1.4

Motivation and challenges

Motivation: Climate change

Current excess from burning fossil fuels: add 26 Small contribution in percentage terms, but the input-output balance is no longer maintained.

1769

• Fig. 1.4

Motivation and challenges

Motivation: Climate change

Current excess from burning fossil fuels: add 26 Small contribution in percentage terms, but the input-output balance is no longer maintained. The earth is a complex dynamical system – the exact details of the impact are not certain.

1769

• Fig. 1.4

Motivation and challenges

Motivation: Climate change

Current excess from burning fossil fuels: add 26 Small contribution in percentage terms, but the input-output balance is no longer maintained. The earth is a complex dynamical system – the exact details of the impact are not certain.

1769

• Fig. 1.4

IPCC Report: ‘severe and pervasive’ impacts of climate change will be felt

• everywhere. March 31, 2014

The Intergovernmental Panel on Climate Change (IPCC) today issued a report saying the effects of climate change are already occurring on all continents and across the oceans. The world, in many cases, is ill-prepared for risks from a changing climate.

Motivation and challenges

Motivation: Climate change

More from UN

Motivation and challenges

Motivation: Climate change

Responsibility for pollution: The width of each rectangle is proportional to population, and the height to CO2 contribution per person

(average emission rate over the period 1880–2004).

Link to p. 14 of monograph by MacKay

Motivation and challenges

Challenges to sustainability

The red stack in figure 18.1 adds up to 195 kWh per day per person (based on a typical British citizen). The green stack adds up to about 180 kWh/d/p. A close race!

• Fig. 18.1

Motivation and challenges

Challenges to sustainability

The red stack in figure 18.1 adds up to 195 kWh per day per person (based on a typical British citizen). The green stack adds up to about 180 kWh/d/p. A close race! Not taking into account all economic, social, and environmental constraints.

• Fig. 18.1

Motivation and challenges

Challenges

There is not enough of it! Please also remember: we get smarter every year!

http://c1cleantechnicacom.wpengine.netdna-cdn.com/files/2014/02/solar-cell-efficiency-records.jpg

Motivation and challenges

Challenges of renewable power generation

Impact of wind and solar on net-load at CAISO Ramp limitations cause price-spikes

Price spike due to high net-load ramping need when solar production ramped out Negative prices due to high mid-day solar production

1200 15 2 4 19 17 21 23 27 25 800 1000 600 400 200

• 200

GW GW Toal Load

Load and Net-load Renewable Generation

Total Wind Net-load: Toal Load, less Wind and Solar $/MWh 24 hrs 24 hrs Peak ramp Peak

Peak ramp Peak

Total Solar

Motivation and challenges

Challenges of renewable power generation

Impact of wind and solar on net-load at CAISO Ramp limitations cause price-spikes

Price spike due to high net-load ramping need when solar production ramped out Negative prices due to high mid-day solar production

1200 15 2 4 19 17 21 23 27 25 800 1000 600 400 200

• 200

GW GW Toal Load

Load and Net-load Renewable Generation

Total Wind Net-load: Toal Load, less Wind and Solar $/MWh 24 hrs 24 hrs Peak ramp Peak

Peak ramp Peak

Total Solar

Jan 01 Jan 02 Jan 03 Jan 04 Jan 05 Jan 06 GW 1 2 3 4

GW (t) = Wind generation in BPA, Jan 2015

Ramps

Motivation and challenges

Challenges: regulation

2 4 6 8 2 4 6 8 0.8

• 0.8

1

• 1

0.8

• 0.8

1

• 1

Sun Mon Tue October 20-25 October 27 - November 1 Hydro Wed Thur Fri Sun Mon Tue Wed Thur Fri

Generation and Laod GW

GW GW

Regulation GW

Thermal Wind Load Generation Regulation

Error Signal in Feedback Loop

Motivation and challenges

Comparison: Flight control

How do we fly a plane through a storm?

Brains

Brawn

Motivation and challenges

Comparison: Flight control

How do we fly a plane through a storm?

Brains

Brawn

What Good Are These?

Motivation and challenges

Comparison: Flight control

How do we operate the grid in a storm? Balancing Authority Ancillary Services Grid

Measurements: Voltage Frequency Phase

Σ −

Brains Brawn What Good Are These?

Motivation and challenges

Challenges of renewable power generation

Increasing needs for ancillary services

20 40 60 80 100 120 140 160 t/hour Reference (from Balancing Authority)

Balancing Authority Ancillary Services Grid Voltage Frequency Phase

Σ −

In the past, provided by the generators - high costs!

Motivation and challenges

Generator Outage

Three or Four Phases to Recovery

Frequency Power deviation 50Hz

Fault Fault+20sec Fault+30mins

Inertia Primary (Governor) Secondary Tertiary Primary (Governor) Secondary Tertiary

Motivation and challenges

Generator Outage

Three or Four Phases to Recovery

Frequency Power deviation 50Hz Fault Fault+20sec Fault+30mins Inertia Primary (Governor) Secondary Tertiary Primary (Governor) Secondary Tertiary

• Primary Frequency Reserve

Based on a local control loop using local system frequency measurements Stabilize frequency in case of major outages of either loads or suppliers. The time scale for activating primary frequency reserve is in the area of 10-30 seconds. Traditionally governor response

Motivation and challenges

Generator Outage

Three or Four Phases to Recovery

Frequency Power deviation 50Hz Fault Fault+20sec Fault+30mins Inertia Primary (Governor) Secondary Tertiary Primary (Governor) Secondary Tertiary

• Primary Frequency Reserve
• Secondary Frequency Reserve

Activated by a BA signal. Restore power balance in a control area, to take part in stabilizing the frequency, and restore the primary reserve. The time scale for activation of secondary reserve is in the magnitude of minutes.

Motivation and challenges

Generator Outage

Three or Four Phases to Recovery

Frequency Power deviation 50Hz Fault Fault+20sec Fault+30mins Inertia Primary (Governor) Secondary Tertiary Primary (Governor) Secondary Tertiary

• Primary Frequency Reserve
• Secondary Frequency Reserve

Activated by a BA signal. Restore power balance in a control area, to take part in stabilizing the frequency, and restore the primary reserve. The time scale for activation of secondary reserve is in the magnitude of minutes. Example: Automatic Generation Control (AGC).

Motivation and challenges

Generator Outage

Three or Four Phases to Recovery

Frequency Power deviation 50Hz Fault Fault+20sec Fault+30mins Inertia Primary (Governor) Secondary Tertiary Primary (Governor) Secondary Tertiary

• Primary Frequency Reserve
• Secondary Frequency Reserve
• Tertiary Frequency Reserve

Activated by a BA signal. Upon activation, the provider of the reserve will change the planned operation such that the necessary upward or downward regulation is

• achieved. The purpose of tertiary reserve is to resolve persistent balance or

congestion problems and in this way restore the secondary and primary frequency reserve.

Motivation and challenges

Generator Outage

Three or Four Phases to Recovery

Frequency Power deviation 50Hz Fault Fault+20sec Fault+30mins Inertia Primary (Governor) Secondary Tertiary Primary (Governor) Secondary Tertiary

• Primary Frequency Reserve
• Secondary Frequency Reserve
• Tertiary Frequency Reserve

Crisis-aversion and reallocation of resources: just one role for the grid operator.

Motivation and challenges

Ancillary Services

... to compensate for energy imbalances ... Kirby 2013

2 4 6 8 2 4 6 8 0.8

• 0.8

1

• 1

0.8

• 0.8

1

• 1

Sun Mon Tue October 20-25 October 27 - November 1 Hydro Wed Thur Fri Sun Mon Tue Wed Thur Fri

Generation and Laod GW

GW GW

Balancing GW

Thermal Wind Load Generation Balancing Reserves Deployed

~ Error Signal in Feedback Loop

Motivation and challenges

Ancillary Services

... to compensate for energy imbalances ... Kirby 2013

2 4 6 8 2 4 6 8 0.8

• 0.8

1

• 1

0.8

• 0.8

1

• 1

Sun Mon Tue October 20-25 October 27 - November 1 Hydro Wed Thur Fri Sun Mon Tue Wed Thur Fri

Generation and Laod GW

GW GW

Balancing GW

Thermal Wind Load Generation Balancing Reserves Deployed

~ Error Signal in Feedback Loop

Time-scale of power deviations are similar to secondary reserves following a fault

Motivation and challenges

Ancillary Services

... to compensate for energy imbalances ... Kirby 2013

2 4 6 8 2 4 6 8 0.8

• 0.8

1

• 1

0.8

• 0.8

1

• 1

Sun Mon Tue October 20-25 October 27 - November 1 Hydro Wed Thur Fri Sun Mon Tue Wed Thur Fri

Generation and Laod GW

GW GW

Balancing GW

Thermal Wind Load Generation Balancing Reserves Deployed

~ Error Signal in Feedback Loop

Time-scale of power deviations are similar to secondary reserves following a fault The Balancing Reserves at BPA are a sum of many error signals in the grid

Motivation and challenges

NERC Coordinating Councils

FRCC Florida Reliability Coordina ng Council MRO Midwest Reliability Organiza on NPCC Northeast Power Coordina ng Council RFC ReliabilityFirst Corpora on SERC SERC Reliability Corpora on SPP RE Southwest Power Pool Regional En ty TRE Texas Reliability En ty WECC Western Electric Coordina ng Council

Interconnec on Requirements for Variable Genera on - September 2012

NERC Regional Entitites From Operating Manual – North American Energy Standards Board (NAESB)

Motivation and challenges

NERC Coordinating Councils

FRCC Florida Reliability Coordina ng Council MRO Midwest Reliability Organiza on NPCC Northeast Power Coordina ng Council RFC ReliabilityFirst Corpora on SERC SERC Reliability Corpora on SPP RE Southwest Power Pool Regional En ty TRE Texas Reliability En ty WECC Western Electric Coordina ng Council

Interconnec on Requirements for Variable Genera on - September 2012

NERC Regional Entitites [FERC] granted NERC the legal authority to enforce reliability standards with all U.S. users, owners, and operators of the bulk power system, and made compliance with those standards mandatory and enforceable. Reliability Standards are also mandatory and enforceable in Ontario and New Brunswick, and NERC is seeking to achieve comparable results in the other Canadian provinces. NERC will seek recognition in Mexico once necessary legislation is adopted.

Motivation and challenges

NERC Coordinating Councils

FRCC Florida Reliability Coordina ng Council MRO Midwest Reliability Organiza on NPCC Northeast Power Coordina ng Council RFC ReliabilityFirst Corpora on SERC SERC Reliability Corpora on SPP RE Southwest Power Pool Regional En ty TRE Texas Reliability En ty WECC Western Electric Coordina ng Council

Interconnec on Requirements for Variable Genera on - September 2012

NERC Regional Entitites [FERC] granted NERC the legal authority to enforce reliability standards with all U.S. users, owners, and operators of the bulk power system, and made compliance with those standards mandatory and enforceable. Reliability Standards are also mandatory and enforceable in Ontario and New Brunswick, and NERC is seeking to achieve comparable results in the other Canadian provinces. NERC will seek recognition in Mexico once necessary legislation is adopted.

NERC = North American Electric Reliability Corporation

Motivation and challenges

NERC Coordinating Councils

Don’t forget FERC’s ISOs & RTOs!

Motivation and challenges

Secondary Control

Balancing Authority

Figure 2 North American Balancing Authorities and Regions

Motivation and challenges

Secondary Control

Balancing Authority

Figure 2 North American Balancing Authorities and Regions

Transmission lines that join two areas are known as tie-lines.

Motivation and challenges

Secondary Control

Balancing Authority

Transmission lines that join two areas are known as tie-lines. The net power out of an area is the sum of the flow on its tie-lines.

Motivation and challenges

Secondary Control

Balancing Authority

Transmission lines that join two areas are known as tie-lines. The net power out of an area is the sum of the flow on its tie-lines. The flow out of an area is equal to

1 T

P

2 T

P

Area 1 G G G G G G G Area 2

total gen - total load - total losses = tie-line flow Remember Economic Dispatch?

Motivation and challenges

Primary and Secondary Control Loops

Turbine Governor Valves or Gates

Turbine

Generator

The Grid

Loads Transmission lines Other generators Speed Balancing Authority

(ACE calculation)

Automatic Generation Control (AGC) (Secondary control loop) Setpoint Calculation Power, Frequency Tie-Line Powers (Internal turbine control loop) (Primary control loop) Setpoint from power generation schedule and tertiary control

ω identical for synchronous

machines in steady state System Frequency: ω Turbine Controller Load Shedding (emergency control)

ETH Dynamics 2012 Regulation Signals

Motivation and challenges

Secondary Control

Area Control Error

Area Control Error: combination of: Deviation of frequency from nominal, and

Motivation and challenges

Secondary Control

Area Control Error

Area Control Error: combination of: Deviation of frequency from nominal, and the difference between the actual tie-line flow, and the scheduled flow (from economic dispatch).

Motivation and challenges

Secondary Control

Area Control Error

Area Control Error: combination of: Deviation of frequency from nominal, and the difference between the actual tie-line flow, and the scheduled flow (from economic dispatch). ACE = Pactual − Pscheduled + B∆ω B is the bias

Motivation and challenges

Secondary Control

Area Control Error

Area Control Error: combination of: Deviation of frequency from nominal, and the difference between the actual tie-line flow, and the scheduled flow (from economic dispatch). ACE = Pactual − Pscheduled + B∆ω B is the bias Provides a measure of whether an area is producing more or less than it should to satisfy schedules and to contribute to controlling frequency. AGC: control signal designed to bring ACE to zero.

Motivation and challenges

Secondary Control

Area Control Error

ETH Dynamics 2012

1 T

P

2 T

P

1

1 1

Tie line power for Area 1 Sum over all tie lines

=

j T T j

P P

Area 1 AGC2 AGC1

1

ω1 ω2

G G G G G G G Area 2 ACE

2

ACE

Motivation and challenges

Secondary Control

Regulation signals

Regulation: on-line generation, responsive load and storage ... helps to maintain interconnection frequency, manage differences between actual and scheduled power flows between balancing areas, and match generation to load within the balancing area. Automatic Generation Control (AGC): commands are typically sent about every four seconds.

PJM RegD Measured

Power

(kW)

1

• 1

10 20 30 40 Time (minute)

Motivation and challenges

Secondary Control

Regulation signals

Regulation: on-line generation, responsive load and storage ... helps to maintain interconnection frequency, manage differences between actual and scheduled power flows between balancing areas, and match generation to load within the balancing area. Automatic Generation Control (AGC): commands are typically sent about every four seconds.

PJM RegD Measured

Power

(kW)

1

• 1

10 20 30 40 Time (minute)

Value is now recognized: FERC Orders

Motivation and challenges

Secondary Control

Regulation signals

Regulation: on-line generation, responsive load and storage ... helps to maintain interconnection frequency, manage differences between actual and scheduled power flows between balancing areas, and match generation to load within the balancing area. Automatic Generation Control (AGC): commands are typically sent about every four seconds.

PJM RegD Measured

Power

(kW)

1

• 1

10 20 30 40 Time (minute)

Value is now recognized: FERC Orders PJM regulation metrics: scores for correlation, precision, and performance. Faster and more accurate response is paid more.

Motivation and challenges

Example: ACE and regulation signals at ERCOT

Power (MW)

13:00 15:00 17:00 19:00 21:00 23:00 01:00 03:00 05:00 07:00 09:00 11:00 1000 800 600 400 200

• 200
• 400
• 600
• 800

Area Control Error Deployed Regulation

Motivation and challenges

Refinement of our Control Loop

K

+ I

KP 1 s

Grid disturbance ACE Measured

GRID BA

+

• ref = 0

ρ(t)

ACE = Pactual − Pscheduled + B∆ω

Motivation and challenges

Refinement of our Control Loop

K

+ I

KP 1 s

Grid disturbance ACE Measured

GRID BA

+

• ref = 0

ρ(t)

ACE = Pactual − Pscheduled + B∆ω

The signal ACE is the sum of two errors.

Motivation and challenges

Refinement of our Control Loop

K

+ I

KP 1 s

Grid disturbance ACE Measured

GRID BA

+

• ref = 0

ρ(t)

ACE = Pactual − Pscheduled + B∆ω

The signal ACE is the sum of two errors. The BA signal ρ(t) is obtained from filtering ACE through Gc ≡ PI-compensator.

Motivation and challenges

Does ACE Work?

Generators are not always good actuators

• 15
• 10
• 5

5 10 15 20 25 30 35 6:00 7:00 8:00 9:00 6:00 7:00 8:00 9:00 Regulation (MW)

Gen 'A' Actual Regulation Gen 'A' Requested Regulation

• 20
• 15
• 10
• 5

5 10 15 20 : Regulation (MW)

Gen 'B' Actual Regulation Gen 'B' Requested Regulation

• Fig. 10. Coal-fired generators do not follow regulation signals precisely....

Some do better than others

Demand Dispatch

Tracking Grid Signal with Residential Loads

Example: 20 pools, 20 kW max load

Each pool consumes 1kW when operating 12 hour cleaning cycle each 24 hours Power Deviation:

20 40 60 80 100 120 140 160 10

• 10

t/hour

kW 20 pools

Input ζt

• 3

3

Output deviation Reference

Nearly Perfect Service from Pools

Meyn, Barooah, B., Chen, Ehren 2015 [IEEE TAC]

using an extension/reinterpretation of Todorov 2007 [NIPS] (linearly solvable MDPs)

Demand Dispatch

Tracking Grid Signal with Residential Loads

Example: 300,000 pools, 300 MW max load

Each pool consumes 1kW when operating 12 hour cleaning cycle each 24 hours Power Deviation:

20 40 60 80 100 120 140 160 t/hour

Output deviation Reference

Input ζt

• 3

3 −100 −50 50 100

MW 300,000 pools

Nearly Perfect Service from Pools

Meyn, Barooah, B., Chen, Ehren 2015 [IEEE TAC]

using an extension/reinterpretation of Todorov 2007 [NIPS] (linearly solvable MDPs)

Demand Dispatch

Demand Dispatch

Frequency Decomposition

Demand Dispatch: Power consumption from loads varies automatically and continuously to provide service to the grid, without impacting QoS to the consumer Frequency decomposition: Each class

• f flexible loads allocated to its own

bandwidth of service, based on QoS constraints and costs

Power Grid Control

Water Pump Batteries Coal Gas Turbine

BP BP BP C BP BP Voltage Frequency Phase

H C

Σ − Actuator feedback loop

A

LOAD

Today: PJM regulation signal: R = RegA + RegD

One Day at CAISO 2020

Net Load Curve Low pass Mid pass High pass

The duck is a sum of a smooth energy signal, and two zero-energy services GW

• 5

5 10 15 20 25 12am 12am 3am 6am 9am 12pm 3pm 6pm 9pm

Demand Dispatch

Demand Dispatch

Frequency Decomposition

Demand Dispatch: Power consumption from loads varies automatically and continuously to provide service to the grid, without impacting QoS to the consumer Frequency decomposition: Each class

• f flexible loads allocated to its own

bandwidth of service, based on QoS constraints and costs

Power Grid Control

Water Pump Batteries Coal Gas Turbine

BP BP BP C BP BP Voltage Frequency Phase

H C

Σ − Actuator feedback loop

A

LOAD

Today: PJM regulation signal: R = RegA + RegD

Jan 01 Jan 02 Jan 03 Jan 04 Jan 05 Jan 06 GW 1 2 3 4

G Goal:

W (t) = Wind generation in BPA, Jan 2015 Ramps

GW (t) + Gr(t) ≡ 4GW

Demand Dispatch

Demand Dispatch

Frequency Decomposition

Demand Dispatch: Power consumption from loads varies automatically and continuously to provide service to the grid, without impacting QoS to the consumer Frequency decomposition: Each class

• f flexible loads allocated to its own

bandwidth of service, based on QoS constraints and costs

Power Grid Control

Water Pump Batteries Coal Gas Turbine

BP BP BP C BP BP Voltage Frequency Phase

H C

Σ − Actuator feedback loop

A

LOAD

Today: PJM regulation signal: R = RegA + RegD

Jan 01 Jan 02 Jan 03 Jan 04 Jan 05 Jan 06 GW 1 2 3 4

Goal: GW (t) + Gr(t) ≡ 4GW

• btained from

generation? Gr(t) Gr(t)

Scary

Demand Dispatch

Demand Dispatch

Frequency Decomposition

Demand Dispatch: Power consumption from loads varies automatically and continuously to provide service to the grid, without impacting QoS to the consumer Frequency decomposition: Each class

• f flexible loads allocated to its own

bandwidth of service, based on QoS constraints and costs

Power Grid Control

Water Pump Batteries Coal Gas Turbine

BP BP BP C BP BP Voltage Frequency Phase

H C

Σ − Actuator feedback loop

A

LOAD

Today: PJM regulation signal: R = RegA + RegD

Jan 01 Jan 02 Jan 03 Jan 04 Jan 05 Jan 06 GW 1 2 3 4

Gr(t) G1 G2 G

Traditional generation DD: Chillers & Pool Pumps DD: HVAC Fans

3

Gr = G1 + G2 + G3

Demand Dispatch

Demand Dispatch

Frequency Decomposition

Demand Dispatch: Power consumption from loads varies automatically and continuously to provide service to the grid, without impacting QoS to the consumer Frequency decomposition: Each class

• f flexible loads allocated to its own

bandwidth of service, based on QoS constraints and costs

Power Grid Control

Water Pump Batteries Coal Gas Turbine

BP BP BP C BP BP Voltage Frequency Phase

H C

Σ − Actuator feedback loop

A

LOAD

Today: PJM regulation signal: R = RegA + RegD

5 10 15

Low pass Mid pass High pass Higher pass

One Day at CAISO: 4GW loss at 6am

Net Load Curve

12am 12am 3am 6am 9am 12pm 3pm 6pm 9pm GW

Zero energy reg signals

Demand Dispatch

Demand Dispatch

Responsive Regulation and desired QoS – A partial list of the needs of the grid operator, and the consumer

High quality Ancillary Service? Customer QoS constraints satisfied?

Demand Dispatch

Demand Dispatch

Responsive Regulation and desired QoS – A partial list of the needs of the grid operator, and the consumer

High quality Ancillary Service? Customer QoS constraints satisfied? Cost effective?

Includes installation cost, communication cost, maintenance, and environmental.

Demand Dispatch

Demand Dispatch

Responsive Regulation and desired QoS – A partial list of the needs of the grid operator, and the consumer

High quality Ancillary Service? Customer QoS constraints satisfied? Cost effective?

Includes installation cost, communication cost, maintenance, and environmental.

Reliable?

Will AS be available each day? (may vary with time, but capacity must be predictable)

Demand Dispatch

Demand Dispatch

Responsive Regulation and desired QoS – A partial list of the needs of the grid operator, and the consumer

High quality Ancillary Service? Customer QoS constraints satisfied? Cost effective?

Includes installation cost, communication cost, maintenance, and environmental.

Reliable?

Will AS be available each day? (may vary with time, but capacity must be predictable)

Is the incentive to the consumer reliable?

If a consumer receives a $50 payment for one month, and only $1 the next, will there be an explanation that is clear to the consumer?

Demand Dispatch

Demand Dispatch

Gr Gr = G1 + G2 + G3 G1 G2 G3 ?

Demand Dispatch

Demand Dispatch

Gr Gr = G1 + G2 + G3 G1 G2 G Traditional generation

3

Demand Dispatch

Demand Dispatch

Gr Gr = G1 + G2 + G3 G1 G2 G Traditional generation Water pumping (e.g. pool pumps) Fans in commercial HVAC

3

Demand Dispatch: Power consumption from loads varies automatically and continuously to provide service to the grid, without impacting QoS to the consumer

Demand Dispatch

Demand Dispatch for Virtual Energy Storage

Responsive Regulation and desired QoS

Demand Dispatch: Power consumption from loads varies automatically and continuously to provide service to the grid, without impacting QoS to the consumer High quality Ancillary Service Reliable Cost effective Customer QoS constraints satisfied Virtual energy storage: achieve these goals simultaneously through distributed control

Control Architecture

Control Goals and Architecture

Macro control

High-level control layer: BA or a load aggregator. The balancing challenges are of many different categories and time-scales: Automatic Generation Control (AGC); time scales of seconds to 20 minutes. Balancing reserves. In the Bonneville Power Authority, the balancing reserves include both AGC and balancing on timescales of many hours. Balancing on a slower time-scale is achieved through real time markets in some other regions of the U.S. Contingencies (e.g., a generator outage) Peak shaving Smoothing ramps from solar or wind generation

Control Architecture

Control Goals and Architecture

Local Control: decision rules designed to respect needs of load and grid

Local feedback loop Local Control Load i

ζt Y i

t

U i

t

Xi

t

Grid signal Local decision Power deviation

• Min. communication: each load monitors its state and a regulation signal

from the grid. Aggregate must be controllable: randomized policies for finite-state loads.

Control Architecture

Randomized Policies

Local feedback loop Local Control Load i

ζt Y i

t

U i

t

Xi

t

Grid signal Local decision Power deviation

Local control architecture For the ith load: Y i

t power, U i t load setpoint, Xi t local state.

Signal ζt is from the grid operator – common to all loads of a certain class.

Control Architecture

Randomized Policies

Local feedback loop Local Control Load i

ζt Y i

t

U i

t

Xi

t

Grid signal Local decision Power deviation

Local control architecture For the ith load: Y i

t power, U i t load setpoint, Xi t local state.

Signal ζt is from the grid operator – common to all loads of a certain class. Policy: Decision rule that maps (ζt, Xi

t) to the input U i t.

Randomized Policy: Decision rule also depends on rand (an intelligent coin-flip)

Control Architecture

Load Model

Controlled Markovian Dynamics

...

Load 1

BA

Reference (MW)

Load 2 Load N

ζ r

+

Gc

Power Consumption (MW)

Discrete time: ith load Xi(t) evolves on finite state space X Each load is subject to common controlled Markovian dynamics. Signal ζ = {ζt} is broadcast to all loads Controlled transition matrix {Pζ : ζ ∈ R}: P{Xi

t+1 = x′ | Xi t = x, ζt = ζ} = Pζ(x, x′)

Questions

• How to analyze aggregate of similar loads?
• Local control design?

Control Architecture

Randomized Policies

• How to analyze aggregate of similar loads?
• How to design Pζ?

−600 −400 −200 200 400 600

Day 1 Day 2 Day 3 Day 4 Day 5 Day 6 Day 7

MW

Mean-field model

How to analyze aggregate?

Mean field model

N loads running independently, each under the command ζ. Empirical Distributions: µN

t (x) = 1

N

N

• i=1

I{Xi(t) = x}, x ∈ X

U(x) power consumption in state x,

yN

t = 1

N

N

• i=1

U(Xi

t) =

• x

µN

t (x)U(x)

Mean-field model:

via Law of Large Numbers for martingales

µt+1 = µtPζt, yt = µt, U ζt = ft(y0, . . . , yt) by design

Example: pool pumps

Example: pool pumps

How Pools Can Help Regulate The Grid

1,5KW 400V

Needs of a single pool ⊲ Filtration system circulates and cleans: Average pool pump uses 1.3kW and runs 6-12 hours per day, 7 days per week ⊲ Pool owners are oblivious, until they see frogs and algae ⊲ Pool owners do not trust anyone: Privacy is a big concern

Example: pool pumps

One Million Pools in Florida

Pools Service the Grid Today

On Call1: Utility controls residential pool pumps and other loads

1Florida Power and Light, Florida’s largest utility.

www.fpl.com/residential/energysaving/programs/oncall.shtml

Example: pool pumps

One Million Pools in Florida

Pools Service the Grid Today

On Call1: Utility controls residential pool pumps and other loads Contract for services: no price signals involved Used only in times of emergency — Activated only 3-4 times a year

1Florida Power and Light, Florida’s largest utility.

www.fpl.com/residential/energysaving/programs/oncall.shtml

Example: pool pumps

One Million Pools in Florida

Pools Service the Grid Today

On Call1: Utility controls residential pool pumps and other loads Contract for services: no price signals involved Used only in times of emergency — Activated only 3-4 times a year Opportunity: FP&L already has their hand on the switch of nearly one million pools!

1Florida Power and Light, Florida’s largest utility.

www.fpl.com/residential/energysaving/programs/oncall.shtml

Example: pool pumps

One Million Pools in Florida

Pools Service the Grid Today

On Call1: Utility controls residential pool pumps and other loads Contract for services: no price signals involved Used only in times of emergency — Activated only 3-4 times a year Opportunity: FP&L already has their hand on the switch of nearly one million pools! Surely pools can provide much more service to the grid

1Florida Power and Light, Florida’s largest utility.

www.fpl.com/residential/energysaving/programs/oncall.shtml

Example: pool pumps

An Intelligent Pool

Local Control Architecture

1 2

. . .

On Off 1 2

. . .

I −1 I I I −1

Local control architecture For the ith load: Y i

t = power: 1kW when running,

Example: pool pumps

An Intelligent Pool

Local Control Architecture

1 2

. . .

On Off 1 2

. . .

I −1 I I I −1

Local control architecture For the ith load: Y i

t = power: 1kW when running,

U i

t = 1 or 0 (pool pump is running or not)

Example: pool pumps

An Intelligent Pool

Local Control Architecture

1 2

. . .

On Off 1 2

. . .

I −1 I I I −1

Local control architecture For the ith load: Y i

t = power: 1kW when running,

U i

t = 1 or 0 (pool pump is running or not)

Xi

t local state: (U i t, Ii t), with Ii t the time in current power state.

Example: pool pumps

An Intelligent Pool

Local Control Architecture

1 2

. . .

On Off 1 2

. . .

I −1 I I I −1

Local control architecture For the ith load: Y i

t = power: 1kW when running,

U i

t = 1 or 0 (pool pump is running or not)

Xi

t local state: (U i t, Ii t), with Ii t the time in current power state.

Randomized Policy: Decision rule that maps (ζt, Xi

t, randt i) to the input U i t.

Example: pool pumps

An Intelligent Pool

Local Control Architecture

1 2

. . .

On Off 1 2

. . .

I −1 I I I −1

Local control architecture For the ith load: Y i

t = power: 1kW when running,

U i

t = 1 or 0 (pool pump is running or not)

Xi

t local state: (U i t, Ii t), with Ii t the time in current power state.

Randomized Policy: Decision rule that maps (ζt, Xi

t, randt i) to the input U i t.

Randomized Policy: As ζ increases, probability of turning on increases:

24 hours 12 0.5 1

Pool with nominal 12 hour cleaning cycle

T

ζ =

• 4

ζ =

• 2

ζ = 4 ζ = 2 ζ =

(ζ)

Example: pool pumps

Pools in Florida Supply G2 – BPA regulation signal∗

Stochastic simulation using N = 106 pools

Reference Output deviation (MW)

−300 −200 −100 100 200 300 20 40 60 80 100 120 140 160 t/hour 20 40 60 80 100 120 140 160

PI control: ζt = 19et + 1.4eI

t ,

et = rt − yt and eI

t = t k=0 ek

Each pool pump turns on/off with probability depending on 1) its internal state, and 2) the BPA reg signal

∗transmission.bpa.gov/Business/Operations/Wind/reserves.aspx

Example: pool pumps

Mean Field Model

Linearized Dynamics

Mean-field model: µt+1 = µtPζt, yt = µt, U ζt = ft(y0, . . . , yt) Linear state space model: Φt+1 = AΦt + Bζt γt = CΦt Interpretations: |ζt| is small, and π denotes invariant measure for P0.

• Φt ∈ R|X|,

a column vector with Φt(x) ≈ µt(x) − π(x), x ∈ X

• γt ≈ yt − y0; deviation from nominal steady-state
• A = P T

0 , C = U T, and input dynamics linearized:

B

T = d

dζ πPζ

• ζ=0

Example: pool pumps

Transfer function

The linearization at a particular value ζ is the state space model with transfer function, Gζ(z) = C[Iz − A]−1B (1) in which A = P T

ζ , Ci =

Uζ(xi) for each i, and

Bi =

• x

πζ(x)Eζ(x, xi), 1 ≤ i ≤ d (2) Eζ =

d dζ Pζ

Uζ = U − ¯

Uζ, with ¯ Uζ = πζ(U).

Example: pool pumps

Tracking Grid Signal with Residential Loads

Example: 20 pools, 20 kW max load

Each pool consumes 1kW when operating 12 hour cleaning cycle each 24 hours Power Deviation:

20 40 60 80 100 120 140 160 10

• 10

t/hour

kW 20 pools

Input ζt

• 3

3

Output deviation Reference

Nearly Perfect Service from Pools

Meyn et al. 2013 [CDC], Meyn et al. 2015 [IEEE TAC]

Example: pool pumps

Tracking Grid Signal with Residential Loads

Example: 300,000 pools, 300 MW max load

Each pool consumes 1kW when operating 12 hour cleaning cycle each 24 hours Power Deviation:

20 40 60 80 100 120 140 160 t/hour

Output deviation Reference

Input ζt

• 3

3 −100 −50 50 100

MW 300,000 pools

Nearly Perfect Service from Pools

Meyn et al. 2013 [CDC], Meyn et al. 2015 [IEEE TAC]

Example: pool pumps

Range of services provided by pools

Example: 10,000 pools, 10 MW max load Reference Power Deviation

20 40 60 80 100 120 140 160

MW

• 4
• 2

2

• 15

15

t/hour

12 hr/day cycle

ζ

Local Control Design

Local Control Design

Local Design

Goal: Construct a family of transition matrices {Pζ : ζ ∈ R}

Myopic Design: P myop

ζ

(x, x′) := P0(x, x′) exp

• ζU(x′) − Λζ(x)
• with Λζ(x) := log
• x′ P0(x, x′) exp
• ζU(x′)
• the normalizing constant.

Local Control Design

Local Design

Goal: Construct a family of transition matrices {Pζ : ζ ∈ R}

Myopic Design: P myop

ζ

(x, x′) := P0(x, x′) exp

• ζU(x′) − Λζ(x)
• with Λζ(x) := log
• x′ P0(x, x′) exp
• ζU(x′)
• the normalizing constant.

Exponential family design: Pζ(x, x′) := P0(x, x′) exp

• hζ(x, x′) − Λhζ(x)
• with

hζ(x, x′) = ζH0(x, x′). The choice of H0 will typically correspond to the linearization of a more advanced design around the value ζ = 0 (or some other fixed value of ζ).

Local Control Design

Local Design

Goal: Construct a family of transition matrices {Pζ : ζ ∈ R}

Individual Perspective Design Consider a finite-time-horizon optimization problem: For a given terminal time T, let p0 denote the pmf on strings of length T, p0(x1, . . . , xT ) =

T −1

• i=0

P0(xi, xi+1) , where x0 ∈ X is assumed to be given. The scalar ζ ∈ R is interpreted as a weighting parameter in the following definition of total welfare. For any pmf p, WT (p) = ζEp T

• t=1

U(Xt)

• − D(pp0)

where the expectation is with respect to p, and D denotes relative entropy: D(pp0) :=

• x1,...,xT

log p(x1, . . . , xT ) p0(x1, . . . , xT )

• p(x1, . . . , xT )

Local Control Design

Local Design

Goal: Construct a family of transition matrices {Pζ : ζ ∈ R}

Myopic design is an optimizer for the horizon T = 1, P myop

ζ

(x0, ·) ∈ arg max

p

W1(p).

Local Control Design

Local Design

Goal: Construct a family of transition matrices {Pζ : ζ ∈ R}

Myopic design is an optimizer for the horizon T = 1, P myop

ζ

(x0, ·) ∈ arg max

p

W1(p). The infinite-horizon mean welfare is denoted, η∗

ζ = lim T →∞

1 T WT (p∗

T )

Local Control Design

Local Design

Goal: Construct a family of transition matrices {Pζ : ζ ∈ R}

Myopic design is an optimizer for the horizon T = 1, P myop

ζ

(x0, ·) ∈ arg max

p

W1(p). The infinite-horizon mean welfare is denoted, η∗

ζ = lim T →∞

1 T WT (p∗

T )

Explicit construction via eigenvector problem: Pζ(x, y) = 1 λ v(y) v(x) ˆ Pζ(x, y) , x, y ∈ X, where ˆ Pζv = λv, ˆ Pζ(x, y) = exp(ζU(x))P0(x, y) for λ Perron-Frobenious eigenvalue.

Extension/reinterpretation of [Todorov 2007] + [Kontoyiannis & Meyn 200X]

Local Control Design

Local Design

Extending local control design to include exogenous disturbances

State space for a load model: X = Xu × Xn. Components Xn are not subject to direct control (e.g. impact of the weather on the climate of a building).

Local Control Design

Local Design

Extending local control design to include exogenous disturbances

State space for a load model: X = Xu × Xn. Components Xn are not subject to direct control (e.g. impact of the weather on the climate of a building). Conditional-independence structure of the local transition matrix P(x, x′) = R(x, x′

u)Q0(x, x′ n),

x′ = (x′

u, x′ n)

Q0 models uncontroled load dynamics and exogenous disturbances.

Local Control Design

Local Design

Goal: Construct a family of transition matrices {Pζ : ζ ∈ R}

Nominal model A Markovian model for an individual load, based on its typical behavior. Finite state space X = {x1, . . . , xd}; Transition matrix P0, with unique invariant pmf π0.

Local Control Design

Local Design

Goal: Construct a family of transition matrices {Pζ : ζ ∈ R}

Nominal model A Markovian model for an individual load, based on its typical behavior. Finite state space X = {x1, . . . , xd}; Transition matrix P0, with unique invariant pmf π0. Common structure for design The family of transition matrices used for distributed control is of the form: Pζ(x, x′) := P0(x, x′) exp

• hζ(x, x′) − Λhζ(x)
• with hζ continuously differentiable in ζ, and the normalizing constant

Λhζ(x) := log

• x′

P0(x, x′) exp

• hζ(x, x′)

Local Control Design

Local Design

Goal: Construct a family of transition matrices {Pζ : ζ ∈ R}

Nominal model A Markovian model for an individual load, based on its typical behavior. Finite state space X = {x1, . . . , xd}; Transition matrix P0, with unique invariant pmf π0. Common structure for design The family of transition matrices used for distributed control is of the form: Pζ(x, x′) := P0(x, x′) exp

• hζ(x, x′) − Λhζ(x)
• with hζ continuously differentiable in ζ, and the normalizing constant

Λhζ(x) := log

• x′

P0(x, x′) exp

• hζ(x, x′)
• Assumption: for all x ∈ X, x′ = (x′

u, x′ n) ∈ X, hζ(x, x′) = hζ(x, x′ u).

Local Control Design

Local Design

Goal: Construct a family of transition matrices {Pζ : ζ ∈ R}

Construction of the family of functions {hζ : ζ ∈ R}

Local Control Design

Local Design

Goal: Construct a family of transition matrices {Pζ : ζ ∈ R}

Construction of the family of functions {hζ : ζ ∈ R} Step 1: The specification of a function H that takes as input a transition matrix. H = H(P) is a real-valued function on X × X.

Local Control Design

Local Design

Goal: Construct a family of transition matrices {Pζ : ζ ∈ R}

Construction of the family of functions {hζ : ζ ∈ R} Step 1: The specification of a function H that takes as input a transition matrix. H = H(P) is a real-valued function on X × X. Step 2: The families {Pζ} and {hζ} are defined by the solution to the ODE:

d dζ hζ = H(Pζ),

ζ ∈ R, in which Pζ is determined by hζ through: Pζ(x, x′) := P0(x, x′) exp

• hζ(x, x′) − Λhζ(x)
• The boundary condition: h0 ≡ 0.

Local Control Design

Local Design

Extending local control design to include exogenous disturbances

For any function H◦ : X → R, one can define H(x, x′

u) =

• x′

n

Q0(x, x′

n)H◦(x′ u, x′ n)

(3)

Local Control Design

Local Design

Extending local control design to include exogenous disturbances

For any function H◦ : X → R, one can define H(x, x′

u) =

• x′

n

Q0(x, x′

n)H◦(x′ u, x′ n)

(3) Then functions {hζ} satisfy hζ(x, x′

u) =

• x′

n

Q0(x, x′

n)h◦ ζ(x′ u, x′ n),

for some h◦

ζ : X → R. Moreover, these functions solve the d-dimensional ODE, d dζ h◦ ζ = H◦(Pζ),

ζ ∈ R, with boundary condition h◦

0 ≡ 0.

Local Control Design

Individual Perspective Design

Local welfare function: Wζ(x, P) = ζU(x) − D(PP0),

where D denotes relative entropy: D(PP0) =

x′ P(x, x′) log

P (x,x′)

P0(x,x′)

.

Markov Decision Process: lim supT →∞

1 T

T

t=1 E[Wζ(Xt, P)]

Average reward optimization equation (AROE): max

P

• Wζ(x, P) +
• x′

P(x, x′)h∗

ζ(x′)

• = h∗

ζ(x) + η∗ ζ

where P(x, x′) = R(x, x′

u)Q0(x, x′ n),

x′ = (x′

u, x′ n)

Local Control Design

Individual Perspective Design

ODE method for IPD design:

Local Control Design

Individual Perspective Design

ODE method for IPD design: Family {Pζ}: Pζ(x, x′) := P0(x, x′) exp

• hζ(x, x′) − Λhζ(x)

Local Control Design

Individual Perspective Design

ODE method for IPD design: Family {Pζ}: Pζ(x, x′) := P0(x, x′) exp

• hζ(x, x′) − Λhζ(x)
• Functions {hζ}: hζ(x, x′

u) = x′

n Q0(x, x′

n)h◦ ζ(x′ u, x′ n),

for h◦

ζ : X → R solutions of the d-dimensional ODE, d dζ h◦ ζ = H◦(Pζ),

ζ ∈ R, with boundary condition h◦

0 ≡ 0.

Local Control Design

Individual Perspective Design

ODE method for IPD design: Family {Pζ}: Pζ(x, x′) := P0(x, x′) exp

• hζ(x, x′) − Λhζ(x)
• Functions {hζ}: hζ(x, x′

u) = x′

n Q0(x, x′

n)h◦ ζ(x′ u, x′ n),

for h◦

ζ : X → R solutions of the d-dimensional ODE, d dζ h◦ ζ = H◦(Pζ),

ζ ∈ R, with boundary condition h◦

0 ≡ 0.

H◦

ζ (x) = d dζ h◦ ζ(x) = x′[Zζ(x, x′) − Zζ(x◦, x′)]U(x′),

x ∈ X,

where Zζ = [I − Pζ + 1 ⊗ πζ]−1 = ∞

n=0[Pζ − 1 ⊗ πζ]n is the fundamental matrix.

Local Control Design

Individual Perspective Design

Linearized dynamics

Unit on Unit off Θmin Θ1 Θm-1 Θ2 Θmax

... ... Magnitude (dB)

• 50
• 45
• 40
• 35
• 30
• 25
• 20

10-4 10-3 10-2 10-1 Phase (deg)

• 180
• 135
• 90
• 45

45 90 Frequency (rad/s) IPD

• 3

2.9 10

ζ = -6 ζ = -3 ζ = 0 ζ = 3 ζ = 6

Bode plots for IPD: Linearizations at five values of ζ

Local Control Design

Individual Perspective Design

Linearized dynamics

Unit on Unit off Θmin Θ1 Θm-1 Θ2 Θmax

... ... Magnitude (dB)

• 50
• 45
• 40
• 35
• 30
• 25
• 20

10-4 10-3 10-2 10-1 Phase (deg)

• 180
• 135
• 90
• 45

45 90 Frequency (rad/s) IPD

• 3

2.9 10

ζ = -6 ζ = -3 ζ = 0 ζ = 3 ζ = 6

Bode plots for IPD: Linearizations at five values of ζ Proof of positive real condition for reversible load dynamics.

Busic & Meyn [CDC’14] Passive Dynamics in Mean Field Control

Local Control Design

System Perspective Design

Strictly positive real by design

Goal: The transfer function of the delay-free linearized aggregate model is passive:

∞

• t=0

utyt+1 ≥ 0, ∀{ut}.

Local Control Design

System Perspective Design

Strictly positive real by design

Goal: The transfer function of the delay-free linearized aggregate model is passive:

∞

• t=0

utyt+1 ≥ 0, ∀{ut}.

Recall: The linearization at a particular value ζ is the state space model with transfer function, Gζ(z) = C[Iz − A]−1B in which A = P T

ζ , Ci =

Uζ(xi) for each i, and

Bi =

• x

πζ(x)Eζ(x, xi), 1 ≤ i ≤ d Eζ =

d dζ Pζ

Uζ = U − ¯

Uζ, with ¯ Uζ = πζ(U).

Local Control Design

System Perspective Design

Strictly positive real by design

Goal: The transfer function of the delay-free linearized aggregate model is passive:

∞

• t=0

utyt+1 ≥ 0, ∀{ut}.

Recall: The linearization at a particular value ζ is the state space model with transfer function, Gζ(z) = C[Iz − A]−1B in which A = P T

ζ , Ci =

Uζ(xi) for each i, and

Bi =

• x

πζ(x)Eζ(x, xi), 1 ≤ i ≤ d Eζ =

d dζ Pζ

Uζ = U − ¯

Uζ, with ¯ Uζ = πζ(U).

Sufficient condition: positive real. A discrete-time transfer function F is positive-real if it is stable (all poles are strictly within the unit disk), and F(ejθ) + F(e−jθ) ≥ 0, θ ∈ R.

Local Control Design

System Perspective Design

Strictly positive real by design

SPD design: P ▽ = P P, with P adjoint of P in L2(π):

P (x, x′) = π(x′)

π(x) P(x′, x), x, x′ ∈ X.

Local Control Design

System Perspective Design

Strictly positive real by design

SPD design: P ▽ = P P, with P adjoint of P in L2(π):

P (x, x′) = π(x′)

π(x) P(x′, x), x, x′ ∈ X.

H◦(x) =

x′[Z▽(x, x′) − Z▽(x◦, x′)]U(x′) x ∈ X

where Z▽ = [I − P ▽ + 1 ⊗ π]−1 the fundamental matrix for P ▽

Local Control Design

System Perspective Design

Strictly positive real by design

SPD design: P ▽ = P P, with P adjoint of P in L2(π):

P (x, x′) = π(x′)

π(x) P(x′, x), x, x′ ∈ X.

H◦(x) =

x′[Z▽(x, x′) − Z▽(x◦, x′)]U(x′) x ∈ X

where Z▽ = [I − P ▽ + 1 ⊗ π]−1 the fundamental matrix for P ▽

• Thm. (SPD design) If P ▽

0 = P 0 P0 is irreducible, and P0 = R0, then the linearized

state-space model at any constant value ζ satisfies G+

ζ (ejθ) + G+ ζ (e−jθ) ≥ σ2 ζ,

θ ∈ R where σ2

ζ is the variance of U under πζ and G+(z) := zG(z).

Local Control Design

Exponential family

Alternative to solving an ODE

For a function H◦

e : X → R, define for each x, x′ u and ζ,

hζ(x, x′

u) = ζHe(x′ u | x)

with He(x′

u | x) :=

• x′

n

Q0(x, x′

n)H◦ e (x′ u, x′ n)

Local Control Design

Exponential family

Alternative to solving an ODE

For a function H◦

e : X → R, define for each x, x′ u and ζ,

hζ(x, x′

u) = ζHe(x′ u | x)

with He(x′

u | x) :=

• x′

n

Q0(x, x′

n)H◦ e (x′ u, x′ n)

Myopic design: H◦

e = U.

Local Control Design

Exponential family

Alternative to solving an ODE

For a function H◦

e : X → R, define for each x, x′ u and ζ,

hζ(x, x′

u) = ζHe(x′ u | x)

with He(x′

u | x) :=

• x′

n

Q0(x, x′

n)H◦ e (x′ u, x′ n)

Myopic design: H◦

e = U.

Linear approximations to the IPD or SPD solutions, with H◦

e = H◦(P0).

Local Control Design

Myopic Design

Linearized dynamics

Unit on Unit off Θmin Θ1 Θm-1 Θ2 Θmax

... ... Magnitude (dB)

• 50
• 45
• 40
• 35
• 30
• 25
• 20
• 15
• 10

10-4 10-3 10-2 10-1 Phase (deg)

• 180
• 135
• 90
• 45

45 90 Myopic

ζ = -6 ζ = -3 ζ = 0 ζ = 3 ζ = 6

• 3

2.9 10 Frequency (rad/s)

Bode plots for myopic design: Linearizations at five values of ζ

Local Control Design

Myopic Design

Linearized dynamics

Magnitude (dB)

• 50
• 45
• 40
• 35
• 30
• 25
• 20
• 15
• 10

10-4 10-3 10-2 10-1 Phase (deg)

• 180
• 135
• 90
• 45

45 90 Myopic

ζ = -6 ζ = -3 ζ = 0 ζ = 3 ζ = 6

• 3

2.9 10 Frequency (rad/s) Magnitude (dB)

• 50
• 45
• 40
• 35
• 30
• 25
• 20

10-4 10-3 10-2 10-1 Phase (deg)

• 180
• 135
• 90
• 45

45 90 Frequency (rad/s) IPD

• 3

2.9 10

ζ = -6 ζ = -3 ζ = 0 ζ = 3 ζ = 6

Local Control Design

Example: Thermostatically Controlled Loads

refrigerators, water heaters, air-conditioning . . . TCLs are already equipped with primitive “local intelligence” based on a deadband (or hysteresis interval) The state process for a TCL at time t: X(t) = (Xu(t), Xn(t)) = (m(t), Θ(t)) , where m(t) ∈ {0, 1} denotes the power mode (“1” indicating the unit is on), and Θ(t) the inside temperature of the load Exogenous disturbances: ambient temperature, and usage

Local Control Design

Example: Thermostatically Controlled Loads

The standard ODE model of a water heater is the first-order linear system, d dtΘ(t) = −λ[Θ(t) − Θa(t)] + γm(t) − α[Θ(t) − Θin(t)]f(t) , Θ(t) temperature of the water in the tank Θin(t) temperature of the cold water entering the tank f(t) flow rate of hot water from the WH m(t) power mode of the WH (“on” indicated by m(t) = 1). Deterministic deadband control: Θ(t) ∈ [Θ−, Θ+] Nominal model for local control design: based on the specification of two CDFs for the temperature at which the load turns on or turns off

F (θ) Θ− Θ+ θ⊕ 1 F ⊕(θ θ ) Θ− Θ+ θ 1 ̺

Local Control Design

Example: Thermostatically Controlled Loads

Discrete-time control. At time instance k, if the water heater is on (i.e., m(k) = 1), then it turns off with probability, p⊖(k + 1) = [F ⊖(Θ(k + 1)) − F ⊖(Θ(k))]+ 1 − F ⊖(Θ(k)) where [x]+ := max(0, x) for x ∈ R; Similarly, if the load is off, then it turns on with probability p⊕(k + 1) = [F ⊕(Θ(k)) − F ⊕(Θ(k + 1))]+ F ⊕(Θ(k)) The nominal behavior of the power mode can be expressed P{m(k) = 1 | θ(k − 1), θ(k), m(k − 1) = 0} = p⊕(k) P{m(k) = 0 | θ(k − 1), θ(k), m(k − 1) = 1} = p⊖(k)

Local Control Design

Example: Thermostatically Controlled Loads

Myopic design - exponential tilting of these distributions: p⊕

ζ (k) := P{m(k) = 1 | θ(k − 1), θ(k), m(k − 1) = 0, ζ(k − 1) = ζ}

= p⊕(k)eζ p⊕(k)eζ + 1 − p⊕(k) p⊖

ζ (k) = P{m(k) = 0 | θ(k − 1), θ(k), m(k − 1) = 1, ζ(k − 1) = ζ}

= p⊖(k) p⊖(k) + (1 − p⊖(k))eζ If p⊕

0 (k) > 0, then the probability p⊕ ζ (k) is strictly increasing in ζ, approaching 1

as ζ → ∞; it approaches 0 as ζ → −∞, if p⊕

0 (k) < 1.

Local Control Design

Example: Thermostatically Controlled Loads

System identification

d dtΘ(t) = −λ[Θ(t) − Θa(t)] + γm(t) − α[Θ(t) − Θin(t)]f(t) ,

Θ(t) temperature of the water in the tank Θin(t) temperature of the cold water entering the tank f(t) flow rate of hot water from the WH m(t) power mode of the WH (“on” indicated by m(t) = 1).

• Temp. Ranges

ODE Pars.

• Loc. Control

Θ+ ∈ [118, 122] F λ ∈ [8, 12.5] × 10−6 Ts = 15 sec Θ− ∈ [108, 112] F γ ∈ [2.6, 2.8] × 10−2 κ = 4 Θa ∈ [68, 72] F α ∈ [6.5, 6.7] × 10−2 ̺ = 0.8 Θin ∈ [68, 72] F Pon = 4.5 kW θ0 = Θ−

Heterogeneous population: 100 000 WHs simulated by uniform sampling of the values in the table Usage data from Oakridge National Laboratory (35WHs over 50 days)

Local Control Design

Tracking performance

and the controlled dynamics for an individual load

100,000 water-heaters When on, individual load consumes 4, 5 kW With no usage, approx. 2% duty cycle, avg. power consumption 10MW.

80 100 120 140 5 10 15 20 5 10 15 20 80 100 120 140 80 100 120 140 50 100

• 50

50

MW MW MW

• 10

10

Nominal power consumption Tracking Tracking Typical Load Response

temp (F) temp (F) temp (F)

rt ≡ 0

No reg:

|rt| ≤ 40 MW |rt| ≤ 10 MW

Load On Load On Load On (hrs)

t

(hrs)

t

BPA Reference: Power Deviation

rt

Local Control Design

Tracking performance

Potential for contingency reserves and ramping

ζ

• 8
• 6
• 4
• 2

2 4 6 8

Power deviation (MW)

• 6
• 5
• 4
• 3
• 2
• 1

1 2

• 8
• 6
• 4
• 2

2 4 6 8

• 6
• 4
• 2

2

Tracking two sawtooth waves with 100,000 water heaters: average power consumption 8MW

5 10 15 20 5 10 15 20

Reference Power Deviation ζ ζ (hrs)

t

Local Control Design

Tracking performance

and the controlled dynamics for an individual load

Heterogeneous setting: 40 000 loads per experiment; 20 different load types in each case Lower plots show the on/off state for a typical load

Stochastic Output Mean-field Model BPA balancing reserves (filtered/scaled)

Open Loop Tracking (MW) Power state

• 15
• 10
• 5

5 10 1 15

• 6
• 4
• 2

2 4 6

• 4
• 2

2 4 6

Refrigerators Fast Electric Water Heaters Slow Electric Water Heaters

24 hrs 24 hrs 6 hrs

Nominal Demand Dispatch

Local Control Design

Example: fleet of batteries

[B. Hashmi, Meyn ACC’17] State: x = (m, s), where m ∈ {ch, dis, id} denotes charging mode, and s ∈ [0, 1] denotes the SoC. The power delivery at state x depends only on charging mode:

U(ch, s) = Uch < 0, U(id, s) = 0, U(dis, s) = Udis > 0.

Local Control Design

Example: fleet of batteries

[B. Hashmi, Meyn ACC’17] State: x = (m, s), where m ∈ {ch, dis, id} denotes charging mode, and s ∈ [0, 1] denotes the SoC. The power delivery at state x depends only on charging mode:

U(ch, s) = Uch < 0, U(id, s) = 0, U(dis, s) = Udis > 0.

Nominal model: Xi

t = (M i t, Si t) denote the state of ith battery at time t.

Si

t+1 = Si t + hδch, if M i t = ch, Si t+1 = Si t − hδdis, if M i t = dis,

Si

t+1 = Si t, if M i t = id, where h is the time step length, and δch and δdis charging

and discharging rates.

The dynamics of the first component are governed by a “two coin-flip” randomized policy. For example, in state (ch, s), the battery changes its mode to idle with probability (1 − pch(s)) × pid(s)/(pid(s) + pdis(s))

0.5 1

pid idle pch charging pdis discharging

1

SoC

Local Control Design

Example: fleet of batteries

1000 batteries, tracking PJM RegD test signal:

Local Control Design

Unmodeled dynamics

[Chen, B., Meyn CDC’15, IEEE TAC’17] Setting: 0.1% sampling, and

1

Heterogeneous population of loads

2

Load i overrides when QoS is out of bounds

0.5 −10 −5 5 10

MW

100 120 110 130

• pt out %

N = 300,000 N = 30,000

100 120 110 130

Closed-loop tracking

−100 −50 50 100 0.5

Output deviation Reference

t/hour t/hour

PI control: ζt = kP et + kIeI

t ,

et = rt − yt, eI

t = t s=0 es

Conclusions and Future Directions

Control Architecture

Frequency Allocation for Demand Dispatch

10-2 10-1 100 101 Frequency (rad/s) 10-5 10-4 10-3 Frequency (rad/s) Magnitude (dB)

• 15
• 10
• 5

5 10 15 20 Phase (deg)

• 90
• 45

45 G r i d T r a n s f e r F u nc t i

• n

Uncertainty Here Fans in Commercial Buildings Residential Water Heaters Refrigerators Water Pumping Pool Pumps Chiller Tanks

Bandwidth centered around its natural cycle

Reference (from Bonneville Power Authority)

10,000 pools

Output deviation

−300 −200 −100 100 200 300

Tracking BPA Regulation Signal (MW)

20 40 60 80 100 120 140 160 t/hour 20 40 60 80 100 120 140 160

Conclusions and Future Directions

Conclusions

Virtual storage from flexible loads

Approach: creating Virtual Energy Storage through direct control of flexible loads

• helping the grid while respecting user QoS

Conclusions and Future Directions

Conclusions

Virtual storage from flexible loads

Approach: creating Virtual Energy Storage through direct control of flexible loads

• helping the grid while respecting user QoS

Challenges: − Stability properties for IPD and myopic design? − Information Architecture: ζt = f(?) Different needs for communication, state estimation and forecast. − Capacity estimation (time varying) − Network constraints − Resource optimization & learning Integrating VES with traditional generation and batteries. − Economic issues Contract design, aggregators, markets . . .

Conclusions and Future Directions

Conclusions

Thank You!

Conclusions and Future Directions

References: this talk

• A. Buˇ

si´ c and S. Meyn. Distributed randomized control for demand dispatch. 55th IEEE Conference on Decision and Control, 2016.

• A. Buˇ

si´ c and S. Meyn. Ordinary Differential Equation Methods For Markov Decision Processes and Application to Kullback-Leibler Control Cost. arXiv:1605.04591v2. Oct 2016.

• S. Meyn, P. Barooah, A. Buˇ

si´ c, Y. Chen, and J. Ehren. Ancillary Service to the Grid Using Intelligent Deferrable Loads. IEEE Trans. Automat. Contr., 60(11): 2847-2862, 2015.

• P. Barooah, A. Buˇ

si´ c, and S. Meyn. Spectral Decomposition of Demand-Side Flexibility for Reliable Ancillary Services in a Smart Grid. 48th Annual Hawaii International Conference on System Sciences (HICSS). 2015.

• A. Buˇ

si´ c and S. Meyn. Passive dynamics in mean field control. 53rd IEEE Conf. on Decision and Control (CDC) 2014.

Conclusions and Future Directions

References: related

Demand dispatch:

• Y. Chen, A. Buˇ

si´ c, and S. Meyn. Individual risk in mean-field control models for decentralized control, with application to automated demand response. 53rd IEEE Conf. on Decision and Control (CDC), 2014.

• Y. Chen, A. Buˇ

si´ c, and S. Meyn. State Estimation and Mean Field Control with Application to Demand

• Dispatch. 54rd IEEE Conference on Decision and Control (CDC) 2015.
• J. L. Mathieu. Modeling, Analysis, and Control of Demand Response Resources. PhD thesis, Berkeley,

2012.

• J. L. Mathieu, S. Koch, D. S. Callaway, State Estimation and Control of Electric Loads to Manage

Real-Time Energy Imbalance, IEEE Transactions on Power Systems, 28(1):430-440, 2013.

Markov processes:

• I. Kontoyiannis and S. P. Meyn. Spectral theory and limit theorems for geometrically ergodic Markov
• processes. Ann. Appl. Probab., 13:304–362, 2003.
• I. Kontoyiannis and S. P. Meyn. Large deviations asymptotics and the spectral theory of multiplicatively

regular Markov processes. Electron. J. Probab., 10(3):61–123 (electronic), 2005.

• E. Todorov. Linearly-solvable Markov decision problems. In B. Sch¨
• lkopf, J. Platt, and T. Hoffman,

editors, Advances in Neural Information Processing Systems, (19) 1369–1376. MIT Press, Cambridge, MA, 2007.

Conclusions and Future Directions

Mean Field Model

Linearized Dynamics

Mean-field model: µt+1 = µtPζt, yt = µt, U ζt = ft(y0, . . . , yt) Linear state space model: Φt+1 = AΦt + Bζt γt = CΦt Interpretations: |ζt| is small, and π denotes invariant measure for P0.

• Φt ∈ R|X|,

a column vector with Φt(x) ≈ µt(x) − π(x), x ∈ X

• γt ≈ yt − y0; deviation from nominal steady-state
• A = P T

0 , C = U T, and input dynamics linearized:

B

T = d

dζ πPζ

• ζ=0