Demand response based desynchronization of thermostatically controlled loads with hardware implementation

Smart buildings are the need of the day with increasing demand-supply ratios and considerable deficiency to generate them. In any modern nonindustrial infrastructure, these demands mainly comprise thermostatically controlled loads (TCLs), which can be maneuvered. TCLs, such as air-conditioners, heaters, and refrigerators, are ubiquitous and their operating times can be controlled to achieve desired aggregate power. This power aggregation, in turn, helps achieve load management targets and thereby serves as an ancillary service (AS) to the power grid. In this work, a distributed averaging protocol is used to achieve the desired power aggregate set by the utility using steady-state desynchronization. The results are verified using a computer program for a homogeneous and heterogeneous population of TCLs. Furthermore, a load following scenario has been implemented using the utility as a reference. Apart from providing a significant AS to the power grid, the proposed idea also helps reduce the amplitude of power system oscillations imparted by the TCLs. Hardware-based results are obtained to verify its implementation feasibility in real-time.


I. INTRODUCTION
Sustainability is the process of maintaining change in a balanced environment, in which the exploitation of resources, the direction of investments, the orientation of technological development, and institutional changes are all in harmony and enhance both current and future potential to meet human needs and aspirations. With the advent of civilization, to serve electrical and transportation needs, fossil fuel reserves have been overutilized. Excessive use of fossil fuels and increasing human settlements have brought imbalance to the environment. Thus, researchers have moved towards building settlements that can serve individual requirements, thereby meaning no harm to the environment utilizing the least natural resources, and, thus, are sustainable. To support the central power grid, ancillary services (ASs), a bidirectional power flow architecture, are employed.
ASs support the transmission of electric power from seller to purchaser, given the obligation of control areas and transmitting utilities within the control area, maintaining reliable operation of the interconnected transmission system. 1 Load following, reactive power voltage regulation, system protection services, loss compensation, system control, load dispatch services, and energy imbalance services are some ASs. Load following is a service program that monitors real-time consumption about scheduled loads to consume what is estimated. In contrast, load dispatch services take real-time action in the case of any discrepancies. System protective services are energy reserve that compensates for any discrepancies in the scheduled loads and dispatches where the loss compensation scheme guides its use. The energy imbalance programs take all these actions using system control service. Reactive power voltage regulation refers to services that mitigate perturbation in the voltage online. As observed, these services either fall into demand-side or supply-side management, but with increasing demand to supply ratio, it is advisable to optimize loads to match available resources.
ASs have gained attention with increasing penetration of renewable energy sources (RESs) in the course of replacing fossil ARTICLE scitation.org/journal/adv fuels. Although renewable energy sources (RESs) are abundantly available in nature, they suffer from intermittencies. These intermittencies are significant enough to render a power grid unstable and hence require additional regulating circuits to avoid direct injection. [2][3][4][5] One of the ways of injection is to store it in a battery and then provide it to the grid using a battery management system (BMS) with a stable discharge rate. Although reliable, this strategy bears a high capital expenditure for installation and hence is commercially less feasible. [6][7][8][9] Therefore, research has moved in the direction of devising methods to regulate the loads with techno-commercial feasibility. Demand Response (DR) programs have been developed to ally with several handshaking criteria and thereby maintain power grid stability. DR programs allow consumers to adjust their electricity consumption in response to energy prices or incentive payments. ASs form a major part of these DR programs, and optimal use of these services can come to an advantage in terms of adjusting loads as per utility demands. As per a survey in the U.S., 10 cooling and heating loads account for the largest annual residential electricity consumption. In addition, the scenario remains the same globally, thereby requiring optimization of these loads to attain required power aggregation. Fundamentally, all heating/cooling loads [or thermostatically controlled loads (TCLs)] comprise a temperature transducer, a thermostat, and other assemblies such as compressor and condenser. It can be noted, though, that TCLs follow a scheme whereby they switch between ON and OFF states to attain required setpoint temperatures. ON states consume full capacity power of the TCL and OFF states imply dormancy. Thus, adjusting these operation states can help achieve the required power aggregate and thus provide load-following ancillary service.
In the literature, 11,12 field experiments are conducted using domestic household refrigerators to quantify the flexibility of TCLs by advancing or delaying operation times, thereby achieving the required aggregate power. To avoid synchronization, parameter randomization of TCLs by compartmentalization of compressor cycles is also an effective scheme. 9 An automatic generation control signal (AGC -a demand-supply balancing technique) tracking 13 and a twolevel scheduling protocol 14 achieve flexible TCL operating cycles in an intraday electricity market, making it a valuable asset to DR programs. However, with increasing heterogeneity in system parameters, the performance worsens in all the schemes mentioned above. On the other hand, a multiobjective model predictive control (MPC) can also be employed for residential heating (or using heat as an energy carrier) with heat pumps. 15,16 Again, the number of parameters used makes it difficult to be implemented for a large population of TCLs.
Aggregated models for TCLs can carry information about the overall power consumption of TCLs and can be framed using discrete state-space dynamics. For instance, a centralized control scheme on an aggregated model can provide effective frequency regulation. 16,17 A Markovian probabilistic framework 18 or a load shedding strategy 19 can also be used to achieve similar objectives. Concepts from optimization theory related to the tracking device, convex polytopes, and binary multiswarm, and particle optimization, are also used to control the desynchronization of TCLs. [20][21][22] Batch reinforcement learning based on Monte Carlo methods is used to allow TCLs to participate in day-ahead markets, thereby achieving DR objectives. 23 A nanogrid solution to store thermal energy in TCLs for later use has also been experimented. 24 This approach avoids wastage by regulating maximum power point tracking (MPPT) of a photovoltaic panel. In the literature, 25 a partial differential equation (PDE) based Fokker-Plank model for demand-side management is also proposed and, as known, it requires highly sophisticated tools for solving. Additionally, an event-triggered control of TCLs on a hierarchical architecture is implemented over a network for optimization of communication cost. 26 All of the methods mentioned above are probabilistic and discretized, and they involve considerable computation cost.
The complexity of the population of TCLs increases with the number of TCLs and with the reference aggregate power set by the utility. The undesired synchronization of TCL duty cycles by the control policies limits the time scale and capacity of ASs that a population can provide. By sequentially changing the set point temperature, safe protocols have been developed to minimize the unwanted power oscillations. 27,28 Other ways of providing power regulation include variable deadband strategy, 27 injecting a delay of M-minutes, 27 and deploying state stacking technique based on priority measures. 29,30 These measures are defined as the distance from the switching boundaries to reduction in aggregate power. TCLs in the ON stack closer to the switching boundary are turned ON first, for an increase in aggregate power; TCLs in the OFF stack closer to the switching boundary are turned OFF first and so on. On similar lines, a phase response curve (PRC) based technique effectively helps to manoeuvre an ensemble of TCLs by using a control input to modulate duty cycles as well as induce delay/advance, thereby changing the set point temperature. 31 The switching dynamics of TCLs are similar to oscillators, and the phase-dependent model has been realized to achieve the desired objective. A continuum model for studying the synchronization dynamics of TCLs to understand and mitigate potential risks of a decentralized response provider offers deeper insights into TCL dynamics. 32 All these methods are effective but involve high computation costs and challenges for practical implementation.
This study extends the Kuramoto phase oscillator model 33,38 used to control the phase delay/advance of the operating cycles of the TCLs, which helps to achieve permissible load requirements set by the utility. We provide an extensive study on the possibilities of practical implementation of the model discussed in the previous work. 33 Based on our findings, we propose a distributed averaging protocol to achieve similar (and in a better way) objectives obtained using the Kuramoto model. 33 Hardware-based outcomes, in conjunction with computer simulations, provide validity to the proposed idea.
This paper is organized as follows: in Sec. II, we lay down fundamentals to develop a platform to propose the mentioned idea. The motivation behind the proposed paper has been discussed briefly in Sec. III. A novel model for the control of TCLs, namely, a distributed averaging protocol, is discussed in Sec. IV. In Sec. V, we provide software and hardware outcomes using light-emitting diodes (LEDs) as acting TCLs, which simulate similar characteristics, thereby validating the proposed idea. Sec. VI presents a competitive benchmark, where we compare the proposed model with those available in the literature and provide its competencies using peer-to-peer reviews.

