# A Stability Analysis of Thermostatically Controlled Loads for Power System Frequency Control

Authors: Ellen Webborn; Robert S. MacKay
Journal: Complexity
 (2017)

---

## Abstract

Thermostatically controlled loads (TCLs) are a flexible demand resource with the potential to play a significant role in supporting electricity grid operation. We model a large number of identical TCLs acting autonomously according to a deterministic control scheme to provide frequency response as a population of coupled oscillators. We perform stability analysis to explore the danger of the TCL temperature cycles synchronising: an emergent phenomenon often found in populations of coupled oscillators and predicted in this type of demand response scheme. We take identical TCLs as it can be assumed to be the worst case. We find that the uniform equilibrium is stable and the fully synchronised periodic cycle is unstable, suggesting that synchronisation might not be as serious a danger as feared. Then detailed simulations are performed to study the effects of a population of frequency-sensitive TCLs acting under real system conditions using historic system data. The potential reduction in frequency response services required from other providers is determined, for both homogeneous and heterogeneous populations. For homogeneous populations, we find significant synchronisation, but very minimal diversity removes the synchronisation effects. In summary, we combine dynamical systems stability analysis with large-scale simulations to offer new insights into TCL switching behaviour.


---

## Body

## 1. Introduction

To operate the electricity grid reliably and securely requires controlling a number of factors, one of which is the electricity grid frequency. The AC frequency continually varies close to its nominal value (50 Hz in Europe) and is kept there by the System Operator (SO) in order to respect regulations and prevent network instabilities and blackouts. This is mainly done by employing flexible power generators, such as gas turbines, to vary their output in response to imbalance between supply and demand. This type of service is often known as frequency response, or frequency regulation. With the arrival of large numbers of wind and solar farms, smart meters, and thousands of domestic solar panels, uncontrolled and largely invisible to the System Operator, new challenges and opportunities have arisen. Rather than relying on a few large (typically fossil-fuelled) power plants to provide system balancing, there is the potential, and perhaps a need, for consumers to play a role. However, for a large number of very small players to participate would require a modelling approach capable of incorporating the complexities of such a system and foreseeing emergent phenomena that may arise as a result of the interactions involved. In this article we explore the potential for certain types of domestic appliances to provide frequency response through simple deterministic rules and consider the possibility of (potentially harmful) demand synchronisation.Thermostatically controlled loads (TCLs) are well-suited for the provision of demand-side response (DSR), due to their simple temperature set point operating rules and ubiquity in society. Examples include fridges, freezers, air-conditioners, hot-water tanks, heat pumps, and swimming pool pumps. Research into the possibility of using TCLs for grid balancing services began in the 1980s with key papers such as [1–4]. Despite the existence of technology for creating frequency-sensitive TCLs for nearly 40 years, implementation remains limited to a relatively small number of trials [5–9]. There are a number of reasons for the absence of large roll-outs of highly distributed DSR schemes [10, 11]. Historically, control paradigms from both technical and economic perspectives have been established for service provision from a (relatively) small number of large power plants. Understandably, the critical nature of electricity grid operation and security deters potentially risky changes and experimentation, and so a great deal of motivation is required for shifts away from traditional approaches. Service availability and reliability can be improved by splitting a service between a multitude of providers, compared to a single unit which will become completely unavailable in the event of a fault or scheduled repair [9, 12]. Yet there is also inherent complexity and potentially reduced predictability in procuring services from thousands of very small demand-side resources if they act autonomously, which is undoubtedly an obstacle to be overcome. Effects on consumers and their appliances will also be of concern to potential participants. Finally, it will be crucial for the success of any scheme, to adequately address the requirements for minimum participation numbers and develop the right business models that ensure fair rewards and effective incentives.The changing energy landscape in the 21st century has brought a new focus to the use of TCLs for electricity grid support and a wealth of literature on the topic [5, 9–11, 13–37]. Work has varied in nature from mathematical frameworks to numerical simulations and real-world trials. Most of the theory can be applied to any type of TCL, and simulations have covered many different possible TCL technologies. A variety of control schemes have been proposed, many of which are discussed and compared in [10, 19, 35]. There are two main classes into which these types of TCL control schemes can be divided:centralised anddecentralised control (with a spectrum in between). Their key features and comparative advantages and disadvantages are summarised in Table 1. It is widely accepted that if millions of TCLs could be used for frequency response they could potentially provide a valuable resource for the system. However, if each device needed constant communication with a central controller, sending data about its temperature and switching history and receiving operation instructions, then the economics and security risks would severely outweigh the benefits of the service. Public perception of the service is also vital for the implementation of any control scheme that involves appliances in people’s homes. For these reasons we choose to focus on decentralised control for our research. A better understanding, however, of the potential undesirable side-effects of decentralised control, is required before any control strategy could be put in place.Table 1
Comparison of centralised and decentralised TCL control strategies. Informed by, for example, [19, 38].
Centralised control
Decentralised control
Key features
(i) TCLs instructed by a central controller
(i) Autonomous local control
(ii) 2-way communication in all TCLs
(ii) Control scheme established once, may be updated periodically
Advantages
(i) Highly controllable
(i) No communications infrastructure required
(ii) Reasonably predictable
(ii) No security risks
(iii) Very fast response possible
Disadvantages
(i) Establishing and maintaining a secure communications network are very expensive(ii) Response time limited by communication speed (iii) Vast amounts of data to manage (iv) Data and appliance control security risks (v) Negative public perceptions of external control of home appliances
(i) Response is less predictable than with centralised control(ii) Synchronisation and instability effects possible and not yet fully understood (iii) Errors and noise in local frequency measurements more likelyThe challenges of implementing a decentralised control scheme largely centre around the propensity for TCL temperature cycling to become synchronised. A number of simulations in the literature indicate TCL synchronisation following a frequency disturbance, for example, [13, 21, 22, 25, 30–32]. As a result, various control schemes have been proposed that aim to prevent such behaviour. A popular choice is some form of stochastic temperature set point control [11, 13, 33, 34, 39, 40]. For example, [11] simulates a heterogeneous population of electric heaters in the Nordic power grid, with deterministic control to respond to sudden frequency fluctuations followed by randomised switch times as the system returns to normal. Domestic fridges are modelled in [13] as Markov-jump linear systems where the on/off switching is governed by transition probability rates rather than temperature set points. These rates are determined by choosing the desired population average temperature or duty cycle, and the temperature probability density is steered towards a desired distribution. While stochastic control schemes do offer solutions to the potential for demand synchronisation, they are typically accompanied by two disadvantages. Firstly, naive stochastic switching can involve a TCL switching twice (or more) within a short time frame, which is detrimental to serving its purpose, and can cause wear on the device (though non-Markovian approaches can avoid this). Secondly, learning that home appliances will randomly switch on and off could cause negative public perceptions of a DSR scheme. We believe further study of potential deterministic schemes would be beneficial before putting attention on stochastic schemes, particularly given the natural diversity in TCL populations that could prevent synchronisation phenomena.An alternative to direct instructions and frequency signals is the use of a price signal that TCLs could respond to. The advantages of price signals are that it is possible to measure the financial benefits to consumers of DSR participation, and individual consumers could potentially make their own choices about the value they place on service disruption at, say, given times of day. However, current price signals typically change on half-hourly or at least several-minute time scales, which makes them ill-suited for dynamic frequency response. Reviews on the use of price signals for demand response can be found in [41–43]. A related approach proposed in [44] involves each device calculating the price directly from the grid frequency, and the authors argue that a well-chosen design for the controller and frequency-price coupling may prevent possible oscillations and instabilities. However, the drawbacks of price-based demand response, as remarked upon in [27], are the potential for synchronisation and the exposure of customers to price volatility, which could prevent sufficient participation for success.In [13] Angeli and Kountouriotis offer theoretical arguments for the long-term tendency of the system towards TCL synchronisation. It is reasoned that any “small periodic ripples in power system frequency will gradually entrain oscillations of refrigerators that have similar frequencies of oscillation, thus reinforcing the frequency ripple and eventually leading to an even larger number of entrained refrigerators.” We concur with the mathematical reasoning presented. However, we believe that further reasoning and inquiry are required for a more complete understanding of this phenomenon.A population of TCLs responding to frequency through preset deterministic rules with no other communications or control can be thought of as a system of coupled oscillators. Synchronisation phenomena emerging from systems of coupled oscillators have been studied in many contexts, from neural signals in the brain to flashing fireflies [45]. The Kuramoto framework was developed [46, 47] which elegantly describes basic features of such systems and allows for stability analysis. Inspired by results for the Kuramoto model, in this article we explore the stability of a system of TCLs and grid frequency using techniques from dynamical systems and agent-based modelling.We study the potential for a population of TCLs to support grid frequency and explore the possibilities for frequency-sensitive TCLs to cause instabilities due to cycle synchronisation. We use techniques from dynamical systems stability analysis along with simulations that incorporate data from the British power grid. In Section2 we introduce our model for a population of TCLs and the electricity grid frequency and explain our choices of parameter values. In Section 3.1 we present a stability analysis of the nominal frequency in the presence of a uniformly distributed population of TCLs. Section 3.2 solves the behaviour of a fully synchronised TCL population and analyses the stability of the population under a split into two groups. Section 4.1 describes our simulations of a large population of TCLs using the model in Section 3.1. Section 4.2 describes our simulations of a population of TCLs acting on the system in the presence of other frequency response providers and naturally occurring frequency fluctuations, including simulations of a heterogeneous population. Section 5 contains our final discussion and conclusions.


## 2. Modelling TCLs and Electricity Grid Frequency

The modelling is kept appliance-neutral where possible, but it is set up for cooling devices such as fridges, (fridge-)freezers, and air-conditioners and would need to be altered in minor ways for other appliances such as heat pumps and hot-water tanks. We make the following assumptions:(i)
Electricity grid frequency is the same everywhere on the network and there are no inter-area oscillations [48] (therefore all machines on the power grid are assumed to rotate in synchrony).(ii)
All TCLs sense frequency deviations with negligible measurement delay or measurement error.(iii)
All system parameters remain constant over time.(iv)
Fridges and freezers are not affected by the fridge/freezer door being opened, or by the addition or removal of food (in effect we assume this never occurs).(v)
TCLs have continuous thermostat control (in temperature and time) and can therefore sense and implement temperature/set point changes with infinite precision.(vi)
TCLs consume constant power when on and zero power when off and are controlled only by the rules outlined in the model.These assumptions allow us to create a tractable model for analytic study. Assumptions (iii) and (iv) are probably the easiest and most natural to relax first and could be relaxed by adding time-dependent forcing effects. For most of the paper we consider a population of identical TCLs, but our formulation can be extended easily to an inhomogeneous population and we believe the effects of sufficient diversity will be stabilising, as supported by our final simulation.


### 2.1. Individual TCLs

For the temperature cycling of a TCL we adopt the linear model and notation presented in [13]. Let the temperature of a TCL at time t be denoted by T(t), the cooling/heating coefficient by α, and the asymptotic temperatures that the TCL would reach if left on and off indefinitely by TON and TOFF, respectively. Then(1)T˙t=αTON-TtwhentheTCLisonαTOFF-TtwhentheTCLisoff.A (cooling) TCL will switch off when the temperature reaches its lower temperature set point T- and switch on when it reaches its upper temperature set point T+. We choose to make these set points sensitive to system frequency deviations away from 50 Hz, denoted f(t) (i.e.,f(t)= Frequency(t)-50Hz). Insufficient generation to meet demand causes f<0 and so we need the TCLs to reduce their power consumption to bring f back to zero. We implement this by increasing the temperature set points so that the TCLs switch off sooner/stay off for longer. Oversupply of electricity to the grid causes f>0, and so in this case we decrease the temperature set points to increase overall power consumption. Thus we define our frequency-sensitive temperature set points,(2a)T-ft≔T-0-β-ftlowerswitchoffsetpoint(2b)T+ft≔T+0-β+ftupperswitchonsetpoint,where β-, β+ are positive constants that determine the sensitivity of the lower and upper temperature set points to frequency deviations. T-0 and T+0 are the uncoupled (a fridge is “uncoupled” from the grid frequency if β-=β+=0) temperature set points, which we typically take to be 2∘C and 7∘C, respectively. This framework is very similar to that suggested in [30], although we allow the upper and lower temperature set points to have different sensitivities to the frequency (β- and β+).We can solve (1) for the temperature of a TCL at time t. If a TCL has temperature T0 at time t0 and does not switch on or off before time t then the temperature T(t) is given by(3a)Tt=T0-TONe-αt-t0+TONwhenon(3b)Tt=T0-TOFFe-αt-t0+TOFFwhenoff.We can rearrange (3a) and (3b) and solve for the on and off durations τON and τOFF, respectively, assuming constant grid frequency:(4a)τONf=1αlog⁡T+f-TONT-f-TON(4b)τOFFf=1αlog⁡TOFF-T-fTOFF-T+f.These variables will be useful when we consider the equilibrium of the system, in which the temperature set points become fixed. In the traditional case when TCLs are uncoupled from the grid (or the special case f≡0) their “natural” on and off cycle durations, τON0 and τOFF0, are given by(5a)τON0=1αlogT+0-TONT-0-TON(5b)τOFF0=1αlogTOFF-T-0TOFF-T+0.In order for the TCLs to operate properly they need to cycle on and off, and so we require that(6)TON<T-ft<T+ft<TOFF∀t.We also need a TCL to respond “appropriately” to a change in frequency, that is to say, for the average power consumption over one cycle to increase when the frequency increases and decrease when the frequency decreases. It is shown in Appendix A that a sufficient condition to ensure this is(7)β+β-∈TOFF-T+TOFF-T-,T+-TONT--TONwhich is a nonempty interval (notably containing 1).


### 2.2. Electricity Grid Frequency

A simplified equation for the frequencyF of a power system can be determined by Newton’s 2nd Law of Motion or the derived equation for energy. If we let f≔F-F0, where F0 is the nominal grid frequency (50 Hz in Europe) and linearise about F0, then we obtain [26](8)Mdfdt+Dft=ΔPg-ΔPland for brevity we introduce new variables along with explicit consideration of TCL power consumption(9)dfdtt=cΔP-ρtPc-γft,whereM≔4π2IF0 stands for 2π times nominal angular momentum of the rotating masses in the system,I stands for total inertia of the rotating masses of the system,D stands for damping factor representing the natural frequency dependence of the load alongside the damping provided by synchronous generator damper windings,ΔPg stands for change in total active power generation, compared to a reference level,ΔPl stands for change in total active power load, compared to a reference level,c≔1/M stands for inverse nominal angular momentum, introduced for brevity,ΔP stands for “surplus power generation for the TCLs;” total system active power generation minus total system active power load, excluding TCL power consumption,ρ stands for proportion of TCLs switched on,Pc stands for power consumed by TCL population when all switched on,γ≔D/M is a variable introduced for brevity.We make the simplifying assumption in Sections2 and 3 that the “surplus” power generation on the system for TCL consumption ΔP is a constant. We use the “∗” notation to denote equilibrium values. In equilibrium (10)cΔP-ρ∗Pc-γf∗=0,hence(11)f∗=cγΔP-ρ∗Pc,and therefore we can rewrite our equation for f˙ in terms of deviations from equilibrium values:(12)f~˙t=cPcρ∗-ρt-γf~,where(13)f~≔f-f∗.


### 2.3. Parameter Choices

We take as reference the Great Britain (GB) electricity system. This covers England, Scotland, and Wales. In 2015 approximately 10.4 m households in the UK, which also includes Northern Ireland, owned a fridge and 19.1 m households owned a fridge-freezer [49]. In the same year approximately 2.8% of the population lived in Northern Ireland [50]. If we assume that the average number of people per household is the same in Northern Ireland and in GB and an even distribution of fridge and fridge-freezer ownership, then approximately 10.1 m and 18.6 m households in GB owned a fridge and fridge-freezer, respectively. If using TCLs for frequency response became standard practice, that would mean that a very large number of appliances could participate in frequency response. We model the case of 1 million fridges participating in frequency response, which corresponds to roughly 10% of fridges in GB. We take the power consumed by an individual fridge when switched on, p, to be 70 W, as assumed in [32, 33]. This means that we let p=7×10-5 MW and the total power consumption if all fridges were switched on, Pc=7×10-5×106=70 MW. Using our approximation for f˙(t) [51], c=50/(2Ek). Our GB system data (discussed later) gives an approximate average value for total stored kinetic energy Ek; Ek=2.5×105 MVAs (note that MVAs = MJ), and so c=1×10-4. We let ρ∗ vary between 0 and 1 by changing ΔP. Our parameters are summarised in Table 2, and throughout this paper take these values unless stated otherwise.Table 2
Parameter values assumed, unless stated otherwise.
Parameter
Value
Units
TOFF
20
°C
TON
−26
°C
T-0
2
°C
T+0
7
°C
α
1.808
×
1
0
-
4
s-1
β+
2.4
°C·Hz−1
β-
2.4
°C·Hz−1
c
1
×
1
0
-
4
Hz(MVAs)−1
γ
0
s−1
p
7
×
1
0
-
5
MW
Pc
70
MW
ΔP
0.3355
MW


