Human in the loop heterogeneous modelling of thermostatically 1 controlled loads for Demand Side Management studies.

Demand Response (DR) is a Smart Grid technology aiming to provide demand regulation for dynamic pricing and ancillary services to the grid. Thermostatically controlled loads (TCLs) are among those with the highest potential for DR. Some of the challenges in modelling TCLs is the various factors that affect their duty cycle, mainly human behaviour and external conditions, as well as heterogeneity of TCLs (load parameters). These add an element of stochasticity, with detrimental impact on the aggregated level. Most models developed so far use Wiener processes to represent this behaviour, which in aggregated models, such as those based on Coupled Fokker-Planck Equations (CFPE), have a negligible effect as “ white noise ” . One of the main challenges is modelling the effect of external factors on the state of TCLs' aggregated population and their impact in heterogeneity during operation. Here we show the impor- tance of those factors as well as their detrimental effect in heterogeneity using cold loads as a case study. A bottom up detailed model has been developed starting from thermal modelling to include these factors, real world data was used as input for realistic results. Based on those we found that the duty cycle of some TCLs in the population can change signi ﬁ cantly and thus the state of the TCLs' population as a whole. Subsequently, the accuracy of aggregation models assuming relative homogeneity and based on small stochasticity (i.e. Wiener process with typical variance 0.01) is questionable. We anticipate similar realistic models to be used for real world applications and aggregation methods based on them, especially for cold loads and similar TCLs, where external factors and heterogeneity in time are signi ﬁ cant. DR control frameworks for TCLs should also be designed with that behaviour in mind and the developed bottom up model can be used to evaluate their accuracy.