A. Behavior of TCLs
A TCL switches ON until it reaches the upper limit of the set threshold temperature, post which it turns OFF until reaching the lower limit of the threshold. Usually, TCLs employ a thermostat that measures the surrounding temperature and thereby generates switching signals to regulate the TCL. Typically, a hybrid model is employed to express the dynamics of TCL mathematically. Incorporating the thermal resistance, capacitance, energy transfer rate, and switching signal, the hybrid model describes the temperature dynamics. It can be stated as follows:Ṫ where Ta and P are the ambient temperature and average energy transfer rate of a TCL, respectively. The minimum and maximum temperature of TCLs are Tmin = (Ts − δ/2) and Tmax = (Ts + δ/2), respectively, with Ts and δ being the thermostat set point temperature and deadband, respectively. s(t) is the switching signal of the TCLs. Figure 1(a) depicts the temperature dynamics, where the TCLs oscillate between 19.5 ○ C and 20.5 ○ C, while Fig. 1(b) is phase portrait of the TCL obtained using hybrid model (1) and data from Table I. Thus, for N TCLs, the aggregate power consumption is given by where ηj is the coefficient of performance of jth TCL. For simplicity of calculation, we consider ηj = η = 1 for all TCLs.
Remark 1. All the parameters used henceforth are taken from Table I unless otherwise specified.

