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.
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 and decentralised 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.
The 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 Section 2 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 be denoted by , the cooling/heating coefficient by , and the asymptotic temperatures that the TCL would reach if left on and off indefinitely by and , respectively. ThenA (cooling) TCL will switch off when the temperature reaches its lower temperature set point and switch on when it reaches its upper temperature set point . We choose to make these set points sensitive to system frequency deviations away from 50 Hz, denoted (i.e., FrequencyHz). Insufficient generation to meet demand causes and so we need the TCLs to reduce their power consumption to bring 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 , 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,where , are positive constants that determine the sensitivity of the lower and upper temperature set points to frequency deviations. and are the uncoupled (a fridge is “uncoupled” from the grid frequency if ) temperature set points, which we typically take to be C and 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 . If a TCL has temperature at time and does not switch on or off before time then the temperature is given byWe can rearrange (3a) and (3b) and solve for the on and off durations and , respectively, assuming constant grid frequency: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 ) their “natural” on and off cycle durations, and , are given byIn order for the TCLs to operate properly they need to cycle on and off, and so we require thatWe 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 iswhich is a nonempty interval (notably containing ).
2.2. Electricity Grid Frequency
A simplified equation for the frequency of a power system can be determined by Newton’s 2nd Law of Motion or the derived equation for energy. If we let , where is the nominal grid frequency (50 Hz in Europe) and linearise about , then we obtain [26]and for brevity we introduce new variables along with explicit consideration of TCL power consumptionwhere stands for times nominal angular momentum of the rotating masses in the system, stands for total inertia of the rotating masses of the system, stands for damping factor representing the natural frequency dependence of the load alongside the damping provided by synchronous generator damper windings, stands for change in total active power generation, compared to a reference level, stands for change in total active power load, compared to a reference level, stands for inverse nominal angular momentum, introduced for brevity, 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, stands for power consumed by TCL population when all switched on, is a variable introduced for brevity.
We make the simplifying assumption in Sections 2 and 3 that the “surplus” power generation on the system for TCL consumption is a constant. We use the “” notation to denote equilibrium values. In equilibrium henceand therefore we can rewrite our equation for in terms of deviations from equilibrium values:where
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 % of fridges in GB. We take the power consumed by an individual fridge when switched on, , to be 70 W, as assumed in [32, 33]. This means that we let MW and the total power consumption if all fridges were switched on, MW. Using our approximation for [51], . Our GB system data (discussed later) gives an approximate average value for total stored kinetic energy ; MVAs (note that MVAs = MJ), and so . We let vary between 0 and 1 by changing . Our parameters are summarised in Table 2, and throughout this paper take these values unless stated otherwise.
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 , in such a way that dictates both the temperature and the state of a TCL. The switched off TCLs are mapped to the interval and the switched on TCLs are mapped to . Then we define the position of a TCL at time with temperature and state on or off byNote that the model implicitly assumes that the temperature set points never change fast enough to leave a TCL outside of the interval . 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 denote the fraction of TCLs that lie between and at time . Then is nonnegative, with period length 2 in , and satisfies the normalisationfor all . The evolution of is governed by the continuity equation [53]where is the velocity of a TCL in -space, . Differentiating (14) giveswhereNote that, for satisfying (7), and . Under a constant grid frequency, and are constants. In equilibrium we have , and therefore (16) implies for some constant . Since , from (17a) and (17b) we haveThen for all , respectively,and is determined by the normalisation criterion (15),The proportion of TCLs switched on, , is given byIn equilibrium (12), 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 by a small amount and to consider the evolution of the perturbation. By (15) the perturbation satisfiesWe writeso that (16) becomesTaking the first-order approximation yieldsRearranging (27) for and substituting (17a) and (17b) for giveHence where we have definedWe have a time-invariant linear system (29), and so it is natural to look for solutions for which the time dependence of our variables and is ; is called an eigenvalue of the system. Defining and renaming to , (29) becomesWe introduce an integrating factor so that on the open intervals we can find an expression for :At the discontinuities and ,We can use (34) to find expressions for and and substitute these into (35b). After substitution for (or ) using (35a) and rearrangement we arrive atwhereRewriting our equation for the rate of change of grid frequency near equilibrium (9) asand setting give Integrating (34) over in (the switched on TCLs), setting the resulting expression equal to the right hand side of (38), and substituting our expression in (36a) for establish the following implicit equation for : where we have defined , which reflects the strength of the effect of the TCLs on grid frequency.
When (no effect of the TCLs on the grid frequency) the eigenvalue equation (39) reduces to , so the eigenvalues are and for (the roots of ). It can also be seen from (39) that for all there is an eigenvalue . 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 to vary from its value derived from the table, by . There is an infinite sequence of eigenvalues going upwards and their reflections in the real axis. Increasing 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.
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 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 subscript denote the th switch on and th switch off event. Without loss of generality, suppose that after our initial start time the next switch event is the population switching on. This implies that, for all , . Figure 2 illustrates the notation. Hence the amount of time the population spends switched on following the th switch on event is given bywhere are the frequencies at the th switch on and off times. The amount of time spent switched off following the th switch off is given byAssuming, as for the numerical analysis in Section 3.1, that the system has no damping, we set in (12) for . In a synchronised population, at time all TCLs are either on () or off (). Then we can define constants such thatHence the values of at the switch off and on times are given by the piecewise-linear functionsAfter substituting for the switching times using (40a) and (40b) and rearranging, these becomeNow since each side of (43a) and (43b) are functions of only one of the variables, we can explicitly name them as suchEach of the four functions is increasing and therefore invertible, and so we can writeand thereforewhich is a mapping from the frequency at one switch on event to the frequency at the next. The mapping is a contraction iff Note thatSimilarlyTherefore a sufficient condition for the mapping to be a contraction is thatwhich is a nonempty interval (containing ), so long as and for all . It is worth recalling our earlier condition on the values of (7) which also imposed that belong to an open interval containing .
3.2.2. Solving for the Periodic Solution
The contraction property of the mapping (45c) under the above conditions implies that there is an attracting fixed point so long as , and hence a periodic solution for the synchronised population. We now seek to solve for this periodic solution. Denote by and 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 and . Therefore the temperature of the population will cycle between upper and lower set points, given by and , respectively. Equations (42a) and (42b) show us thatThe temperature evolution equations (3a) and (3b) allow us to express the switch on and switch off temperatures as follows:which after substituting for using (51a) and rearranging becomeNow we have two equations in terms of , and , which we can combine into one equation and eliminate ,We can also express in terms of by summing (51a) and (51b) to giveand so (53b) and (55) form a pair of coupled equations for and , which can be solved numerically. Figure 3 shows one temperature cycle for the single group under different choices for . Denote by the value of when . As gets further away from the solutions drift further from the uncoupled temperature range (2–7°C). The cycle lengths are symmetric about but the TCLs consume more power per cycle as increases.
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 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 . As before, we denote the switched on duration in this solution by and the switched off duration by . Suppose that the population is split into two groups A and B, such that proportion belongs to group A, and proportion belongs to group B. Suppose also that group B switches on at time , followed soon after by group A switching on, at time . Then after a time period of length similar to group B switches off, which is again followed shortly after by group A switching off. After a time period similar to 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 . 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 and , the temperature difference when B switches on the first and second times, respectively. In addition, let and , the time difference between the two groups switching on the first and second times, respectively. Further notation is shown in Figure 4.
In order to calculate and we need to calculate the switch times and temperatures of the two groups at each switch event leading up to . 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 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 and group B switched off at time 0,In addition, by the temperature evolution equations,Equating (56a) and (56b) and introducing our new notation giveIf we write and take and small, then and linearising in giveswhereMore 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 and linearising about the single group solution, we find in Appendix B that whereSo is a left eigenvector of the linearised map in the space of with eigenvalue . We can plot against for various to see whether (in which case the two groups merge into one) or whether (they move apart). The results are shown in Figure 5. We find that the second eigenvalue is within the interval for any and for our parameter values (taken from Table 2 with the exception of 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 .
(a) Full interval
(b) Enlargement
By solving for the dividing case 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 . 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.
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 ) 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 Section 3.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 ) 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 Section 3.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 Section 2, 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 “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 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 divided by the number of agents, . 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 distribution we can alter the number of TCLs switched on or off from the equilibrium proportions and , 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 , and we refer to the normalised wave peak amplitude (normalised by dividing by ). This normalisation means that when we plot , 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.
Figure 8 shows the effects of these perturbations on the initial conditions in each case, plotting 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.
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 Section 3.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 Section 2.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.
Figure 11 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.
4.2.2. Inputs
As shown in Figure 11, 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 () 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.
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,” , we need to calculate it from the (measured) demand data input, .
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 and total imbalance. By underlying imbalance, , 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, , includes both the underlying imbalance and, additionally, what we shall refer to as dynamic imbalance.
There are two sources of dynamic imbalance: generator response (frequency response provided by power generators as the frequency changes) and demand 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, , 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 , which we take to be 1 second. In contrast, demand response, , is assumed to occur instantaneously and is defined as the measured system demand minus the demand at 50 Hz, (see “calculating demand at 50 Hz”). Therefore by (62)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 bywhere 50 is the nominal frequency 50 Hz and is total stored kinetic energy in MVAs. Sincewe are able to findwhich for simulation time step size givesWe take s, so , the generator response time lag. Generator response is calculated using historic response holdings and the frequency 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 Figure 11). 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 frequency using the equations from “calculating the underlying imbalance” and beginning with the approximation and since we get Note that we let , the original frequency value at time 0.(3) Calculate the new fridge conditions by updating their temperature set points with the new frequency calculated in step , 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), MW, as in previous simulations, and the other with 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 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.
Figure 13 shows cumulative response savings over time (introduced above) for the 36 data sample simulations in each case. For both values of , 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 ( 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 MW, with the addition of a very small amount of diversity (less than 0.25% in relative terms) to the parameters.
(a) MW
(b) MW
We 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 take 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 , , , , and . This results in , the duty cycle approximately , and the cycle period in minutes . 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.
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 Section 3.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 goes to infinity of 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.
Appendix
A. Derivation of (7)
Sufficient condition (7) for a fridge to change its power consumption as required by the system frequency is derived as follows. We would like to know how the typical power consumption over one cycle changes as system frequency changes. Typical power consumption per TCL per cycle, , is given by where is the instantaneous power consumption of a TCL when switched on (assumed to be independent of time in its on-phase). From (2a), (2b), (4a), and (4b),ThereforeSince , is strictly positive if and only ifwhich is the case if and only ifTherefore a sufficient condition for the derivative of with respect to to be positive is that both terms in the preceding equation should be strictly positive. Since , , and the squared terms are always strictly positive this leaves us with two sufficient criteria:The upper bound is greater than 1 and the lower bound is less than 1 if .
B. Derivation of (61): The First Eigenvalue in the Two-Group Case
The temperature evolution equations tell us thatThe temperature set point equations provide us withEquations (B.1a) and (B.2a) determine in terms of , (B.1a) and (B.1b) determine and so forth, and hence are determined by . To analyse the linear stability of the fixed point of this map corresponding to the one group solution (Section 3.2.2) we find that depends only on and some differences of switching times, by eliminating the temperatures at the intermediary switch times.Defining and linearising about the single group solution using the notation from Figure 4 (signed displacement from the single group switch times) giveSince (B.1a) and (B.2a) are two equations for , (B.1d) and (B.2b) are two equations for , (B.1e) and (B.2c) are two equations for , and (B.1h) and (B.2d) are two equations for , we haveWe used (B.5a) already to determine in terms of to first-order (59). Denote , then to find to first order we linearise (B.5c) to obtainNow we can find a substitution for and and we can use this substitution to arrive atWith our expression for and for using (59), we arrive at where
C. Derivation of the Second Eigenvalue in the Two-Group Case
In Appendix B we found one of the eigenvalues (61) for the case of two groups of TCLs whose cycling was close to the single group solution and claimed that the other eigenvalue was insignificant for determining the stability of the system. Here we derive bounds on the second eigenvalue to prove this claim.
The temperature cycles of groups A and B are initially very close to the single group temperature cycle and therefore, for , we write Our approach is to seek a map such that Taking linear approximations as in Section 3.2.3, (B.5a) approximates to whereWe can use this expression for and take a first-order approximation of (B.1a) to find an expression for in terms of :Similarly, using (B.1b) we find thatWe repeat this process to find expressions for each time interval (recall that is the first interval between switch times) for and each subsequent and . We use (B.5b) to find whereEquations (B.1c) and (B.1d) thus yieldUsing (B.5c) we find whereEquations (B.1e) and (B.1f) yieldFinally (B.5d) gives whereand (B.1g) and (B.1h) giveFor each we can write and sowhereWe introduce the following simplifying notation before defining each matrix . LetThenWe already know one of the eigenvalues of the system (), and the second eigenvalue is given by . Note that , andWe can rewrite (see (61)) in our new notation aswhich allows us to writeDenote the second eigenvalue by , and useSince , , and , which for our choice of parameter values gives Therefore the absolute value of the second eigenvalue is (significantly) less than 1, and so the stability of the system is determined by the first eigenvalue .
Conflicts of Interest
Ellen Webborn was formerly an employee at National Grid. National Grid contributed code and data (but no financial support) to this project. The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the EPSRC through Grant no. EP/I01358X/1. The authors are very grateful to National Grid for data, code, and many engaging discussions, in particular to Roy Cheung, Matthew Roberts, and William Ramsay. They are also very grateful to Michael Waterson and Lisa Flatley for their valued input and to numerous colleagues at the Centre for Complexity Science at the University of Warwick for all of their support.