Introduction
As more intermittent renewable and distributed energy comes onto the grid, old generation resources are retired and costly new plant expenditures are avoided. Transmission congestion in Medium Voltage and Low Voltage becomes a larger problem 1], and therefore new Demand Response (DR) programs are emerging. Economic optimization under dynamic pricing, as well as provision of ancillary services and regulation through DR has been in the spotlight of Smart Grid research [2,3]. Balancing demand and supply has predominantly been achieved by supply regulation, but in some cases by demand regulation too. Theoretically any regulation measure from the "supply side" has an equivalent countermeasure from the "demand side". In practice, generation units have different characteristics than loads, and DR consequently has to be engaged as such, thus examine DR's suitability for ancillary services and form regulations specifically for DR. It is imperative to note that DR comes from aggregation of large populations of loads with different characteristics. It would be challenging to track the state of each one in real time, but estimating the state of the population probabilistically is possible and can be accurate [4]. Moreover, based on the type of loads aggregated, there are different characteristics of the population; accordingly, constraints and potential for ancillary services of each population. The problem of aggregation and scalability is one of the major challenges DR faces, whose potential solution could be clustering, a topic discussed in this paper and will be examined more in depth in future work. DR programs (mainly from industrial units) provide ancillary services [6] or regulation [7] and can replace wholly or partial conventional units; commercial and residential DR services are under development and have relatively low enrolment to date [5]. Basic types of DR services can be seen in Fig. 1.
Not only can DR substitute conventional units for ancillary services, but they also bring several advantages including lower cost, no CO 2 emissions, instantaneous response [6], ability to reduce load peaks [8] and synergetic potential with Renewable Energy Sources (RES) intermittence [9]. Additionally, unique characteristics of DR can change the way balancing services work today. DR comes from loads, which are ubiquitous throughout the grid, giving the unique ability of spatial actions. Large distributed populations also have inherently higher reliability than conventionally units, in case an individual conventional unit becomes unavailable there is a great impact in overall availability. In addition, contingencies can be countered through DR, through grid relief, whilst a generation unit cannot provide such service. Modern approaches include "third party" aggregators, such as the OpenADR, originally developed at Lawrence Berkeley National Laboratory (LBNL), Berkley, CA [10]. The OpenADR Alliance formed in 2010, is an example of successful aggregated DR implementation with over 130 member companies (as of May 2016). Unfortunately, though, current implementations are unsophisticated regarding the impact to consumers. As the community has stated, with new AMI being deployed access to DR is inexpensive and what separates us from realising the full potential of DR is development of proper load models and their control strategies [11]. This paper focuses on examining the importance of external factors and creating a realistic model using cold loads as a case study.
A lot of effort has been focused on TCLs, due to the fact that they consume a large portion of the total electrical energy, in some buildings 40%e50% [7,12], while cold loads contribute 25%e36% of domestic consumption in UK [13]. TCLs are the most ubiquitous type of load, and especially in the case of cold loads they are plugged in 24/7, with the highest availability time wise [13]. On the other hand, when electric loads are used for space heating and cooling (i.e. air conditioners), they are usually one of the main sources of peaks in demand. Prediction of load dynamics is essential for DR optimization, both non-dispatchable and dispatchable, hence for power systems cost and reliability. For load forecasting purposes and energy efficiency, regression methods [14,15], can successfully be applied. Yet, prove inadequate for precise control of TCLs and similar loads (refrigeration, electric space and/or water heating, air conditioners, heat pumps, electric vehicles etc.), which is essential for ancillary services such as provision of secondary frequency control or short term operating reserve [2]. DR actions alter demand out of its "ordinary" state, thus rendering past load measurements at best marginally useful, even more so in the case of dispatchable DR for emergencies (i.e. reserve, frequency control) [16]. Additionally, load behaviour and rebound effects at recovery time are not captured by forecasting models [16]. The model developed in this work is applicable to all TCLs, the focus is on DR for ancillary services, though it might also be used alongside regression models for faulty operation detection or energy efficiency. We choose cold loads for simulations and case studies as one of the most suitable type of TCLs.
Older DR concepts predefine aggregated load trajectory by suggesting load interruption [7,17], which causes consumer discomfort and risks withdrawal from participation in DR. Even those who do, such as fuzzy logic-based approaches for ancillary services [18] or peak load shaving and operating cost reduction [19], do not take into account the realistic limits of TCLs due to thermal characteristics. In addition, rebound effects were not taken into consideration.
Therefore, more sophisticated models have been favoured for TCLs, where thermal models are used extensively. Totu et al. [20] use a stochastic hybrid system for individual loads and a system of CFPE to describe the aggregated population, in order to develop a control algorithm for the aggregated population. The same individual modelling is used by Callaway [21], who focuses on finding the exact solution of the CFPE for homogeneous TCLs. In Ref. [22], Koch et al. a different modelling approach is considered, through stored thermal energy, and concludes in similar first order stochastic equations. This TCL model (or variations of it) has been used by the majority, whilst the main effort has been on obtaining more accurate aggregation methods or on control theory, such as Perfumo's et al. model-based control [23], or in Callaway's aggregation through CFPE [6,24], or Mathieu's et al. Markov Chain aggregation and control method [9] to name a few. A comparison of the main methods in summarized in Ref. [25]. It is important to note, that the above mention TCL model originates from Mortensen and Haggerty's work [26]; in general models described by Stochastic Partial Deferential Equations (S-PDE). The purpose of it was to study the "cold load pick up" phenomenon, for specific moments of the day, thus "static" conditions at those given moments. In Malham e and Chong's work [27], it was proved that it is possible to construct a system of CFPE (equation (5) and (6) temperature in short periods and c) the above S-PDE does indeed represent the population and its behaviour accurately. More recently, Callaway went a step ahead [21], giving an exact solution to the CFPE, which provides the eigenvalues governing transient deviations from the steady state temperature distribution. Additionally, showing that small heterogeneity (described as "relative homogeneity" in this paper) within the population does actually improve the behaviour of the linear TCL model and associated control strategies, but for high heterogeneity renders the model inadequate. Those assumptions and points were the motivation for this paper. The main equation (1) has 3 basic parts; the (simplified) time varying solution of the differential equation that governs heat transfer mechanics (Newton's law of cooling), a discrete state varying parameter denoting the switch on/off state of TCLs and a stochastic part, commonly a Wiener process (Gaussian with mean 0, variance s typically around 0.01 Cs À1/2 or similar distributions).
Annotation: q temperature ( C), m switch state {0,1} [21], W Wiener process (Gaussian distribution with variance hs 2 ), P power rate (Watts), C thermal capacitance (kWh/ C) and R thermal resistance ( C/kW), q s set-point temperature, d width of temperature deadband, q min minimum (switch OFF) temperature, q max maximum (switch ON) temperature, h efficiency coefficient, l time step.
The main equations of Malham e and Chong [27], as used by Callaway [21], can be written in continuous time for cooling/cold loads as:  [27]. This Wiener process is mainly used to model "random" factors (i.e. cold loads temperature change, door opening etc.) affecting the operation of TCLs, with mean value 0 and small variance. The argument here is whether a Wiener process is adequate for such modelling and if the above equations, as well as the following (equation (7) and (8)), hold true otherwise.
r≡ À 1 CR ðqðtÞ À q a Þð C=hÞ (8) Equations (7) and (8) express the cooling (c) and heating rates (r) respectively, Malham e & Chong [27] approximated the stationary solution of the CFPE assuming that the above rates are independent of temperature. Callaway [21] instead solved the CFPE eigenvalue problem to find the exact solution. An important observation at this point is that the rates are practically linearly dependent on the temperature difference, also steady ambient (q a ) temperature in time is assumed.
One of the key contributions of this paper is increased insight in the importance of human behaviour in the operation of TCLs and subsequently aggregated behaviour of TCLs' population. Moreover, a new model is developed for that purpose which describes an individual TCL realistically. Our results suggest that clustering TCLs based on appliance parameters (R, C, P) will not result in homogeneous or relatively homogeneous populations, an essential assumption for most aggregation models and the developed control algorithms. Different clustering approaches and/or aggregation models need to be developed. The detailed bottom up MC model developed has high heterogeneity inherently, including the human factor, thus suitable for realistic DR simulations and to access the performance of aggregation models and control algorithms before real world experimental studies.
Section 2 reviews existing load models in literature and differences with real world data, Section 3 details our model and the human factor, Section 4 details Monte Carlo aggregation for heterogeneous populations (case study of cold loads), Section 5 presents our simulation results, and Section 6 concludes.