## 2.1. Individual TCLs

For the temperature cycling of a TCL we adopt the linear model and notation presented in [13]. Let the temperature of a TCL at time t be denoted by T(t), the cooling/heating coefficient by α, and the asymptotic temperatures that the TCL would reach if left on and off indefinitely by TON and TOFF, respectively. Then(1)T˙t=αTON-TtwhentheTCLisonαTOFF-TtwhentheTCLisoff.A (cooling) TCL will switch off when the temperature reaches its lower temperature set point T- and switch on when it reaches its upper temperature set point T+. We choose to make these set points sensitive to system frequency deviations away from 50 Hz, denoted f(t) (i.e.,f(t)= Frequency(t)-50Hz). Insufficient generation to meet demand causes f<0 and so we need the TCLs to reduce their power consumption to bring f back to zero. We implement this by increasing the temperature set points so that the TCLs switch off sooner/stay off for longer. Oversupply of electricity to the grid causes f>0, and so in this case we decrease the temperature set points to increase overall power consumption. Thus we define our frequency-sensitive temperature set points,(2a)T-ft≔T-0-β-ftlowerswitchoffsetpoint(2b)T+ft≔T+0-β+ftupperswitchonsetpoint,where β-, β+ are positive constants that determine the sensitivity of the lower and upper temperature set points to frequency deviations. T-0 and T+0 are the uncoupled (a fridge is “uncoupled” from the grid frequency if β-=β+=0) temperature set points, which we typically take to be 2∘C and 7∘C, respectively. This framework is very similar to that suggested in [30], although we allow the upper and lower temperature set points to have different sensitivities to the frequency (β- and β+).We can solve (1) for the temperature of a TCL at time t. If a TCL has temperature T0 at time t0 and does not switch on or off before time t then the temperature T(t) is given by(3a)Tt=T0-TONe-αt-t0+TONwhenon(3b)Tt=T0-TOFFe-αt-t0+TOFFwhenoff.We can rearrange (3a) and (3b) and solve for the on and off durations τON and τOFF, respectively, assuming constant grid frequency:(4a)τONf=1αlog⁡T+f-TONT-f-TON(4b)τOFFf=1αlog⁡TOFF-T-fTOFF-T+f.These variables will be useful when we consider the equilibrium of the system, in which the temperature set points become fixed. In the traditional case when TCLs are uncoupled from the grid (or the special case f≡0) their “natural” on and off cycle durations, τON0 and τOFF0, are given by(5a)τON0=1αlogT+0-TONT-0-TON(5b)τOFF0=1αlogTOFF-T-0TOFF-T+0.In order for the TCLs to operate properly they need to cycle on and off, and so we require that(6)TON<T-ft<T+ft<TOFF∀t.We also need a TCL to respond “appropriately” to a change in frequency, that is to say, for the average power consumption over one cycle to increase when the frequency increases and decrease when the frequency decreases. It is shown in Appendix A that a sufficient condition to ensure this is(7)β+β-∈TOFF-T+TOFF-T-,T+-TONT--TONwhich is a nonempty interval (notably containing 1).


## 2.2. Electricity Grid Frequency

A simplified equation for the frequencyF of a power system can be determined by Newton’s 2nd Law of Motion or the derived equation for energy. If we let f≔F-F0, where F0 is the nominal grid frequency (50 Hz in Europe) and linearise about F0, then we obtain [26](8)Mdfdt+Dft=ΔPg-ΔPland for brevity we introduce new variables along with explicit consideration of TCL power consumption(9)dfdtt=cΔP-ρtPc-γft,whereM≔4π2IF0 stands for 2π times nominal angular momentum of the rotating masses in the system,I stands for total inertia of the rotating masses of the system,D stands for damping factor representing the natural frequency dependence of the load alongside the damping provided by synchronous generator damper windings,ΔPg stands for change in total active power generation, compared to a reference level,ΔPl stands for change in total active power load, compared to a reference level,c≔1/M stands for inverse nominal angular momentum, introduced for brevity,ΔP stands for “surplus power generation for the TCLs;” total system active power generation minus total system active power load, excluding TCL power consumption,ρ stands for proportion of TCLs switched on,Pc stands for power consumed by TCL population when all switched on,γ≔D/M is a variable introduced for brevity.We make the simplifying assumption in Sections2 and 3 that the “surplus” power generation on the system for TCL consumption ΔP is a constant. We use the “∗” notation to denote equilibrium values. In equilibrium (10)cΔP-ρ∗Pc-γf∗=0,hence(11)f∗=cγΔP-ρ∗Pc,and therefore we can rewrite our equation for f˙ in terms of deviations from equilibrium values:(12)f~˙t=cPcρ∗-ρt-γf~,where(13)f~≔f-f∗.


## 2.3. Parameter Choices

We take as reference the Great Britain (GB) electricity system. This covers England, Scotland, and Wales. In 2015 approximately 10.4 m households in the UK, which also includes Northern Ireland, owned a fridge and 19.1 m households owned a fridge-freezer [49]. In the same year approximately 2.8% of the population lived in Northern Ireland [50]. If we assume that the average number of people per household is the same in Northern Ireland and in GB and an even distribution of fridge and fridge-freezer ownership, then approximately 10.1 m and 18.6 m households in GB owned a fridge and fridge-freezer, respectively. If using TCLs for frequency response became standard practice, that would mean that a very large number of appliances could participate in frequency response. We model the case of 1 million fridges participating in frequency response, which corresponds to roughly 10% of fridges in GB. We take the power consumed by an individual fridge when switched on, p, to be 70 W, as assumed in [32, 33]. This means that we let p=7×10-5 MW and the total power consumption if all fridges were switched on, Pc=7×10-5×106=70 MW. Using our approximation for f˙(t) [51], c=50/(2Ek). Our GB system data (discussed later) gives an approximate average value for total stored kinetic energy Ek; Ek=2.5×105 MVAs (note that MVAs = MJ), and so c=1×10-4. We let ρ∗ vary between 0 and 1 by changing ΔP. Our parameters are summarised in Table 2, and throughout this paper take these values unless stated otherwise.Table 2
Parameter values assumed, unless stated otherwise.
Parameter
Value
Units
TOFF
20
°C
TON
−26
°C
T-0
2
°C
T+0
7
°C
α
1.808
×
1
0
-
4
s-1
β+
2.4
°C·Hz−1
β-
2.4
°C·Hz−1
c
1
×
1
0
-
4
Hz(MVAs)−1
γ
0
s−1
p
7
×
1
0
-
5
MW
Pc
70
MW
ΔP
0.3355
MW


## 3. Stability Analysis

Concerns that frequency-responsive TCLs controlled by deterministic rules will exhibit herding behaviour and create frequency oscillations have been raised in various previous works, either by predictions or examples from simulations [13, 21, 22, 25, 30–32]. The simplicity of our model allows for a rigorous mathematical treatment of the stability of a population of TCLs responding according to the scheme introduced above. In the first part of this section we model a TCL population as a continuum on the temperature cycle and linearise about the equilibrium discussed in Section 2.2. In the second part of this section we consider the opposite extreme for a TCL distribution (one or two Dirac delta distributions), solving for the behaviour of a fully synchronised population of TCLs and studying the dynamics of two synchronised groups.


### 3.1. Uniform Distribution of TCLs

We begin by studying the stability of a population of TCLs uniformly distributed in phase (meaning the time since last switch on). This means that under constant temperature set point conditions the TCLs would switch on at a constant rate and switch off at a (possibly different) constant rate (note that since TCLs heat (or cool) at different rates depending on their current temperature, uniformly distributing the TCLs within each part of the cycle does not correspond to uniformly distributing the population over the temperature scale). In the context of the Kuramoto model this is usually referred to as the “incoherent solution,” for example [47, 52]. Just as in Strogatz and Mirollo’s treatment of the Kuramoto model [52], we model the infinite-N limit of a population of TCLs as a continuum of TCLs distributed over an interval with periodic boundary conditions.In order to obtain a tractable model, comparable to the Kuramoto model, three key challenges must be addressed. Firstly, the TCL temperature cycling is described by the piecewise-smooth nonlinear function (see (3a) and (3b)), with nondifferentiability at each temperature set point. Secondly, these set points are continuously changing with grid frequency, and so any map to a periodic regime must be sufficiently flexible to accommodate this. Finally, in order to know a TCL’s rate of change of temperature at any time, one needs to know both its current temperature and its current (on/off) state. We therefore propose a new modelling framework to overcome these challenges and permit stability analysis for the model.We map each TCL with temperature and on/off state to a pointθ on the interval [-1,1), in such a way that θ dictates both the temperature and the state of a TCL. The switched off TCLs are mapped to the interval [-1,0) and the switched on TCLs are mapped to [0,1). Then we define the position θ(t) of a TCL at time t with temperature T(t) and state on or off by(14)θt=θONt=1ατONftlog⁡T+ft-TONTt-TON∈0,1ifonθOFFt=1ατOFFftlog⁡TOFF-T+ftTOFF-Tt∈-1,0ifoff.Note that the model implicitly assumes that the temperature set points never change fast enough to leave a TCL outside of the interval [T-(f(t)),T+(f(t))]. Since in this paper we use this model for only linear stability analysis about the equilibrium, we consider this to be a reasonable assumption. Our choice of θ means that uniformly distributing a population of TCLs over each part of the temperature cycle (as discussed above) corresponds to a uniform distribution of on and off TCLs in their respective halves of θ-space.As in [52], we consider the population density in θ-space; let u(θ,t)dθ denote the fraction of TCLs that lie between θ and θ+dθ at time t. Then u is nonnegative, with period length 2 in θ, and satisfies the normalisation(15)∫-1+1uθ,tdθ=1for all t. The evolution of u is governed by the continuity equation [53](16)∂u∂t+∂∂θuv=0,where v is the velocity of a TCL in θ-space, v(θ,t)≔θ˙(t). Differentiating (14) gives(17a)vONθ,t=1τONft1+1αϕONftθ-β+T+ft-TONf˙t(17b)vOFFθ,t=1τOFFft1+1αϕOFFftθ+β+TOFF-T+ftf˙t,where(18a)ϕONf≔β+T+f-TON-β-T-f-TON(18b)ϕOFFf≔β+TOFF-T+f-β-TOFF-T-f.Note that, for β+/β- satisfying (7), ϕONf<0 and ϕOFFf>0. Under a constant grid frequency, θ˙ON and θ˙OFF are constants. In equilibrium u∗ we have u˙∗=0, and therefore (16) implies (19)uON∗θ=k0vON∗θ;uOFF∗θ=k0vOFF∗θ,for some constant k0. Since f˙∗=0, from (17a) and (17b) we have(20)vON∗=1τON∗;vOFF∗=1τOFF∗.Then for all θ∈[-1,0),[0,1), respectively,(21)uON∗θ=k0τON∗uOFF∗θ=k0τOFF∗and k0 is determined by the normalisation criterion (15),(22)k0=1τON∗+τOFF∗.The proportion of TCLs switched on, ρ(t), is given by(23)ρt=∫01uθ,tdθ.In equilibrium ρ(t)=ρ∗ (12), (24)ρ∗=∫01u∗θ,tdθ=τON∗τON∗+τOFF∗.We introduce the notation “•” to imply that an equation holds for the variable with either of two values, “on” or “off.” Our approach is to perturb the system about the equilibrium (u∗,f∗) by a small amount τ•∗η(θ,t) and to consider the evolution of the perturbation. By (15) the perturbation satisfies(25)∫-1+1ηθ,tdθ=0.We write(26)u•=k0+ηθ,tτ•∗(27)v•=1τ•∗1+wθ,tso that (16) becomes(28)τ•∗∂∂tη+∂∂θk0+k0w+η+ηw=0.Taking the first-order approximation yields(29)τ•∗∂∂tη+k0∂∂θw+∂∂θη=0.Rearranging (27) for w and substituting (17a) and (17b) for v• give(30a)wON=1αϕON∗θ-β+T+∗-TONf˙t-δτONtτON∗wOFF=1αϕOFF∗θ+β+TOFF-T+∗f˙t-δτOFFtτOFF∗,δτONt=-ϕON∗f~tα;δτOFFt=-ϕOFF∗f~tα.Hence (31)∂∂θw•t=1αϕ•∗f˙+wtθ=0+-wtθ=0-δθ+wtθ=-1-wtθ=1δθ-1∂∂θw•t=1αϕ•∗f˙-ν0δθ+ν1δθ-1f˙t+μαδθ-1-δθf~,where we have defined(32a)ν0≔β+T+∗-TON+β+TOFF-T+∗>0(32b)ν1≔β-T-∗-TON+β-TOFF-T-∗>0(32c)μ≔ϕOFF∗τOFF∗-ϕON∗τON∗>0ifβ+β-satisfies7withT±∗.We have a time-invariant linear system (29), and so it is natural to look for solutions for which the time dependence of our variables f~ and η is eλt; λ∈C is called an eigenvalue of the system. Defining k≔k0/α and renaming f~ to f, (29) becomes(33)τ•∗λη+∂η∂θ+kϕ•-ν0δθ+ν1δθ-1λf+kμδθ-1-δθf=0.We introduce an integrating factor so that on the open intervals (-1,0)∪(0,1) we can find an expression for η(θ):(34)∂∂θeλτ•∗θη+kϕ•λfeλτ•∗θ=0∴ηθ=η•0+kfϕ•∗τ•∗e-λτ•∗θ-kfϕ•∗τ•∗.At the discontinuities θ=0 and θ=±1,(35a)ηON0-ηOFF0=kfλν0+μ(35b)ηOFF-1-ηON1=-kfλν1+μ.We can use (34) to find expressions for η(-1) and η(1) and substitute these into (35b). After substitution for ηOFF(0) (or ηON(0)) using (35a) and rearrangement we arrive at(36a)ηON0gλ=-kfϕON∗τON∗gλ+λν1-ν0eλτOFF∗(36b)ηOFF0gλ=-kfϕOFFτOFF∗gλ+λν1-ν0e-λτON∗,where(36c)gλ=eλτOFF∗-e-λτON∗.Rewriting our equation for the rate of change of grid frequency near equilibrium (9) as(37)f˙t=-γft-cPcτON∗∫01ηθ,tdθand setting f˙=λf give (38)∫01ηθ,tdθ=-λ+γfcPcτON∗.Integrating (34) over [0,1) in θ (the switched on TCLs), setting the resulting expression equal to the right hand side of (38), and substituting our expression in (36a) for ηON(0) establish the following implicit equation for λ: (39)λ+γ-ZϕON∗gλ=Zν1-ν0eλτOFF∗1-e-λτON∗,where we have defined Z≔kcPc, which reflects the strength of the effect of the TCLs on grid frequency.WhenZ=0 (no effect of the TCLs on the grid frequency) the eigenvalue equation (39) reduces to (λ+γ)g(λ)=0, so the eigenvalues are λ=-γ and λ=2nπi/τON∗+τOFF∗ for n∈Z (the roots of g(λ)=0). It can also be seen from (39) that for all Z there is an eigenvalue λ=0. It corresponds to conservation of the number of TCLs. This eigenvalue 0 is removed by the normalisation condition (25). The real and imaginary parts of λ that solve (39) can be solved for numerically, using, for example, [54]. Figure 1 shows numerical solutions for the first five eigenvalues above (or on) the real axis for the parameter values given in Table 2 in Section 2.3 and allowing Z to vary from its value Z0 derived from the table, by Z=hZ0. There is an infinite sequence of eigenvalues going upwards and their reflections in the real axis. Increasing Z from zero by powers of 10 is seen to decrease the real part of the eigenvalues from zero and therefore the system is stable to small perturbations.Figure 1
Numerical solutions for the first five eigenvalues above the real axis (there is an infinite sequence going further along the imaginary axis, and they are reflected in the real axis). We use multiplierh to increase Z, and the real part of each eigenvalue we have followed decreases from 0 as Z increases from 0.This is a surprising result because intuitively identical TCLs are vulnerable to synchronisation which would cause instabilities on the system, which is the general view in the literature as discussed previously. The result is not due to the damping constantγ, because we chose γ=0 so as not to mask the effect of the TCLs. What the analysis does not tell us is how small any perturbations would have to be for a population of TCLs to have a stabilising effect on grid frequency. It might be that a larger perturbation than valid for linearisation leads to instability. In Section 4 we study the effects of different sized perturbations using simulations and indeed find growth of synchronisation. In the next section we consider the behaviour of a population of TCLs under the opposite type of perturbation; namely, all TCLs synchronised into one or two groups.