A. Coupled pendulums
To further understand the effect of TCL switching oscillations on aggregate power consumption, we correlate the periodic behavior of TCL switching with a setup of spring coupled pendulums. It is well known that the Kuramoto model can be well understood using a system of coupled oscillators and various modes of oscillations can be analyzed. Consider two pendulums coupled with a spring, as shown in Fig. 2, where mi are masses, ϕ 1 and ϕ 2 the phases, and K is the spring constant. For the sake of simplicity, assume both the pendulums to be identical in terms of system parameters and

ARTICLE
scitation.org/journal/adv dissipative losses to be zero. As known, the phases of coupled pendulums exhibit "in-phase" (phase difference = 0 radians), "anti-phase" (phase difference = π radians), and "quasi in-phase" (0 < phase difference < π radians) type oscillatory behavior. 34 Let all the phases measured in the anticlockwise direction be positive and assume a constant external force that drives these oscillations. These external forces can be impulsive or persistent in time. It can be seen analogous to the addition of an external delay to maintain the required steady-state behavior (see Fig. 2). Now, all external forces (delay = 0 radians) that keep aggregate switching amplitudes (a mathematical sum of phases) to its maximum can be considered to be oscillating in "in-phase" modes. All forces (delay = π radians) that keep aggregate switching amplitudes to its minimum can be considered to be oscillating in "antiphase" modes. All the other states oscillate between these two extrema can be termed "quasi in-phase" modes. We portray "inphase" and "anti-phase" oscillation modes pictorially, as shown in Fig. 2. Furthermore, using the concept of clustering, we show how these modes come to benefit in achieving the required aggregate power.