Real world data and TCL loads' population analysis
The original model, is suitable for "cold load pick-up", which is basically the response (demand) of TCLs after brown outs and/or black outs. The demand at that period may be assumed as the demand likely to occur in a similar given moment (date, day period, weather conditions), with a small "noise" factor for stochasticity. Even though this model is suitable to represent a "snap shot" of TCLs' population for a given moment (with knowledge of similar conditions, based on historic data), it doesn't reflect the dynamic nature of TCLs' population throughout the day, as well as the main driving factors behind this behaviour. As seen in some of these studies, the aggregated profile has significant differences with actual demand profiles. The best showcase being the cold loads, which are plugged in 24 h a day, where their aggregated demand is displayed as an almost straight line with a small fluctuation (noise). In reality, throughout the day, cold loads' aggregated demand is not almost linear, but rather seems to follow the trend of the rest demand, lower overnight and higher during the day, showing peaks around the same time as the total demand, in both residential and commercial cases. El-F erik & Malham e [16] proposed an identification algorithm to calculate and update the parameters of the model, thus the model is able to cope with real world changes, with sampling intervals 15.198 min. or less. Assumptions of relatively steady ambient temperature (q a ), small s (practically no significant external disturbance) and homogeneity are required for this model.
In Fig. 2 the daily mean cold load from real world measurements is shown [28,29]. There are significant deviations in daily cold loads, yet the majority of the test cases used in the literature do not take that into account; e.g. in the Markov Chain aggregation used in Ref. [9], or the CFPE aggregation methods [21], which would normally change their statistic properties dynamically. The changes in consumption can be attributed to various factors [30], mainly ambient temperature changes and human interaction. Ambient temperature is the most important factor for cold loads [30], their contribution to the average domestic load is at least 36% in summer and 25% in winter, but there are big variations between households [13]. As observed in Fig. 2, early morning hours show a significant increase in mean fridge consumption, which can be attributed to human interaction (door opening) and heating when residents wake up.
The important note here is that only part of the population causes this increase, thus in reality that part has on average a higher increase than the mean displayed. It also varies between individuals. This is a good indication of the cold loads' population dynamics and how even identical appliances (homogeneous population characteristics) will behave differently due to external factors in time (heterogeneous population in operation, duty cycle). This is mainly due to human behaviour and the effect of human interaction on cold loads power consumption was studied in Refs. [30,31], where all concluded to it been the most important factor in energy consumption. The overall deviation in domestic demand (and behaviour) was studied in Ref. [32], which was based on a detailed model of domestic demand using real data from Time of Use Survey (TUS) [33], in order to keep human behaviour as realistic as possible. Fig. 3 shows two clusters (classification based) where there is an obvious difference in demand during working hours, as expected. This shows the clear difference in behaviour and human presence in households, as seen around working hours (9am e 7pm) there is distinct difference in demand around 150e200 W. The lowest demand for the second cluster (employed occupant) is attribute mainly to passive consumption of cold loads, with minimal to no interaction, whilst the first cluster shows human presence; interaction with appliances and most likely different household temperature (heating). As such cold loads' operation of even identical appliances (parameters R, C, P of Malham e and Chong's model [27]) between these groups differ due to different ambient temperature and human interaction.
The aggregated power in (3) is not completely accurate, since it assumes that units on idle state consume no power. In reality most TCLs have a small consumption during idle state. Power demand is written more accurately as: For a large enough number of loads N (Kolmogorov's law): However, this can be assumed for homogenous or almost homogenous populations only. An operational characteristic of TCLs is the duty cycle (D), defined as the runtime ratio that a unit is on (t on ) within the period (T) of a cycle (D ¼ t on /T). Even though each TCL is either fully on or fully off, the probability to be on in random moment is equal to its duty cycle. The phase of the cycling of a TCL relative to another TCL is random, so for a large population the aggregate of all the individuals' demand tends to be independent of time. Given no external interactions or changes, the quasiequilibrium condition is called the natural diversity of a cycling load. Thus based on these assumptions, Equation (9) can be written as: where D and p are the mean duty cycle and mean power rating respectively. The second part of (9) is the passive consumption of TCLs. In a population where duty cycle changes in time, p idle;i can be perceived as a constant base load and by introducing (13), (9) or (12) become (14) where P idle the sum of P idle;i : That base consumption is not affected by DR and thus can be disregarded in DR models as long as P i is used instead of P on;i . Lastly a similar approach can be taken for a dynamic heterogeneous population. The first step is to create clusters of TCLs in time. Assume a number Z of clusters C j ðtÞ; whose population (z j ) and representative means are a function of time (15), then total consumption is approximated by (16), the larger the Z and the smaller the time steps the more accurate, but we also want to keep small number of clusters for computation efficiency but also larger population ðm j Þ; thus higher statistical accuracy.
As we can see, the population's demand is approximately linearly dependent on the duty cycle. Upon synchronization (e.g. a DR signal or power outage) this is no longer the case, synchronization causes units to operate simultaneously in phase and natural diversity will be re-established slowly in time due to heterogeneity and random factors or through corrective actions [21]. This synchronization causes the rebound effect and its magnitude is based on the number of synchronized units.

Non-linear consumption factors
The aim of this paper is to study the cause of the aggregated behaviour of TCLs and its effect on aggregated demand. Using main driving factors, create pragmatic models, which can be used to simulate real world scenarios of DSM and DR applications. A realistic model was created to incorporate the most important dynamics while also focusing on simplicity when possible, then validated against real world data from DECC [28] and Smart-A [29]. It was based on thorough and extensive research in experimental data of TCLs' model parameters [5,21], mainly on cold loads as a case study [8,20,34,35], thermal properties [30,31], human behaviour [30,31], and simulations. Additionally, on Laguerre's et al. thorough research which focuses on cold loads [36e38]. The model was created such that it may be used for either cold/cooling loads or heating loads. The governing Stochastic PDE (Partial Differential Equations) used can be applied to all normal thermal loads, with only the input changing. Cold loads were used for the base case study, being the best example but also the most common among TCLs, with the highest demand and potential for DR. Important TCLs consumption factors and human in the loop: a) Ambient temperature, a function of human behaviour and weather b) Appliance characteristics and operation c) Human behaviour, preferences and socio-economic factors