### 3.2. Synchronised Groups of TCLs

In the previous section we studied the stability of a uniformly distributed (continuum) population of TCLs at the 50 Hz equilibrium and found it to be stable almost everywhere in parameter space. In this section we consider the opposite extreme of possible TCL distributions, the Dirac delta distribution. That is to say, we explore the behaviour of a fully synchronised population of TCLs, all switching on and off at the same time, all with the same temperature, and (again) identical parameters. This is equivalent to a single TCL with the power consumption of the whole population.


#### 3.2.1. Mapping the Switch Times

We begin by constructing a map from one (whole population) switch on event to the next. We show that under certain conditions such a mapping is a contraction. Let the subscriptn denote the nth switch on and nth switch off event. Without loss of generality, suppose that after our initial start time t0 the next switch event is the population switching on. This implies that, for all n∈N, tnOFF>tnON. Figure 2 illustrates the notation. Hence the amount of time the population spends switched on following the nth switch on event is given by(40a)tnOFF-tnON=1αlog⁡T+0-β+fnON-TONT-0-β-fnOFF-TON,where fnON,fnOFF are the frequencies at the nth switch on and off times. The amount of time spent switched off following the nth switch off is given by(40b)tn+1ON-tnOFF=1αlog⁡TOFF-T-0+β-fnOFFTOFF-T+0+β+fn+1ON.Assuming, as for the numerical analysis in Section 3.1, that the system has no damping, we set γ=0 in (12) for f˙(t). In a synchronised population, at time t all TCLs are either on (ρ(t)=1) or off (ρ(t)=0). Then we can define constants cON,cOFF>0 such that(41)f˙=-cON≔cPcρ∗-1whenthepopulationison+cOFF≔cPcρ∗whenthepopulationisoff.Hence the values of f at the switch off and on times are given by the piecewise-linear functions(42a)fnOFF=fnON-cONtnOFF-tnON(42b)fn+1ON=fnOFF+cOFFtn+1ON-tnOFF.After substituting for the switching times using (40a) and (40b) and rearranging, these become(43a)fnOFF-cONαlog⁡T-0-β-fnOFF-TON=fnON-cONαlog⁡T+0-β+fnON-TON(43b)fn+1ON+cOFFαlog⁡TOFF-T+0+β+fn+1ON=fnOFF+cOFFαlog⁡TOFF-T-0+β-fnOFF.Now since each side of (43a) and (43b) are functions of only one of the fn• variables, we can explicitly name them as such(44a)ϕON-fnOFF≔fnOFF-cONαlog⁡T-0-β-fnOFF-TON(44b)ϕON+fnON≔fnON-cONαlog⁡T+0-β+fnON-TON(44c)ϕOFF-fnOFF≔fnOFF+cOFFαlog⁡TOFF-T-0+β-fnOFF(44d)ϕOFF+fn+1ON≔fn+1ON+cOFFαlog⁡TOFF-T+0+β+fn+1ON.Each of the four ϕ functions is increasing and therefore invertible, and so we can write(45a)fnOFF=ϕON--1ϕON+fnON(45b)fn+1ON=ϕOFF+-1ϕOFF-fnOFFand therefore(45c)fn+1ON=ϕOFF+-1ϕOFF-ϕON--1ϕON+fnON,which is a mapping from the frequency at one switch on event to the frequency at the next. The mapping is a contraction iff (46)ϕOFF+-1ϕOFF-ϕON--1ϕON+′<1(47)iffϕOFF-′ϕOFF+′ϕON+′ϕON-′<1evaluatedattheappropriateplaces.Note that(48)ϕOFF-′ϕOFF+′=1+β-cOFF/αTOFF-Tn+1+β+cOFF/αTOFF-Tn+1+<1iffβ+β->TOFF-Tn+1+TOFF-Tn-.Similarly(49)ϕON+′ϕON-′=1+β+cON/αTn+-TON1+β-cON/αTn--TON<1iffβ+β-<Tn+-TONTn--TON.Therefore a sufficient condition for the mapping to be a contraction is that(50)β+β-∈TOFF-Tn+1+TOFF-Tn-,Tn+-TONTn--TONwhich is a nonempty interval (containing {1}), so long as TON<Tn-<Tn+<TOFF and Tn-<Tn+1+ for all n. It is worth recalling our earlier condition on the values of β± (7) which also imposed that β+/β- belong to an open interval containing {1}.Figure 2
Illustration of thenth and n+1th switching events of the fully synchronised population and the frequency-sensitive temperature set points T-(f(t)) and T+(f(t)).


#### 3.2.2. Solving for the Periodic Solution

The contraction property of the mappingfnON↦fn+1ON (45c) under the above conditions implies that there is an attracting fixed point so long as TON<Tn-<Tn+<TOFF, and hence a periodic solution for the synchronised population. We now seek to solve for this periodic solution. Denote by lON and lOFF the amount of time spent on and off during one (periodic) cycle, respectively. Since power consumption for the population is constant during each on/off phase, the frequency moves linearly between upper and lower values which we denote by f+ and f-. Therefore the temperature of the population will cycle between upper and lower set points, given by T+0-β+f+ and T-0-β-f-, respectively. Equations (42a) and (42b) show us that(51a)f-=f+-cONlON(51b)f+=f-+cOFFlOFF.The temperature evolution equations (3a) and (3b) allow us to express the switch on and switch off temperatures as follows:(52a)T+0-β+f+=T-0-β-f--TOFFe-αlOFF+TOFF(52b)T-0-β-f-=T+0-β+f+-TONe-αlON+TONwhich after substituting for f- using (51a) and rearranging become(53a)f+β-e-αlOFF-β+=T-0-TOFF+β-cONlONe-αlOFF+TOFF-T+0,(53b)f+β+e-αlON-β-=T+0-TONe-αlON+TON-T-0-β-cONlON.Now we have two equations in terms of f+, lON and lOFF, which we can combine into one equation and eliminate f+,(54)β+e-αlON-β-T-0-TOFF+β-cONlONe-αlOFF+TOFF-T+0=β-e-αlOFF-β+T+0-TONe-αlON+TON-T-0-β-cONlON.We can also express lOFF in terms of lON by summing (51a) and (51b) to give(55)cONlON=cOFFlOFFor,equivalently,1-ρ∗lON-ρ∗lOFF=0and so (53b) and (55) form a pair of coupled equations for lON and lOFF, which can be solved numerically. Figure 3 shows one temperature cycle for the single group under different choices for ρ∗. Denote by ρ0 the value of ρ when f=0. As ρ∗ gets further away from ρ0 the solutions drift further from the uncoupled temperature range (2–7°C). The cycle lengths are symmetric about ρ∗=1/2 but the TCLs consume more power per cycle as ρ∗ increases.Figure 3
One cycle of the single group solution for different values ofρ∗ when ρ0≈0.3355. They include values that lead to unrealistic results for real fridges but are there to illustrate the effect.To begin some analysis of the stability of this fully synchronised solution we address the question: given a population split into two synchronised groups, will the groups merge into one fully synchronised population, or will they remain distinct forever?


#### 3.2.3. Two-Group Dynamics

Suppose we have a population of frequency-sensitive TCLs that are split into two synchronised groups. We would like to understand the dynamics of the switch times, and we ask whether, given sufficient time, the groups will merge, or whether they will remain distinct, possibly settling down to separated periodic solutions. In particular, we consider the initial difference between the switch on timesΔt to be very small and the switch on temperatures very close to the single group periodic solution from the previous subsection.LetΓ denote the single group periodic solution, which cycles periodically through temperature space with temperature TΓ(t). As before, we denote the switched on duration in this solution by lON and the switched off duration by lOFF. Suppose that the population is split into two groups A and B, such that proportion σ belongs to group A, and proportion 1-σ belongs to group B. Suppose also that group B switches on at time t=0, followed soon after by group A switching on, at time t1>0. Then after a time period of length similar to lON group B switches off, which is again followed shortly after by group A switching off. After a time period similar to lOFF each of the groups then switches back on. We shall assume that the switching order does not change, since if they do swap, we need only repeat this process with σ replaced by 1-σ. Simulations show that the switching order will not continue to change indefinitely.We would like to compare the temperature cycles of these two groups with the single group periodic solutionΓ. Without loss of generality suppose that group B initially switches on at the same time as a fully synchronised population solution. We compare the cycling of the groups A and B using the following measures, along with all those shown in Figure 4. Let ΔT≔TB(0)-TA(0) and ΔT′≔TB(t4)-TA(t4), the temperature difference when B switches on the first and second times, respectively. In addition, let Δt≔t1-0=t1=ϵ1 and Δt′≔t5-t4, the time difference between the two groups switching on the first and second times, respectively. Further notation is shown in Figure 4.Figure 4
Linearisation about the single group solution. (a) Single group solutionTΓ(t). (b) Temperature cycling of groups A and B close to the single group solution.In order to calculateΔT′ and Δt′ we need to calculate the switch times and temperatures of the two groups at each switch event leading up to t5. Solving for the switch times and temperatures when there are two groups is a little more complicated than for the fully synchronised case. It requires solving the temperature set point equations using the system conditions at the previous switch and the equation for f˙ which now takes one of four values depending on which combination of groups is switched on (both, neither, A only, or B only). We begin by making the simplifying assumption β-=β+≔β. Now since group A is switching on at time t1 and group B switched off at time 0,(56a)TAt1=T+0-βft1ft1=f0-cPc1-σ-ρ∗t1f0=1βt+0-TB0∴TAt1=TB0+βcPc1-σ-ρ∗t1.In addition, by the temperature evolution equations,(56b)TAt1=TA0-TOFFe-αt1+TOFF.Equating (56a) and (56b) and introducing our new notation give(57)βcPc1-σ-ρ∗Δt+ΔT=TA0-TOFFe-αΔt-1.If we write TA(0)=TΓ(0)+δTA(0) and take δTA(0) and Δt small, then (58)ΔT=TΓ0+δTA0-TOFFe-αΔt-1-βcPc1-σ-ρ∗Δtand linearising in Δt gives(59)ΔT≈ξΔt,where(60)ξ≔αTOFF-TΓ0-βcPc1-σ-ρ∗.More generally, at each switch event we have the temperature evolution equations that describe the temperature of each group as a function of their temperature at the previous switch (such as (56b)) and an additional equation for the temperature of the switching group, using the temperature set point equations (such as (56a)). Writing ΔT′=TB(t4)-TA(t4) and linearising about the single group solution, we find in Appendix B that ΔT′=λΔT where(61)λ≔1-αTOFF-TONαTOFF-TΓ0-βcPc1-σ-ρ∗1-αTOFF-TONαTΓlON-TON+βcPcσ-ρ∗e-αL.So [-1,+1] is a left eigenvector of the linearised map in the space of δTAδTB with eigenvalue λ. We can plot λ against σ for various ρ∗ to see whether |λ|<1 (in which case the two groups merge into one) or whether |λ|>1 (they move apart). The results are shown in Figure 5. We find that the second eigenvalue is within the interval (-1,+1) for any σ and ρ∗ for our parameter values (taken from Table 2 with the exception of β=0.1 which has been reduced to limit the rate of change of the frequency and ensure model validity), see Appendix C for details, and therefore the stability is governed by λ.Figure 5
Solutions forλ (see (61)) (solid lines) for different values of ρ∗. Dashed lines show reflection in σ=1/2 to show the effect of reversing the switching order of the groups. Blue: ρ∗=0.1, yellow: ρ∗=0.2, green: ρ∗=0.3, and red: ρ∗=0.4. Black line shows the boundary of stability (stable below, unstable above). The results are identical when ρ∗ is replaced by 1-ρ∗. (b) shows an enlargement centred at σ=1/2, showing that either switching order of the groups leads to λ2>1 on a small interval of σ. In this case the groups never merge and in all other cases they will.
(a)
Full interval
(b)
EnlargementBy solving for the dividing caseλ=1 we can create a bifurcation diagram in terms of the parameters σ and ρ∗ to show where the single group solution is attracting and repelling. Figure 6 sketches the solution, along with the solution for the case when the switching order of the groups is reversed, found by replacing σ with 1-σ. If in either case (group A switching first or group B switching first) the solution is attracting, then the two groups will merge together into the one group solution. However, if both cases have unstable dynamics then the solutions will never merge. Simulations show that in this parameter region the two groups will settle down to a fixed phase distance apart. If the solution is attracting for one switching order and repelling for the other, we find that the typical behaviour is for a small separation in the unstable direction to grow until the phase difference becomes almost a whole cycle, when they merge. Figure 7 illustrates how the cycles of the two groups can change over time relative to one another, depending on which of the three regions in the bifurcation diagram their parameters belong to.Figure 6
Bifurcation diagram for the stability of the single group solution to splitting in two. Stable if the parameters lie in the yellow or blue regions (the groups will ultimately merge); unstable in the green parameter region (the groups will never merge). Boundary lines are solutions to (61) as a function of σ (or 1-σ to capture switching order reversal) and ρ∗.Figure 7
Illustration of the three types of cycling behaviour of two groups relative to one another, based on simulations. Arrows indicate the occurrence of many cycles and the central illustrations are snapshots of the cycling behaviour between the start and the final behaviour. Synchronisation occurs in cases (a) and (b), while in case (c) each group tends towards a fixed phase difference apart.What these results show is that when a population is split into two groups, if they are sufficiently similar in size then they will remain apart, effectively trying to counteract one another and balance the frequency fluctuations. Conversely, if one of the groups is significantly larger (“significantly” here depends on the size ofρ∗, and may be very small if ρ∗≈0.5) than the other then it will have too strong an effect on the frequency and “pull” on the smaller group’s cycle. The closer the proportion switched on in equilibrium is to the proportion switched off (i.e., the closer it is to 0.5), the more similar the groups have to be in size to remain distinct.With more than two groups of TCLs the modelling becomes far more complicated, since there are now far more possibilities to be considered for the switching order of the groups. Simulations have shown that for three groups it is possible for all three cycles to settle down to a fixed, separated pattern. This occurs if the groups are very similar in size, just as in the two-group case. Once one group is too large (or too small), the groups collapse into two, before synchronising completely. Taking the number of groups to infinity is equivalent to modelling a continuum of TCLs as in Section3.1. Taking all groups of equal size and uniformly distributed in phase θ is equivalent to the continuum population equilibrium studied earlier. From above we found analytically that small perturbations to the population distribution should relax back to the uniform distribution, that is, that the equilibrium was stable. Now we find that if the population is discretised into 2 (and hypothetically N) groups then so long as they are of close to equal sizes, they will attempt to settle the frequency back to its nominal value by “spreading out” their cycles.In reality we will never have a continuum of TCLs, and they may exhibit nonlinear dynamics not captured by our analysis. This motivates our use of simulations to gain further insights into how a large population of TCLs would behave according to our switching rules and how the grid frequency would be affected.


## 3.1. Uniform Distribution of TCLs