B. Boolean phase oscillator model for TCLs
In the previous work, 33 the author makes use of the Boolean phase oscillator model to incorporate the system dynamics and thereby control the aggregate power consumed (Pagg). Visualizing the TCLs as oscillators, 33 the author employs the Kuramoto framework to desynchronize TCLs. Beginning with the conventional Kuramoto equation, 33 to accommodate the digital switching signals, uses the Heaviside function, thereby imparting the Kuramoto framework the capability to couple digital signals. The standard Kuramoto equation can be expressed as follows: where ωi, ϕi, ε, and N are the natural frequency, phase of ith oscillator, gain and number of oscillators, respectively, ∀i ∈ Z + . The inherent coupling dynamics of (3) enables the oscillators to establish communication between all TCLs. As mentioned earlier, the author using the Heaviside function modifies (3) 33 and the resultant equation can be given aṡ where Θ is the Heaviside function given by and K is the coupling between oscillators. The phase separation αij in (4) induces a phase separation between the i th and j th oscillator, or in this case TCLs as shown in Fig. 3(a). The following equations characterize the temperature dynamics and stabilization of aggregate where TON,i is the ON time of the ith switching signal (mapped to its equivalent phase in radians), ωi is the frequency of switching, and K is the coupling constant. A bias, (si ,0 ), helps the generated signal map the duty cycle in the Kuramoto framework.

ARTICLE scitation.org/journal/adv
In the previous work, 33 the idea is further supported by simulation results. It was thus conclusively stated that the power oscillations and aggregation are minimum when the phase delay (αij) is 2π/N.
To delineate this theory, the author 33 takes a case of N = 4 homogeneous TCLs, with the natural frequencies as ωi = 0.55 Hz. Furthermore, the TCLs desynchronize at their respective phase differences with their common frequency being ω FFT = 0.298 Hz, where ω FFT is calculated using fast Fourier transform (FFT) of the switching signal. Finally, the Kuramoto model desynchronizes TCLs at a common frequency lesser than their mean frequency. Figure 3(c) shows the power vs time plot for N = 4 TCLs.
The author then takes a case of N = 100 heterogeneous TCLs with their natural frequencies in the range of ±5% of its nominal value. In Fig. 4, it is observed that the load fluctuations are minimized to within ±2% after applying the Kuramoto based model.

Remark 2.
Frequencies calculated in the previous work 33 are scaled down by hours, i.e., ωi is calculated as 1/T, where T is in hours. However, in this paper, we consider time in seconds.

C. Shortcomings of Boolean phase oscillator model
Although the Boolean phase oscillator model 33 provides motivating results, some shortcomings in the model proposed can be listed as follows: 1. It must be noted that the natural frequency ωi used in (6) can be easily calculated using mathematical manipulations as follows: However, it is observed that as the number of TCLs increases, the value of ωi needs correction. For instance, N = 16 requires ωi = 8.81 Hz and N = 100 needs ωi = 62.18 Hz, to achieve ω FFT = 0.298 Hz. Moreover, as the number of TCLs increases, ωi computed from (7) requires an additional correction factor e to give desired results. Also, since the system under consideration consists of mechanical devices such as compressors and switches, such high frequencies are undesirable. Additionally, the TCLs are frequency sensitive devices, and drastic changes in the frequency must be avoided to prolong the product life. 28 In a practical scenario, where system parameters (e.g., N, αij) change dynamically, computation of ωi in real-time would seem difficult. Thus, rendering the Boolean phase oscillator model is difficult to be practically realized. 2. K is a parameter that is computed using dynamical equations governing the system parameters. 35 In the case of TCLs, the value of K can be calculated from the induced time delays in a programmable frame. 33 However, the interdependence of K and ωi in (6) adds complexity to its accurate computation.