Ambient temperature
Thermal loads used by the industrial, commercial and residential sector make use of heat transfer mechanisms; advection, conduction, convection and radiation. For relatively small temperature changes, such as in the case of TCLs (domestic and commercial), Newton's law of cooling applies. Therefore, heat transfer and consecutively the thermal load is practically linearly dependent on the temperature difference, dðqðtÞÞ dt ¼ Àl$½qðtÞ À q a ðtÞ þ mðtÞ$Q g (17) where m(t) defines the operational state of a TCL (1 for ON state and 0 for OFF/idle), l constant (Newton's heat transfer coefficient), Q g is the heat/cooling transferred during operation (ON state). As expected TCLs' energy consumption depends mainly on the ambient temperature and any change affects them directly, as seen in (7) & (8), cooling and heating rates are linearly dependent of Dq (q set -q a ). Refrigeration demand during winter is about 2/3 of the one in summer, while there is also a deviation between daytime and nigh time [28,34]. Previous work on TCL modelling and simulations has assumed steady or relatively steady ambient temperature around a set value (i.e. [20e27]). In Ref. [39] variations of 10 C in all house types during heating periods were found, which highlights the significance of q a in heterogeneity and deviation of heating practices. Zehir et al. [8], used real data to study cold loads' demand side management, where room temperature was also monitored, showing changes of 3 C in a few minutes. Mean temperature during the day can be seen in Fig. 4, based on measurements over nine months for a UK city-wide survey of over 500 homes [39].
These are used as a reference of maximum, minimum and average temperature in the model developed. As discussed previously periods where there are ambient temperature changes for a part of the population are of significance to the accuracy of the load models. An important note here is how equations (7)e (10) and (17) are linked. The thermal load or rather loss of it, is expressed by (17), which defines the cooling and heating duration (also rates in a population) and is linearly dependent of DT. Cooling and heating rates are also directly linked to duty cycle (D ¼ c/r þ c) [40]. Equation (18) is a natural outcome, since when the thermal demand changes power rating doesn't, what does is duty cycle, proportionally to the change of the thermal load. This is rather useful since it applies to aggregated populations too, thus knowledge of D(t) (t on , t idle ) through P(t) can be used to estimate the state of a TCL, or with multiple reading across the whole population, a principle used in Ref. [16]. In Masjuki et al.'s measurements [30], energy consumption increased from 0.56 kWh/day to 1.12 kWh/day when the ambient temperature was raised from 16 C to 31 C, a 100% overall increase or equivalently around 37.3 Wh/day for a 1 C increase in temperature. Hasanuzzaman et al. [31] has reported an increase in consumption from 1.2 kWh/day to 1.7 kWh/day, when ambient temperature was raised from 18 C to 30 C, 41.66% increase or 46 Wh/day per 1 C. Therefore, for changes in short periods, such effects on state (D(t)) should be modelled.

TCL characteristics and operation
During the operation of a TCL around the set-point value q set , there are q min and q max values which are the switching (ON/OFF) points. Excluding external interference, thermal losses are defined by the thermal characteristics of the appliance, namely thermal capacitance C, resistance R, efficiency h and also the difference of ambient q a and set-point temperature q set (Dq as discussed above).
Consequently, set-point temperature is of equal importance to ambient temperature. Set point actuation control algorithms are based on the above fact (thermal load f Dq) [6,21,23]. The effect, as expected, is essential the same as ambient q a changes [30].

Human behaviour, preferences and socio-economic factors
Human behaviour, preferences and socio-economic factors vary across consumers and TCL end function. They are multi-variable dependent and stochastic, posing the biggest challenge to model, especially because they introduce heterogeneity in operation even among identical appliances, thus making a homogeneous population behave heterogeneously in time, directly affecting the accuracy of models and control actions. Yet for large populations, statistical approaches are fit for such tasks, especially given the fact humans are "creatures of habit", thus proper examination can lead to appropriate clustering and control frameworks. Each of these factors affect each load type differently and thus the examination breaks down per load type: space heating & cooling, water heating, cold loads.
Space heating/cooling loads operation is essentially a combination of weather and human behaviour. Human behaviour is rather complicated, being a function of three connecting variables; time, comfort zones (conditions, preferences) and socio-economic factors. For instance, a household of employed individuals, working office hours 9e5, is statistically less likely to use heating/cooling during office hours (individuals not present) on a given working day [32], Fig. 3 shows the difference in demand. Once a household is active (individuals present), heating/cooling might be used; it depends primarily on current conditions and individuals' preference of "comfort zones" but also socio-economic factors. De Cian et al. [41] examined the interaction between income, temperature and energy demand, where an income interaction model was created, examining the income/temperature elasticity of electricity demand. Additionally, in Kane T. et al.'s [39] work (real world measurements), the rebound effect was greater than expected, which is attributed to the above socio-economic behaviour. Thus geographical clustering for thermal loads, such as in Ref. [23] where relative homogeneity is assumed (air conditioners) and similar operation characteristics (q a ) is inadequate, especially for TCLs such as colds loads. Water heating is similar, but less reliant on weather conditions or comfort zones, rather based on preferences and habits instead.
The major electric TCLs consumers, cold loads, have one external affecting factor, human behaviour. It can be broken down to door opening and loading the compartments. Experimental tests, using ISO standards [30], with door opening of 12 s at a 90 angle, report an increase in consumption from 0.85 kWh/day to 1.42 kWh/day, for 75 such events, 7.6 Wh (or 0.894%) increase per event. In practise, during door opening the insulation (and therefore overall heat transfer coefficient) change drastically, warm air mixes with cool air inside, heat transfer occurs through convection. A note at this point is that such events are dependent on Dq; it is the same Dq during normal operation and during the event, with the change of overall heat transfer coefficient (h), thus the proportional increase (0.894%) can be assumed relatively constant given not high deviation of conditions. This can be expressed as: Experimental studies on refrigeration have actually showed that new load can have the greatest impact in consumption, especially in short term [30,31,37]. Zehir et al. [8] mention that these effects need to be considered during simulations, yet they are hard to model. In Masjuki et al.'s measurements [30], energy consumption increased from 0.96 kWh/day to almost 2.3 kWh/day with 18 kg of water added (room temperature), though not linearly. Until about 9 kg the increase was linear, rate of about 37.5 Wh/kg, an increase of 3.9% per kg. Hasanuzzaman et al. [31] has reported an increase in consumption from 1.2 kWh/day to 1.9 kWh/day, with 12 kg of water added (room temperature), about 58.3 Wh/kg, increase of 4.83% per kg. Taking these findings into account, as well as the fact that those are not randomly distributed during the day (Fig. 5), it is obvious that a Wiener process (mean 0) with small variance (s) doesn't properly model such effects or their impact in consumption and the TCLs populations' dynamic nature in time.
Moreover, (7), (8) do not include any of the above (q a , human interaction), since ambient temperature is assumed almost constant and the Wiener process isn't included due to having 0 mean and practically insignificant deviation. CFPE should have included this dynamic behaviour (heterogeneity) in time, thus some error is expected.