We begin by studying the stability of a population of TCLs uniformly distributed in phase (meaning the time since last switch on). This means that under constant temperature set point conditions the TCLs would switch on at a constant rate and switch off at a (possibly different) constant rate (note that since TCLs heat (or cool) at different rates depending on their current temperature, uniformly distributing the TCLs within each part of the cycle does not correspond to uniformly distributing the population over the temperature scale). In the context of the Kuramoto model this is usually referred to as the “incoherent solution,” for example [47, 52]. Just as in Strogatz and Mirollo’s treatment of the Kuramoto model [52], we model the infinite-N limit of a population of TCLs as a continuum of TCLs distributed over an interval with periodic boundary conditions.In order to obtain a tractable model, comparable to the Kuramoto model, three key challenges must be addressed. Firstly, the TCL temperature cycling is described by the piecewise-smooth nonlinear function (see (3a) and (3b)), with nondifferentiability at each temperature set point. Secondly, these set points are continuously changing with grid frequency, and so any map to a periodic regime must be sufficiently flexible to accommodate this. Finally, in order to know a TCL’s rate of change of temperature at any time, one needs to know both its current temperature and its current (on/off) state. We therefore propose a new modelling framework to overcome these challenges and permit stability analysis for the model.We map each TCL with temperature and on/off state to a pointθ on the interval [-1,1), in such a way that θ dictates both the temperature and the state of a TCL. The switched off TCLs are mapped to the interval [-1,0) and the switched on TCLs are mapped to [0,1). Then we define the position θ(t) of a TCL at time t with temperature T(t) and state on or off by(14)θt=θONt=1ατONftlog⁡T+ft-TONTt-TON∈0,1ifonθOFFt=1ατOFFftlog⁡TOFF-T+ftTOFF-Tt∈-1,0ifoff.Note that the model implicitly assumes that the temperature set points never change fast enough to leave a TCL outside of the interval [T-(f(t)),T+(f(t))]. Since in this paper we use this model for only linear stability analysis about the equilibrium, we consider this to be a reasonable assumption. Our choice of θ means that uniformly distributing a population of TCLs over each part of the temperature cycle (as discussed above) corresponds to a uniform distribution of on and off TCLs in their respective halves of θ-space.As in [52], we consider the population density in θ-space; let u(θ,t)dθ denote the fraction of TCLs that lie between θ and θ+dθ at time t. Then u is nonnegative, with period length 2 in θ, and satisfies the normalisation(15)∫-1+1uθ,tdθ=1for all t. The evolution of u is governed by the continuity equation [53](16)∂u∂t+∂∂θuv=0,where v is the velocity of a TCL in θ-space, v(θ,t)≔θ˙(t). Differentiating (14) gives(17a)vONθ,t=1τONft1+1αϕONftθ-β+T+ft-TONf˙t(17b)vOFFθ,t=1τOFFft1+1αϕOFFftθ+β+TOFF-T+ftf˙t,where(18a)ϕONf≔β+T+f-TON-β-T-f-TON(18b)ϕOFFf≔β+TOFF-T+f-β-TOFF-T-f.Note that, for β+/β- satisfying (7), ϕONf<0 and ϕOFFf>0. Under a constant grid frequency, θ˙ON and θ˙OFF are constants. In equilibrium u∗ we have u˙∗=0, and therefore (16) implies (19)uON∗θ=k0vON∗θ;uOFF∗θ=k0vOFF∗θ,for some constant k0. Since f˙∗=0, from (17a) and (17b) we have(20)vON∗=1τON∗;vOFF∗=1τOFF∗.Then for all θ∈[-1,0),[0,1), respectively,(21)uON∗θ=k0τON∗uOFF∗θ=k0τOFF∗and k0 is determined by the normalisation criterion (15),(22)k0=1τON∗+τOFF∗.The proportion of TCLs switched on, ρ(t), is given by(23)ρt=∫01uθ,tdθ.In equilibrium ρ(t)=ρ∗ (12), (24)ρ∗=∫01u∗θ,tdθ=τON∗τON∗+τOFF∗.We introduce the notation “•” to imply that an equation holds for the variable with either of two values, “on” or “off.” Our approach is to perturb the system about the equilibrium (u∗,f∗) by a small amount τ•∗η(θ,t) and to consider the evolution of the perturbation. By (15) the perturbation satisfies(25)∫-1+1ηθ,tdθ=0.We write(26)u•=k0+ηθ,tτ•∗(27)v•=1τ•∗1+wθ,tso that (16) becomes(28)τ•∗∂∂tη+∂∂θk0+k0w+η+ηw=0.Taking the first-order approximation yields(29)τ•∗∂∂tη+k0∂∂θw+∂∂θη=0.Rearranging (27) for w and substituting (17a) and (17b) for v• give(30a)wON=1αϕON∗θ-β+T+∗-TONf˙t-δτONtτON∗wOFF=1αϕOFF∗θ+β+TOFF-T+∗f˙t-δτOFFtτOFF∗,δτONt=-ϕON∗f~tα;δτOFFt=-ϕOFF∗f~tα.Hence (31)∂∂θw•t=1αϕ•∗f˙+wtθ=0+-wtθ=0-δθ+wtθ=-1-wtθ=1δθ-1∂∂θw•t=1αϕ•∗f˙-ν0δθ+ν1δθ-1f˙t+μαδθ-1-δθf~,where we have defined(32a)ν0≔β+T+∗-TON+β+TOFF-T+∗>0(32b)ν1≔β-T-∗-TON+β-TOFF-T-∗>0(32c)μ≔ϕOFF∗τOFF∗-ϕON∗τON∗>0ifβ+β-satisfies7withT±∗.We have a time-invariant linear system (29), and so it is natural to look for solutions for which the time dependence of our variables f~ and η is eλt; λ∈C is called an eigenvalue of the system. Defining k≔k0/α and renaming f~ to f, (29) becomes(33)τ•∗λη+∂η∂θ+kϕ•-ν0δθ+ν1δθ-1λf+kμδθ-1-δθf=0.We introduce an integrating factor so that on the open intervals (-1,0)∪(0,1) we can find an expression for η(θ):(34)∂∂θeλτ•∗θη+kϕ•λfeλτ•∗θ=0∴ηθ=η•0+kfϕ•∗τ•∗e-λτ•∗θ-kfϕ•∗τ•∗.At the discontinuities θ=0 and θ=±1,(35a)ηON0-ηOFF0=kfλν0+μ(35b)ηOFF-1-ηON1=-kfλν1+μ.We can use (34) to find expressions for η(-1) and η(1) and substitute these into (35b). After substitution for ηOFF(0) (or ηON(0)) using (35a) and rearrangement we arrive at(36a)ηON0gλ=-kfϕON∗τON∗gλ+λν1-ν0eλτOFF∗(36b)ηOFF0gλ=-kfϕOFFτOFF∗gλ+λν1-ν0e-λτON∗,where(36c)gλ=eλτOFF∗-e-λτON∗.Rewriting our equation for the rate of change of grid frequency near equilibrium (9) as(37)f˙t=-γft-cPcτON∗∫01ηθ,tdθand setting f˙=λf give (38)∫01ηθ,tdθ=-λ+γfcPcτON∗.Integrating (34) over [0,1) in θ (the switched on TCLs), setting the resulting expression equal to the right hand side of (38), and substituting our expression in (36a) for ηON(0) establish the following implicit equation for λ: (39)λ+γ-ZϕON∗gλ=Zν1-ν0eλτOFF∗1-e-λτON∗,where we have defined Z≔kcPc, which reflects the strength of the effect of the TCLs on grid frequency.WhenZ=0 (no effect of the TCLs on the grid frequency) the eigenvalue equation (39) reduces to (λ+γ)g(λ)=0, so the eigenvalues are λ=-γ and λ=2nπi/τON∗+τOFF∗ for n∈Z (the roots of g(λ)=0). It can also be seen from (39) that for all Z there is an eigenvalue λ=0. It corresponds to conservation of the number of TCLs. This eigenvalue 0 is removed by the normalisation condition (25). The real and imaginary parts of λ that solve (39) can be solved for numerically, using, for example, [54]. Figure 1 shows numerical solutions for the first five eigenvalues above (or on) the real axis for the parameter values given in Table 2 in Section 2.3 and allowing Z to vary from its value Z0 derived from the table, by Z=hZ0. There is an infinite sequence of eigenvalues going upwards and their reflections in the real axis. Increasing Z from zero by powers of 10 is seen to decrease the real part of the eigenvalues from zero and therefore the system is stable to small perturbations.Figure 1
Numerical solutions for the first five eigenvalues above the real axis (there is an infinite sequence going further along the imaginary axis, and they are reflected in the real axis). We use multiplierh to increase Z, and the real part of each eigenvalue we have followed decreases from 0 as Z increases from 0.This is a surprising result because intuitively identical TCLs are vulnerable to synchronisation which would cause instabilities on the system, which is the general view in the literature as discussed previously. The result is not due to the damping constantγ, because we chose γ=0 so as not to mask the effect of the TCLs. What the analysis does not tell us is how small any perturbations would have to be for a population of TCLs to have a stabilising effect on grid frequency. It might be that a larger perturbation than valid for linearisation leads to instability. In Section 4 we study the effects of different sized perturbations using simulations and indeed find growth of synchronisation. In the next section we consider the behaviour of a population of TCLs under the opposite type of perturbation; namely, all TCLs synchronised into one or two groups.


## 3.2. Synchronised Groups of TCLs

In the previous section we studied the stability of a uniformly distributed (continuum) population of TCLs at the 50 Hz equilibrium and found it to be stable almost everywhere in parameter space. In this section we consider the opposite extreme of possible TCL distributions, the Dirac delta distribution. That is to say, we explore the behaviour of a fully synchronised population of TCLs, all switching on and off at the same time, all with the same temperature, and (again) identical parameters. This is equivalent to a single TCL with the power consumption of the whole population.


### 3.2.1. Mapping the Switch Times

We begin by constructing a map from one (whole population) switch on event to the next. We show that under certain conditions such a mapping is a contraction. Let the subscriptn denote the nth switch on and nth switch off event. Without loss of generality, suppose that after our initial start time t0 the next switch event is the population switching on. This implies that, for all n∈N, tnOFF>tnON. Figure 2 illustrates the notation. Hence the amount of time the population spends switched on following the nth switch on event is given by(40a)tnOFF-tnON=1αlog⁡T+0-β+fnON-TONT-0-β-fnOFF-TON,where fnON,fnOFF are the frequencies at the nth switch on and off times. The amount of time spent switched off following the nth switch off is given by(40b)tn+1ON-tnOFF=1αlog⁡TOFF-T-0+β-fnOFFTOFF-T+0+β+fn+1ON.Assuming, as for the numerical analysis in Section 3.1, that the system has no damping, we set γ=0 in (12) for f˙(t). In a synchronised population, at time t all TCLs are either on (ρ(t)=1) or off (ρ(t)=0). Then we can define constants cON,cOFF>0 such that(41)f˙=-cON≔cPcρ∗-1whenthepopulationison+cOFF≔cPcρ∗whenthepopulationisoff.Hence the values of f at the switch off and on times are given by the piecewise-linear functions(42a)fnOFF=fnON-cONtnOFF-tnON(42b)fn+1ON=fnOFF+cOFFtn+1ON-tnOFF.After substituting for the switching times using (40a) and (40b) and rearranging, these become(43a)fnOFF-cONαlog⁡T-0-β-fnOFF-TON=fnON-cONαlog⁡T+0-β+fnON-TON(43b)fn+1ON+cOFFαlog⁡TOFF-T+0+β+fn+1ON=fnOFF+cOFFαlog⁡TOFF-T-0+β-fnOFF.Now since each side of (43a) and (43b) are functions of only one of the fn• variables, we can explicitly name them as such(44a)ϕON-fnOFF≔fnOFF-cONαlog⁡T-0-β-fnOFF-TON(44b)ϕON+fnON≔fnON-cONαlog⁡T+0-β+fnON-TON(44c)ϕOFF-fnOFF≔fnOFF+cOFFαlog⁡TOFF-T-0+β-fnOFF(44d)ϕOFF+fn+1ON≔fn+1ON+cOFFαlog⁡TOFF-T+0+β+fn+1ON.Each of the four ϕ functions is increasing and therefore invertible, and so we can write(45a)fnOFF=ϕON--1ϕON+fnON(45b)fn+1ON=ϕOFF+-1ϕOFF-fnOFFand therefore(45c)fn+1ON=ϕOFF+-1ϕOFF-ϕON--1ϕON+fnON,which is a mapping from the frequency at one switch on event to the frequency at the next. The mapping is a contraction iff (46)ϕOFF+-1ϕOFF-ϕON--1ϕON+′<1(47)iffϕOFF-′ϕOFF+′ϕON+′ϕON-′<1evaluatedattheappropriateplaces.Note that(48)ϕOFF-′ϕOFF+′=1+β-cOFF/αTOFF-Tn+1+β+cOFF/αTOFF-Tn+1+<1iffβ+β->TOFF-Tn+1+TOFF-Tn-.Similarly(49)ϕON+′ϕON-′=1+β+cON/αTn+-TON1+β-cON/αTn--TON<1iffβ+β-<Tn+-TONTn--TON.Therefore a sufficient condition for the mapping to be a contraction is that(50)β+β-∈TOFF-Tn+1+TOFF-Tn-,Tn+-TONTn--TONwhich is a nonempty interval (containing {1}), so long as TON<Tn-<Tn+<TOFF and Tn-<Tn+1+ for all n. It is worth recalling our earlier condition on the values of β± (7) which also imposed that β+/β- belong to an open interval containing {1}.Figure 2
Illustration of thenth and n+1th switching events of the fully synchronised population and the frequency-sensitive temperature set points T-(f(t)) and T+(f(t)).


### 3.2.2. Solving for the Periodic Solution

The contraction property of the mappingfnON↦fn+1ON (45c) under the above conditions implies that there is an attracting fixed point so long as TON<Tn-<Tn+<TOFF, and hence a periodic solution for the synchronised population. We now seek to solve for this periodic solution. Denote by lON and lOFF the amount of time spent on and off during one (periodic) cycle, respectively. Since power consumption for the population is constant during each on/off phase, the frequency moves linearly between upper and lower values which we denote by f+ and f-. Therefore the temperature of the population will cycle between upper and lower set points, given by T+0-β+f+ and T-0-β-f-, respectively. Equations (42a) and (42b) show us that(51a)f-=f+-cONlON(51b)f+=f-+cOFFlOFF.The temperature evolution equations (3a) and (3b) allow us to express the switch on and switch off temperatures as follows:(52a)T+0-β+f+=T-0-β-f--TOFFe-αlOFF+TOFF(52b)T-0-β-f-=T+0-β+f+-TONe-αlON+TONwhich after substituting for f- using (51a) and rearranging become(53a)f+β-e-αlOFF-β+=T-0-TOFF+β-cONlONe-αlOFF+TOFF-T+0,(53b)f+β+e-αlON-β-=T+0-TONe-αlON+TON-T-0-β-cONlON.Now we have two equations in terms of f+, lON and lOFF, which we can combine into one equation and eliminate f+,(54)β+e-αlON-β-T-0-TOFF+β-cONlONe-αlOFF+TOFF-T+0=β-e-αlOFF-β+T+0-TONe-αlON+TON-T-0-β-cONlON.We can also express lOFF in terms of lON by summing (51a) and (51b) to give(55)cONlON=cOFFlOFFor,equivalently,1-ρ∗lON-ρ∗lOFF=0and so (53b) and (55) form a pair of coupled equations for lON and lOFF, which can be solved numerically. Figure 3 shows one temperature cycle for the single group under different choices for ρ∗. Denote by ρ0 the value of ρ when f=0. As ρ∗ gets further away from ρ0 the solutions drift further from the uncoupled temperature range (2–7°C). The cycle lengths are symmetric about ρ∗=1/2 but the TCLs consume more power per cycle as ρ∗ increases.Figure 3
One cycle of the single group solution for different values ofρ∗ when ρ0≈0.3355. They include values that lead to unrealistic results for real fridges but are there to illustrate the effect.To begin some analysis of the stability of this fully synchronised solution we address the question: given a population split into two synchronised groups, will the groups merge into one fully synchronised population, or will they remain distinct forever?


### 3.2.3. Two-Group Dynamics