IV. DISTRIBUTED AVERAGING MODEL FOR THE TCLs
Consider an undirected graph G = {V, E} with a set of nodes N in V and edge set E. Each edge of the graph E ∈ G is an unordered pair of distinct nodes. A real scalar quantity xi(t) is associated with node i at any given time t. The distributed averaging consensus algorithm computes the average 1 N ∑ N i=1 xi(0) iteratively, thus facilitating a local communication link between the nodes. The nodes are updated continuously based on their states as well as the states of all other neighboring nodes.
A widely accepted and used linear iterative algorithm can be stated as follows:ẋ for i ∈ {1, . . ., n} and t ∈ R + and Wij being the weight associated with each edge Eij. Since the nodes, in the presented case, are undirected, the weights associated with them are symmetric, i.e., Wij = Wji.
where wij is the weight of edge Eij and diag( N ) is the diagonal matrix with N , N×N denoting a vector and matrix of ones, respectively. Thus, the matrix W is the communication matrix that describes the strength of the link between two nodes. Analogous to frequency synchronization seen in the Kuramoto framework (6), we extend the distributed theory to the system of TCLs. The nodes are similarized to TCLs and the scalar quantity to the frequency of switching signal as follows: where W is the tuning parameter of synchronization, left to the control designer. (10) synchronizes the frequencies randomly chosen from distribution to their mean value, as shown in Fig. 5. Thus, (6) can be reframed as follows: It should be noted that (11) bifurcates the onus of controlling the frequency and angular separation.
For the matter of comparison, a similar ensemble of heterogeneous TCLs, as in Fig. 4 (natural frequencies in the same range as of ±5% of ωi = 0.271 Hz), are taken. Additionally, the value of weight (W) is chosen as 0.06. By using (11), the system desynchronizes, and Fig. 6 depicts the results obtained. Similar results can be noted by comparing the two plots, however, as opposed to the Kuramoto based model in Fig. 6, high frequency switching and computational hindrance have been avoided.

V. SOFTWARE AND HARDWARE IMPLEMENTATION
This section presents the application of the distributed averaging protocol in the form of software simulations and later describes its practical implementation.

A. Case study
To make the two cases put forth as realistic as possible, real time data of a TCL have been used to simulate the cases. For this, a TCL of a certain make was observed and the data listed in Tables II and III were collected. This TCL had a power rating of 1.66 kW and a capacity of 1.5 tons. The values R = 9.8 ○ C/kW and C = 0.0746 kW/ ○ C were computed using (1), Tables II and III. 1. Case 1: N = 1000 For N = 1000, heterogeneous TCLs are defined having a duty cycle in the range of 0.422-0.482, δ as 3.0 ○ C and f between 0.0029 Hz and 0.0033 Hz, respectively. The TCLs have been set to work at a set point temperature of 27 ○ C.
Following the distributed averaging model, TCLs settle at their respective phase difference, i.e., αij = α = 2π/N, which in this case is α = π/500. As depicted in Fig. 7, the power fluctuation for N = 1000 heterogeneous TCLs dampens over time and moves towards the steady state. The aggregate power at which Following the distributed averaging model, TCLs settle at their respective phase difference, i.e., αi = α = 2π/N, which in this case is α = π/5000. As depicted in Fig. 8, the power fluctuation for N = 10 000 heterogeneous TCLs dampens over time and moves towards the steady state. The aggregate power at which the TCLs settle is 8.4453 MW, and fluctuations in power are ±2.5%. In Fig. 9, phase portraits of 100 TCLs chosen randomly from the population are shown. It is clear that the portrait dissimilarity is minimal, although the set point temperatures remain unaffected.

B. Load following
Next, the authors show how a load-following case can be emulated, providing a valuable AS to the utility. In the previous work, 33 a function mapping RMS (root mean square) aggregate power of a population of TCLs with the effective delay induced was proposed. Using the same idea, we choose different utility demands and the corresponding delay of αij is calculated. The effective power so obtained is normalized using total RMS power capacities, i.e., where Prms,agg is the RMS of aggregate power, Prms,α is the RMS of aggregate power at the steady-state at αij = α, and Pnorm is the normalized power of all the TCLs along a given period, respectively. These results can be seen from Fig. 10(a), where utility demands for random reduction/increase in the total allowable aggregate power and the corresponding change in the delay induced causes controlled desynchronization and hence the required load following scenario is achieved. As seen from Fig. 10(b) (N = 100, P = 14 kW), it can be inferred that the resultant aggregate power is a significant AS to the centralized power grid.

