A Phase Model Approach for Thermostatically Controlled Load Demand Response
AA Phase Model Approach for Thermostatically Controlled Load Demand Response
Walter Bomela , ∗ Anatoly Zlotnik , † and Jr-Shin Li ‡ Department of Electrical and Systems Engineering,Washington University, Saint Louis, MO 63130, USA Theoretical Division, Los Alamos National Laboratory,Los Alamos, NM 87545, USA
A significant portion of electricity consumed worldwide is used to power thermostatically con-trolled loads (TCLs) such as air conditioners, refrigerators, and water heaters. Because the short-term timing of operation of such systems is inconsequential as long as their long-run average powerconsumption is maintained, they are increasingly used in demand response (DR) programs to bal-ance supply and demand on the power grid. Here, we present an ab initio phase model for generalTCLs, and use the concept to develop a continuous oscillator model of a TCL and compute its phaseresponse to changes in temperature and applied power. This yields a simple control system modelthat can be used to evaluate control policies for modulating the power consumption of aggregatedloads with parameter heterogeneity and stochastic drift. We demonstrate this concept by comparingsimulations of ensembles of heterogeneous loads using the continuous state model and an establishedhybrid state model. The developed phase model approach is a novel means of evaluating DR pro-vision using TCLs, and is instrumental in estimating the capacity of ancillary services or DR ondifferent time scales. We further propose a novel phase response based open-loop control policy thateffectively modulates the aggregate power of a heterogeneous TCL population while maintainingload diversity and minimizing power overshoots. This is demonstrated by low-error tracking of aregulation signal by filtering it into frequency bands and using TCL sub-ensembles with duty cyclesin corresponding ranges. Control policies that can maintain a uniform distribution of power con-sumption by aggregated heterogeneous loads will enable distribution system management (DSM)approaches that maintain stability as well as power quality, and further allow more integration ofrenewable energy sources.
I. INTRODUCTION
Significant efforts have been made in recent years to re-duce the greenhouse gas emissions caused by fossil-fueledpower plants throughout the world. In the United States,many states have been adopting or increasing their re-newable energy generation portfolio standards [1]. Theincreasing penetration of such renewable energy sources(RESs) affects power quality on electric distribution sys-tems and complicates load balancing in power systems[2]. These issues compel the development of new ap-proaches to regulate the inherently fluctuating and un-controllable power outputs of RESs [3, 4].Demand response (DR) programs enable electricityusers to adjust their consumption in response to energyprices or incentive payments [5, 6], and thus provide sig-nificant capability to balance supply and demand on thepower grid. Indeed, demand-side management technol-ogy has been under development since the late 1970’sand early 1980’s [7], with consideration of diverse objec-tives that include peak shaving, valley filling, and strate-gic energy conservation [8]. With the advent of smartgrid technologies and the integration of RESs in distri-bution systems around the world, system operators will ∗ [email protected] † [email protected] ‡ [email protected] now have to rely on DR and ancillary services (AS) tobalance supply and demand of energy more than everbefore. It is therefore imperative to develop technologiesthat can exploit the resources available for DR to thefullest extent.A comprehensive review on the current challenges andbarriers to a full deployment of demand response pro-grams is given in [9]. One important observation madein the review is that, despite receiving information abouttheir energy consumptions, the majority of the partici-pants in a DR study continued with their everyday rou-tines and habits. This kind of consumer behavior willdiminish the potential of DR. In order to maximize DRcapabilities, it may be better to remove human decisionsfrom the loop. This could be achieved by using a trans-active control paradigm, which is a popular direct loadcontrol approach in which customers have agreed to al-low the power company to control some of their flexibleloads [6]. This control approach reacts faster to marketprice fluctuations and enables TCLs to participate in thereal-time retail electricity market [5].The integration of RESs into the power grid reducesthe inertia of the system, which poses a challenge forautomatic generation control (AGC), which is used tocontrol frequency by balancing supply and demand. Sev-eral studies have examined ways to utilize electric vehi-cles (EVs) and battery energy storage systems (BESSs)to provide frequency regulation reserves [10, 11]. Forinstance, [12] developed a techno-economic evaluationframework to quantify the challenges and assess bene- a r X i v : . [ n li n . AO ] M a r fits to primary frequency control by using EV batteries.Although BESSs can provide primary frequency control,the cost of deployment in large capacities remains pro-hibitively high [13]. Hence, the focus has been turnedto TCL loads such as refrigerators [14, 15], heating ven-tilation and cooling (HVAC), and electric water heaters(EWH) that have been shown to be suitable for providingancillary services to the grid [16, 17].In this paper, we consider the direct control of TCLsthrough thermostat set-point changes for DR provision.Although we consider AC systems, the modeling andcontrol approach presented herein is applicable to anygeneral TCL. In the recent past, various modeling andcontrol methods for TCLs were developed, and exper-iments were conducted to demonstrate the ability ofTCLs to provide ancillary services. Field experimentsconducted using domestic refrigerators in [15] and [14]aimed at quantifying the flexibility of household TCLs,and the computational resource constraints for the con-trol of large TCL populations. The frequency controllerin [14] switches the TCLs OFF one by one based onthe ability of each refrigerator to stay OFF longer. In[15], the same authors introduced a delay in their con-trol scheme to improve the performance of the controller,which reduced the power overshoot. They noted thatwithdrawing a large number of loads abruptly producedinstability and caused the loads to synchronize their dutycycles. In [13], such synchronization of TLCs is avoidedby randomizing the parameters in the control scheme.This approach consisted of dividing the compressor cycleinto four different states based on the compartment tem-perature. Their control architecture was based on theutility sending a regulation signal to Cooperative HomeEnergy Management (CoHEM) systems located at dis-tribution transformers, with each CoHEM then sendingdifferent control signals to home energy management sys-tems that control the refrigerators. The investigation in[18] combined three different control protocols to improvethe accuracy in tracking an AGC signal. The controllerswitched between a temperature priority based control,a sliding mode control, and a two-stage regulation con-trol. A two-level scheduling method intended to facilitatethe scheduling of flexible TCLs in the intra-day electric-ity market was proposed in [16]. However, the perfor-mance of this method deteriorates with parameter het-erogeneity, and worsens as forecast uncertainty increases.Control approaches based on model predictive controltechniques have also been considered. In [19], a multi-objective model predictive control strategy for residentialheating with heat pumps was proposed. This approachtakes into account the users’ energy cost, the environ-mental impact, and expansion of electricity generationcapacity. In addition, it considers detailed models forheat pump and thermal energy storage, and accounts forthe feedback effect of individual controllers on electric-ity generation. This level of detail will make scaling theapproach to a large TCL population problematic.One very important question that was posed in [20], and that still needs to be accurately answered, is howlarge is the potential of DR and on which time scale canDR be the most effective? This is not a particularly easyquestion to answer. In this paper we propose a novelprocedure for evaluating the maximum capacity and theappropriate time scale for DR provision by a collectionof TCLs. Based on phase reduction theory [21], this ap-proach can evaluate the capacity of a TCL population indifferent bandwidths. In [22], a theoretical upper boundon the capacity of TCLs to provide ancillary services asa function of frequency was proposed. The capacity-bandwidth constraints was derived based on standardlinear dynamical system models. The proposed boundindicates that for a given TCL population with a fixedtime constant, the capacity for AS provision is inverselyproportional to the frequency of the regulation signalto track. Our method supplements this observation byfurther showing that the capacity does not monotoni-cally decreases with frequency. Indeed, while the capac-ity monotonically decreases for a regulation signal withfrequency below the mean natural frequency of a givenTCL population, it then increases and decreases again insome frequency bandwidths higher than the TCLs’ meanfrequency. This is particularly significant because thistheory can be used to analyze different TCL populationsand classify them based on their capacities to provide ASin different bandwidths (time scales).Unwanted synchronization of TCL duty cycles by thecontrol policy is most likely the main factor that lim-its the time scale and capacity of ancillary services thata given population can provide. This phenomenon hasmotivated the research and development of control poli-cies such as the safe control protocols [23, 24], whichaim to minimize unwanted power oscillations in responseto pulse-like changes of the set-point temperature. Oneprotocol allows the TCL to stay in its current state untilthe temperature hits one of the transition points, thenstarts following a new pair of deadband limits [23], andthe other adds a delay of M-minutes before changing thestatus of the TCL [24]. Other control strategies with sim-ilar goals consist of turning ON/OFF some TCLs in theOFF/ON stack based on their priority measure [25], [26].The priority measure is defined as the distance from theswitching boundary and hence, if the TCLs’ aggregatepower needs to be reduced, the TCLs in the ON stackcloser to the switching boundary will be turned OFF first.Conversely, if the aggregate power needs to be increasedto follow an increase in the generated power, the TCLsin the OFF stack closer to the switching boundary willbe turned ON first.In an effort to further the understanding of the oscil-latory behavior of the aggregate power consumption, thedamping of oscillations in aggregate power was charac-terized as a function of parameter heterogeneity by ex-ploiting the similarities that exist between a populationof mass-springs systems and an ensemble of TCLs [27].Concurrently, the dependence of the mixing rate of thepopulation on the model parameters was characterizedin [28].To summarize, the work in this paper provides an novelapproach for evaluating DR capacity and time scale us-ing a novel phase model representation of TCLs. Theevaluation is conducted by computing the entrainmentregions also known as Arnold tongues [21]. To obtain thephase model, we started by first extending the commonhybrid-state model of TCL dynamics to a neuroscience-inspired representation in the form of a continuous two-dimensional system, then applying a widely-used phasereduction computation [29]. The phase model is a one-dimensional system with scalar state, and its simplicityhas made it one of the most popular models for studyingoscillatory systems, including power grids [30] and neuraloscillators [31], with particular advantages for control de-sign in the presence of parameter uncertainty [32]. Thesecond novel contribution of this work is the proposedphase response based open-loop control policy, which ex-ploits the natural diversity that exists amongst the phaseresponse functions of TCLs. This approach is promis-ing for developing hierarchical grid control architecturesthat maintain power quality on distribution systems withlarge penetration of RESs while providing load balancingservices to regional transmission operators (RTOs). Tothe best of our knowledge, this paper is the first to con-sider the modeling and control of TCLs using the phasemodel reduction approach.The rest of the paper is structured as follows. An abinitio scalar phase model of a basic hybrid system rep-resentation of a single TCL is derived in Section II B.Then, a continuous two-state TCL model and its phase-reduced representation are derived in Sections II C andII D. Monte Carlo simulations of an example TCL de-scribed using the two models are then compared in Sec-tion III. After this model validation, in Section IV weanalyze the synchronization properties of TCL ensemblesin order to enable the evaluation of DR capacity and timescales, and explore a phase response based open-loop con-trol policy whose efficacy for controlling the aggregatepower of TCL ensembles is demonstrated by tracking aregulation signal from the Balancing Authority (BA) de-composed into power spectrum components. A discus-sion of applications and promising future research direc-tions is given in Section V.
II. THERMOSTATICALLY CONTROLLEDLOAD MODELSA. One-Dimensional Hybrid Model
The dynamics of the internal temperature θ ( t ) of ahouse equipped with an air-conditioning system is oftendescribed using a simple hybrid-state model [33]. Themodel describes how an AC unit regulates the averagetemperature by means of a thermostat and a relay withstate s ( t ) ∈ { , } [34]. The hybrid state model describ-ing the evolution of the internal temperature θ ( t ) is given by ˙ θ ( t ) = − RC [ θ ( t ) − θ a + s ( t ) P R ] , (1) s ( t ) = θ ( t ) < θ min θ ( t ) > θ max s ( t ) otherwise , (2)where θ a represents the ambient temperature, P is theaverage energy transfer rate of the TCL in the ON state, C and R are the thermal capacitance and resistance ofthe house, respectively. The minimum and maximumtemperature of the TCL are θ min = θ s − δ/ θ max = θ s + δ/
2, respectively, with θ s the thermostat temperaturesetpoint and δ the deadband. Consider a population of N TCLs with temperature states θ i ( t ) for i ∈ { , · · · , N } that evolve according to (1). Then the aggregate power P agg ( t ) drawn by all N TCLs is [35] P agg ( t ) = N (cid:88) i =1 η i s i ( t ) P i , (3)where η i is the coefficient of performance of each TCLunit. B. Ab Initio Deterministic TCL Phase Model
We formulate a basic model of a TCL as a deterministicswitched oscillating system, whose state is given by thetemperature θ ∈ [ θ min , θ max ] and the switching function s ∈ { OFF, ON } . Suppose without loss of generality thatthe system is an air conditioner that cycles through aduty cycle to maintain a temperature θ s = ( θ max − θ min )while remaining within the deadband δ = θ max − θ min .We denote the beginning of the cycle as the point ( θ, s ) =( θ max , ON) just after the unit has turned on. The state ofthe unit changes according to ˙ θ = r − < s = ON, then it switches to s = OFF when the statereaches θ = θ min . Then, the state of the unit changesaccording to ˙ θ = r + > θ = θ max , andthe unit turns on again. Let us call the length of timewhen the unit is ON T ON , and the length of time whenthe unit is OFF T OFF , so that the period of oscillation is T = T ON + T OFF . onoff FIG. 1. Phase model of a switched oscillating system.
Our goal is to map this behavior to a phase model,as illustrated in Fig. 1. Such models are desirable be-cause the homogeneous dynamics are linear and scalar.The state is represented by a scalar phase φ , whichadvances linearly with time according to a frequency ω = 2 π/T . For an unforced system, this yields a sim-ple solution φ = ωt (mod2 π ). Thus, we map the statepoint ( θ max , ON) where the unit turns ON to the phasepoint φ = φ ON ≡
0, and the state point ( θ min , OFF)where the unit turns OFF to the phase point φ = φ OFF .This yields a continuous representation of the switchedsystem.To complete the picture, we first determine the phase φ OFF ∈ [0 , π ) when the unit switches OFF. If r + = r − , then it is straightforward to show that φ OFF = π .However, most units will have different values of r + and r − , which will also depend on other factors such as theambient temperature. Suppose then that r + (cid:54) = r − , andlet us denote the rate ratio, which is equivalent to theduty ratio, by γ = − r − r + = T OFF T ON . (4)It is straightforward to show that these ratios are equiv-alent. Integrating the s = ON dynamics ˙ θ = r − from t = 0 to t = T ON yields r − T ON = θ min − θ max = − δ, (5)and integrating the s = OFF dynamics ˙ θ = r + from t = T ON to t = T = T ON + T OFF yields r + T OFF = θ max − θ min = δ. (6)Then, equation (4) follows directly from (5) and (6), andwe then have γ >
0, and T OFF = γT ON . We can thencompute T ON and T OFF according to T ON = T − T OFF = T − γT ON = 11 + γ T = r + r + − r − T,T
OFF = T − T ON = T −
11 + γ T = γ γ T = − r − r + − r − T. It follows that the switch-OFF phase is given by φ OFF = 2 πT T ON = 2 π r + r + − r − = 2 π T ON T ON + T OFF . It is straightforward to show that the parameters in thehybrid-state model given in equations (1)-(2) are relatedto the phase model parameters by T OFF = RC ln (cid:18) θ a − θ min θ a − θ max (cid:19) , T ON = RC ln (cid:18) θ max − θ a + P Rθ min − θ a + P R (cid:19) , with the period given by T = RC ln (cid:18) ( θ a − θ min )( θ max − θ a + P R )( θ a − θ max )( θ min − θ a + P R ) (cid:19) . (7)Hence, we can express the switch-off phase by φ OFF = 2 π ln (cid:16) θ max − θ a + P Rθ min − θ a + P R (cid:17) ln (cid:16) ( θ a − θ min )( θ max − θ a + P R )( θ a − θ max )( θ min − θ a + P R ) (cid:17) . Note that to control the average power of a TCL, oneneeds to modulate the duty cycle of the power utilization.This is equivalent to controlling the phase φ OFF of thephase model illustrated in Fig. 1 assuming for instancethat φ ON = 0 is the reference phase. C. Continuous Representation of SwitchingDynamics
We now present a two-dimensional continuous state-space model that approximates the dynamics of the in-ternal temperature of a TCL unit and the thermostatswitching. This is an important step that precedes thederivation of the complete phase model, which will re-quire the computation of the phase response curve (PRC)of the continuous oscillator model. Subsequently, we willexamine the phase response of the temperature dynamicsto control action applied to the TCL.Consider the evolution of the temperature θ ( t ) (see Fig.2a) described by (1), and the corresponding phase por-trait in Fig. 2c, simulated with the parameters providedin Table 1. Observe that the unperturbed behavior of aTCL is similar to that of an oscillator with a stable limitcycle. Time [h] ( t ) [ o C ] s (t) Time [h] s ( t ) OFF ON s(t) ( t ) [ o C ] Phase portrait, s = 20 o C maxmax min min (b)(c)(a) FIG. 2. Simulation of the hybrid model (1)-(2) with theparameters in Table 1 and the deadband δ = 1 . ◦ C. (a)Time evolution of the temperature θ ( t ) around the setpoint θ s = 20 ◦ C. (b) Thermostat switching function s ( t ). (c) Phaseportrait. TABLE I. Nominal TCLs parameter values [36]Parameter Meaning Value θ s temperature setpoint 20 ◦ C θ a ambient temperature 32 ◦ C δ thermostat deadband 0 . ◦ C R thermal resistance 2 ◦ C/kW C thermal capacitance 10 kWh/ ◦ C P energy transfer rate 14 kW η coefficient of performance 2 . The hybrid-state nature of the system described by(1) and (2) is due to the thermostat switching function s ( t ) that transitions between 1 and 0 states. There-fore, modeling the behavior of a TCL using continuousstates requires a continuous approximation of the switch-ing function s ( t ) and its dynamics. Ideally, the evolutionof the temperature and the continuous switching func-tion could be represented using a system of two coupleddifferential equations that has a stable limit cycle sim-ilar to the one shown in Fig. 2c. Our motivation forthe proposed model is the Van der Pol oscillator, whichis a simple model of the limit cycle observed in circuitswith vacuum tubes. More so, a similar phase portraitis observed in the FitzHugh-Nagumo model, which is asimple mathematical description of the firing dynamicsof a neuron [29]. Inspired by these examples, we proposethe continuous-state TCL model given by˙ x ( t ) = µ (cid:18) ( δ σ ) x − x θ − θ s (cid:19) , ˙ θ ( t ) = − RC ( θ − θ a + ¯ s ( t ) P R ) , (8)where x ( t ) is the state variable of the switching function, θ ( t ) the internal temperature, and ¯ s ( t ) is an approxima-tion of the ideal switching function s ( t ). The parameter σ was introduced to compensate for the reduction of theeffective deadband in (8). The constant µ is a damp-ing parameter that controls the oscillation frequency fora fixed time constant τ = RC as well as the shape ofthe phase portrait. Just like for the Van der Pol oscilla-tor, as µ gets smaller the limit cycle takes the shape ofa circle and for large values of µ the shape of the limitcycle appears more rectangular. Given the rectangularshape of the limit cycle of the hybrid model in (1)-(2),we chose a large value of µ = 100 such that the limit cy-cle of (8) is similar to the one in Fig. 2c while oscillatingapproximately at the same frequency as (1).Once the value of µ is fixed, the small difference in thedeadband that translates into frequency deviation can becompensated for by the parameter σ < δ . The parameter σ can quickly be determined through numerical simula-tions as follows. Knowing that σ ∈ [0 , δ ), one can sample n values of σ over this interval, simulate the dynamics(8) for each value of σ , and compare the oscillation fre- quencies at each value to the nominal frequency of thehybrid state model in (1) to determine the most appro-priate value of σ . This is only done once using the TCLwith a natural frequency corresponding to the averagefrequency of the population of TCLs considered. For theparameters in the Table 1, σ = 0 . s ( t )from the variable x ( t ) that is as similar as possible tothe switching function s ( t ) in the hybrid model, we usean approximation of the Heaviside step function. Suchstep function approximations have been used in a similarway in the Morris-Lecar neuron model, whose behavioris similar to the Hodgkin-Huxley spiking neuron model[31]. The resulting switching function approximation isgiven by ¯ s ( t ) = 12 (1 + tanh( kx )) , (9)where the parameter k ∈ [5 ,
10] controls the sharpnessof the switching action. As shown in Fig. 3d, the phaseportrait is similar to that of the hybrid-state model shownin Fig. 2c.
Time [h] -2-1012 x ( t ) , s ( t ) x(t) s bar (t) Time [h] [ ° C ] s (t) -1 0 1 x(t) [ ° C ] s(t) [ ° C ] (d)(c) (b)(a) FIG. 3. Simulation of the system in (8). (a) Evolution oftemperature θ ( t ) around the setpoint θ s = 25 ◦ C. (b) Theswitching state variable x ( t ), and the switching function ¯ s ( t ).(c, d) Limit cycles for both θ ( t ) vs. x ( t ) in (c) and θ ( t ) vs. ¯ s ( t ) for the new switching function (9) in (d). The parametersused are in Table 1, but with C = 2 kWh / ◦ C and δ = 1 . ◦ C. The model given in (8) reproduces the dynamical be-havior of the hybrid-state model in (1). However, ourgoal is to design an ensemble control of a TCL popula-tion such that the aggregate power closely tracks a givenreference power. Therefore, we modify the model in (8)by substituting for the state variable x ( t ) with the in-stantaneous power y ( t ). The derivation of this model isprovided in Appendix A. The resulting model is givenby˙ y ( t ) = µk (cid:18) ( δ σ )¯ y −
13 ¯ y + θ − θ s (cid:19) (1 − ηP y ) y, ˙ θ ( t ) = − CR ( θ − θ a + ηy ( t ) R ) , (10)where ¯ y ( t ) = − k ln( Pηy − P divided by the co-efficient of performance η . The aggregate electric powerof an ensemble of N TCLs described by (10) is given by P agg ( t ) = N (cid:88) i =1 y i ( t ) . (11)The evolution of the TCL’s internal temperature in (10),its instantaneous power and the limit cycle are shown inFig. 4. Time [h] ( t ) [ o C ] s (t) Time [h] y ( t ) [ k W ] Power y(t) [kW] ( t ) [ o C ] Phase portrait, s = 20 o C (b)(a)(c) FIG. 4. Simulation of the system in (10). (a) Evolution ofthe temperature θ ( t ) around the setpoint θ s = 20 ◦ C. (b) Evo-lution of the electric power y ( t ) drawn by the TCL. (c) TCLphase portrait. The parameters used are given in Table 1. D. TCL Phase Model
In Section II B, we have shown that the unforced cy-cling dynamics of a TCL can be represented by a phasemodel ˙ φ = ω , where ω is the natural frequency. How-ever, in the presence of an external control input u ( t ),the phase model takes the form dφdt ( t ) = ω + (cid:15)Z ( φ ) u ( t ) , (12)where Z ( φ ) is the phase sensitivity function, also knownas phase response curve (PRC) and (cid:15) is the intensity of the perturbation or control input u ( t ) [29]. For (8) and(10), the PRCs are vectors Z ( φ ) = ( Z s ( φ ) , Z θ ( φ )) and Z ( φ ) = ( Z y ( φ ) , Z θ ( φ )), respectively, where Z s ( φ ) , Z y ( φ )and Z θ ( φ ) are phase sensitivity functions of the switch-ing variable x ( t ), the instantaneous power y ( t ), and thetemperature θ ( t ), respectively (see Fig. 5). Because weare only interested in controlling one state variable ofthe TCL e.g., the switching s ( t ) or the power P ( t ), thePRC in (12) will be taken as one of the scalar functions Z ( φ ) = Z s ( t ) or Z ( φ ) = Z y ( t ), and the scalar controlinput will be a temperature signal to offset the set-point.If on the other hand we were to consider the phase modelof the temperature θ ( t ), the PRC will be Z ( φ ) = Z θ ( t )and the corresponding control will be the input power.The natural frequency is computed as ω = 2 π/T us-ing the period T in (7) and it is the same for all threevariables s ( t ) , P ( t ) and θ ( t ). The PRC itself must becomputed numerically, using for example the method ofthe adjoint which requires the computation of the Jaco-bian of the dynamics of the system described by a con-tinuous function. This method consists of linearizing thesystem around its periodic orbit Γ( t ) and solving the ad-joint equations with a backward integration [29, 31]. Astandard software package used for the computation ofthe PRC is the XPPAUT [37]. The ability to computethe PRC using the method of the adjoint was one of thereasons for the derivation of the continuous models de-scribed in Section II C. Another important reason forsuch modeling is the ability to access the switching vari-able x ( t ) in (8) or the power variable y ( t ) in (10) in orderto characterize their dependence on a stimulus u ( t ). Thephase sensitivity quantifies how much the phase of thesevariables changes in response to impulse I ( t ). In Ap-pendix B, we review some basics of the phase reductiontheory as well as different techniques for determining thephase sensitivity curve of an oscillatory system.From the phase model (12), one can see that the con-trol action is applied through the PRC, consequently theoscillator phase will either be advanced or delayed. Thismeans that the time at which the TCL turns ON or OFFcan either be delayed or advanced by ∆ T when an ap-propriate control input is applied. The implication isthat the duty cycle of the input power to the TCL canbe modulated, and therefore its instantaneous or averagepower can be increased or decreased accordingly.The phase sensitivity functions for the temperatureand power evolution described by (10) are shown in Fig.5 for various values of the thermal resistance R . Thisshows that the parameter heterogeneity of the TCL en-sembles will cause each unit to respond differently to acommon control. The PRC of (8) looks similar to that of(10) except for the amplitude of Z s ( φ ) which is smallerthan Z y ( φ ). [rad] -10-50510 Z y () [r a d k W - ] [rad] -505 Z () [r a d o C - ] R = 1.8R = 2.0R = 2.2 [rad] P () [ k W ] [rad] () [ o C ] (c)(a) (d)(b) FIG. 5. Phase response curves of the power and temperaturefor a range of thermal resistance R ∈ [1 . , . φ . (c) and (d) show their phase sensitivityfunctions Z y ( φ ) and Z θ ( φ ), respectively, for the dynamics in(10). The parameters used are shown in Table 1. III. SIMULATION COMPARISON OF TCLMODELS
In this section, we provide simulation results that areintended to show how well the continuous model (10)proposed in Section II C approximates the hybrid modeldynamics in (1)-(2). We first compare the phase modelin (12) to the hybrid-state model in (1)-(2) by plottingthe time evolutions of the phase φ and the temperaturetogether in Fig. 6a. The phase φ OF F indicates wherethe TCL turns OFF after being ON for a time T ON . Theinput powers for both models are shown Fig. 6b, withthe control u ( t ) = 0.Second, using the data in Table 1, we simulate the ag-gregate power of 1,000 heterogeneous TCLs for a periodof 30 hours. The initial values of temperature and theparameters C , R and P were randomly distributed uni-formly within ±
5% of their nominal values. The resultsin Fig. 6c show good agreement between the transient os-cillations of both the hybrid and continuous models, andthe stationary variation about the long-run mean whichoccurs after 10 hours appears similar as well.Now that we have shown that the phase model (12)captures the cycling dynamics of a TCL with sufficientaccuracy, in the next section we will analyze the syn-chronization properties of phase model representationsof TCLs, and then derive a PRC-based control policy.
IV. PRC-BASED CONTROL POLICY
Direct control of TCLs by the Balancing Authority canenable regulation of power consumed by such loads on adistribution subsystem over faster time scales than the
Hybrid model Phase model
Hybrid model Phase model
Time [h] P o w e r [ M W ] Hybrid model Continuous model
OFF T ON (a)(b)(c) FIG. 6. Comparison between the Hybrid model and its de-rived phase model representation. (a) Time evolution ofthe temperature (blue, solid) and the phase evolution (red,dashed). (b) Time evolution of the power consumed by theTCL for both models. Note that the TCL turns OFF exactlyafter being ON for time t = T ON and the switching phaseis φ OFF . (c) Aggregate power of 1000 heterogeneous TCLsdescribed by the hybrid and continuous models (1) and (10),respectively. Simulated with the parameters in Table 1. price response approach allows. Among possible con-trol architectures proposed in the literature, one mayfind centralized, hierarchical, and distributed controllers[38]. Implementing a centralized feedback control pol-icy to modulate the power consumption of an ensem-ble of TCLs is challenging. It requires the measurementand transmission of the states of each TCL to a cen-tralized controller, resulting in extra cost for measure-ment equipment and communication channels that needto be established, which has also raised some privacyconcerns [23]. Other issues such as temporary synchro-nization arise when a step change in the control input u ( t ) is applied to a large number of TCLs. This causesthe aggregate power to strongly oscillate for up to a fewhours before settling to the desired steady state value[24, 34, 39]. In this section we analyze and implementa PRC-based control policy that effectively regulates theaggregate power of an ensemble of TCLs on a relativelyfast time scale while preventing undesirable power fluctu-ation that can be caused by temporary synchronization. A. Control of a Single TCL
Phase reduction theory has been used in various scien-tific areas including physics, neural engineering and biol-ogy. This powerful technique has enabled the analysis ofsynchronization properties of limit cycle oscillators [29].One can find various applications in the literature suchas entrainment of chemical oscillator [40], optimal en-trainment of neural oscillator ensembles [32] and phaseadvance or delay in circadian oscillators where light isused as a control input [41]. Most applications of thephase reduction theory consist of either controlling thephase or the frequency of oscillators e.g., controlling thespiking time of neurons [42], and different control tech-niques have been developed for that purpose [32, 41].In the present application, in addition to regulatingthe switching frequency of a TCL, we wish to also mod-ulate the duty cycle of the input power such that theaverage power consumption during a period of time T iseither increased or decreased based on a regulation sig-nal ξ ( t ). By appropriately switching the TCLs ON andOFF as illustrated in Fig. 7, one can modulate the aggre-gate power of an ensemble over a short time scale with-out impacting the average temperature of the individualunits over the long run. The controls used to producethe results in Fig. 7 are of the froms u ( t ) = Z + ( φ ) ξ ( t )or u ( t ) = − Z − ( φ ) ξ ( t ) for increasing or decreasing thepower consumption, respectively, where Z ± ( φ ) representthe positive and negative parts of the switching PRC( Z s ( φ )), and the regulation signal ξ ( t ) = ∓ . ◦ C wasused to increase and decrease the power by forcing theTCL to switch either ON (Fig. 7b) or OFF (Fig. 7c)early. set (t) (a) (b) (c) FIG. 7. Illustration of the phase advances induced by thecontrol input u ( t ). Column (a) shows the unperturbed sys-tem. Columns (b) and (c) show the controlled systems. In(b), the control switches the TCL ON before its normal turnON time by temporarily decreasing the set-point. The oppo-site happens in (c). The first and second rows of plots showthe temperatures and their corresponding phases. The thirdand fourth rows show the corresponding power and controlwaveforms. The phase advances induced by ξ ( t ) in the phase model(12) are depicted in the second row of Fig. 7. By switch- ing a large number of TCLs ON, the aggregate powergiven by (3) will instantaneously increase, and converselyswitching them OFF will decrease the aggregate power.It is equally possible to delay the phase which results inthe TCL staying ON or OFF longer than it would natu-rally. An example is illustrated in Fig. 8. Set Temp.(t) P (d)(c) (b)(a) FIG. 8. Illustration of a phase delay induced by the control u ( t ) = − Z − ξ ( t ), with ξ ( t ) = − . ◦ C . (a) Evolution of thephase φ ( t ) showing a slow down (delay) when the control isapplied. (b) Evolution of the temperature showing that theTCL stays ON longer, hence the temperature went below thelower limit of 19 . ◦ C. (c) and (d) show the control input(shift of the set-point temperature) and the TCL power overthat period of time, respectively.
B. Control of an Ensemble of TCLs
One of the main challenges encountered when design-ing an open-loop control policy for a large populationof TCLs is the temporary synchronization that can becaused by a sudden change of the set-point. Although thechanges can be small (0.1-0.5 ◦ C) and barely noticeableby the customer, they can induce large power fluctuations[36, 43]. Alternatively, such policies can track the aggre-gate power reference closely when it is relatively slowlyvarying and smooth. The thermostat in this case mustbe assumed to be adjustable with infinitesimally fine pre-cision, although this requirement may be relaxed in ap-plication. To address the synchronization issues, varioussolutions have been proposed. For instance, a feedbackcontrol for TCLs with finite thermostat resolution wasexplored in [34]. Among the open-loop control policies,[44] developed a switching-fraction broadcast signal thatdetermines the number of loads to be switched ON orOFF by solving an optimization problem with a multi-linear objective.The open-loop control policy proposed herein ad-dresses a similar issue as encountered in the applicationof the safe methods of [23, 24] or the priority stack frame-work [25, 26]. Specifically, we aim to prevent excessiveaggregate power fluctuations caused by synchronizationof a large number of TCLs in the ensemble that in turnarise from sudden changes of the control signal (e.g., set-point temperature). The synchronization is prevented byletting the phase sensitivity function dictate the responseof a unit to a quickly changing control signal. Note thatbecause of TCL heterogeneity, the PRCs of all the unitsin a population are different. The PRC-based controlarchitecture is presented in Fig. 9. It is assumed thateach TCL can measure its state variables (temperatureand switching status) and has knowledge of its own phasesensitivity function Z ( φ ). Hence the controller can gener-ate the corresponding control u ( t ) in response to a globalregulation signal ξ ( t ) from the BA. Aggregate powerξ(t) ... TCL 1Local Control 1Balancing Authority Power GridTCL NLocal Control N u(t)Local feedback loop P(t)P(t)Local feedback loopu(t)
FIG. 9. Block diagram of the control architecture. The BAregulation signal ξ ( t ) is sent to all the TCLs in the popula-tion. This signal requests that each TCLs changes its powerconsumption by a given fraction that, if aggregated, will com-pensate for demand on the power grid. The local controllerreceives ξ ( t ) as well as the states of the TCL, then determinesthe appropriate control input u ( t ). The open-loop control structure illustrated in Fig. 9presents several practical advantages. First of all, it con-siderably limits the information complexity required atthe Balancing Authority level in the sense that no feed-back is required to form the reference signal ξ ( t ). Thisalso limits the computation and communication coststhat would have occurred if feedback was needed. Bysuppressing the need for a two-way communication, thiscontrol policy also addresses privacy concerns [23]. Theassumptions that we make on capabilities of the TCL andcommunication with the BA can be summarized as 1) theBA has knowledge of the power utilization of the TCLensemble and its capacity to service demand response, 2)in some instances the BA can totally or partially estimatethe aggregate power of the ensemble being controlled, 3)each TCL has knowledge of its own PRC and can mea-sure its internal states, and 4) each TCL is equipped witha control unit and a thermostat that has fine resolutionon the deadband.As a proof of concept, we propose a PRC-based inte-gral controller, and evaluate its performance against atraditional integral controller and the direct control. Bydirect control, we refer to situations with feedback andno local controller, in which the signal ξ ( t ) directly con-trol the TCL ensemble. We further show its efficacy by tracking a real Area Control Error (ACE) signal takenfrom the Bonneville Power Administration (BPA) web-site [45]. The proposed controller is of the form˙ u ( t ) = I [ I sgn( − s ( t )) Z ( φ ( t )) ξ ( t ) − u ( t )] , (13)where sgn is the signum function that extracts the signof 1 / − s ( t ), I and I are control gains. Discretizingthe differential form of the control in (13), we arrive at u ( t ) = u k for t ∈ [ t k , t k +1 ] where u k +1 = u k + I h [ I sgn( − s k ) Z ( φ k ) ξ k − u k ] , (14)where h = t k +1 − t k is the time step, and we have substi-tuted dependence on the discrete time t k by the subscript k for simplicity. The control equation (14) describes whateach local controller in Fig. 9 is doing when the regu-lation signal ξ ( t ) is received. To understand how thiscontrol policy is able to track a regulation signal withoutexcessive synchronization of TCL dynamics, observe that ξ ( t ) enters the control signal through a product with thePRC ( Z ( φ ) = Z y ( φ )), which is different from one TCLto another (see Fig. 5). This implies that each TCLwill respond differently according to its own parameters.Without loss of generality, we may assume that at time t = t k , u k = 0 and ξ k >
0, requesting that the TCLsreduce their power consumptions by increasing their set-points by a small fraction of the nominal value. Giventhe shape of each PRC, the requested change will notbe instantaneous for all TCLs. Depending on where theTCLs are in their cycles, some will switch their statusright away while others will do so with a delay that isfunction of the PRC. The gains I and I control the re-sponse time of the controller and the steady state error,respectively. The choice of an integral controller of theform (14) was motivated by the desire to minimize higherorder harmonics present in the PRCs. C. Analysis of Temporary Synchronization in TCLEnsembles
The potential demand response (DR) service that anensemble of TCLs could be used to provide to a powergrid is limited by many factors such as the specified limitsfor maintaining customer quality of service, the number N of TCLs in the population (the power capacity of aTCL ensemble increases with N ), and crucially the fre-quency bandwidth in which DR can be extracted. Tem-porary synchronization that causes undesirable fluctua-tion of the aggregate power is a consequence of the band-width limit.In this section, we provide some useful tools that elu-cidate the synchronization behavior of TCL ensembles.Observe what happens when a population of TCLs ismade to track a step change in the reference power (Fig.10). Before the step change is applied at time t = 5 h , theinitial conditions and the finite number of TCLs in theensemble cause the aggregate power to oscillate with a0frequency ω ≈ /N (cid:80) Ni =1 ω i , where the ω i ’s are the nat-ural frequencies of the TCLs in the population. Usingthe wavelet transform, we compute the power spectrumof the aggregate power as a function of time. It appearsthat the TCL ensemble naturally has damped oscillationswith a mean frequency ω = 5 .
69 rad/h and that a stepchange at t = 5 h amplifies these oscillations (Fig. 10c).Observe that by superimposing a decaying sinusoidal sig-nal whose sign is opposed to the step change, the poweror amplitude of the oscillations greatly decreases (Fig.10b-d), which implies that it is possible to design a con-trol that can change the power usage of the ensembleover a short time scale while minimizing the unwantedoscillations. P a gg r . [ M W ] PaggP ref
Magnitude Scalogram
Time [h] F r e qu e n c y [r a d / h ] PaggP ref
Magnitude Scalogram
Time [h] M a gn it ud e (b) (d)(c) (a) FIG. 10. Response of a heterogeneous TCL ensemble to astep change in the reference power. (a)-(c) Aggregate powerresponse to a step change and its wavelet transform magni-tude scalogram. (b)-(d) Aggregate power response to a stepchange with an opposing decaying sinusoidal of the same fre-quency as the induced oscillations. Comparing (c) and (d)reveals that the introduction of a decaying sinusoid with signopposite the induced oscillation reduces the power content ofthe undesired oscillations considerably.
It is crucial to note that the step at t = 5 h behaves likean impulse stimulus whose power content extends overall the frequencies and the population is strongly excitedby the frequency closer to its natural mean frequency.It appears as if all the TCLs are now oscillating withthe same frequency ω . This phenomenon is known inthe study of rhythmic systems as frequency entrainment[40, 46, 47]. The phase model representation of nonlinearoscillators becomes highly valuable in this case becausethe phase sensitivity function Z ( φ ) can be used to providethe theoretical linearized limits of the entrainment regioncommonly referred to as Arnold tongue (see Fig. 11) [48–50]. For more details see Appendix C.The shape and width of Arnold tongues depend onthe PRC and the control input waveform. For entrain-ment purposes, the control input waveform can be de-signed to increase the width of the Arnold tongue result- / P R M S FIG. 11. Theoretical Arnold tongues for a sinusoidal inputwith frequency Ω applied to one TCL with natural frequency ω . The entrainment ratios are indicated as N:M, with Ω = N/Mω . For this example, the TCL natural frequency is ω =7 .
95 rad/h, the parameters in Table 1 were used with differentthermal capacitance C = 1 . ◦ C and the deadband δ =1 . ◦ C, respectively. ing in a maximum entrainment range [50] or fast entrain-ment [51]. While maximizing the width of the Arnoldtongue is good for entrainment, for control of TCL en-sembles it is to be avoided. In Fig. 12 we have generatedArnold tongues for three different sets of TCL ensem-bles with different mean frequencies ω by measuring theRoot Mean Square Error (RMSE) of the aggregate powertracking with respect to a sinusoidal reference signal withfrequency Ω and amplitude A . The tracking error wasmeasured at different power level A and frequency Ω.As it can be seen from Fig. 12c, the Arnold tonguescorresponding to the 1:1 entrainment result in the highesttracking error. The 2:1, 3:1 and 4:1 entrainment regionsalso have a relatively high tracking error. This is be-cause temporary synchronization happens and there aremore TCLs turning ON or OFF at the same time than itis needed for tracking a reference. Hence by computingthe Arnold tongue one can determine the regions in thepower vs. frequency (P,Ω) space where the TCLs can pro-vide ancillary service with minimal oscillatory responseor better accuracy. Unlike the bound provided in [22],which suggests that the tracking capacity decreases lin-early with frequency, the results in this paper show thatthis is not entirely the case. Fig. 12c for instance showsthat indeed the tracking capacity decreases linearly withthe input frequency and reaches its minimum at Ω = ω but, it increases and decreases again forming bell shapesin the intervals [ ω , ω ], [2 ω , ω ] and [3 ω , ω ]. Thisimplies that it is possible to extract responsive regula-tion from a TCL ensemble at time-scales well beyond theaverage of natural frequencies of the TCLs by applyingan appropriate control policy that minimizes the area ofthe Arnold tongue. An example of the Arnold tongueof a TCL ensemble controlled with a proportional and1 Ω [rad/h] P o w e r [ M W ] Ω [rad/h] P o w e r [ M W ] Ω [rad/h] P o w e r [ M W ] (c)(a) (b) FIG. 12. Arnold tongues of three different TCL ensemblesof 10,000 units obtained by numerical simulations withouta controller. The contour maps represent the RMS errors intracking a sinusoidal input with frequency Ω. (a) The averagefrequency of the ensemble was ω = 4 . C = 10 kWh/ ◦ C and deadband δ = 0 . ◦ C.(b) The average frequency of the ensemble was ω = 17 . C = 2 . ◦ C anddeadband δ = 0 . ◦ C. (c) The average frequency of the ensem-ble was ω = 7 . C = 1 . ◦ C and deadband δ = 1 . ◦ C. The white lines ( −− ) cor-respond to the theoretical Arnold tongues with whole numberratios in Fig. 11. The power is normalized to 1 MW. The peakpower actually corresponds to 15 MW. integral (PI) controller is shown in Fig. 13a, and for thePRC-based controller the Arnold tongue is shown in Fig.13b. Ω [rad/h] P o w e r [ M W ] Ω [rad/h] P o w e r [ M W ] (b)(a) FIG. 13. Arnold tongues of TCL ensembles of 10,000 unitswith controllers. (a) With PI controller. (b) With PRC-based controller. The parameter in Table 1 were used withthe thermal capacitance C = 1 . ◦ C and the deadband δ = 1 . ◦ C. The power is normalized to 1. The peak powercorresponds to 15 MW.
The PI controller reduces the width of the 1:1 Arnoldtongue, and the PRC-based controller further reduces itand significantly minimizes its intensity, but the Arnoldtongue corresponding to the 1:2 entrainment becomes ap-parent. This once more confirms that it is possible to ob-tain significantly more demand response capacity from a given TCL population by using an appropriate controller(or control waveform).
D. Tracking of a Power Regulation Signal Based onSpectral Decomposition
In this section we demonstrate the tracking of an ACEsignal by using the Arnold tongues in Fig. 12 to deter-mine the appropriate spectral decomposition to apply tothe ACE signal. We identify each TCL ensemble by itsmean natural frequency ω . The ACE signal to track isshown in Fig. 14a and its wavelet transform is shownin Fig. 14b, in which we can see that the dominantfrequency is ω = 0 .
59 rad/h and the second dominantfrequency is around ω = 2 π rad/h. Magnitude Scalogram
Time [h] F r e qu e n c y (r a d / h ) M a gn it ud e Time [h] -300-200-1000100200300 P o w e r [ M W ] ACE signalFiltered ACE (b)(a)
FIG. 14. Area control error (ACE) signal and its powerspectrum. The power spectrum was computed using wavelettransform. (a) ACE signal with its lowpass filtered content.The cutoff frequency was 4 rad/h. (b) ACE frequency contentvs. time.
The ACE signal contains different frequency compo-nents (low and high frequencies) that can be decomposedinto different bands so that a certain population of TCLswith a given mean natural frequency is able to accuratelytrack the reference power. In Fig. 15, we show differentgroups of TCLs tracking the ACE signal that has beenfiltered in specific frequency bands and scaled such thatit can be tracked by the ensemble. Each TCL group iscomposed of 10 ,
000 heterogeneous units. For the TCLswith ω = 17 . = 17.1 rad/s, banpass [0,4] rad/s = 17.1 rad/s, banpass [0,8] rad/s = 17.1 rad/s, banpass [23,25] rad/s P [ M W ] = 4.3 rad/s, banpass [0,2] rad/s P [ M W ] = 7.9 rad/s, banpass [0,3] rad/s P [ M W ] = 7.9 rad/s, banpass [0,4] rad/s Time [h] -40-2002040 P [ M W ] Total power
ACE ref
Pagg (g) (f)(c) (e)(b)(a) (d)
FIG. 15. Tracking of the ACE by three different groups ofTCL populations with 10,000 units each. The groups areidentified by their mean natural frequencies ω correspondingto specific sets of parameters C and δ . The rest of the pa-rameters are the same as in Table 1. (a) Power tracking ofan ensemble of TCL with ω = 4 . ω c = 2 rad/h. (b) Power tracking ofan ensemble of TCL with ω = 7 . ω c = 3 rad/h. (c) Power tracking ofan ensemble of TCL with ω = 7 . ω c = 4 rad/h. (d) Power tracking ofan ensemble of TCL with ω = 17 . ω c = 4 rad/h. (e) Power tracking ofan ensemble of TCL with ω = 17 . ω c = 8 rad/h. (f) Power tracking ofan ensemble of TCL with ω = 17 . ω l = 23 and ω h = 25 rad/h. (g)The total power P agg is obtained as a sum of the power in (a),(b), (e) and (f). The ACE ref signal is the sum of the filteredACE signal used in (a), (b), (e) and (f) which almost recov-ers all the power content of the scaled down original signal inFig. 14. Comparing P agg and ACE shows that we are ableto recover most of the power spectrum of the ACE referencepower. In Fig. 16 we show the relative percent error of eachgroup tracking different frequency bands. The relative er-ror is computed as err ( t ) = ( P ref ( t ) − P agg ( t )) /P ref ( t ) × P base , which is the av- erage power consumption of 10,000 TCLs. The trackingRMSE (normalized by the average aggregate power) foreach case in Fig. 15 are given in Table II where the label-ing of TCL groups (a) to (g) corresponds to the labelingused in Fig. 15 and 16. The RMSE is computed asRMSE % = (cid:115) T (cid:82) T ( P ref ( t ) − P agg ( t )) dtP base × . e rr o r [ % ] (a) e rr o r [ % ] (b) e rr o r [ % ] (c) e rr o r [ % ] (d) e rr o r [ % ] (e) e rr o r [ % ] (f) Time [h] -20-10010 e rr o r [ % ] (g) hourly average Time [h] -0.500.5 c on t r o l s [ ° C ] (t) u i (t) (g) (h)(e) (f)(d)(c)(a) (b) FIG. 16. Relative tracking errors and sample TCL controlinputs. The legends (a) to (g) indicate the relative trackingerrors in percentage of the corresponding powers in Fig. 15(a) to (g). The hourly average error in (g) represents thecomputed average of the error using a sliding window of onehour. This shows that on an hourly average the relative erroris less than 10%. (h) Shows the global reference signal ξ ( t )sent by the BPA that requests of each TCL to change theirpower consumptions such the ACE signal can be tracked. Af-ter going through the local controllers (see Fig. 9), the TCLscontrol signals u i ( t ) are generated.TABLE II. Results on the tracking of the ACE regulationsignal TCL groups (a) (b) (c) (d) (e) (f) (g) P max [MW] 10.3 11.2 10.5 9.85 16.5 8.72 39.5RMSE % 2.31 1.92 3.09 1.38 2.92 3.39 5.89 The performance of the PRC-based controller is note-worthy given that communication between the BA andTCLs is open-loop. For comparison, consider for example3the tracking performances in [26, 36, 52]. The minimumvariance control law (MVC) in [36] achieves a trackingrelative error of less than 5% with RMSE < .
8% from 3 . P agg as feedbackin our controller, its peak relative error is reduced to ap-proximately 4% from approximately 8%. The control isthen of the form u k +1 = u k + I h [ I sgn( 12 − s k ) Z ( φ k )( ξ k − ξ P aggk ) − u k ] , where ξ P aggk is an appropriately scaled regulation signalthat depends on the aggregate power. In [26] four differ-ent methods were evaluated, two of which use state feed-back along with the aggregate power, and the other twoonly use the aggregate power in the construction of thecontrol signal. The two approaches that use state feed-back achieved the RMSE of 1 .
18% and 8 . ≤ . V. CONCLUSIONS AND FUTUREDIRECTIONS
We have presented a new modeling paradigm for ther-mostatically controlled loads using neuroscience-inspiredoscillator modeling and phase model reduction. The sim-plicity of the model opens the door to a rich variety ofanalysis tools and control strategies that could improvethe way in which TCLs are operated to provide demandresponse on a power distribution system. Successful con-trol policies should be able to provide demand responseon time-scales of interest while satisfying the TCLs’ oper-ating constraints and customer comfort. Specifically, thefrequency at which a load is made to switch ON and OFFin a given time period should be kept within a reasonablelimit to avoid excessive wear on equipment.We have developed a methodology for open-loop mod-ulation of TCL power consumption by using phase mod-els of the response to modulation of power, which only ad-vances or delays the phase of the TCL power cycle whenit is in a particular state, e.g., close to turning ON orturning OFF, as seen in Figure 5c. We have shown thattracking of aggregate power set-points requested by a sys-tem operator could be achieved using a PRC-based con-trol policy that prevents excessive synchronization andunwanted power fluctuations. Moreover, we have usedthe concept of the Arnold tongue to characterize howwell such unwanted effects are prevented. This approachcould be used to quantify the performance of control poli-cies that do not rely on communicating feedback of theTCL states to the balancing authority.The use of simplified phase models leads to controlpolicies that are robust to heterogeneity in system pa-rameters. Such models also simplify the potential for-mulation of optimization problems that can be solvedto estimate the maximum actuation of aggregated powerconsumption by TCLs given state constraints and dis-tribution subsystem size or structure. Solution of suchan optimization problem would quantify the demand re-sponse that can be obtained from a given collection ofTCLs while maintaining distribution system stability andpower quality. The application of phase models to evalu-ate the ability of random ON/OFF switching policies toprovide demand response on different time-scales is alsopromising.
Appendix A: TCL Model with Power andTemperature as Variables
Consider the model in (8) with the switching functiondefined in (9) together with the aggregate power equa-tion (3). Let y ( t ) = 1 /ηs ( t ) P be the instantaneous powerdrawn by a single TCL. The time derivative of the in-stantaneous power is given by ˙ y ( t ) = 1 /η ˙ s ( t ) P . Theswitching function in (9) can also be written using an ex-ponential function as s ( t ) = (1 + exp( − kx )) − . We can4then write x ( t ) in terms of the switching function as x ( t ) = − k ln (cid:18) − ss (cid:19) . (A1)Taking the time derivative of x ( t ) yields − k dxdt = ddt (ln(1 − s ) − ln( s )) . Let u = 1 − s , by applying the chain rule we get − k dxdt = ddt ln( u ) − ddt ln( s ) , = ddu ln( u ) dudt − dds ln( s ) dsdt , = 1 u dudt − s dsdt . After substituting for u and dudt = − dsdt , we obtain˙ x = ˙ sks (1 − s ) . (A2)We may now write (A2) in terms of the instantaneouspower y ( t ) as ˙ x = ˙ yk (1 − ηP y ) . (A3)Substituting for (A3) and (A1) in (8) with s ( t ) = ηP y , weobtain˙ y ( t ) = µk ( δ y −
13 ¯ y + θ − θ s )(1 − ηP y ) y, ˙ θ ( t ) = − CR ( θ − θ a + ηRy ) , (A4)where ¯ y ( t ) = − k ln( Pηy − Appendix B: Phase Model and Phase ResponseCurve
The phase reduction theory is a powerful tool forstudying multi-dimensional rhythmic systems that are re-duced to a scalar differential equation that is much easierto analyze and control. The autonomous oscillatory sys-tem is then described by its phase variable φ rotating ona circle S . This is represented by the phase equation˙ φ ( t ) = ω . In neuroscience the origin of the phase φ is de-fined as the time since the last spike of a neuron [21], andin our work we define it as the phase corresponding tothe time the TCL turns ON. When an oscillator receivesa pulse of strength A and duration ∆ T , the magnitude ofthe induced phase shift is given by PRC( φ ) = φ new − φ [21, 29].For completeness, we summarize the derivation of thephase model here. More details can be found in [21, 53].Consider a smooth dynamical system described by˙ x = f ( x, u ) , where the state variable x ( t ) ∈ R n and the input u ( t ) ∈U ⊆ R m . Suppose that the unforced system ˙ x = f ( x, ⊂ R n withperiod T . The limit cycle is then described by a non-constant periodic trajectory γ ( t ) = γ ( t + T ) ∈ χ, ∀ t ≥ δ ˙ x = A ( t ) δx ( t ) + B ( t ) u ( t ) , where A ( t ) = ∂f∂x ( γ ( t ) ,
0) and B ( t ) = ∂f∂u ( γ ( t ) ,
0) are T -periodic. Given that the limit cycle Γ is a one-dimensional closed curve [21], the position of any point x ∈ Γ can be uniquely described by a scalar phase φ ∈ S = [0 , π ) [29]. Let’s introduce the phase functionΘ( x ) that maps each point x on the limit cycle to itsphase φ = Θ( x ). The phase variable φ : R ≥ → S isdefined for each trajectory on the limit cycle as φ ( t ) = γ ( t + ω − φ ), and is periodic because of the periodicityof γ ( t ).From the linearized model and the asymptotic phasevariable, one can derive the phase-reduced model in aneighborhood of the limit cycle Γ for sufficiently smallinputs [21, 41]. By differentiating φ ( t ) with respect totime in the neighborhood of γ ( t ) using the chain rule,one obtains d Θ( x ) dt = ∂ Θ ∂x ( γ ( t )) · ddt x ( t ) + ∂ Θ ∂x ( γ ( t )) · B ( t ) u ( t ) , = ∂ Θ ∂x ( γ ( t )) · f ( x ) + ∂ Θ ∂x ( γ ( t )) · B ( t ) u ( t ) , = ω + Z ( φ ) · B ( t ) u ( t ) , where we have used to the fact that ∂ Θ ∂x ( γ ( t )) · f ( x ) = ω and Z ( φ ) = ∂ Θ ∂x ( γ ( t )) is the phase sensitivity functionalso referred to as infinitesimal PRC. The input matrixfunction B ( t ) depends of the differential equations de-scribing the system, for example B ( t ) = ( − µ,
0) for thesystem in (8) where the control is a perturbation of theset-point temperature θ s ( t ).In the following we review the main three methodsthat are commonly used to compute the phase sensitiv-ity function. These techniques are explained with greatdetails and illustrations in [21]. • Winfree’s ApproachIn a sufficiently small neighborhood of the limitcycle, the PRC scales linearly with respect to thestrength of the pulse. Hence one can writePRC( φ, A ) ≈ Z ( φ ) A, where Z ( φ ) = ∂ PRC( φ, A ) /∂A at A = 0 is thelinear response or sensitivity function that quanti-fies the small change in the instantaneous frequencycaused by the weak stimulus that was applied. Nowassume that we apply a sufficiently small stimu-lus (cid:15)p ( t ) and that the perturbed trajectory remainsnear the limit cycle attractor at all time. Replacingthe continuous input function (cid:15)p ( t ) with the equiv-alent pulse train of strength A = (cid:15)p ( t n ) h , where h t n is the timing of the n th pulse, one can write thePoincar´e phase map as φ ( t n +1 ) = { φ ( t n ) + Z ( φ ( t n )) (cid:15)p ( t n ) + h } mod T , in the form φ ( t n + h ) − φ ( t n ) h = 1 + Z ( φ ( t n ) (cid:15)p ( t n ) , which is a discrete version of˙ φ = 1 + (cid:15)Z ( φ ) · p ( t ) , (B1)in the limit h →
0. Note that the phase model(B1) is valid for any arbitrary input function p ( t ).To summarize, Winfree’s approach consists of mea-suring the phase shift induced by a pulse train todetermine the PRC. • Kuramoto ApproachKuramoto considered the unperturbed oscillatorwith φ ( x ) denoting the phases of points near itslimit cycle attractor. Differentiating φ ( x ) using thechain rule yields dφ ( x ) dt = grad φ · dxdt = grad φ · f ( x ) , where grad φ is the gradient of φ ( x ) with respect tothe state vector of the oscillator x ∈ R n . However,given that on the limit cycle the flow of the vectorfield f ( x ) is exactly in the direction of the periodicorbit so that dφ ( x ) dt = 1, we obtain the importantequality grad φ · f ( x ) = 1 . (B2)By applying the chain rule to the perturbed system dφ ( x ) dt = grad φ · dxdt , = grad φ · { f ( x ) + (cid:15)p ( t ) } , = grad φ · f ( x ) + (cid:15) grad φ · p ( t ) , and using (B2), the phase model is obtained as˙ φ = 1 + (cid:15) grad φ · p ( t ) . (B3)Kuramoto phase model (B3) and Winfree’s model(B1) are equivalent. Hence we have Z ( φ ) = grad φ . • Malkin’s ApproachHere we formally state Malkin’s theorem as in [21].Suppose the unperturbed oscillator has an expo-nentially stable limit cycle of period T . Its phaseevolution is described by˙ φ = 1 + (cid:15)Q ( φ ) · p ( t ) , (B4) where Q is a T -periodic function that is the solutionto the linear “adjoint” equation˙ Q = −{ Df ( x ( t )) T } Q, with Q (0) · f ( x (0)) = 1 , (B5)where Df ( x ( t )) T is the transposed Jacobian of theflow f at the point x ( t ) on the limit cycle, and thenormalization condition can be replaced by Q ( t ) · f ( x ( t )) = 1 , ∀ t. The phase models (B4) and (B3) or (B1) are equivalent,hence one can see that Z ( φ ) = grad φ ( x ) = Q ( φ ) . The method of the adjoint was used in this paper to nu-merically compute the phase sensitivity function Z ( φ ).Examples of computer codes for this method can befound in [21, 29]. Appendix C: Entrainment Region
The PRC defines the synchronization properties of anoscillator and the synchronized states as fixed points ofthe corresponding Poincar´e phase map [21]. The phe-nomenon of entrainment by weak forcing of limit-cycleoscillators can be modeled by˙ φ = ω + AZ ( φ ) v (Ω t ) , (C1)where ω and Ω are the natural frequencies of the oscilla-tor, respectively, and the forcing input is u ( t ) = Av (Ω t ),where v is 2 π -periodic with unit energy [47]. The re-gion of existence of a synchronized state is called Arnoldtongue [49, 54]. This region of phase-locked states on(Ω , A )-plane shrinks as the intensity A of the stimulusapproaches 0, with Ω the frequency of the stimulus. Ingeneral, m : n entrainment occurs when ω/ Ω ≈ n/m with positive relative prime integers n and m . This im-plies that the oscillator rotates exactly m times while theexternal forcing oscillates n times. Let ∆ = ω − mn Ω andby formal averaging we can write (C1) as dψdt = ∆ + A Γ m/n ( ψ ) , (C2)where ψ = φ − mn Ω t is a slow varying phase variable, andΓ m/n ( ψ ) is the interaction function determined by Z and v as Γ m/n ( ψ ) = 1 T ext (cid:90) T ext Z (cid:16) ψ + mn Ω t (cid:17) v (Ω t ) dt, = 12 π (cid:90) π Z ( ψ + mθ ) v ( nθ ) dψ, = 12 π (cid:104) Z ( ψ + mθ ) , v ( nθ ) (cid:105) , where T ext = π Ω is the period of the external forcing and θ = Ω tn ∈ [0 , π ). Without loss of generality, let A = 16and consider the 1 : 1 entrainment i.e., n = m = 1 andwrite Γ m/n ( ψ ) as simply Γ( ψ ). It can be shown thatwhen the conditionmin Γ( ψ ) < − ∆ < max Γ( ψ ) , is satisfied, (C2) has at least two fixed points at which dψ ( t ) /dt = 0 holds, and one of them is stable [29, 47].The interval ∆ of phase locking for a fixed input strength A decreases as A →
0. So, for different values of A andentrainment ratios n/m one obtains the Arnold tonguesas shown in Fig. 11. ACKNOWLEDGMENT
The authors wish to thank Scott Backhaus, MichaelChertkov, and Sean Meyn for valuable discussions. Thiswork was carried out under the auspices of the NationalNuclear Security Administration of the U.S. Departmentof Energy at Los Alamos National Laboratory underContract No. DE-AC52-06NA25396 with the support ofthe Advanced Grid Modeling Research Program in theU.S. Department of Energy Office of Electricity. [1] E. C. Kara, M. D. Tabone, J. S. MacDonald, D. S. Call-away, and S. Kiliccote, in
Proceedings of the 1st ACMConference on Embedded Systems for Energy-EfficientBuildings , BuildSys ’14 (ACM, New York, NY, USA,2014) pp. 140–147.[2] J. Yan, Y. Zhai, P. Wijayatunga, A. M. Mohamed, andP. E. Campana, Applied Energy , 241 (2017).[3] E. C. Kara, M. Berg´es, and G. Hug, IEEE Transactionson Smart Grid , 2560 (2015).[4] S. P. Meyn, P. Barooah, A. Buˇsi´c, Y. Chen, andJ. Ehren, IEEE Transactions on Automatic Control ,2847 (2015).[5] S. Behboodi, D. P. Chassin, N. Djilali,and C. Crawford, Applied Energy (2017),https://doi.org/10.1016/j.apenergy.2017.07.058.[6] C. Sijie and L. Chen-Ching, Journal of Modern PowerSystems and Clean Energy , 10 (2017).[7] R. F. Bischke and R. A. Sella, IEEE Transactions onPower Apparatus and Systems PAS-104 , 1290 (1985).[8] R. M. Delgado, Proceedings of the IEEE , 1471 (1985).[9] S. Nolan and M. OMalley, Applied Energy , 1 (2015).[10] P. Siano, Renewable and Sustainable Energy Reviews ,461 (2014).[11] I. Pavi´c, T. Capuder, and I. Kuzle, Applied Energy ,724 (2016).[12] F. Teng, Y. Mu, H. Jia, J. Wu, P. Zeng, and G. Strbac,Applied Energy , 353 (2017).[13] A. Malik and J. Ravishankar, Applied Energy (2017),https://doi.org/10.1016/j.apenergy.2017.08.160.[14] V. Lakshmanan, M. Marinelli, J. Hu, and H. W. Bind-ner, Applied Energy , 470 (2016).[15] V. Lakshmanan, M. Marinelli, A. M. Kosek, P. B. Nrgrd,and H. W. Bindner, Energy , 705 (2016).[16] Y. Zhou, C. Wang, J. Wu, J. Wang, M. Cheng, andG. Li, Applied Energy , 456 (2017).[17] A. Ghaffari, S. Moura, and M. Krsti´c, Journal of Dy-namic Systems, Measurement, and Control , 1213(2015).[18] K. Ma, C. Yuan, J. Yang, Z. Liu, and X. Guan, Energies , 953 (2017).[19] B. Baeten, F. Rogiers, and L. Helsen, Applied Energy , 184 (2017).[20] A. Aryandoust and J. Lilliestam, Applied Energy ,749 (2017). [21] E. M. Izhikevich, Dynamical systems in neuroscience (MIT press, 2007).[22] P. Barooah, A. Buˇsi´c, and S. Meyn, in
System Sciences(HICSS), 2015 48th Hawaii International Conference on (IEEE, 2015) pp. 2700–2709.[23] N. A. Sinitsyn, S. Kundu, and S. Backhaus, Energy Con-version and Management , 297 (2013).[24] N. Mehta, N. A. Sinitsyn, S. Backhaus, and B. C. Lesieu-tre, Energy Conversion and Management , 784 (2014).[25] H. Hao, B. M. Sanandaji, K. Poolla, and T. L. Vincent,IEEE Transactions on Power Systems , 189 (2015).[26] E. Vrettos, S. Koch, and G. Andersson, in InnovativeSmart Grid Technologies (ISGT Europe), 2012 3rd IEEEPES International Conference and Exhibition on (IEEE,2012) pp. 1–8.[27] D. Docimo and H. K. Fathy, Journal of Dynamic Sys-tems, Measurement, and Control , 101009 (2017).[28] M. Chertkov and V. Chernyak, arXiv preprintarXiv:1701.04939 (2017).[29] H. Nakao, Contemporary Physics , 188 (2016),http://dx.doi.org/10.1080/00107514.2015.1094987.[30] F. D¨orfler and F. Bullo, SIAM Journal onControl and Optimization , 1616 (2012),http://dx.doi.org/10.1137/110851584.[31] B. Ermentrout, Neural Computation , 979 (1996).[32] A. Zlotnik and J.-S. Li, Journal of Neural Engineering ,046015 (2012).[33] S. Ihara and F. C. Schweppe, IEEE Transactions onPower Apparatus and Systems PAS-100 , 4142 (1981).[34] C. Perfumo, E. Kofman, J. H. Braslavsky, and J. K.Ward, Energy Conversion and Management , 36(2012).[35] S. Bashash and H. K. Fathy, in Proceedings of the 2011American Control Conference (2011) pp. 4546–4553.[36] D. S. Callaway, Energy Conversion and Management ,1389 (2009).[37] B. Ermentrout, Simulating, analyzing, and animating dy-namical systems: a guide to XPPAUT for researchersand students (SIAM, 2002).[38] D. S. Callaway and I. A. Hiskens, Proceedings of theIEEE , 184 (2011).[39] S. Kundu, N. Sinitsyn, S. Backhaus, and I. Hiskens, 17thPower System Computation Conference (2011). [40] T. Harada, H.-A. Tanaka, M. J. Hankins, and I. Z. Kiss,Physical review letters , 088301 (2010).[41] D. Efimov, P. Sacr´e, and R. Sepulchre, in Decision andControl, 2009 held jointly with the 2009 28th ChineseControl Conference. CDC/CCC 2009. Proceedings of the48th IEEE Conference on (IEEE, 2009) pp. 7692–7697.[42] I. Dasanayake and J.-S. Li, Physical Review E , 061916(2011).[43] W. Zhang, J. Lian, C.-Y. Chang, and K. Kalsi, IEEEtransactions on power systems , 4655 (2013).[44] L. C. Totu and R. Wisniewski, IFAC Proceedings Vol-umes , 9956 (2014).[45] B. P. Administration, “Bpa area control error (ace) an-nual reports,”.[46] A. Zlotnik, R. Nagao, and I. Z. K. J.-S. Li, Nature com-munications (2016).[47] H.-A. Tanaka, Physica D: Nonlinear Phenomena , 1(2014). [48] S. Shirasaka, W. Kurebayashi, and H. Nakao, PhysicalReview E , 012212 (2017).[49] A. Granada, R. Hennig, B. Ronacher, A. Kramer, andH. Herzel, Methods in enzymology , 1 (2009).[50] H.-A. Tanaka, I. Nishikawa, J. Kurths, Y. Chen, andI. Z. Kiss, EPL (Europhysics Letters) , 50007 (2015).[51] A. Zlotnik and J.-S. Li, SIAM Journal on Applied Dy-namical Systems , 1654 (2014).[52] J. L. M. Stephan Koch and D. S. Callaway, 17th PowerSystem Computation Conference (2011).[53] D. Efimov, P. Sacr´e, and R. Sepulchre, in Decision andControl, 2009 held jointly with the 2009 28th ChineseControl Conference. CDC/CCC 2009. Proceedings of the48th IEEE Conference on (2009) pp. 7692–7697.[54] M. J. Schaus and J. Moehlis, in