Suppose we have a population of frequency-sensitive TCLs that are split into two synchronised groups. We would like to understand the dynamics of the switch times, and we ask whether, given sufficient time, the groups will merge, or whether they will remain distinct, possibly settling down to separated periodic solutions. In particular, we consider the initial difference between the switch on timesΔt to be very small and the switch on temperatures very close to the single group periodic solution from the previous subsection.LetΓ denote the single group periodic solution, which cycles periodically through temperature space with temperature TΓ(t). As before, we denote the switched on duration in this solution by lON and the switched off duration by lOFF. Suppose that the population is split into two groups A and B, such that proportion σ belongs to group A, and proportion 1-σ belongs to group B. Suppose also that group B switches on at time t=0, followed soon after by group A switching on, at time t1>0. Then after a time period of length similar to lON group B switches off, which is again followed shortly after by group A switching off. After a time period similar to lOFF each of the groups then switches back on. We shall assume that the switching order does not change, since if they do swap, we need only repeat this process with σ replaced by 1-σ. Simulations show that the switching order will not continue to change indefinitely.We would like to compare the temperature cycles of these two groups with the single group periodic solutionΓ. Without loss of generality suppose that group B initially switches on at the same time as a fully synchronised population solution. We compare the cycling of the groups A and B using the following measures, along with all those shown in Figure 4. Let ΔT≔TB(0)-TA(0) and ΔT′≔TB(t4)-TA(t4), the temperature difference when B switches on the first and second times, respectively. In addition, let Δt≔t1-0=t1=ϵ1 and Δt′≔t5-t4, the time difference between the two groups switching on the first and second times, respectively. Further notation is shown in Figure 4.Figure 4
Linearisation about the single group solution. (a) Single group solutionTΓ(t). (b) Temperature cycling of groups A and B close to the single group solution.In order to calculateΔT′ and Δt′ we need to calculate the switch times and temperatures of the two groups at each switch event leading up to t5. Solving for the switch times and temperatures when there are two groups is a little more complicated than for the fully synchronised case. It requires solving the temperature set point equations using the system conditions at the previous switch and the equation for f˙ which now takes one of four values depending on which combination of groups is switched on (both, neither, A only, or B only). We begin by making the simplifying assumption β-=β+≔β. Now since group A is switching on at time t1 and group B switched off at time 0,(56a)TAt1=T+0-βft1ft1=f0-cPc1-σ-ρ∗t1f0=1βt+0-TB0∴TAt1=TB0+βcPc1-σ-ρ∗t1.In addition, by the temperature evolution equations,(56b)TAt1=TA0-TOFFe-αt1+TOFF.Equating (56a) and (56b) and introducing our new notation give(57)βcPc1-σ-ρ∗Δt+ΔT=TA0-TOFFe-αΔt-1.If we write TA(0)=TΓ(0)+δTA(0) and take δTA(0) and Δt small, then (58)ΔT=TΓ0+δTA0-TOFFe-αΔt-1-βcPc1-σ-ρ∗Δtand linearising in Δt gives(59)ΔT≈ξΔt,where(60)ξ≔αTOFF-TΓ0-βcPc1-σ-ρ∗.More generally, at each switch event we have the temperature evolution equations that describe the temperature of each group as a function of their temperature at the previous switch (such as (56b)) and an additional equation for the temperature of the switching group, using the temperature set point equations (such as (56a)). Writing ΔT′=TB(t4)-TA(t4) and linearising about the single group solution, we find in Appendix B that ΔT′=λΔT where(61)λ≔1-αTOFF-TONαTOFF-TΓ0-βcPc1-σ-ρ∗1-αTOFF-TONαTΓlON-TON+βcPcσ-ρ∗e-αL.So [-1,+1] is a left eigenvector of the linearised map in the space of δTAδTB with eigenvalue λ. We can plot λ against σ for various ρ∗ to see whether |λ|<1 (in which case the two groups merge into one) or whether |λ|>1 (they move apart). The results are shown in Figure 5. We find that the second eigenvalue is within the interval (-1,+1) for any σ and ρ∗ for our parameter values (taken from Table 2 with the exception of β=0.1 which has been reduced to limit the rate of change of the frequency and ensure model validity), see Appendix C for details, and therefore the stability is governed by λ.Figure 5
Solutions forλ (see (61)) (solid lines) for different values of ρ∗. Dashed lines show reflection in σ=1/2 to show the effect of reversing the switching order of the groups. Blue: ρ∗=0.1, yellow: ρ∗=0.2, green: ρ∗=0.3, and red: ρ∗=0.4. Black line shows the boundary of stability (stable below, unstable above). The results are identical when ρ∗ is replaced by 1-ρ∗. (b) shows an enlargement centred at σ=1/2, showing that either switching order of the groups leads to λ2>1 on a small interval of σ. In this case the groups never merge and in all other cases they will.
(a)
Full interval
(b)
EnlargementBy solving for the dividing caseλ=1 we can create a bifurcation diagram in terms of the parameters σ and ρ∗ to show where the single group solution is attracting and repelling. Figure 6 sketches the solution, along with the solution for the case when the switching order of the groups is reversed, found by replacing σ with 1-σ. If in either case (group A switching first or group B switching first) the solution is attracting, then the two groups will merge together into the one group solution. However, if both cases have unstable dynamics then the solutions will never merge. Simulations show that in this parameter region the two groups will settle down to a fixed phase distance apart. If the solution is attracting for one switching order and repelling for the other, we find that the typical behaviour is for a small separation in the unstable direction to grow until the phase difference becomes almost a whole cycle, when they merge. Figure 7 illustrates how the cycles of the two groups can change over time relative to one another, depending on which of the three regions in the bifurcation diagram their parameters belong to.Figure 6
Bifurcation diagram for the stability of the single group solution to splitting in two. Stable if the parameters lie in the yellow or blue regions (the groups will ultimately merge); unstable in the green parameter region (the groups will never merge). Boundary lines are solutions to (61) as a function of σ (or 1-σ to capture switching order reversal) and ρ∗.Figure 7
Illustration of the three types of cycling behaviour of two groups relative to one another, based on simulations. Arrows indicate the occurrence of many cycles and the central illustrations are snapshots of the cycling behaviour between the start and the final behaviour. Synchronisation occurs in cases (a) and (b), while in case (c) each group tends towards a fixed phase difference apart.What these results show is that when a population is split into two groups, if they are sufficiently similar in size then they will remain apart, effectively trying to counteract one another and balance the frequency fluctuations. Conversely, if one of the groups is significantly larger (“significantly” here depends on the size ofρ∗, and may be very small if ρ∗≈0.5) than the other then it will have too strong an effect on the frequency and “pull” on the smaller group’s cycle. The closer the proportion switched on in equilibrium is to the proportion switched off (i.e., the closer it is to 0.5), the more similar the groups have to be in size to remain distinct.With more than two groups of TCLs the modelling becomes far more complicated, since there are now far more possibilities to be considered for the switching order of the groups. Simulations have shown that for three groups it is possible for all three cycles to settle down to a fixed, separated pattern. This occurs if the groups are very similar in size, just as in the two-group case. Once one group is too large (or too small), the groups collapse into two, before synchronising completely. Taking the number of groups to infinity is equivalent to modelling a continuum of TCLs as in Section3.1. Taking all groups of equal size and uniformly distributed in phase θ is equivalent to the continuum population equilibrium studied earlier. From above we found analytically that small perturbations to the population distribution should relax back to the uniform distribution, that is, that the equilibrium was stable. Now we find that if the population is discretised into 2 (and hypothetically N) groups then so long as they are of close to equal sizes, they will attempt to settle the frequency back to its nominal value by “spreading out” their cycles.In reality we will never have a continuum of TCLs, and they may exhibit nonlinear dynamics not captured by our analysis. This motivates our use of simulations to gain further insights into how a large population of TCLs would behave according to our switching rules and how the grid frequency would be affected.


## 3.2.1. Mapping the Switch Times

We begin by constructing a map from one (whole population) switch on event to the next. We show that under certain conditions such a mapping is a contraction. Let the subscriptn denote the nth switch on and nth switch off event. Without loss of generality, suppose that after our initial start time t0 the next switch event is the population switching on. This implies that, for all n∈N, tnOFF>tnON. Figure 2 illustrates the notation. Hence the amount of time the population spends switched on following the nth switch on event is given by(40a)tnOFF-tnON=1αlog⁡T+0-β+fnON-TONT-0-β-fnOFF-TON,where fnON,fnOFF are the frequencies at the nth switch on and off times. The amount of time spent switched off following the nth switch off is given by(40b)tn+1ON-tnOFF=1αlog⁡TOFF-T-0+β-fnOFFTOFF-T+0+β+fn+1ON.Assuming, as for the numerical analysis in Section 3.1, that the system has no damping, we set γ=0 in (12) for f˙(t). In a synchronised population, at time t all TCLs are either on (ρ(t)=1) or off (ρ(t)=0). Then we can define constants cON,cOFF>0 such that(41)f˙=-cON≔cPcρ∗-1whenthepopulationison+cOFF≔cPcρ∗whenthepopulationisoff.Hence the values of f at the switch off and on times are given by the piecewise-linear functions(42a)fnOFF=fnON-cONtnOFF-tnON(42b)fn+1ON=fnOFF+cOFFtn+1ON-tnOFF.After substituting for the switching times using (40a) and (40b) and rearranging, these become(43a)fnOFF-cONαlog⁡T-0-β-fnOFF-TON=fnON-cONαlog⁡T+0-β+fnON-TON(43b)fn+1ON+cOFFαlog⁡TOFF-T+0+β+fn+1ON=fnOFF+cOFFαlog⁡TOFF-T-0+β-fnOFF.Now since each side of (43a) and (43b) are functions of only one of the fn• variables, we can explicitly name them as such(44a)ϕON-fnOFF≔fnOFF-cONαlog⁡T-0-β-fnOFF-TON(44b)ϕON+fnON≔fnON-cONαlog⁡T+0-β+fnON-TON(44c)ϕOFF-fnOFF≔fnOFF+cOFFαlog⁡TOFF-T-0+β-fnOFF(44d)ϕOFF+fn+1ON≔fn+1ON+cOFFαlog⁡TOFF-T+0+β+fn+1ON.Each of the four ϕ functions is increasing and therefore invertible, and so we can write(45a)fnOFF=ϕON--1ϕON+fnON(45b)fn+1ON=ϕOFF+-1ϕOFF-fnOFFand therefore(45c)fn+1ON=ϕOFF+-1ϕOFF-ϕON--1ϕON+fnON,which is a mapping from the frequency at one switch on event to the frequency at the next. The mapping is a contraction iff (46)ϕOFF+-1ϕOFF-ϕON--1ϕON+′<1(47)iffϕOFF-′ϕOFF+′ϕON+′ϕON-′<1evaluatedattheappropriateplaces.Note that(48)ϕOFF-′ϕOFF+′=1+β-cOFF/αTOFF-Tn+1+β+cOFF/αTOFF-Tn+1+<1iffβ+β->TOFF-Tn+1+TOFF-Tn-.Similarly(49)ϕON+′ϕON-′=1+β+cON/αTn+-TON1+β-cON/αTn--TON<1iffβ+β-<Tn+-TONTn--TON.Therefore a sufficient condition for the mapping to be a contraction is that(50)β+β-∈TOFF-Tn+1+TOFF-Tn-,Tn+-TONTn--TONwhich is a nonempty interval (containing {1}), so long as TON<Tn-<Tn+<TOFF and Tn-<Tn+1+ for all n. It is worth recalling our earlier condition on the values of β± (7) which also imposed that β+/β- belong to an open interval containing {1}.Figure 2
Illustration of thenth and n+1th switching events of the fully synchronised population and the frequency-sensitive temperature set points T-(f(t)) and T+(f(t)).


## 3.2.2. Solving for the Periodic Solution

The contraction property of the mappingfnON↦fn+1ON (45c) under the above conditions implies that there is an attracting fixed point so long as TON<Tn-<Tn+<TOFF, and hence a periodic solution for the synchronised population. We now seek to solve for this periodic solution. Denote by lON and lOFF the amount of time spent on and off during one (periodic) cycle, respectively. Since power consumption for the population is constant during each on/off phase, the frequency moves linearly between upper and lower values which we denote by f+ and f-. Therefore the temperature of the population will cycle between upper and lower set points, given by T+0-β+f+ and T-0-β-f-, respectively. Equations (42a) and (42b) show us that(51a)f-=f+-cONlON(51b)f+=f-+cOFFlOFF.The temperature evolution equations (3a) and (3b) allow us to express the switch on and switch off temperatures as follows:(52a)T+0-β+f+=T-0-β-f--TOFFe-αlOFF+TOFF(52b)T-0-β-f-=T+0-β+f+-TONe-αlON+TONwhich after substituting for f- using (51a) and rearranging become(53a)f+β-e-αlOFF-β+=T-0-TOFF+β-cONlONe-αlOFF+TOFF-T+0,(53b)f+β+e-αlON-β-=T+0-TONe-αlON+TON-T-0-β-cONlON.Now we have two equations in terms of f+, lON and lOFF, which we can combine into one equation and eliminate f+,(54)β+e-αlON-β-T-0-TOFF+β-cONlONe-αlOFF+TOFF-T+0=β-e-αlOFF-β+T+0-TONe-αlON+TON-T-0-β-cONlON.We can also express lOFF in terms of lON by summing (51a) and (51b) to give(55)cONlON=cOFFlOFFor,equivalently,1-ρ∗lON-ρ∗lOFF=0and so (53b) and (55) form a pair of coupled equations for lON and lOFF, which can be solved numerically. Figure 3 shows one temperature cycle for the single group under different choices for ρ∗. Denote by ρ0 the value of ρ when f=0. As ρ∗ gets further away from ρ0 the solutions drift further from the uncoupled temperature range (2–7°C). The cycle lengths are symmetric about ρ∗=1/2 but the TCLs consume more power per cycle as ρ∗ increases.Figure 3
One cycle of the single group solution for different values ofρ∗ when ρ0≈0.3355. They include values that lead to unrealistic results for real fridges but are there to illustrate the effect.To begin some analysis of the stability of this fully synchronised solution we address the question: given a population split into two synchronised groups, will the groups merge into one fully synchronised population, or will they remain distinct forever?


## 3.2.3. Two-Group Dynamics

Suppose we have a population of frequency-sensitive TCLs that are split into two synchronised groups. We would like to understand the dynamics of the switch times, and we ask whether, given sufficient time, the groups will merge, or whether they will remain distinct, possibly settling down to separated periodic solutions. In particular, we consider the initial difference between the switch on timesΔt to be very small and the switch on temperatures very close to the single group periodic solution from the previous subsection.LetΓ denote the single group periodic solution, which cycles periodically through temperature space with temperature TΓ(t). As before, we denote the switched on duration in this solution by lON and the switched off duration by lOFF. Suppose that the population is split into two groups A and B, such that proportion σ belongs to group A, and proportion 1-σ belongs to group B. Suppose also that group B switches on at time t=0, followed soon after by group A switching on, at time t1>0. Then after a time period of length similar to lON group B switches off, which is again followed shortly after by group A switching off. After a time period similar to lOFF each of the groups then switches back on. We shall assume that the switching order does not change, since if they do swap, we need only repeat this process with σ replaced by 1-σ. Simulations show that the switching order will not continue to change indefinitely.We would like to compare the temperature cycles of these two groups with the single group periodic solutionΓ. Without loss of generality suppose that group B initially switches on at the same time as a fully synchronised population solution. We compare the cycling of the groups A and B using the following measures, along with all those shown in Figure 4. Let ΔT≔TB(0)-TA(0) and ΔT′≔TB(t4)-TA(t4), the temperature difference when B switches on the first and second times, respectively. In addition, let Δt≔t1-0=t1=ϵ1 and Δt′≔t5-t4, the time difference between the two groups switching on the first and second times, respectively. Further notation is shown in Figure 4.Figure 4
Linearisation about the single group solution. (a) Single group solutionTΓ(t). (b) Temperature cycling of groups A and B close to the single group solution.In order to calculateΔT′ and Δt′ we need to calculate the switch times and temperatures of the two groups at each switch event leading up to t5. Solving for the switch times and temperatures when there are two groups is a little more complicated than for the fully synchronised case. It requires solving the temperature set point equations using the system conditions at the previous switch and the equation for f˙ which now takes one of four values depending on which combination of groups is switched on (both, neither, A only, or B only). We begin by making the simplifying assumption β-=β+≔β. Now since group A is switching on at time t1 and group B switched off at time 0,(56a)TAt1=T+0-βft1ft1=f0-cPc1-σ-ρ∗t1f0=1βt+0-TB0∴TAt1=TB0+βcPc1-σ-ρ∗t1.In addition, by the temperature evolution equations,(56b)TAt1=TA0-TOFFe-αt1+TOFF.Equating (56a) and (56b) and introducing our new notation give(57)βcPc1-σ-ρ∗Δt+ΔT=TA0-TOFFe-αΔt-1.If we write TA(0)=TΓ(0)+δTA(0) and take δTA(0) and Δt small, then (58)ΔT=TΓ0+δTA0-TOFFe-αΔt-1-βcPc1-σ-ρ∗Δtand linearising in Δt gives(59)ΔT≈ξΔt,where(60)ξ≔αTOFF-TΓ0-βcPc1-σ-ρ∗.More generally, at each switch event we have the temperature evolution equations that describe the temperature of each group as a function of their temperature at the previous switch (such as (56b)) and an additional equation for the temperature of the switching group, using the temperature set point equations (such as (56a)). Writing ΔT′=TB(t4)-TA(t4) and linearising about the single group solution, we find in Appendix B that ΔT′=λΔT where(61)λ≔1-αTOFF-TONαTOFF-TΓ0-βcPc1-σ-ρ∗1-αTOFF-TONαTΓlON-TON+βcPcσ-ρ∗e-αL.So [-1,+1] is a left eigenvector of the linearised map in the space of δTAδTB with eigenvalue λ. We can plot λ against σ for various ρ∗ to see whether |λ|<1 (in which case the two groups merge into one) or whether |λ|>1 (they move apart). The results are shown in Figure 5. We find that the second eigenvalue is within the interval (-1,+1) for any σ and ρ∗ for our parameter values (taken from Table 2 with the exception of β=0.1 which has been reduced to limit the rate of change of the frequency and ensure model validity), see Appendix C for details, and therefore the stability is governed by λ.Figure 5
Solutions forλ (see (61)) (solid lines) for different values of ρ∗. Dashed lines show reflection in σ=1/2 to show the effect of reversing the switching order of the groups. Blue: ρ∗=0.1, yellow: ρ∗=0.2, green: ρ∗=0.3, and red: ρ∗=0.4. Black line shows the boundary of stability (stable below, unstable above). The results are identical when ρ∗ is replaced by 1-ρ∗. (b) shows an enlargement centred at σ=1/2, showing that either switching order of the groups leads to λ2>1 on a small interval of σ. In this case the groups never merge and in all other cases they will.
(a)
Full interval
(b)
EnlargementBy solving for the dividing caseλ=1 we can create a bifurcation diagram in terms of the parameters σ and ρ∗ to show where the single group solution is attracting and repelling. Figure 6 sketches the solution, along with the solution for the case when the switching order of the groups is reversed, found by replacing σ with 1-σ. If in either case (group A switching first or group B switching first) the solution is attracting, then the two groups will merge together into the one group solution. However, if both cases have unstable dynamics then the solutions will never merge. Simulations show that in this parameter region the two groups will settle down to a fixed phase distance apart. If the solution is attracting for one switching order and repelling for the other, we find that the typical behaviour is for a small separation in the unstable direction to grow until the phase difference becomes almost a whole cycle, when they merge. Figure 7 illustrates how the cycles of the two groups can change over time relative to one another, depending on which of the three regions in the bifurcation diagram their parameters belong to.Figure 6
Bifurcation diagram for the stability of the single group solution to splitting in two. Stable if the parameters lie in the yellow or blue regions (the groups will ultimately merge); unstable in the green parameter region (the groups will never merge). Boundary lines are solutions to (61) as a function of σ (or 1-σ to capture switching order reversal) and ρ∗.Figure 7
Illustration of the three types of cycling behaviour of two groups relative to one another, based on simulations. Arrows indicate the occurrence of many cycles and the central illustrations are snapshots of the cycling behaviour between the start and the final behaviour. Synchronisation occurs in cases (a) and (b), while in case (c) each group tends towards a fixed phase difference apart.What these results show is that when a population is split into two groups, if they are sufficiently similar in size then they will remain apart, effectively trying to counteract one another and balance the frequency fluctuations. Conversely, if one of the groups is significantly larger (“significantly” here depends on the size ofρ∗, and may be very small if ρ∗≈0.5) than the other then it will have too strong an effect on the frequency and “pull” on the smaller group’s cycle. The closer the proportion switched on in equilibrium is to the proportion switched off (i.e., the closer it is to 0.5), the more similar the groups have to be in size to remain distinct.With more than two groups of TCLs the modelling becomes far more complicated, since there are now far more possibilities to be considered for the switching order of the groups. Simulations have shown that for three groups it is possible for all three cycles to settle down to a fixed, separated pattern. This occurs if the groups are very similar in size, just as in the two-group case. Once one group is too large (or too small), the groups collapse into two, before synchronising completely. Taking the number of groups to infinity is equivalent to modelling a continuum of TCLs as in Section3.1. Taking all groups of equal size and uniformly distributed in phase θ is equivalent to the continuum population equilibrium studied earlier. From above we found analytically that small perturbations to the population distribution should relax back to the uniform distribution, that is, that the equilibrium was stable. Now we find that if the population is discretised into 2 (and hypothetically N) groups then so long as they are of close to equal sizes, they will attempt to settle the frequency back to its nominal value by “spreading out” their cycles.In reality we will never have a continuum of TCLs, and they may exhibit nonlinear dynamics not captured by our analysis. This motivates our use of simulations to gain further insights into how a large population of TCLs would behave according to our switching rules and how the grid frequency would be affected.