Population size and stochasticity
When dealing with stochastic approaches, population size is detrimental for statistical accuracy. Especially for realistic studies with heterogeneous populations, the number should be large enough to be able to model real world dynamics and have statistical validity. At the same time, we want to have the smallest possible acceptable size for computational efficiency. Relatively small population sizes might be also more appropriate to model and study commercial units, which have considerably higher demand than domestic loads but are also relatively limited in numbers (i.e. supermarkets refrigerators). In such a case bottom up models (such as Monte Carlo) are appropriate, others such as CFPE, Euler-Maruyama approximation etc. yield statistical errors, since the distribution of the population is Binomial.
At this point we will examine the model's population size. For highly heterogeneous populations, subject to external factors, it is hard to define a minimum acceptable population size. Yet, for the case of homogeneous, free of external factors TCL population, minimum population size can be calculated simply and efficiently through confidence intervals (CI). As previously mentioned, estimating available power in time is directly linked to estimation of the state of the population in that given time. It is thus essential to have a metric of the state in time, as well as a metric of its validity. The parameter that describes the state of an individual TCL (on/off probability) is its duty cycle D 2 (0,1). Based on above assumptions (homogeneous population, no external factors), this falls under a Bernoulli process, thus the Confidence Interval is directly linked to the duty cycle and the size of the population.
The probability of each value x of a Binomial distributed random variable X is defined through its probability mass function: Expected value: The above also holds true for heterogeneous populations when clustered in M relatively homogeneous clusters, where D is replaced by D i ; i 2 (1, M). Taking into consideration the operational heterogeneity (D function of time for each individual), then clusters cannot be defined in a classic way, but rather dynamically in time blocks. In which case D(t) can be written for time blocks b ¼ f(t), as D(b) and the above becomes: The larger the population the closer to the expected value in a random moment, thus "random/uniform" distribution can be assumed with a small enough error. Confidence Intervals can be used as a metric; adjusted Wald, Wilson-Score, and exact Clopper-Pearson methods were considered. For large populations (n¼> 1000) and .1 < p < .9, any of those 3 methods give practically the same results. For smaller populations the exact method is preferable.
The above values of duty cycle were selected since in reality cold loads' duty cycle fluctuates usually within this region. We can see that for populations (clusters) around 1,000, stochasticity is significant and their expected state wouldn't be as statistically accurate. The fluctuation introduced in the model due to being a binomial distribution will affect the accuracy of assessing the impact of the factors described above. For a population size of 100,000 it is well within limits and preferable, yet computationally slow. A model with size of 10,000 has acceptable accuracy to study TCLs' dynamics, whilst remaining computationally efficient. In a heterogeneous population, the above values are expected to change slightly. This can also serve as a metric for clusters' size (minimum requirements for statistical accuracy).