C. Effects in a power grid
It should be noted that apart from the control of power aggregation, the proposed algorithm benefits the grid by reducing the power system fluctuations. To compare/visualize the same, we take a small population of N = 100 TCLs. To exhibit random behavior (i.e., depicting as -is characteristics without applying control), the frequency and phases of TCLs are chosen randomly from a distribution. Next, the same set of heterogeneous TCLs is acted upon by the distributed averaging protocol achieving effective desynchronization. As shown in Fig. 11(a), the system without application of the distributed averaging strategy induced higher oscillations in the grid as compared to after application of the proposed scheme. In Fig. 11(b), the quantum of reduction is observed to be ∼40% (560 kW) of a random case, showing a significant reduction [P red (%)], where Similarly, a case of N = 10 000 TCLs is also shown in Fig. 11(c). It can be seen that a reduction of ∼2% can be noticed, which translates to 2.8 MW. Hence, apart from providing load following objectives, the proposed scheme can also help reduce the amplitude of power system oscillations.

D. Hardware implementation
To validate the presented theory, the authors have implemented the model on a hardware system. The authors use the circuit diagram as shown in Fig. 12, which practically looks as seen in Fig. 13.
Light-emitting diodes (LEDs) provide close resemblance to a switching device, analogous to TCLs. Hence, for the demonstration purpose, LEDs have been used in place of actual TCLs. Since the number of LEDs was limited to four, a microcontroller with necessary computational capabilities and decent clock frequency was sufficient.

Hardware
The components used for the implementation include 4 LEDs, microcontroller, a breadboard, 5 V DC supply, wires for connection, and probes to display the output on the oscilloscope. The digital storage oscilloscope (DSO) used has the following specifications:

Implementation
The setup for the practical implementation of the distributed averaging model is shown in Fig. 13. Based on the distributed averaging model, as described by (11), the logic is coded and burned on Arduino. The Arduino generates switching signals for the LEDs. The output of the LEDs is observed on the oscilloscope.
As it can be seen, the output shown in Fig. 14 is similar to that of Fig. 6.  period, the LEDs have no phase difference between them. Figure 15 shows the synchronized state of the LEDs. The synchronized state of the TCLs is a highly undesirable condition in a power system. Furthermore, succeeding t = 12 s, the LEDs induce a phase difference of α = π/2 between them. As seen from Fig. 16, the output of the LEDs is minimized. Additionally, the fluctuations are also minimized.
The results of the hardware implementation were in conjunction with the results of computer simulations. Thus, the results of the proposed theory are validated.

VI. COMPETITIVE BENCHMARK
In this section, the authors compare the results obtained from the proposed model with the models mentioned in the literature. From Fig. 7, Pagg is 750 kW (45.18% of maximum power consumption). The fluctuations are quantified using root mean square error (RMSE%) given as follows: For a heterogeneous ensemble of TCLs, the relative error computed is within 2.49% and the RMSE% is 4.8685%. Juxtaposing these results alongside with the minimum variance law (MVC), 36 the latter achieves relative error less than 5%. Table IV shows the efficacy of a distributed averaging model against several state-of-the-art methods. The relative errors obtained  37 the RMSE% of the presented model is marginally higher. As compared to the Boolean phase oscillator model, 33 the relative errors are similar. Thus, the distributed averaging model achieves competent results while avoiding computational complexities.

VII. CONCLUSION
The major highlight of this work is the application of a distributed averaging protocol for the desynchronization of a population of TCLs. The system of TCLs is visualized as oscillators and their power fluctuations have been minimized. A significant cornerstone of the proposed work is that user comfort is unimpaired. The parameters governing the dynamics of the system are then articulated as a set of differential equations.
Apart from power fluctuations, aggregation, and user comfort, another breakthrough is the simplicity of the presented model and reduced computation. Load following helps the utility by improving the usability of the generated power, thereby supporting in loadfrequency control and maintaining voltage profile. A reduction in power fluctuation and aggregation would ease off the demand pressure from the industries directly or indirectly associated with power generation, distribution, protection, and maintenance. The results of both simulations, as well as hardware, support the findings of the presented works. The results obtained are motivating and show great potential to be implemented as a plug-n-play device and be made available in the market.