## 4. Simulations

### 4.1. Perturbations of a Uniform Distribution of TCLs

In Section3.1 we analysed the stability of a large population of TCLs uniformly distributed in each part of the on/off cycle. In this section we simulate a large population of fridges with initial conditions close to the equilibrium distribution (the uniform distribution) and compare the results with our analytical work.The model is the same as presented in Section2, and unless stated otherwise, the parameter values are as in Table 2. In Section 3.1 we modelled our population as a continuum. For our simulations we split the fridge population into 104 “agents” (groups of fridges) that are each represented by a temperature and state, and who operate according to the switching rules and temperature progression equations in Section 2.1. These 104 agents are representative of the million fridges we assume are participating in our DSR scheme (i.e., operating in frequency-sensitive mode), since one million (or more) individuals would require very large amounts of computing time and memory. The power consumption of each agent is taken to be the total possible population consumption Pc divided by the number of agents, 104. Each time step is taken to be 1 s, and at each time step each agent updates its temperature and based on the frequency at the previous time step may switch on or off. The exact switch time is approximated using linear interpolation between the current and previous time step, and the new temperature is adjusted accordingly.To perturb the TCL distributionu(θ) we can alter the number of TCLs switched on or off from the equilibrium proportions ρ∗ and (1-ρ∗), respectively, and we can perturb the uniform distributions within each on/off half of the θ interval. We choose to perturb the distributions by the addition of a sine wave to u∗, and we refer to the normalised wave peak amplitude Δu (normalised by dividing by u∗). This normalisation means that when we plot u(θ,0)/u∗, the zero perturbation case is 1 for all θ both on and off and the results are more clear. Table 3 shows eight combinations of choices for these perturbation parameters. All other parameters are as stated in Table 2.Table 3
Parameter values for plots in Figures8, 9, and 10.
Plot number
ρ
o
n
(
0
)
Δ
u
a(i)
ρ
∗
0
a(ii)
ρ
∗
0.1
a(iii)
ρ
∗
0.25
a(iv)
ρ
∗
0.5
b(i)
1.5
ρ
∗
0
b(ii)
1.5
ρ
∗
0.1
b(iii)
1.5
ρ
∗
0.25
b(iv)
1.5
ρ
∗
0.5Figure8 shows the effects of these perturbations on the initial conditions in each case, plotting u(θ,0)/u∗ against θ. Figure 9 shows the final fridge distributions after ten days. The unperturbed case a(i) has remained uniform, while the peaks of the perturbation cases have all grown by varying amounts. In cases a(ii)–a(iv) (no perturbation to the proportion switched on) the final distributions exhibit increasing levels of synchronisation, but the clustering is far less than in cases b(i)–b(iv) which see the population synchronised into seven or fewer groups. The effects of this synchronisation on the electricity grid frequency can be seen in Figure 10.Figure 8
Initial fridge distributions inθ-space, with labels matching those in Table 3. Pink indicates switched off fridges, and blue indicates switched on. Distributions scaled by 1/u∗ and histograms formed of 100 bins. (a) has no perturbation to the proportion of fridges switched on, and (b) has increasing perturbation (going downwards) to the number of fridges switched on. All involve sinusoidal distribution perturbations.Figure 9
Final fridge distributions after 10 days inθ-space, with properties as given in Table 3 and initial distributions as shown in Figure 8. Perturbations have grown (except for the zero perturbation case a(i)).Figure 10
Electricity grid frequency over 10 days (values plotted once per 5 minutes), with fridge distributions as described in Table3 and Figures 8 and 9. The perturbed systems (all but a(i)) see an initial reduction in oscillation amplitude followed by oscillation growth.Interestingly, in each case with perturbations, the frequency oscillations initially die down to close to 50 Hz. This means that to begin with the fridges are controlling the frequency oscillations caused by their initial condition perturbations. This aligns with our analysis from Section3.1, in which we found that the uniform distribution of a continuum population is stable to small perturbations. What that analysis was unable to capture was the long-term effects of frequency sensitivity. In each case the frequency oscillations grow after less than a day, becoming very large in several cases. Before the large spikes in b(iii) we see that the frequency oscillations shrink down. This shows the inherently volatile nature of the system and potentially explains why the oscillations in b(iv) are ultimately less severe. It could be that these lower oscillations will shortly become much larger. In either case, the size of most of the final oscillations would be too large for the system to cope without frequency response from other providers.These simulations reveal that while a homogeneous population of TCLs will act to dampen system perturbations, their behaviour to support the electricity grid will, given sufficient time, lead to further oscillations. The larger the perturbations are, the sooner these detrimental effects will occur.


### 4.2. Simulating TCLs on the GB Electricity Grid

Our model and simulations have thus far reduced the complexity of the problem by assuming that, apart from the TCL population and the grid frequency, all other network conditions remain constant. This was necessary for our model to be tractable and to ensure that any results from the simulations were attributable to the frequency-sensitive TCL population. An important next step is to consider the TCL population in the context of a real system. In collaboration with the GB System Operator National Grid, we are able to model the GB system with real data from 36 separate 10-day periods during 2015-2016 and simulate what would have happened if a frequency-sensitive fridge population had been active. We consider how the distribution of TCLs changes over this period, and the reduction in the amount of response that other providers needed to supply because of the contribution from the fridges.


#### 4.2.1. Methods

We simulate a population of TCLs (specifically fridges) that respond to the grid frequency according to the rules in Section2.1. We use various historic data from National Grid to model real system conditions and simulate the effects of a frequency-sensitive fridge population acting on the GB system. By considering the population in the context of real data including response provision from other sources such as power generators, we are able to get a better understanding of the potential impact of the fridges compared to, say, modelling them in isolation responding to a one-off frequency event.Figure11 gives an overview of the simulation process. Rhombi indicate inputs and outputs; rectangles indicate methods used in the simulation. Methods are applied working downwards, except for the dashed arrows which create an iterative loop.Figure 11
Simulation methodology diagram. Rhombi indicate input or calculated data/simulated data, and rectangles indicate methods/calculations. Events occur from top to bottom with the exception of the dashed arrows which form the iterative loop.


#### 4.2.2. Inputs

As shown in Figure11, there are four types of data input, in addition to the fridge population initial conditions. We use 36 consecutive ten-day data samples from the period July 2015–June 2016.Kinetic energy data is an estimate for the total kinetic energy in MVAs (megavolt-ampere seconds) [55]. Values are calculated by summing the kinetic energy of all running synchronised generators (a generator-specific constant provided to the System Operator by each power generator) with an estimate of kinetic energy from demand. The kinetic energy data provided (confidentially from National Grid) is per settlement period (settlement periods split the day into 48 half hour units starting on the hour and half hour) and repeats each value for the full 30 minutes (rather than interpolating). Typical kinetic energy values are in the range 20000–40000 MVAs.Demand data consists of per-second metered demand from National Grid. This is a sum of the power leaving the electricity transmission system, including any power exports through the interconnectors. Half-hourly demand data is accessible via National Grid’s “Data Explorer” [56].Historic frequency data consists of per-second system frequency data in hertz. Frequency measurements are taken in multiple locations to ensure reliable data availability in the event of any metering faults. The frequency data provided by National Grid has undergone a cleaning process that takes advantage of the multiple readings. It is available via National Grid’s “Enhanced Frequency Response” [57].Response holdings are the amount of frequency response delivery in MW (as a function of grid frequency) that National Grid expects each second. Response holdings are positive (or negative) for “low (high) frequency response delivery” when the frequency is below (above) 50 Hz, respectively. For each time step (1 second), 9 different values for response holding are listed. These take the form of primary, secondary, and high response.Primary response values are given for trigger points at 49.9 Hz, 49.5 Hz, and 49.2 Hz. This means that at these frequencies the power response provided through various types of primary response service are the historic response holding values given, subject to a 1 second reaction delay. We assume that the response increases linearly from 0 between 49.985 Hz and 49.8 Hz and likewise linearly between all other frequency trigger values. Below 49.2 Hz the response is assumed to be the constant 49.2 Hz response value. The starting frequency trigger value of 49.985 Hz is used to take into account the Grid Code deadband of (50±0.015) Hz, within which response is not required. Secondary response values are given for frequency trigger points 49.8 Hz and 49.5 Hz, and response is modelled in the same was as for primary response, only with an 11 s response delay. High response values have trigger points 50.2 Hz and 50.5 Hz. Just as for primary response, the time lag is 1 s and again, response is modelled as linear interpolation through these points, starting at the edge of the deadband at 50.015 Hz and remaining constant beyond 50.5 Hz. Figure 12 illustrates an example of how response holding data (Table 4) are interpreted in the model. Values given are indicative only of possible values.Table 4
Illustrative historic response holding data behind Figure12.
Primary
Secondary
High
Frequency trigger (Hz)
49.2
49.5
49.8
49.5
49.8
50.2
50.5
Response (MW)
850
800
430
950
500
−350
−680Figure 12
Representative historic response data with interpolation method for primary response (solid line below 50 Hz), secondary response (dashed line), and high response (solid line above 50 Hz). Zero response in the deadband (50±0.015) Hz.Fridge conditions are the initial on/off state and initial temperature of each fridge in the population. For the simulations presented here we take the zero perturbation case a(i) in Table 3 from the previous section.


#### 4.2.3. Calculating the Demand at 50 Hz

Deviations in grid frequency away from 50 Hz affect the total system demand. We make the assumption that demand increases linearly by approximately 2.5% of its value at 50 Hz for every 1 Hz increase in frequency above 50 Hz (and decreases by the same amount as frequency decreases below 50 Hz). In order to know the demand at the nominal frequency, “demand at 50 Hz,”Demω0(t), we need to calculate it from the (measured) demand data input, D(t).(62)Dt=Demω0t1+0.025ft-50(63)Demω0t=Dt1+0.025ft-50.


#### 4.2.4. Calculating the Underlying Imbalance

In order to calculate the effects of the fridge population on the system frequency, we first need to calculate the underlying supply-demand imbalance (in MW) that caused the original system frequency deviations away from 50 Hz. At this point it is necessary to distinguish between two important, similar-sounding terms:underlying imbalance andtotal imbalance. By underlying imbalance, Imbunder(t), we mean the generation-demand imbalance that occurs independently of the system frequency. This may be due to, for example, fluctuations in wind or solar power generation or discrepancies between the total predicted system demand and the actual real-time demand. In contrast, total imbalance, Imbtot(f,t), includes both the underlying imbalance and, additionally, what we shall refer to asdynamic imbalance.There are two sources of dynamic imbalance:generator response (frequency response provided by power generators as the frequency changes) anddemand response (the automatic change in demand as frequency changes). Note that in this context “demand response” is completely different to demand-side response services, which, given their current low penetration of the response market, we exclude from our simulations. Generator response, Genresp, consists of the actual response delivered by generators, calculated as described above from the response holdings and the historic system frequency. Generator response is assumed to have a small time lag δt, which we take to be 1 second. In contrast, demand response, Demresp, is assumed to occur instantaneously and is defined as the measured system demand D(t) minus the demand at 50 Hz, Dω0(t) (see “calculating demand at 50 Hz”). Therefore by (62)(64)Demrespft,t=0.025Dω0tft-50.Both sources of dynamic imbalance will change when we introduce the population of responsive fridges (because of their impact on the frequency) and will therefore need to be recalculated.We use a linear approximation for the rate of change of frequency [51], which in our notation is given by(65)dfdt=50Imbtott2Ekt,where 50 is the nominal frequency 50 Hz and Ek=F0M/2 is total stored kinetic energy in MVAs. Since(66)Imbtott=Imbundert+Genrespft-δt,t-Demrespft,twe are able to find(67)Imbundert=Ekt25dfdt-Genrespft-δt,t+Demrespft,twhich for simulation time step size Δt gives(68)Imbundert=Ektft-ft-Δt25Δt-Genrespft-δt,t+Demrespft,t.We take Δt=1 s, so Δt=δt, the generator response time lag. Generator response is calculated using historic response holdings and the frequency t-δt seconds ago along with some constraints on the generator ramp rates.


#### 4.2.5. Iterative Loop

Once the underlying imbalance has been calculated for all time steps it can be used along with the response holdings and fridge conditions to begin a loop formed of three calculation steps that iterates over all time steps (see the “iterative loop” in Figure11). The steps are as follows:(1)
Calculate the frequency response delivery from the fridge population and from the dynamic response providers based on the previous frequency value (the first iteration takes the first historic frequency value, after which the “new frequency” values are used). For the fridge population this requires summing the switched on fridges multiplied by their individual power consumption and subtracting the power consumption of the population if the fridges were not frequency-sensitive. Response from the dynamic response providers is described above.(2)
Calculate the new frequencyf∗(t) using the equations from “calculating the underlying imbalance” and beginning with the approximation (69)f∗t=f∗t-Δt+Δtdf∗dtt=f∗t-Δt+Δt25Imbtot∗tEktand since(70)Imbtot∗t=Imbundert+Genrespf∗t-δt,t-0.025Demω0f∗t-50we get(71)f∗t=f∗t-Δt+25Δt/EktImbundert+Genrespf∗t-δt,t+1.25Demω01+0.625δt/EktDemω0.Note that we letf∗(0)=f(0), the original frequency value at time 0.(3)
Calculate the new fridge conditions by updating their temperature set points with the new frequency f∗ calculated in step 2, according to (2a) and (2b). Each fridge temperature is evolved one time step according to (3a) or (3b). If a switch on or off should have occurred during the time step then the exact time of switch is estimated and the temperature is recalculated from the switch time to the end of the time step using linear interpolation.


#### 4.2.6. Outputs