Physically-based model of a single TCL
A TCL can be described by _ qðtÞ ¼ Àl$½qðtÞ À q a ðtÞ þ mðtÞ$q g þ uðtÞ$q e (21) Note that (21) is different to the classic one, as q a (t) is not constant, and u(t), q e represent interactions due to human behaviour, which in theory can be either positive or negative, most commonly in practise though, for cold loads positive, for water heating negative (use of water), for space heating either and space cooling positive. There is a definite difference between m(t) and u(t) since u(t) is stochastic.
We introduce the parameter q g ¼ ð1=lÞq H and q g ¼ Àð1=lÞq C for heating and cooling accordingly, to simplify the above equations. q g can be perceived as the temperature gain (or loss for cooling), to which the system tends to during the on state (and would eventually approach for a prolong t on '[t on ). At that point dðqðtÞÞ dt $ 0/qðt 0 on Þyq g þ q a : Similarly, q e ¼ ð1=lÞq e . Equations (22a), (22c) can be written as:

Discrete solution
Solving (21) for discrete time steps (n to t) and using the same notations as above for q g ; we get: 3.1.3. At this point we will examine a fit for q a ðtÞ as a function of time For space heating/cooling (heat pumps, air conditioners, electric space heating, gas heating), ambient temperature is perceived as the outdoors, which has relatively small variations from hour to hour and very small on a minute basis (the simulation time step). In reality the change in external temperature depends on many factors and the best fit between small time steps is polynomial of small order or linear. For cold loads and water heaters, ambient temperature is indoors, which can have more sudden changes, such as early morning or when returning home and heat is turned on. The temperature increase in this case, based on experimental data seems to also follow polynomial or linear (Fig. 4 [39]). We will solve (23) and (24), difference is minimal thus linear is preferred for simplicity, if there is any case where Dq a is significant (maybe in an industrial environment), the exponential one is a better fit though probably still of small importance. A note here is that the model is developed for a simulation time step of 1min, with temperature readings every 1 h, q a (t) is thus fit to 1min steps through these.
The term ½l À 1 þ e Àl =l defines the importance of a; which is the rate of change of the ambient temperature. This rate also depends on the time step, for larger time steps of q a (i.e. 10min instead of 1min), a has a value 10 times higher. HðnÞ defines the effect of human behaviour, which as examined earlier has a considerable effect in cold loads consumption. It can be modelled to include a Wiener process though not mandatory when Monte Carlo simulation is used, since Monte Carlo introduces stochasticity. In this case a Wiener Process was deemed unnecessary since Monte Carlo introduced higher stochasticity and it would only be extra computational cost.
The model, similarly to (1), is composed of two interconnected subsystems, a linear part characterized by a stochastic continuous state q(t), whose evolution depends on u(t) and a nonlinear part u(t) which changes states deterministically based on q(t). This vector process ½qðtÞ mðtÞ T is a hybrid-state Markov process.
3.1.5. Constant q a ðtÞ ¼ q a ¼ const, continuous and discrete time solutions Practically these solutions are special cases of the above for q a ðtÞ ¼ q a ¼ cons: These can be used for simplicity when q a is relatively constant or can be assumed as such based on the application.

Equivalent thermal models for multi-compartment TCLs
Fridge-freezers are the most common cold loads in households (about 69.7% ownership in 2014), whilst also being the cold load with the largest average consumption per unit [28] and are found in commercial buildings too. ISO 8187, ISO 8561, and ISO 7371 are the relevant standards for testing the energy consumption of household refrigeratorefreezers having two or more compartments. At least one compartment (the fresh food storage compartment) is suitable for storing unfrozen food, and at least one compartment (the food freezer compartment) is suitable for freezing fresh food and for the storage of frozen food at À18 C or lower [30,34]. They may be equipped with one or two compressors. In the case of one compressor, the operating cycle is controlled by both the refrigerator's and freezer's air temperature, while commonly a damper or fan is used to assist heat transfer from the refrigerator to the freezer. In the case of two compressors, each operating cycle is independent. The temperature in each compartment is better regulated but the price is higher, thus less common.
Deterministic equations governing temperature and operation state in each compartment are given by equation (30) for the case of one compressor and by equation (31) The main problem with modelling those as in [24], is the fact that they consist of 2 compartments with different set point temperatures and temperature dead-bands and as consequently different thermal dynamics (Dq compartment-room ). Therefore, the effects of q a change or human interaction are different. In both cases the operating cycles are determined by the compound dynamics of the system, and an equivalent (one compartment) simplified model of the system needs to be developed for computational purposes, to avoid slow simulation but also enable the use of a model as designed above. This becomes clear by looking (30), (31) or (32), (33), which after adding the stochastic elements of cold loads would be much harder to compute.
Note that in the case of two compressors on/idle states might differ, which means that the power consumption of the two differs in time. Yet, for an aggregated model of thousands of units, the percentage of units on is still directly linked to the duty (ή operating) cycle, consequently control of aggregated loads has statistically the same result.
The total thermal load losses of a refrigerator-freezer in static state (minimal to zero disturbance from external factors) is the sum of the freezer's and refrigerator's load, can be calculated as: where _ Q f heat transfer from room to freezer, _ Q f Àr heat transfer from refrigerator to freezer and _ Q r heat transfer from room to refrigerator. The heat transfer can be approximated as _ Q ¼ hAðÀDqÞ; where h is the heat transfer coefficient and A the surface respectively. The above equation can be written as: The equivalent single compartment model can be described by h eqÀroom A eqÀroom ðÀDq eqÀroom Þ; where The objective now is to find the q eq . Where w weight factors of compartment capacity to total capacity.
The two compartments share two dimensions, meaning that capacities are proportional to A; thus: Dq ref Àroom (38) Heat transfer coefficients are the sum of the internal cabinet coefficient (h fr , h ref ) and the appliance's outside surface coefficient with the room (h room ). Studies have shown a difference less than 9% between freezer's and refrigerator's heat transfer coefficients (convective plus radiative), by adding the surface's coefficient the sums' difference is less than 4% (3.7%) [31]. The heat transfer coefficients can be set as: For simplicity h f Àroom ; h rÀroom will be written as h f , h r for the following equations.
Using (36), (38) & (39) on (35) gives: It is important to note that the aim of this equivalent model is to simulate the effects of external factors on energy consumption and the non-linear aggregated demand of cold loads. Typical values are given in Table 3, for different types of fridges. Approximations are used for simplicity and computational efficiency; in reality thermal laws have a different effect. As we can see in Fig. 6, the heat transfer is essentially always from the room to the freezer, directly or indirectly (through refrigerator), thus for normal operation the thermal behaviour is practically closer to the freezer's thermal properties. Yet human interaction with the refrigerator, i.e. door opening, extra load etc., will have an effect (heat transfer) closer to the refrigerator's thermal properties.
An examination of the above can be done with the experimental measurements in [31]. Using (18), in this case duty cycle change  from 0.303 to 0.429, or dðDðtÞÞ Dðt0Þ ¼ 0:4158; dðDðtÞÞ Dðt0Þ z dðPðtÞÞ Pðt0Þ ¼ 0:4166, a difference less than 0.2%. Table 4 shows the relative consumption (P) and temperature (q) increase.

Aggregation methodology and heterogeneity
Even though the simulations of CFP partial differential equations are much faster than the Monte Carlo ones and the accuracy satisfactory (under conditions), Monte Carlo is a bottom-up detailed model, better suitable in assessing the above developed model, where not all CFPE's conditions are met. Additionally, the aggregation process for CFP partial differential equations is more difficult and the complexity is higher, while also losing accuracy in heterogeneous models. The major advantage of Monte Carlo-based techniques is flexibility, as well as the highest possible accuracy in this case of highly heterogeneous populations with human behaviour, but also effective with smaller population sizes.
A note here is that defrost heater power and cycles are ignored for simplicity. Defrost heaters are typically on for less than 5% of the time [8,24], yet their demand during this time is usually considerably higher than that of the cold load's operation (480 W [8]). They *data inferred or assumed based on similar experiments and loads **calculated from consumption and power rating. a Laguerre's experimental data in [37], but not referring exact ambient temperature as in [36,38]. Similar conditions are expected since it's the same climate and environment as the others. b Freezer with separate compressor running its own cycle [36]. Table 6 Cold load l, q g (calculations).

4.1.
Step 1: calculation of basic TCLs' parameters using experimental data In [20] and [25], parameters such as appliance's thermal conductance, thermal mass, power consumption, coefficient of performance were used to determine the values l, q g (Q g ) of (14).
Those are constant, dependent on appliances' characteristics and independent of ambient temperature and human interaction.
Hence simply knowing l, q g is enough to build (14), an accurate representation of a TCL's physical model, a principle used in [16] where an algorithm was developed to identify them through readings of duty cycle (D).
Using values from Table 5, l, q g can be calculated. During a deterministic cycle (without external effects) and constant (q a ) ambient temperature, a TCLs idle cycle is described by For t ¼ t idle : qðt idle Þ ¼ ð1 À e Àlt idle Þq a þ q 0 $e Àlt idle : Case of cooling cycle, qðt idle Þ ¼ q þ ; q 0 ¼ q À ; substituting on (41) and solving for l gives: Similarly; q g ¼ h q À À q þ $e Àl$ton .
1 À e Àl$ton i À q a The values of [8] have significant variations compared to the rest; as stated by the authors, data was taken from a G class top mounted refrigerator, and was chosen as a good example for a low efficient refrigerator with high energy demand (relatively older model). According to Residential Energy Consumption Survey 2009, in USA, this was one of the most common (more than 60% percent) refrigerators in households, thus chosen by the authors for their study on 2009, also mentioning that customers mostly use the same appliance for more than 5 years up to 14 years. Newer models (post 2009) though are more efficient and more similar to those used in [20]. Similarly, for fridges (refrigerator-freezer) some of the above have significantly different conditions than what commonly observed in UK (i.e. 33 C room temperature) or missing data. For the case of UK specifically, values from [20,31,36,38] are more suitable for a realistic model.
These are used via Monte Carlo to construct heterogeneous populations of TCLs (3 basic types of cold loads in this case) with randomized variance in parameters q g , l, q a , and initial conditions q (t), u(t); u(t) based on on/off probability which is equal to the duty cycle D of each. The population's mixture, as a percentage of each cold load type, is based on UK ownership statistics [28]. Specifically, for UK ownership of cold loads per household is about 107% (expressed to overall population), which means that 7% have an extra unit.

Step 2. Human interaction
The effect of human interaction as discussed in section 2 is relative to Dq similarly to ambient thus better expressed in (%). We need to link that to the term H(b) of (20b), which varies for each type of load. To the best of the authors knowledge this is the first realistic human interaction modelling based on experimental data, such as measured effects of door openings or increased consumption per kg of extra load, combined with real data of frequency of human interaction. The exact procedure, in more detail can be seen in Fig. 8. Effectively, human interaction increases t on or reduces t idle . Using the fact that increase in consumption is linked to increase in duty cycle dðDðtÞÞ Dðt0Þ z dðPðtÞÞ Pðt0Þ (18), we can express an increase in consumption by x (%) as an equal increase in duty cycle: D 0 ¼ Dð1 þ xÞ: Simulation results and analysis

Model Flowchart
Step 1 Step 2 Step similarly during ton:H ¼ À q À À q a À q g Á 1 À q À À q a À q g q þ À q a À q g ε ! (50) Using the above and experimental data described in Section 2 [30,31], human behaviour effect per event for each cold load type is calculated, separately for door opening and new load. Then using the frequency of door-openings per day, as described in [29] (see Fig. 5 in Section 2) an hourly PMF (probability mass function) is created, and a Gaussian distribution for events per day with expected mean 40e60 events and new load per day [30]. Those are used to create a Monte Carlo model with minute resolution for events per household throughout the day (Figs. 9e11). Events are based to reflect experimental ones, such as in [30], "door remained open for 12 s at an angle of 90 ", but since one such event might differ in reality (i.e. longer duration than 12 s), non-integer values were used (i.e. reflect the extra effect of an event with longer duration). Monte Carlo events are linked to the term Hb through (36), (37), (39) & (40) to express the effect in cold loads' state, consequently cycles and consumption.

Step 3. Ambient (room) temperature
Temperature fluctuations as discussed in Section 2 affect the state of TCLs directly, especially rather fast changes, such as those occurring in early morning and evening after work. Kane T. et al. [39] examined the variation of indoor temperatures and heating practices in UK dwellings, variation of up to 10 C were reported, we use UK heating practise (probability of human behaviour essentially, as developed in model [39]), reported "preferred" comfort zones and outdoors temperature to create a realistic MC model. The heating PMF is used to create a Monte Carlo model of possible heating on action, which is then compared to probable indoor temperature (based on outdoor temperature) and thus hourly temperature is defined per dwelling. In reality, based on knowledge of weather conditions and historic data (population's behaviour to these conditions) the temperature can be probabilistically estimated. The resulting MC model is then converted to minute resolution. It could also be used in hourly resolution, with the term of (28c) a converted to hourly changes (in this case would be 60 times higher). Temperature changes in dwellings due to use of heating can occur in shorter periods than an hour [8], but for this model more "mild" Human interaction (PMF) Step 2: Human interaction Step 1 Step 2 Step 3 Hb per load type (per event) Events per day (Gaussian) MC (events, minute resolution)  fluctuations were used to describe the general behaviour of the population. Simulated average household temperature and reported average temperature ( [39]) can be seen in Fig. 12, as expected they are very close. Assuming that this report is a good representation of UK's heating practises overall, the cold load demand is also expected relatively close to the actual report in [28].

Simulation results
One important aspect of a model described by (28c) is that it can focus on the state of the TCLs specifically without the power demand of each. Thus, it can give a clear picture of TCLs' state and how they respond to different commands as well as other factors. Demand can be added through (9) to determine the total power demand of the TCLs' population.
The model can be used as classic models (Malham e and Chong), by assuming q a constant and minimal external factors (Wiener process with mean 0 and small variance s). In which case the heterogeneity of the population due to technical parameters can be studied and the effects of stochasticity. Figs. 13 and 14 show the difference between population sizes of 10,000 and 20,000. As expected the stochastic behaviour reduces as the population size increases, the values (fluctuation) are close but higher than those in Tables 1 and 2, due to heterogeneity. Populations of smaller size will have significant elements of stochasticity (i.e. 1000) and thus are not proposed for DR studies of TCLs. Unless the study is specifically aimed for small populations, then one of the most important elements to be studied is stochasticity and control actions need to take it into account. Increasing the population size from 10,000 to 20,000 has a small improvement, in this study a population of 10,000 was deemed satisfactory for examining effects of human interaction and ambient temperature. Assumptions of almost uniform distribution (around 1% fluctuation or less) for probabilistic models in large heterogeneous populations, require population sizes at the order of 100,000. Detailed bottom up models, like this one has the highest fidelity study these behaviour characteristics. Arguably, one of the other important uses of this model is the realistic representation of TCLs under different conditions during the day and their dynamic behaviour in time. Especially when considering the human factor. Fig. 15 displays simulation results for varying ambient (room) temperature and human behaviour, as described in steps 2 and 3. It is also compared to experimental data from DECC and Smart-A as seen in Fig. 16 (converted to hourly means and experimental data scaled for comparison purposes). The difference with the "state of     the art" TCL model can be observed in Fig. 17.
Other important characteristics and phenomena that can be studied in detail with this model are the behaviour in external commands during different conditions (periods of the day), rebound effects and oscillation damping after those, due to partial synchronization ( Figs. 18 and 19). The exact response of a DR actions (or other forms of interruptions such as brownouts), as shown below, are directly affected by TCLs aggregated state at that given period. Additionally, as these cause part of the TCL population to synchronize, rebound effects occur, which are detrimental for DR actions. They can be significant for dispatchable and nondispatchable DR. For instance, when a certain amount of power is required for frequency control, TCLs can provide it by reducing demand output equal to that required amount, but after some time when the rebound effect is introduced an almost equal extra amount of energy is required, which if not countered properly by additional control actions could practically recreate the original problem [20].
It is thus essential for control actions to take such effects into consideration and try to minimize/counter them or at least postpone them, as well as consider the natural damping occurring due to heterogeneity [21]. It is crucial to note that, in non-intrusive actions (thermal limits are maintained), the overall thermal load remains essentially the same. Consecutively, the electrical demand is also the same, but allocated in time, unless the control framework is designed otherwise.
One of the starting arguments of this paper and motivation for this work was whether small heterogeneity can be assumed for TCLs, such as cold loads, and the accuracy of models based on such assumption. Fig. 20 shows the variation of duty cycles during the day for 4 TCLs with relatively similar starting duty cycle (and thus operation), yet there is significant change during the day. This is true for the majority of the population, and the original point, that even in a (relatively) homogeneous population there will be significant heterogeneity in operation in time holds true, as well as that Wiener process with small variance is probably inadequate for representing it. This model can be used for detailed DR studies of TCLS, rebound effects, efficiency of control algorithms being a bottom up MC realistic model. Moreover, it is realistic point of reference to check the feasibility of models with less computational requirements (such as CFPE models) and their accuracy, which will be used in aggregators [21]. It is also important to note that response to control actions varies during the day (i.e. nightmorning hours have relatively shorter cycles than evening), based on the state of the TCLs' population, also a factor in identifying the    available power for DR.

Conclusions and future work
TCLs modelling was examined, and a detailed bottom-up realistic model was developed for DR studies, including human interaction modelling and its effect. Controlling large amounts of loads is inevitably only possible through probabilistic models, but also accurate enough. One of the keys issues for TCLs' applications is the high heterogeneity and the dynamic behaviour in time. Future work will focus on exploring solutions to this problem, potentially clustering frameworks as well as the impact of DR. Other important topics are computationally efficient aggregation models, control frameworks and data requirements for those. This paper focuses on data from experimental studies to determine the most important factors in TCLs' dynamic behaviour. Based on those and the developed model, simulations were carried out for the case of colds loads. The initial hypothesis of high heterogeneity in operation even among relatively homogeneous loads was shown to hold true, thus the accuracy of various existing aggregation models is questionable.
Important points and outcomes can be summarized in 6 points (with DR significance order). First, showing the dynamic nature of TCLs and how their duty cycle (operation) changes during the day significantly, thus the problematic nature of using aggregation models such as CFPE by assumptions of relative heterogeneous groups. This is due to external factors, essentially human behaviour, which were modelled and simulated using real world data to highlight their actual effects in TCLs' consumption, subsequently also in TCLs' population state and thus DR. Additionally a simplified equivalent thermal model of multi-compartment cold loads was developed, which converts the actual second order thermal model to a first order one, thus useable for the S-PDE models developed so far for DR without loss of realistic TCLs characteristics. The bottom up approach used, can be computationally demanding due to its detail, yet that detail allows for high flexibility and the highest possible accuracy to study DR actions, rebound effects, test the performance of control algorithms. Additionally, it can be used to evaluate the accuracy of more computationally efficient aggregation models and simulate smaller populations taking into consideration the magnitude of stochasticity in such cases. Lastly, due to the main equations including external factors, minimum and maximum limits of "normal" operation can be deduced, as such abnormal (faulty) TCL consumption can be tracked, something useful mainly for commercial loads. So far this has mainly been done through regression models, but it can be used in conjunction to improve accuracy or to estimate accurately consumption under different operational conditions for energy efficient.