There are two key outputs for our analysis: firstly, the temperatures and states of each fridge over time and secondly, the frequency response supplied by all other providers on the grid. Since response can be positive or negative depending on the frequency, but both incur payment, we take the absolute value of the response at each time step. We take the cumulative sum of the difference between this response in the presence of TCLs and the original system response and call it “cumulative response savings,” which we measure in MWh. This allows us to find out how much benefit (or detriment) the fridges provided the system and how that changes over time as they respond to frequency perturbations.


#### 4.2.7. Results

We begin by comparing the results of two populations of fridges, one with total power consumption (if every fridge were switched on),Pc=70 MW, as in previous simulations, and the other with Pc=700 MW. The second case is the extreme with 10 million frequency-sensitive fridges, similar to Trovato et al. who modelled 11 million fridges in [33]. Rather than modelling all 1 million or 10 million fridges, we split the population into 104 groups, where the fridges within one group have the same temperature and state. We begin both simulations with the groups uniformly distributed in phase (the unperturbed case) just as in figures a(i) of the previous simulation section. We repeat these simulations on 36 10-day data samples from July 2015–June 2016.Figure13 shows cumulative response savings over time (introduced above) for the 36 data sample simulations in each case. For both values of Pc, in at least a third of the simulations, the fridge population ended up doing more harm than good (negative savings) due to synchronisation. Increasing participation tenfold increased the best results by about a factor of 7 but worsened the worst results by a factor of 15. When there were fewer participants (Pc=70 MW), the results were more erratic over time, which we attribute to there being less response on the system, and so a less smooth frequency trace to respond to. In light of these findings, we present the results from another set of 36 simulations over the year for Pc=700 MW, with the addition of a very small amount of diversity (less than 0.25% in relative terms) to the parameters.Figure 13
Cumulative response savings (MWh), the difference between other providers’ response with and without the frequency-sensitive fridge population (cumulatively) for the 36 data samples over one year for two different participation levels. We sum high frequency response and the absolute value of low frequency response, in MWh. Negative results indicate that the other providers had to compensate for detrimental fridge behaviour. In both cases the population is homogeneous.
(a)
P
c
=
7
0 MW
(b)
P
c
=
700 MWWe find that very small amounts of parameter diversification can eradicate the detrimental fridge behaviour in our simulations. A full presentation of all of our simulations can be found in Webborn Ph.D. thesis, to appear. As an example we takePc=700 and the other parameters to have the same mean as in all of our previous simulations but draw them from normal distributions with nonzero standard deviation. We choose TON~N(-26,0.020), TOFF~N(20,0.0133), T-0~N(2,0.005), (T+0-T-0)~N(5,0.003), and α~N(18.08×105,0.030×105). This results in T+0~N(7,0.005), the duty cycle approximately ~N(33.55%,0.027%), and the cycle period in minutes ~N(41.144,0.108). The results from these simulations are shown in Figure 14. We see that even with this very small amount of parameter diversity, all the simulations show a net and growing benefit to the system over the ten days. Since diversity naturally occurs in any real population, this offers reasonable evidence that TCLs could be a valuable resource for the grid, without the need for stochastic switching or regular communications from a central controller.Figure 14
Cumulative response savings (MWh) (see Figure13) for a heterogeneous fridge population with Pc=700 MW. The introduction of a very small amount of diversity has eradicated the detrimental behaviour in all cases and the fridge population provides a clear benefit to the system.


## 4.1. Perturbations of a Uniform Distribution of TCLs

In Section3.1 we analysed the stability of a large population of TCLs uniformly distributed in each part of the on/off cycle. In this section we simulate a large population of fridges with initial conditions close to the equilibrium distribution (the uniform distribution) and compare the results with our analytical work.The model is the same as presented in Section2, and unless stated otherwise, the parameter values are as in Table 2. In Section 3.1 we modelled our population as a continuum. For our simulations we split the fridge population into 104 “agents” (groups of fridges) that are each represented by a temperature and state, and who operate according to the switching rules and temperature progression equations in Section 2.1. These 104 agents are representative of the million fridges we assume are participating in our DSR scheme (i.e., operating in frequency-sensitive mode), since one million (or more) individuals would require very large amounts of computing time and memory. The power consumption of each agent is taken to be the total possible population consumption Pc divided by the number of agents, 104. Each time step is taken to be 1 s, and at each time step each agent updates its temperature and based on the frequency at the previous time step may switch on or off. The exact switch time is approximated using linear interpolation between the current and previous time step, and the new temperature is adjusted accordingly.To perturb the TCL distributionu(θ) we can alter the number of TCLs switched on or off from the equilibrium proportions ρ∗ and (1-ρ∗), respectively, and we can perturb the uniform distributions within each on/off half of the θ interval. We choose to perturb the distributions by the addition of a sine wave to u∗, and we refer to the normalised wave peak amplitude Δu (normalised by dividing by u∗). This normalisation means that when we plot u(θ,0)/u∗, the zero perturbation case is 1 for all θ both on and off and the results are more clear. Table 3 shows eight combinations of choices for these perturbation parameters. All other parameters are as stated in Table 2.Table 3
Parameter values for plots in Figures8, 9, and 10.
Plot number
ρ
o
n
(
0
)
Δ
u
a(i)
ρ
∗
0
a(ii)
ρ
∗
0.1
a(iii)
ρ
∗
0.25
a(iv)
ρ
∗
0.5
b(i)
1.5
ρ
∗
0
b(ii)
1.5
ρ
∗
0.1
b(iii)
1.5
ρ
∗
0.25
b(iv)
1.5
ρ
∗
0.5Figure8 shows the effects of these perturbations on the initial conditions in each case, plotting u(θ,0)/u∗ against θ. Figure 9 shows the final fridge distributions after ten days. The unperturbed case a(i) has remained uniform, while the peaks of the perturbation cases have all grown by varying amounts. In cases a(ii)–a(iv) (no perturbation to the proportion switched on) the final distributions exhibit increasing levels of synchronisation, but the clustering is far less than in cases b(i)–b(iv) which see the population synchronised into seven or fewer groups. The effects of this synchronisation on the electricity grid frequency can be seen in Figure 10.Figure 8
Initial fridge distributions inθ-space, with labels matching those in Table 3. Pink indicates switched off fridges, and blue indicates switched on. Distributions scaled by 1/u∗ and histograms formed of 100 bins. (a) has no perturbation to the proportion of fridges switched on, and (b) has increasing perturbation (going downwards) to the number of fridges switched on. All involve sinusoidal distribution perturbations.Figure 9
Final fridge distributions after 10 days inθ-space, with properties as given in Table 3 and initial distributions as shown in Figure 8. Perturbations have grown (except for the zero perturbation case a(i)).Figure 10
Electricity grid frequency over 10 days (values plotted once per 5 minutes), with fridge distributions as described in Table3 and Figures 8 and 9. The perturbed systems (all but a(i)) see an initial reduction in oscillation amplitude followed by oscillation growth.Interestingly, in each case with perturbations, the frequency oscillations initially die down to close to 50 Hz. This means that to begin with the fridges are controlling the frequency oscillations caused by their initial condition perturbations. This aligns with our analysis from Section3.1, in which we found that the uniform distribution of a continuum population is stable to small perturbations. What that analysis was unable to capture was the long-term effects of frequency sensitivity. In each case the frequency oscillations grow after less than a day, becoming very large in several cases. Before the large spikes in b(iii) we see that the frequency oscillations shrink down. This shows the inherently volatile nature of the system and potentially explains why the oscillations in b(iv) are ultimately less severe. It could be that these lower oscillations will shortly become much larger. In either case, the size of most of the final oscillations would be too large for the system to cope without frequency response from other providers.These simulations reveal that while a homogeneous population of TCLs will act to dampen system perturbations, their behaviour to support the electricity grid will, given sufficient time, lead to further oscillations. The larger the perturbations are, the sooner these detrimental effects will occur.


## 4.2. Simulating TCLs on the GB Electricity Grid

Our model and simulations have thus far reduced the complexity of the problem by assuming that, apart from the TCL population and the grid frequency, all other network conditions remain constant. This was necessary for our model to be tractable and to ensure that any results from the simulations were attributable to the frequency-sensitive TCL population. An important next step is to consider the TCL population in the context of a real system. In collaboration with the GB System Operator National Grid, we are able to model the GB system with real data from 36 separate 10-day periods during 2015-2016 and simulate what would have happened if a frequency-sensitive fridge population had been active. We consider how the distribution of TCLs changes over this period, and the reduction in the amount of response that other providers needed to supply because of the contribution from the fridges.


### 4.2.1. Methods

We simulate a population of TCLs (specifically fridges) that respond to the grid frequency according to the rules in Section2.1. We use various historic data from National Grid to model real system conditions and simulate the effects of a frequency-sensitive fridge population acting on the GB system. By considering the population in the context of real data including response provision from other sources such as power generators, we are able to get a better understanding of the potential impact of the fridges compared to, say, modelling them in isolation responding to a one-off frequency event.Figure11 gives an overview of the simulation process. Rhombi indicate inputs and outputs; rectangles indicate methods used in the simulation. Methods are applied working downwards, except for the dashed arrows which create an iterative loop.Figure 11
Simulation methodology diagram. Rhombi indicate input or calculated data/simulated data, and rectangles indicate methods/calculations. Events occur from top to bottom with the exception of the dashed arrows which form the iterative loop.


### 4.2.2. Inputs

As shown in Figure11, there are four types of data input, in addition to the fridge population initial conditions. We use 36 consecutive ten-day data samples from the period July 2015–June 2016.Kinetic energy data is an estimate for the total kinetic energy in MVAs (megavolt-ampere seconds) [55]. Values are calculated by summing the kinetic energy of all running synchronised generators (a generator-specific constant provided to the System Operator by each power generator) with an estimate of kinetic energy from demand. The kinetic energy data provided (confidentially from National Grid) is per settlement period (settlement periods split the day into 48 half hour units starting on the hour and half hour) and repeats each value for the full 30 minutes (rather than interpolating). Typical kinetic energy values are in the range 20000–40000 MVAs.Demand data consists of per-second metered demand from National Grid. This is a sum of the power leaving the electricity transmission system, including any power exports through the interconnectors. Half-hourly demand data is accessible via National Grid’s “Data Explorer” [56].Historic frequency data consists of per-second system frequency data in hertz. Frequency measurements are taken in multiple locations to ensure reliable data availability in the event of any metering faults. The frequency data provided by National Grid has undergone a cleaning process that takes advantage of the multiple readings. It is available via National Grid’s “Enhanced Frequency Response” [57].Response holdings are the amount of frequency response delivery in MW (as a function of grid frequency) that National Grid expects each second. Response holdings are positive (or negative) for “low (high) frequency response delivery” when the frequency is below (above) 50 Hz, respectively. For each time step (1 second), 9 different values for response holding are listed. These take the form of primary, secondary, and high response.Primary response values are given for trigger points at 49.9 Hz, 49.5 Hz, and 49.2 Hz. This means that at these frequencies the power response provided through various types of primary response service are the historic response holding values given, subject to a 1 second reaction delay. We assume that the response increases linearly from 0 between 49.985 Hz and 49.8 Hz and likewise linearly between all other frequency trigger values. Below 49.2 Hz the response is assumed to be the constant 49.2 Hz response value. The starting frequency trigger value of 49.985 Hz is used to take into account the Grid Code deadband of (50±0.015) Hz, within which response is not required. Secondary response values are given for frequency trigger points 49.8 Hz and 49.5 Hz, and response is modelled in the same was as for primary response, only with an 11 s response delay. High response values have trigger points 50.2 Hz and 50.5 Hz. Just as for primary response, the time lag is 1 s and again, response is modelled as linear interpolation through these points, starting at the edge of the deadband at 50.015 Hz and remaining constant beyond 50.5 Hz. Figure 12 illustrates an example of how response holding data (Table 4) are interpreted in the model. Values given are indicative only of possible values.Table 4
Illustrative historic response holding data behind Figure12.
Primary
Secondary
High
Frequency trigger (Hz)
49.2
49.5
49.8
49.5
49.8
50.2
50.5
Response (MW)
850
800
430
950
500
−350
−680Figure 12
Representative historic response data with interpolation method for primary response (solid line below 50 Hz), secondary response (dashed line), and high response (solid line above 50 Hz). Zero response in the deadband (50±0.015) Hz.Fridge conditions are the initial on/off state and initial temperature of each fridge in the population. For the simulations presented here we take the zero perturbation case a(i) in Table 3 from the previous section.


### 4.2.3. Calculating the Demand at 50 Hz

Deviations in grid frequency away from 50 Hz affect the total system demand. We make the assumption that demand increases linearly by approximately 2.5% of its value at 50 Hz for every 1 Hz increase in frequency above 50 Hz (and decreases by the same amount as frequency decreases below 50 Hz). In order to know the demand at the nominal frequency, “demand at 50 Hz,”Demω0(t), we need to calculate it from the (measured) demand data input, D(t).(62)Dt=Demω0t1+0.025ft-50(63)Demω0t=Dt1+0.025ft-50.


### 4.2.4. Calculating the Underlying Imbalance

In order to calculate the effects of the fridge population on the system frequency, we first need to calculate the underlying supply-demand imbalance (in MW) that caused the original system frequency deviations away from 50 Hz. At this point it is necessary to distinguish between two important, similar-sounding terms:underlying imbalance andtotal imbalance. By underlying imbalance, Imbunder(t), we mean the generation-demand imbalance that occurs independently of the system frequency. This may be due to, for example, fluctuations in wind or solar power generation or discrepancies between the total predicted system demand and the actual real-time demand. In contrast, total imbalance, Imbtot(f,t), includes both the underlying imbalance and, additionally, what we shall refer to asdynamic imbalance.There are two sources of dynamic imbalance:generator response (frequency response provided by power generators as the frequency changes) anddemand response (the automatic change in demand as frequency changes). Note that in this context “demand response” is completely different to demand-side response services, which, given their current low penetration of the response market, we exclude from our simulations. Generator response, Genresp, consists of the actual response delivered by generators, calculated as described above from the response holdings and the historic system frequency. Generator response is assumed to have a small time lag δt, which we take to be 1 second. In contrast, demand response, Demresp, is assumed to occur instantaneously and is defined as the measured system demand D(t) minus the demand at 50 Hz, Dω0(t) (see “calculating demand at 50 Hz”). Therefore by (62)(64)Demrespft,t=0.025Dω0tft-50.Both sources of dynamic imbalance will change when we introduce the population of responsive fridges (because of their impact on the frequency) and will therefore need to be recalculated.We use a linear approximation for the rate of change of frequency [51], which in our notation is given by(65)dfdt=50Imbtott2Ekt,where 50 is the nominal frequency 50 Hz and Ek=F0M/2 is total stored kinetic energy in MVAs. Since(66)Imbtott=Imbundert+Genrespft-δt,t-Demrespft,twe are able to find(67)Imbundert=Ekt25dfdt-Genrespft-δt,t+Demrespft,twhich for simulation time step size Δt gives(68)Imbundert=Ektft-ft-Δt25Δt-Genrespft-δt,t+Demrespft,t.We take Δt=1 s, so Δt=δt, the generator response time lag. Generator response is calculated using historic response holdings and the frequency t-δt seconds ago along with some constraints on the generator ramp rates.


### 4.2.5. Iterative Loop

Once the underlying imbalance has been calculated for all time steps it can be used along with the response holdings and fridge conditions to begin a loop formed of three calculation steps that iterates over all time steps (see the “iterative loop” in Figure11). The steps are as follows:(1)
Calculate the frequency response delivery from the fridge population and from the dynamic response providers based on the previous frequency value (the first iteration takes the first historic frequency value, after which the “new frequency” values are used). For the fridge population this requires summing the switched on fridges multiplied by their individual power consumption and subtracting the power consumption of the population if the fridges were not frequency-sensitive. Response from the dynamic response providers is described above.(2)
Calculate the new frequencyf∗(t) using the equations from “calculating the underlying imbalance” and beginning with the approximation (69)f∗t=f∗t-Δt+Δtdf∗dtt=f∗t-Δt+Δt25Imbtot∗tEktand since(70)Imbtot∗t=Imbundert+Genrespf∗t-δt,t-0.025Demω0f∗t-50we get(71)f∗t=f∗t-Δt+25Δt/EktImbundert+Genrespf∗t-δt,t+1.25Demω01+0.625δt/EktDemω0.Note that we letf∗(0)=f(0), the original frequency value at time 0.(3)
Calculate the new fridge conditions by updating their temperature set points with the new frequency f∗ calculated in step 2, according to (2a) and (2b). Each fridge temperature is evolved one time step according to (3a) or (3b). If a switch on or off should have occurred during the time step then the exact time of switch is estimated and the temperature is recalculated from the switch time to the end of the time step using linear interpolation.


### 4.2.6. Outputs

There are two key outputs for our analysis: firstly, the temperatures and states of each fridge over time and secondly, the frequency response supplied by all other providers on the grid. Since response can be positive or negative depending on the frequency, but both incur payment, we take the absolute value of the response at each time step. We take the cumulative sum of the difference between this response in the presence of TCLs and the original system response and call it “cumulative response savings,” which we measure in MWh. This allows us to find out how much benefit (or detriment) the fridges provided the system and how that changes over time as they respond to frequency perturbations.


### 4.2.7. Results

We begin by comparing the results of two populations of fridges, one with total power consumption (if every fridge were switched on),Pc=70 MW, as in previous simulations, and the other with Pc=700 MW. The second case is the extreme with 10 million frequency-sensitive fridges, similar to Trovato et al. who modelled 11 million fridges in [33]. Rather than modelling all 1 million or 10 million fridges, we split the population into 104 groups, where the fridges within one group have the same temperature and state. We begin both simulations with the groups uniformly distributed in phase (the unperturbed case) just as in figures a(i) of the previous simulation section. We repeat these simulations on 36 10-day data samples from July 2015–June 2016.Figure13 shows cumulative response savings over time (introduced above) for the 36 data sample simulations in each case. For both values of Pc, in at least a third of the simulations, the fridge population ended up doing more harm than good (negative savings) due to synchronisation. Increasing participation tenfold increased the best results by about a factor of 7 but worsened the worst results by a factor of 15. When there were fewer participants (Pc=70 MW), the results were more erratic over time, which we attribute to there being less response on the system, and so a less smooth frequency trace to respond to. In light of these findings, we present the results from another set of 36 simulations over the year for Pc=700 MW, with the addition of a very small amount of diversity (less than 0.25% in relative terms) to the parameters.Figure 13
Cumulative response savings (MWh), the difference between other providers’ response with and without the frequency-sensitive fridge population (cumulatively) for the 36 data samples over one year for two different participation levels. We sum high frequency response and the absolute value of low frequency response, in MWh. Negative results indicate that the other providers had to compensate for detrimental fridge behaviour. In both cases the population is homogeneous.
(a)
P
c
=
7
0 MW
(b)
P
c
=
700 MWWe find that very small amounts of parameter diversification can eradicate the detrimental fridge behaviour in our simulations. A full presentation of all of our simulations can be found in Webborn Ph.D. thesis, to appear. As an example we takePc=700 and the other parameters to have the same mean as in all of our previous simulations but draw them from normal distributions with nonzero standard deviation. We choose TON~N(-26,0.020), TOFF~N(20,0.0133), T-0~N(2,0.005), (T+0-T-0)~N(5,0.003), and α~N(18.08×105,0.030×105). This results in T+0~N(7,0.005), the duty cycle approximately ~N(33.55%,0.027%), and the cycle period in minutes ~N(41.144,0.108). The results from these simulations are shown in Figure 14. We see that even with this very small amount of parameter diversity, all the simulations show a net and growing benefit to the system over the ten days. Since diversity naturally occurs in any real population, this offers reasonable evidence that TCLs could be a valuable resource for the grid, without the need for stochastic switching or regular communications from a central controller.Figure 14
Cumulative response savings (MWh) (see Figure13) for a heterogeneous fridge population with Pc=700 MW. The introduction of a very small amount of diversity has eradicated the detrimental behaviour in all cases and the fridge population provides a clear benefit to the system.


## 4.2.1. Methods

We simulate a population of TCLs (specifically fridges) that respond to the grid frequency according to the rules in Section2.1. We use various historic data from National Grid to model real system conditions and simulate the effects of a frequency-sensitive fridge population acting on the GB system. By considering the population in the context of real data including response provision from other sources such as power generators, we are able to get a better understanding of the potential impact of the fridges compared to, say, modelling them in isolation responding to a one-off frequency event.Figure11 gives an overview of the simulation process. Rhombi indicate inputs and outputs; rectangles indicate methods used in the simulation. Methods are applied working downwards, except for the dashed arrows which create an iterative loop.Figure 11
Simulation methodology diagram. Rhombi indicate input or calculated data/simulated data, and rectangles indicate methods/calculations. Events occur from top to bottom with the exception of the dashed arrows which form the iterative loop.


## 4.2.2. Inputs

As shown in Figure11, there are four types of data input, in addition to the fridge population initial conditions. We use 36 consecutive ten-day data samples from the period July 2015–June 2016.Kinetic energy data is an estimate for the total kinetic energy in MVAs (megavolt-ampere seconds) [55]. Values are calculated by summing the kinetic energy of all running synchronised generators (a generator-specific constant provided to the System Operator by each power generator) with an estimate of kinetic energy from demand. The kinetic energy data provided (confidentially from National Grid) is per settlement period (settlement periods split the day into 48 half hour units starting on the hour and half hour) and repeats each value for the full 30 minutes (rather than interpolating). Typical kinetic energy values are in the range 20000–40000 MVAs.Demand data consists of per-second metered demand from National Grid. This is a sum of the power leaving the electricity transmission system, including any power exports through the interconnectors. Half-hourly demand data is accessible via National Grid’s “Data Explorer” [56].Historic frequency data consists of per-second system frequency data in hertz. Frequency measurements are taken in multiple locations to ensure reliable data availability in the event of any metering faults. The frequency data provided by National Grid has undergone a cleaning process that takes advantage of the multiple readings. It is available via National Grid’s “Enhanced Frequency Response” [57].Response holdings are the amount of frequency response delivery in MW (as a function of grid frequency) that National Grid expects each second. Response holdings are positive (or negative) for “low (high) frequency response delivery” when the frequency is below (above) 50 Hz, respectively. For each time step (1 second), 9 different values for response holding are listed. These take the form of primary, secondary, and high response.Primary response values are given for trigger points at 49.9 Hz, 49.5 Hz, and 49.2 Hz. This means that at these frequencies the power response provided through various types of primary response service are the historic response holding values given, subject to a 1 second reaction delay. We assume that the response increases linearly from 0 between 49.985 Hz and 49.8 Hz and likewise linearly between all other frequency trigger values. Below 49.2 Hz the response is assumed to be the constant 49.2 Hz response value. The starting frequency trigger value of 49.985 Hz is used to take into account the Grid Code deadband of (50±0.015) Hz, within which response is not required. Secondary response values are given for frequency trigger points 49.8 Hz and 49.5 Hz, and response is modelled in the same was as for primary response, only with an 11 s response delay. High response values have trigger points 50.2 Hz and 50.5 Hz. Just as for primary response, the time lag is 1 s and again, response is modelled as linear interpolation through these points, starting at the edge of the deadband at 50.015 Hz and remaining constant beyond 50.5 Hz. Figure 12 illustrates an example of how response holding data (Table 4) are interpreted in the model. Values given are indicative only of possible values.Table 4
Illustrative historic response holding data behind Figure12.
Primary
Secondary
High
Frequency trigger (Hz)
49.2
49.5
49.8
49.5
49.8
50.2
50.5
Response (MW)
850
800
430
950
500
−350
−680Figure 12
Representative historic response data with interpolation method for primary response (solid line below 50 Hz), secondary response (dashed line), and high response (solid line above 50 Hz). Zero response in the deadband (50±0.015) Hz.Fridge conditions are the initial on/off state and initial temperature of each fridge in the population. For the simulations presented here we take the zero perturbation case a(i) in Table 3 from the previous section.


## 4.2.3. Calculating the Demand at 50 Hz

Deviations in grid frequency away from 50 Hz affect the total system demand. We make the assumption that demand increases linearly by approximately 2.5% of its value at 50 Hz for every 1 Hz increase in frequency above 50 Hz (and decreases by the same amount as frequency decreases below 50 Hz). In order to know the demand at the nominal frequency, “demand at 50 Hz,”Demω0(t), we need to calculate it from the (measured) demand data input, D(t).(62)Dt=Demω0t1+0.025ft-50(63)Demω0t=Dt1+0.025ft-50.


## 4.2.4. Calculating the Underlying Imbalance

In order to calculate the effects of the fridge population on the system frequency, we first need to calculate the underlying supply-demand imbalance (in MW) that caused the original system frequency deviations away from 50 Hz. At this point it is necessary to distinguish between two important, similar-sounding terms:underlying imbalance andtotal imbalance. By underlying imbalance, Imbunder(t), we mean the generation-demand imbalance that occurs independently of the system frequency. This may be due to, for example, fluctuations in wind or solar power generation or discrepancies between the total predicted system demand and the actual real-time demand. In contrast, total imbalance, Imbtot(f,t), includes both the underlying imbalance and, additionally, what we shall refer to asdynamic imbalance.There are two sources of dynamic imbalance:generator response (frequency response provided by power generators as the frequency changes) anddemand response (the automatic change in demand as frequency changes). Note that in this context “demand response” is completely different to demand-side response services, which, given their current low penetration of the response market, we exclude from our simulations. Generator response, Genresp, consists of the actual response delivered by generators, calculated as described above from the response holdings and the historic system frequency. Generator response is assumed to have a small time lag δt, which we take to be 1 second. In contrast, demand response, Demresp, is assumed to occur instantaneously and is defined as the measured system demand D(t) minus the demand at 50 Hz, Dω0(t) (see “calculating demand at 50 Hz”). Therefore by (62)(64)Demrespft,t=0.025Dω0tft-50.Both sources of dynamic imbalance will change when we introduce the population of responsive fridges (because of their impact on the frequency) and will therefore need to be recalculated.We use a linear approximation for the rate of change of frequency [51], which in our notation is given by(65)dfdt=50Imbtott2Ekt,where 50 is the nominal frequency 50 Hz and Ek=F0M/2 is total stored kinetic energy in MVAs. Since(66)Imbtott=Imbundert+Genrespft-δt,t-Demrespft,twe are able to find(67)Imbundert=Ekt25dfdt-Genrespft-δt,t+Demrespft,twhich for simulation time step size Δt gives(68)Imbundert=Ektft-ft-Δt25Δt-Genrespft-δt,t+Demrespft,t.We take Δt=1 s, so Δt=δt, the generator response time lag. Generator response is calculated using historic response holdings and the frequency t-δt seconds ago along with some constraints on the generator ramp rates.


## 4.2.5. Iterative Loop

Once the underlying imbalance has been calculated for all time steps it can be used along with the response holdings and fridge conditions to begin a loop formed of three calculation steps that iterates over all time steps (see the “iterative loop” in Figure11). The steps are as follows:(1)
Calculate the frequency response delivery from the fridge population and from the dynamic response providers based on the previous frequency value (the first iteration takes the first historic frequency value, after which the “new frequency” values are used). For the fridge population this requires summing the switched on fridges multiplied by their individual power consumption and subtracting the power consumption of the population if the fridges were not frequency-sensitive. Response from the dynamic response providers is described above.(2)
Calculate the new frequencyf∗(t) using the equations from “calculating the underlying imbalance” and beginning with the approximation (69)f∗t=f∗t-Δt+Δtdf∗dtt=f∗t-Δt+Δt25Imbtot∗tEktand since(70)Imbtot∗t=Imbundert+Genrespf∗t-δt,t-0.025Demω0f∗t-50we get(71)f∗t=f∗t-Δt+25Δt/EktImbundert+Genrespf∗t-δt,t+1.25Demω01+0.625δt/EktDemω0.Note that we letf∗(0)=f(0), the original frequency value at time 0.(3)
Calculate the new fridge conditions by updating their temperature set points with the new frequency f∗ calculated in step 2, according to (2a) and (2b). Each fridge temperature is evolved one time step according to (3a) or (3b). If a switch on or off should have occurred during the time step then the exact time of switch is estimated and the temperature is recalculated from the switch time to the end of the time step using linear interpolation.


## 4.2.6. Outputs

There are two key outputs for our analysis: firstly, the temperatures and states of each fridge over time and secondly, the frequency response supplied by all other providers on the grid. Since response can be positive or negative depending on the frequency, but both incur payment, we take the absolute value of the response at each time step. We take the cumulative sum of the difference between this response in the presence of TCLs and the original system response and call it “cumulative response savings,” which we measure in MWh. This allows us to find out how much benefit (or detriment) the fridges provided the system and how that changes over time as they respond to frequency perturbations.


## 4.2.7. Results

We begin by comparing the results of two populations of fridges, one with total power consumption (if every fridge were switched on),Pc=70 MW, as in previous simulations, and the other with Pc=700 MW. The second case is the extreme with 10 million frequency-sensitive fridges, similar to Trovato et al. who modelled 11 million fridges in [33]. Rather than modelling all 1 million or 10 million fridges, we split the population into 104 groups, where the fridges within one group have the same temperature and state. We begin both simulations with the groups uniformly distributed in phase (the unperturbed case) just as in figures a(i) of the previous simulation section. We repeat these simulations on 36 10-day data samples from July 2015–June 2016.Figure13 shows cumulative response savings over time (introduced above) for the 36 data sample simulations in each case. For both values of Pc, in at least a third of the simulations, the fridge population ended up doing more harm than good (negative savings) due to synchronisation. Increasing participation tenfold increased the best results by about a factor of 7 but worsened the worst results by a factor of 15. When there were fewer participants (Pc=70 MW), the results were more erratic over time, which we attribute to there being less response on the system, and so a less smooth frequency trace to respond to. In light of these findings, we present the results from another set of 36 simulations over the year for Pc=700 MW, with the addition of a very small amount of diversity (less than 0.25% in relative terms) to the parameters.Figure 13
Cumulative response savings (MWh), the difference between other providers’ response with and without the frequency-sensitive fridge population (cumulatively) for the 36 data samples over one year for two different participation levels. We sum high frequency response and the absolute value of low frequency response, in MWh. Negative results indicate that the other providers had to compensate for detrimental fridge behaviour. In both cases the population is homogeneous.
(a)
P
c
=
7
0 MW
(b)
P
c
=
700 MWWe find that very small amounts of parameter diversification can eradicate the detrimental fridge behaviour in our simulations. A full presentation of all of our simulations can be found in Webborn Ph.D. thesis, to appear. As an example we takePc=700 and the other parameters to have the same mean as in all of our previous simulations but draw them from normal distributions with nonzero standard deviation. We choose TON~N(-26,0.020), TOFF~N(20,0.0133), T-0~N(2,0.005), (T+0-T-0)~N(5,0.003), and α~N(18.08×105,0.030×105). This results in T+0~N(7,0.005), the duty cycle approximately ~N(33.55%,0.027%), and the cycle period in minutes ~N(41.144,0.108). The results from these simulations are shown in Figure 14. We see that even with this very small amount of parameter diversity, all the simulations show a net and growing benefit to the system over the ten days. Since diversity naturally occurs in any real population, this offers reasonable evidence that TCLs could be a valuable resource for the grid, without the need for stochastic switching or regular communications from a central controller.Figure 14
Cumulative response savings (MWh) (see Figure13) for a heterogeneous fridge population with Pc=700 MW. The introduction of a very small amount of diversity has eradicated the detrimental behaviour in all cases and the fridge population provides a clear benefit to the system.


## 5. Discussion and Conclusions

The combination of our mathematical analysis and simulations has improved our understanding of how a population of frequency-sensitive TCLs might act on the electricity grid, in a number of ways. Our analysis in Section3.1 was able to capture the short-term benefits of using identical TCLs for frequency response. However, the simulations of the model were important in showing that, in response to realistic forcing, in the long term, nonlinear dynamics occur that could not be captured with our linearisation. Studying the two-group case indicated that the fully synchronised periodic solution is unstable to splitting it into two synchronised groups of similar size. A continuum of TCLs, as analysed in Section 3.1, is the limit as N goes to infinity of N equally sized groups of TCLs. Therefore the stability found for the uniformly distributed continuum of TCLs can be compared to the case with two groups of similar size remaining separate. If the distribution of TCLs is perturbed too far from uniform, as in the simulations in Section 4, then this is similar to the two-group case in which the groups are not similarly sized.Simulating a population of identical frequency-sensitive TCLs under typical system conditions revealed that for many of the data samples the short-term benefits were outweighed by switching behaviour that requires greater frequency response from the rest of the system than when the TCLs were not frequency-sensitive. In these cases regular communications would need to be sent to the TCLs to desynchronise. However, we find that the addition of a very small amount of parameter diversity, for example, 0.24% variation in the natural cycle period, which is likely to occur naturally in a population, can eradicate these issues. In our example the population was able to reduce the response required from other providers in all periods of the year studied.A number of open questions remain. How would factors such as daily room temperature variations, opening the door, and changing the contents of TCLs like fridge-freezers affect the cycle distribution? How would a population cope during more severe frequency incidents than existed during 2015-16? What would be the effects of modelling the population on a network where frequency variations can be spatially dependent? How much diversity exists in real TCL populations and can the effects of diversity be understood theoretically?In conclusion, this study indicates that a population of fridges might perform a valuable service to the grid without requiring centralised or stochastic control.



---

*Source: 5031505-2017-12-06.xml*