Access the full text.
Sign up today, get DeepDyve free for 14 days.
Loading next page...
/lp/arxiv-cornell-university/super-relaxation-of-space-time-quantized-ensemble-of-energy-loads-to-RhUYlFrF9g
- ISSN
- 0306-2619
- eISSN
- ARCH-3344
- DOI
- 10.1016/j.apenergy.2020.116419
- Publisher site
- See Article on Publisher Site
Abstract
Ensembles of thermostatically controlled loads (TCL) provide a signi cant de- mand response reserve for the system operator to balance power grids. However, this also results in the parasitic synchronization of individual devices within the ensemble leading to long post-demand-response oscillations in the integrated en- ergy consumption of the ensemble. The synchronization is eventually destructed by uctuations, thus leading to the (pre-demand response) steady state; how- ever, this natural desynchronization, or relaxation to a statistically steady state, is too long. A resolution of this problem consists in measuring the ensemble's instantaneous consumption and using it as a feedback to stochastic switching of the ensemble's devices between on- and o - states. A simpli ed continuous-time model showed that carefully tuned nonlinear feedback results in a fast (super-) relaxation of the ensemble energy consumption. Since both state information and control signals are discrete, the actual TCL devices operation is space-time quantized, and this must be considered for realistic TCL ensemble modelling. Here, assuming that states are characterized by indoor temperature (quantify- ing comfort) and air conditioner regime (on, o ), we construct a discrete model Corresponding author Email address: h.ouerdane@skoltech.ru (H. Ouerdane) Preprint submitted to Applied Energy January 5, 2021 arXiv:2008.03118v2 [eess.SY] 4 Jan 2021 based on the probabilistic description of state transitions. We demonstrate that super-relaxation holds in such a more realistic setting, and that while it is stable against randomness in the stochastic matrix of the quantized model, it remains sensitive to the time discretization scheme. Aiming to achieve a balance be- tween super-relaxation and customer's comfort, we analyze the dependence of super-relaxation on details of the space-time quantization, and provide a simple analytical criterion to avoid undesirable oscillations in consumption. Keywords: Demand response; thermostatically controlled loads; energy consumption dynamics 1. Introduction Power grids of today are uncertain with the major sources of uncertainty being uctuations of renewables, especially of wind and solar, and market un- certainty. To deal with the uncertainties, grid operators need new exible and inexpensive resources. Demand response (DR) came up prominently as a way if not to resolve the problem completely, then at least to reduce its consequences [1]. The main idea consists in exploiting the fact that many consumers of elec- tricity, also called loads, can tolerate delays provided that their comfort zone is not violated. The complexity of the power system and the electricity markets dynamics call for the development of practical approaches capable to implement DR while accounting for the increasing penetration of the renewables. Demand side man- agement and automatic meter management systems on the customers' side, and on the supply side remote control of power ows, have been proposed to facil- itate renewables' integration [2]. To mitigate risks associated with imbalances of the supply and demand, balancing demand response on an hourly basis was suggested [3]. A market-based approach with correct price signals, fostering exible, smart power system has been discussed as a way forward to integrate large shares of renewables [4]. Further, while DR presents bene ts and costs [5], which depend on the technology used for control schemes [6], it can reduce 2 the market prices and consumers' cost owing to its contribution to the capacity market [7]. Moderate energy consumers like the service sector may participate in DR if barriers (e.g., restructuring costs and regulations) are mitigated and drivers (e.g., positive public image) enhanced [8]. As power generation a ected by uctuating renewables, may result in transmission line overload, electricity grids' topology must be designed to sustain such overloads [9]. Electric vehi- cles are increasingly suggested as exible grid components that permit stability [10, 11] while partaking in DR. But in case of contingency, distributed DR scheduling can help with frequency regulation [12]. While involving big stable loads, like aluminium smelters, in DR services is a well established practice, there is also great potential in utilizing opportunities in DR which can be o ered by many small loads. To unlock this potential, a statistics-based approach as pioneered in, e.g., [13], is necessary as shown by a large body of works. Models within this approach quickly evolved to account for characteristics such as lifestyle and weather [14], and propose a methodol- ogy for classi cation of elementary component loads and aggregation [15]. Us- ing Fokker-Planck equations [16, 17], large aggregates of loads controlled with thermostats such as heaters and air conditioners, were considered given their important impact on the power system dynamics [18, 19]. Interestingly, as the dynamics of electrical heaters and air conditioners can be characterized as an alternating renewal process, it was shown how consumption data can be used to identify electric load models [20]. Considering aggregated thermostatically con- trolled loads (TCL), a state-queuing model revealed the important in uence of load state diversity on the aggregated pro le dynamics, as synchronization leads to unwanted peak loads [21, 22]. While a direct way to achieve peak shaving is by interruption of power delivery to the loads, other, smarter ways propose to control the loads consumption by variation of the temperature set points thus permitting exibility on fast time scales, and ensemble diversity [23, 24, 25]. To ensure ecacy of load management to lower TCLs' aggregated power consump- tion when needed, by feedback control, a model-based feedback control strategy must aim for high accuracy in the characterization of the aggregate dynamics 3 [26], as control errors may occur with conventional thermostats [27]. To mit- igate the negative e ects of power uctuations in the distribution networks it was shown how to coordinate multiple TCL groups employing a two-stage op- timization model applied in real-time [28]. Recent works also considered DR as a perturbation of the TCL ensemble dynamics driving it away from its steady state [29, 30, 31]. Relaxation after perturbation is of importance to ensure the stability of the power system. This article contributes this later line of work. Several hurdles must be overcome to make the DR contribution of many small loads meaningful. It is not economically viable to expect a small load, e.g. a thermostatically controlled loadlike air-conditioner or heater, to be en- gaged in a sophisticated individual control. Instead, aggregation of many small loads would be a preferred solution [32]. In this scheme the aggregator is an authority receiving DR requests from the system operator and broadcasting the same signal to all their consumers. It is assumed that the consumers obey and perform the requested action, that is switch o or switch on, follows when requested. An unfortunate side e ect of all consumers following the same sig- nal is a parasitic synchronization/oscillations seen long after engagement of the ensemble in the DR [23]. Consumer-speci c uctuations will lead, eventually, through mixing to a decay of oscillation (de-synchronization). However, natural mixing is typically weak, thus leading to long transients, delaying availability of the ensemble for the next DR session. As shown in [29], the randomization of switching, implemented through the broadcast of a Poisson rate of the switch on/o delay, helps to reduce the mixing time while also providing an accept- able \comfort zone" guaranteed to loads. Diversity of loads contributing to the ensemble helps to reduce the mixing time even further [30]. The solution suggested in [29, 30] did not depend on any knowledge of the current system state (temperature and switch on/o status). The next signif- icant step in improving control of the ensemble was made in [31], where the following question was addressed: is it possible to set up a viable aggregation model that would rely only on receiving instantaneous integrated consumption of the entire ensemble as a feedback? Notice that even though the absence of 4 the individual response of a load makes the problem of organizing the aggrega- tor control harder, the ability to receive one signal, integrated over the entire ensemble, makes the approach desirable from the viewpoint of keeping the con- sumption of individual loads private. It was shown in [31] that the question just posed has an armative answer: making nonlinear feedback on the instan- taneous integrated consumption of the ensemble allows to accelerate relaxation (de-synchronization) of the integrated consumption to the steady-state. No- tice that this approach, coined the \mean- eld" control in reference to related methods originating from plasma physics, control, management sciences and applied mathematics [33, 34, 35, 36, 37], has this strong e ect, dubbed super- relaxation [31], only on a specially selected expectation over the ensemble's probability distribution (mean instantaneous consumption) while other expec- tations over instantaneous probability distribution over the ensemble, continue to relax slowly. The model in [31] assumed continuous temperature variation and time but operations of the actual TCL devices are space-time quantized. In this work we develop a space-time-quantized model of the TCL ensemble, which incorpo- rates mean- eld control, that is feedback on instantaneous total consumption, of the switch on/o rates of all the loads of the ensemble. We show that the super-relaxation e ect is also observed in the space-time quantized model, better representing the real-world of energy management than the continuous model studied before. We experiment with the model parameters { the size of space- time quantization steps and degree of the mean- eld control nonlinearity in the Poisson switching on/o rates { to make a recommendation on the choice of parameters achieving a reasonable balance between fast mixing of the ensemble (faster post-DR restoration) and \comfort zone" of the consumers. The article is organized as follows. In Section 2, we derive and describe the basic equations of our space-time quantized model generalizing the space-time continuous model of [31]. Section 3 is devoted to discussion of the numerical results and of the insight they provide. Section 4 is reserved for conclusions and discussion of the path forward. Technical details are presented in Appendices. 5 2. Problem formulation 2.1. Space-time continuous TCL model We start with an overview of the basic elements of the continuous model [31]. Assume that at every moment t each TCL load is characterized by two pa- rameters: (a) consumer's instantaneous temperature x(t) and, (b) binary state, (t) = 0; 1, characterizing the on or o state of a consumer's thermal (heating or cooling) device. The dynamics of each TCL in the phase space, character- ized by the tuple fx(t); (t)g, can be complex as it depends on various factors such as operating power, desired temperature, outside temperature, as well as on the local level of noise and uncertainty associated with details of the consumer's operation regime (e.g. frequency of the doors or windows opening, trac through the consumer space, etc). To manage this complexity, we con- sider the following set of simplifying assumptions (also focusing without loss of generality on air-conditioning, thus cooling, as our enabling example): i/ When the device is switched on, the temperature decreases, and the tem- perature raises when the device is switched o . We assume that the re- laxation of x(t) is linear in both switch on and switch o regimes with the relaxation rates equal to each other by the absolute value. Linearity of x(t) is justi ed in our work by the assumption that the outside tempera- ture x and the minimal possible temperature x that can be achieved by the air conditioner being constantly in the on position, are both far apart from the indoor temperature perceived as comfortable by the consumers. ii/ All TCL devices and their settings are identical, both in terms of their relaxation rates and the temperature extent of the comfort zone. iii/ Stochastic e ects, associated with device-speci c uncertainties, are as- sumed small and are thus neglected. iv/ A TCL does not switch on or o immediately after crossing either of the boundaries of the comfort zone. The switching is delayed according to a 6 Poisson process with rate, r. It is assumed that the operator broadcasts the same r to all the consumers. Considered within these assumptions, the basic model of the continuous-time TCL dynamics in the (x; ) space is described by the following set of equations [13, 14]: ; = "; dx = (1) dt ; = #; #; with probability rdt otherwise (t) and x < x ; (t + dt) = (2) "; with probability rdt otherwise (t) and x > x ; where are the cooling/heating rates, and r is the constant rate of (Poisson) switching (from on to o and vice versa) de ning switching delay after x(t) crosses the threshold temperature, x or x . Then, the two-dimensional proba- # " bility distribution vector, P (xjt), satis es the following Fokker-Planck equation: 0 0 1 1 0 1 1 0 P (xjt) : " @ @ A A @ A @ L P (xjt) = 0; P (xjt) = ; (3) 0 1 P (xjt) 0 1 0 1 1 0 (x x) (x x ) : # " @ A @ A L = @ r ; (4) 0 1 (x x) (x x ) # " where is the Heaviside step function and P ; P are the two components of " # P (xjt), corresponding to the probability distributions for a consumer to be in the switched-on and switched-o states, respectively. This basic model and generalizations were discussed extensively in [29, 30]. To complete the model one needs to describe actions of the aggregator. We assume that the aggregator has instant access to the integrated consumption of the ensemble. This is realistic in the case when all participants of the ensem- ble are collocated geographically within the power distribution system, i.e., all reside in the same power distribution feeder area. Then we measure the inte- grated consumption with a physical device sitting at the sub-station connecting 7 the feeder to the rest of the power system. The instantaneous aggregated con- sumption of the ensemble is U (t) = dxP (xjt), where the integral accounts for the number (proportion) of consumers which are switched on at time t. Assum- ing that all consumers are of the same type, e.g. similar ats or houses with a similar set of devices, we switch to dimensionless characteristics where the power consumed by an individual participant of the ensemble is unity. We then assume that the signal representing the Poisson switching on/o rate, q(t), sent by the aggregator to individual consumers is a functional of U (t): q(t) = C [U (t)], and hence acquires a dynamical character. This type of control is called the "mean- eld control" [33, 34, 35, 36, 37, 38] because it involves feedback (in choosing the rate r) on the global measured quantity, which is instantaneous integrated consumption of the ensemble. A particular form of the Poisson rate dependence on the integrated consumption, C [U ] = r(2U ) , where r is the basic rate intro- duced in Eq. (2), and the parameter s controls the degree of nonlinearity, was considered in [31]. 2.2. Space-time quantized TCL model We now proceed with the details of the space-time-quantized version of the continuous model summarized in Eq. (3). Using the same set of assumptions as in the continuous model, we bin the temperature range and denote the quan- tized space states, . Then, we consider transitions from state to state in discrete time. This space-time-quantized description re ects realistic practice of built-in controls of practical (small and inexpensive) TCLs. In the simplest approach with no feedback (no mean- eld control) the transition probability matrix, describing probability for a load to transition from the state to the state during the discrete period of time t, reads: (0) " # p 0 = p + rp + rp (5) 0 0 0 (0) where p describes the transition from to , associated with cyclic evolution as if it would occur exactly at the thresholds (immediately after entering the # " discomfort zone) and rp and rp are corrections due to the Poisson delay 0 0 8 (0) in switching between the on and o state. The term p also includes di usion which is described by random transitions to neighboring nodes with probability . The physical meaning of can be explained as follows. In a typical real-life situation, x(t) uctuates due to uncertainties, e.g. linked to incidental move- ment of people inside, door and windows opening and closing and related. These uctuations can be simulated by a white noise acting on x(t). The white noise contributes the di usion term in the Fokker-Planck equation. (0) # " The matrices p , p and p can be graphically represented as shown 0 0 0 in Figs. 1, 2, and 3 respectively. Notice that the transition probability matrix p de ned by Eq. (5) is stochastic, i.e. p = 1: (6) Transitions between states are then governed by the following time-space-discrete master equation: (t + 1) = p 0 0 (t) (7) where (t) is a probability mass function, which stands for the probability of a TCL to be in the state at time t. x x x x − ↓ ↑ + 1 − 1 − 2 1 − 2 1 − 2 1 − 2 1 − 2 1 − 2 1 − 2 0 0 0 0 0 0 0 0 0 X X X X X X X X X 1 2 3 4 5 6 7 8 9 1 1 1 1 1 1 1 1 1 1 1 X X X X X X X X X 1 2 3 4 5 6 7 8 9 1 − 2 1 − 2 1 − 2 1 − 2 1 − 2 1 − 2 1 − 2 1 − (0) j Figure 1: Graphical representation of p ; X denotes the state with a particular temperature and regime of the air-conditioner (on, o ) where i is the temperature and j is the operation regime of the air-conditioner. The arrows denote the possible transitions with the associated non-zero probability. This part of the transition matrix governs the transitions without Pois- son switchings in the out-of-comfort zone. The space between x and x is the comfort zone; # " x and x are points where the load must turn to another state. The same-state transitions and the two-step transitions are characterized by a di usion rate . 9 x x x x − ↓ ↑ + 0 0 0 0 0 0 0 0 0 X X X X X X X X X 1 2 3 4 5 6 7 8 9 −1 −1 1 1 1 1 1 1 1 1 1 1 1 X X X X X X X X X 1 2 3 4 5 6 7 8 9 Figure 2: Graphical representation of p , which governs the Poisson switchings from on to o . Negative elements ensure the preservation of the stochastic property of the resulting transition matrix. x x x x − ↓ ↑ + 0 0 0 0 0 0 0 0 0 X X X X X X X X X 1 2 3 4 5 6 7 8 9 1 1 −1 −1 1 1 1 1 1 1 1 1 1 X X X X X X X X X 1 2 3 4 5 6 7 8 9 Figure 3: Graphical representation of p , which governs the Poisson switchings from o to on. Mean- eld control amounts to allowing the switching rate to be dependent on the energy consumption of a device in a particular state averaged over the probability mass function N (t) = (t)U ; (8) where 1; if 2 set of on states U = (9) 0; otherwise: With this de nition, N (t) can also be understood as the fraction of loads 10 switched on at time t, and we may then generalize Eq. (5) as follows (0) " # p 0 (t) = p + q (t)p + q (t)p ; 0 0 0 " # q (t) = f (r[2N (t)] ) ; # " q (t) = f (r[2(1 N (t))] ) ; (10) " " where q (t) and q (t) are the Poisson rates modi ed by the mean- eld control, " # via the function f explicitly de ned further below, and denotes the degree of nonlinearity. The corresponding master equation takes a similar form as Eq. (7): 0 0 (t + 1) = p (t) (t): (11) According to the graphical representation of the matrix p (t) in Fig. 4, each particular load may experience four types of transition while in the out-of- comfort zone: i/ it may remain in the same state with probability ; ii/ it may go one step deeper in the out-of-comfort zone with probability 1 2 q (t) "=# (where for ease of notation " = # means either " or #); iii/ it may go two steps deeper in the out-of-comfort zone with probability ; iv/ it may switch state from on (resp. o ) to o (resp. on) with probability q (t) (resp. q (t)). Each # " particular probability must be non-negative; so the Poisson rates must satisfy q (t) 1 2. Consequently, f (x) is restricted to the [0; 1 2] interval. "=# Acknowledging that many choices are possible, we choose to work with the following form of the saturation function: x; x < 1 2; f (x) = (12) 1 2; otherwise Note that the discrete schemes described by the transition matrices Eqs. (5) and (10) have a proper continuous limit, as the corresponding master equations transform into Fokker-Planck equations discussed in [29, 30, 31] in this limit. (See Appendix Appendix A for details.) To measure the system evolution toward steady state, we use two quantities: (st) H (t) = k (t)k , which is the L distance describing how the probability 1 1 11 (st) mass function (t) goes toward its steady state; and the jN (t) N j, which provides a measure of how the total energy consumption of the ensemble goes toward its steady-state value set by the aggregator. We show below that in the case of the mean- eld control the rate of the two quantities relaxation to the (st) steady state may be dramatically di erent. Speci cally, jN (t) N j may converge to 0 much faster than H (t). x x x x − ↓ ↑ + 1 − − q 1 − 2 − q ↓ ↓ 1 − 2 1 − 2 1 − 2 1 − 2 1 − 2 1 − 2 0 0 0 0 0 0 0 0 0 X X X X X X X X X 1 2 3 4 5 6 7 8 9 q q q q ↓ ↓ ↑ ↑ 1 1 1 1 1 1 1 1 1 1 1 X X X X X X X X X 1 2 3 4 5 6 7 8 9 1 − 2 1 − 2 1 − 2 1 − 2 1 − 2 1 − 2 1 − 2 − q 1 − − q ↑ ↑ (0) # " Figure 4: Full transition matrix p 0 (t) = p + q (t)p + q (t)p , 0 q (t) 1 2, 0 # 0 " 0 # 0 q (t) 1 2. 2.3. Linear analysis of the decay rate Standard eigenvalue analysis of the linear master equation Eq. (7) (with the constant switching rates) shows that there is a unique maximal eigenvalue equal (st) (st) to unity if > 0, so that the steady state is unique. Then, if N = (this is proven in Appendix Appendix B), the nonlinear master equation Eq. (11) also has the same steady-state, which is unique in a small neighborhood. In our numerical simulations, we did not encounter other steady states. To see how fast the system approaches its steady state with the mean- eld control, we proceed with the eigenvalue analysis of Eq. (11). (st) Applying the decomposition, (t) = + (t), and keeping only the linear term we arrive at 0 0 (t + 1) = S (t); X X (0) (st) (st) # " # " S 0 = p + rp + rp + 2 rU 0 p 2 rU 0 p (13); 0 0 0 00 00 00 00 00 00 12 (st) where we have used the equality N = 1=2. The transition matrix S de ned in Eq. (13) can be split in two parts: p, which is the transition matrix of the ensemble without mean- eld control, Eq. (5), and the term V , which can be 0 0 0 treated as a perturbation. The matrix elements S , p , and V read: S 0 = p 0 + V 0; (0) # " p 0 = p + rp + rp ; 0 0 0 X X # (st) " (st) V 0 = 2 rU 0 p 2 rU 0 p : (14) 00 00 00 00 00 00 where r is constant. The spectral decomposition of the transition matrix S yields: (i) (i) (i) S = ; (15) n o n o (i) (i) (i) where is the set of eigenvalues, and and are respectively (i) (j) the right eigenvectors and the left eigenvectors sets such that = . ij The time evolution of the perturbation (t) is then given by: XX (i) (i) (i) (t) = (0): (16) (i) Note that the matrix S is a real matrix, so if is an eigenvalue of S then (i) is also an eigenvalue i.e. all complex eigenvalues are paired. As shown in Appendix Appendix D, there is at least one distinct eigenvalue (0) (0) = 1 whose corresponding eigenvector, , characterizes a mode that does not decay towards the steady state with time, but whose amplitude is always equal to 0, in order to satisfy the normalization condition: = 1. The (i) (i) other modes associated with the right eigenvectors (i 6= 0) decay as . It is convenient to introduce the relaxation constant of the mode as the complex (i) (i) number = log(j j) + i arg . With these notations we rewrite the 13 dynamics of the perturbation as: XX (i) (i) (t) = exp [ Re( )t + iIm( )t] 0 (0): (17) i i The relaxation constants introduced here are the space-time-quantized ana- logues to the Fokker-Planck operator's eigenvalues of the continuous models [29, 30, 31]. From now on we focus only on the exponential rates with the smallest real part, dominating the long-time relaxation. In [30], authors found two classes of modes when analysing the continuous version of the model with mean eld control: one that does not contribute at all to the total energy N and one that (i) (i) (i) does. Following this idea we sort our modes into two families u and v : (i) (i) The family of vectors u for which we have V 0u = 0 () 0 0 (i) U u = 0, where i denotes the label of eigenvectors in the set, which we call the ghost family, fGFg, as vectors of this set do not contribute, after the DR perturbation is applied, neither to the energy consumption U (or equivalently, N ) nor to respective relaxation. However, and unless degeneracy, these modes contribute to other observables, in particular the H { distance between steady-state and the current time state. (i) (i) The family of vectors v for which we have V v 6= 0, which we 0 0 call the signi cant family, fSFg, as it in uences the energy consumption and its relaxation, contributing as well to the relaxation of the entire ensemble. The whole family, fWFg, of eigenvalues , is the union of the ghost family and the signi cant family of the eigenvalues: fWFg = fSFg[fGFg. A similar classi cation of eigenvalues of the Fokker-Planck operator eigenvalues was done for the continuous models [30, 31]. Using the de nitions above we can now introduce the relaxation rate as min fReg and reintroduce the relaxation 2fSFg= rate for the entire ensemble as min fReg. The detailed comparison of 2fWFg= these relaxation rates is performed in the next Section. 14 3. Numerical Results and Discussion 3.1. Relaxation Constants 0.40 0.4 0.35 0.30 0.3 0.25 λ λ 0.20 0.0 2.5 5.0 7.5 10.0 0.2 0.15 0.10 0.1 0.05 λ λ 1 3 λ λ 2 4 0.00 0.0 0 10 20 30 40 0 10 20 30 40 α α Figure 5: Real and imaginary parts, Re( ) and Im( ), of the rst 4 relaxation constants as i i functions of the degree of nonlinearity are shown for r = 0:1. The number of states in the comfort and out-of-comfort zones are n = 12 and n = 18, respectively. out in We compute leading relaxation constants numerically. The behavior of the rst four constants, i.e. those with the smallest real parts, Re( ), of the whole set excluding = 0, is shown in Fig. 5 as functions of the strength of the mean- eld signal, , for two di erent values of the Poisson rate, r. Eigen- values associated with the ghost eigenvectors (even indexes) do not depend on by de nition, so the corresponding relaxation rates are also -independent. Hence, only the eigenvalues of the signi cant family contribute relaxation of the consumption. The jump of Im( ), seen on the right panel, is due to the fact the the imaginary part is de ned modulo 2. The blue stripe on the left panel marks the domain where total consumption of the signi cant ensemble and the whole ensemble show di erent relaxation rates at r = 0:1, i.e. Re( ) < Re( ). The 2 1 inset is a magni ed view of the [0; 10] [0; 0:75] domain in the f ; Re()g plane, Re( ) Im( ) with crosses mark intersection where the relaxation rates of the signi cant and ghost families start to deviate. Such domain does not exist at r = 0:3. 3.2. Super-relaxation in space-time quantized model r = 0.05, α = 10 (st) |N − N | − 1 min{Reλ } {EF} − 3 min{Reλ } {WF} − 5 − 7 − 9 − 11 − 13 0 100 200 300 400 500 600 Figure 6: Illustration of the super-relaxation, with = 0:05, = 10 and r = 0:05 by compar- (st) (0) ison of jN (t) N j and H = k (t) k , re ecting how the whole ensemble relaxes 1 1 to its steady state. The dashed-doted curves are the relaxation rates obtained by spectral decomposition. Both quantities decay as e with di erent . Super-relaxation is essentially the fast relaxation of the total consumption N (t) while the L -distance H (t) is much slower to reach steady state. Since " 1 only the signi cant family fSFg and its set of eigenvalues govern the relaxation of consumption, we may derive a criterion for the super-relaxation. In the general case, the relaxation rate of the ensemble takes the value min Re() as it yields the fastest characteristic decay time, while the relaxation rate for N takes a di erent value: min Re(). Mismatch between the two minima, 2fEFg G, is called the "gap": G = minfRe()g minfRe()g : (18) fSFg fWFg 16 n = 30, n = 18 out ° 1 24 0.8 ° 2 0.6 0.4 ° 3 0.2 1.0 0.8 ° 4 0.6 0.0 10 0 10 20 30 40 50 10 r 0.4 0.2 0.0 Figure 7: The left panel shows the gap G characterizing the super-relaxation regime, for di erent values of r and , number of states in the comfort zone n , number of states out of in the comfort zone n . The blue line marks the frontier between standard/super relaxation out areas. The gap is zero across the whole white area. The right panel displays areas with z coordinates that indicate the relaxation areas for various numbers of states in the out-of- comfort zone varying from n = 8 to n = 24 with the step n = 4, while the number out out out of states in the comfort zones is constant with n = 12 in this example. Plots in both panels in were produced with a di usion coecient = 0:05. The gap determines the relaxation regime: standard or super-relaxation. If G = 0 the system dynamics follows the standard regime; if G > 0 the system under- goes super-relaxation. To better understand peculiarities of the super-relaxation regime in the space-time-quantized framework, we compute the \phase diagram" of the gap in the f ; rg plane, using Eq. (18). An illustrative example of the (st) dynamics with super-relaxation is shown in Fig. 6: jN N j goes to zero faster than H does; we also see that minfRe()g and minfRe()g 1 fSFg fWFg have di erent slopes ( and respectively), which means that the gap G is 1 2 nonzero. The two curves may cross each other thus closing their gap at some point in time. The particular point when the gap is zero (no super-relaxation) depends on both model parameters r and . We analyze this further by calcu- lating the phase diagram of G, showing two possible phases: standard relaxation and super-relaxation in Fig. 7. We denote n the total number of states in up (down) position, and n , out out r the number of states in the out-of-comfort zone in up (down) position. For particular values of n and n , and a small di usion coecient , the typical out behavior of the gap as a function of both r and is shown in Fig. 7. The super- relaxation domain in the ( ; r) plane for di erent characteristics of interest in the out-of-comfort zones is also shown in Fig. 7: for large values of n , the out super-relaxation f ; rg-domain decreases signi cantly and tends asymptotically to a xed shape as shown on the right panel of Fig. 7. Convergence to the xed shape is rather fast due to the fact that the probability mass is localized around the comfort zone, and that the far-lying out-of-comfort zone solutions do not in uence dynamics of the model. These domains overlap: the deep orange domain is partly covered by the other domains, which are smaller in sizes. This result is consistent with the fact that the ensemble mixing, which favors fast relaxation, can hardly be achieved if the number of states in the out-of-comfort zone remains high. We have checked numerically that variation of the di usion coecient , has a rather limited impact on the super-relaxation surface. We also veri ed that in the limit of in nite number of states in the out- of-comfort and comfort zones, value of the converges to the one correspondent to the continuous model, thus implying that the dynamical behavior described in Ref. [31] is recovered in this limit. 3.3. Undamped oscillations in consumption As seen in Fig. (8), reporting experimental observations, at some values of r and , and depending on how time discretization is implemented, undamped oscillations in consumption are observed. The oscillations are also preceded by the period of growth. This is clearly an undesirable phenomenon which needs to be explained. In the following we are discussing results of comparison of the experiments with nonlinear system juxtaposed against the linear stability analysis. The compari- son shows that there exist a range of parameters where the linear analysis shows an instability fully consistent with the amplitude growth observed in the experi- ment. The oscillations are seen in the regime which is beyond the linear stability 18 analysis. Some details and discussions of the phenomenon are discussed in the following. r = 0.2, α = 25 (st) |N − N | 4 1 min{Reλ } {WF} − 1 − 6 − 11 0 50 100 150 200 250 Figure 8: System dynamics with unstable steady state, as observed for n = 30, n = 18, out r = 0:2 and = 25. In this case N and H grow in time according to e as is negative. " 1 Note that at some point, the exponential growth stops and the dynamics stabilizes; from this point on the linear analysis (red-dashed curve) no longer applies. To understand better dynamics of the energy consumption after a DR pertur- bation we ought to monitor for super-relaxation but also for instability, checking the dominant relaxation rate, , i.e. one with the smallest real part. Contrary to the continuous model where only one crossing (correspondent to equal real parts of and ) is observed as we change , in the discrete case multiple 1 2 events of level crossings are possible, e.g. as illustrated in Fig. 5. This means that the entire spectrum of the relaxation constants, f g, need to be considered to resolve which mode dominates the relaxation. Consider the case depicted in Fig. 9 and follow, as varies, the peculiar (1) behavior of the eigenvalue , which is related to the relaxation constant, (1) (1) , according to = log(j j) + i arg( ). Since at each time step the 1 1 (1) amplitude of the corresponding mode is multiplied by , the mode decays in 19 (1) ReΛ 1.0 (1) ImΛ 0.5 0.0 − 0.5 α α α 0 1 2 − 1.0 0 10 20 30 40 (1) Figure 9: Real and imaginary parts of the eigenvalue as functions of the degree of (1) nonlinearity . Here, becomes a real number from = 5:38. The particular values at (1) (1) which = 0 and becomes smaller than 1 are = 19:21 and = 39:45 respectively. 1 2 (1) (1) time if j j < 1, and it grows if j j 1. At suciently small and before reaches the value , , the mode decays. (Notice that there is also 1 1 (1) another special value, , where 0 < < , such that Im( ) is nite at 0 0 1 < and it is zero at 0. Crossing does not have implications on 0 0 how the mode decays.) The aforementioned instability occurs when becomes (1) larger than at which point = 1. It is useful to have a simple, albeit not absolutely precise, criterion which allows to avoid undesirable instability following by oscillations. We suggest a criterion based on estimations of and . As shown in Appendix Appendix 1 2 C, considering a simpli ed version of the dynamical equation, Eq. (C.6), yields 2(n 1) n 2 1 and 1. The estimation results in the following 1 2 r(n 2) r(n 2) out out estimation of the frontier separating stable and unstable regimes (see Appendix Appendix C for details) n 2 n =2 1 out 2 (1 + )r + 1 = 0 (19) n n 20 Summary of the behavior, illustrating frontier (red dashed curve) where the instability occurs, is also shown in Fig. 10 for an exemplary values of the di u- sion coecient and the size of the out-of-comfort zone. We observe that only eigenvalues from the main family can lead to instability since the ghost fam- ily is not a ected by the mean- eld feedback. Even though the (red dashed) boundary correspondent to the criterion (19) is not precise it nevertheless gives a conservative guidance on the range of parameters where the instability can be safely avoided. We conclude emphasizing that the instability is an unfortunate artifact of the discrete regime and it does not occur in the continuous regime discussed in [31]. n = 30, n = 18 out 1.0 0.8 0.6 0.4 0.2 0.0 0 10 20 30 40 Figure 10: Phase diagram showing instability observed for n = 30 and n = 18. The white out zone, where the real part of all relaxation constants is positive, is stable. The green zone is unstable with at least one relaxation constant with negative real part. The red dashed curve shows the analytical estimation of for di erent values of r, Eq. (19), approximately marking the frontier separating the two regimes. The dashed purple curve is the analytical estimation of the value of for di erent r that marks the frontier between dynamical regimes. r 3.4. Remarks on the linearity of x(t) To demonstrate super-relaxation in a model re ecting the discrete nature of an actual TCL ensemble, the relevant parameters to consider are the degree of nonlinearity , the Poisson rate r, and the number of load states within and outside the comfort zone, n and n respectively. Our space-time quantized in out model and hence our numerical results depend on the strong assumption of the linearity of x(t), the impact of which can be quanti ed with a concrete example. Note rst that since we work in the limit where x(t) is linear, the dynamics is invariant under a shift of all temperature parameters by some constant x, i.e. the transformation x ! x + x does not change the dynamics. There is "(#) "(#) thus no need for an explicit knowledge of the temperature threshold values, espe- cially as we aim for generality of our objective: observation of super-relaxation in discrete models pertinent to real-life cases. Now, let us consider the following illustrative case of a hot summer period when air conditioners are operating at high power. The outside temperature is x = 32 C, the minimal possible indoor temperature with a conditioner constantly being in the on position is x = 12 C, and the indoor comfortable temperature zone goes from x = 21 C to x = 23 C, with the mean being # " x = 22 C. The temperature dynamics, when the conditioner is in its o state, dx(t) satis es = (x x(t)), where is some variation rate. One can decompose dt x(t) as follows x(t) = x + x(t), where x(t) is a deviation of x(t) about a dx(t) mean comfortable temperature x . Then, we obtain = (x x x(t)). c + c dt The solution of this equation reads x(t) x(0) = (x x )t x( )d . + c The second term on the r.h.s. is the nonlinear term in t that is neglected when x x " # assuming linear evolution in Eq. (1). x(t) is at most of order (for large switching rates which is our case). Comparing the linear and nonlinear x( )d x x " # terms, we nd . = 1=10 for the realistic parameters (x x )t 2(x x ) + c + c validating our linear dynamics assumption. Similarly for the on state we get x( )d x x " # . = 1=10 1: Note that in this illustrative example (x x )t 2(x x ) c c both cooling and heating rates give 1/10, which shows at the same time that 22 the symmetric rate assumption is also valid. 3.5. In uence of the Poisson rate on the consumers' comfort From the consumers' viewpoint, remaining in the comfort zone is the most important aspect of the problem. Generally speaking, comfort can been seen as based on a subjective perception of indoor microclimate, but one may still use objective parameters and values such as indoor temperature [39, 40], humidity, CO concentration [41, 42], etc. to de ne bounds for comfort in the frame of a model. In our framework, the level of comfort can be de ned as C (t) = (t)C , where C is the weight of a particular state from the consumers' comfort point of view: C = 1 in the comfort zone, and 0 outside. In the regime when the stationary state of the ensemble is slightly perturbed because of DR, the dynamics of C (t) reads X t X X (i) (i) (i) (st) C (t) = C (0) + C ; (20) (st) (st) where C = C . 1.0 0.9 0.8 0.7 0.6 0.5 0.0 0.2 0.4 0.6 0.8 1.0 Figure 11: Dependence of the comfort level C on the Poisson rate r. C In our work we are interested in air conditioners as part of a large ensemble, and the key parameter to consider is temperature. The air conditioners are in our model two-state (on/o ) devices subjected to randomized switching char- acterized by the Poisson rate r. From the power systems viewpoint, the main goal is to recover a sucient level of ensemble mixing to avoid oscillations after DR perturbation. The smaller r is the faster the mixing is, but this obviously implies an increased Poisson time (1=r) and the resulting delayed switching on the level of consumers comfort as a deviation of temperature from thresholds as illustrated in Fig. 11. This particular point was addressed and discussed in [29]. − 2 − 4 − 6 (st) |N − N | − 8 ↑ min{Reλ } {EF} − 10 10 min{Reλ } {WF} (st) |C − C | 0 100 200 300 400 500 600 Figure 12: Temporal dynamics of di erent observable of the ensemble. = 10, r = 0:05. We may also comment on the impact of non-linearity, if it were to be in- troduced, on the consumers' comfort level. One can see that in the general case, the consumers' comfort relaxes to its stationary value obeying the same relaxation constants as other observables of the ensemble. Indeed, as shown in Fig. 12, the consumers comfort relaxes with the same rate as H , the L distance describing how the probability mass function (t) goes toward its steady state. 24 Hence, introduction of non-linearity would not have any impact on the level of consumers comfort in the steady-state, because it does not have any impact on the steady-state. 4. Conclusions and path forward We start the concluding Section of the manuscript with a brief summary of the results reported: • E ect of the super-relaxation, previously observed in the continuous time model, extends to more realistic discrete time models where it becomes a useful practical tool for demand response. • We show that the super-relaxation is stable with respect to variations, uctuations and uncertainty in the operationally sensible range of the model parameters. • We also observe that dynamics of the TCL ensemble is sensitive to some details of the discretization scheme. In particular, for values of the ensem- ble parameters correspondent to large accumulations of nonlinear e ects (including feedback) over a time step the system becomes linearly unstable then resulting in parasitic oscillations. We analyze the instability and pro- vide a simple to implement criteria which allows to avoid the undesirable regime. Discussing the last point in some extra details, it is important to emphasize that emergence of the parasitic instability is a special feature of the discrete time model not observed in the continuous time model. We observed that the super-relaxation in space-time-quantized models entails a more complicated spectral structure than that obtained with the continuous model [31]. We saw that undamped oscillations may arise if the discretization scheme is not cali- brated proper. To uncover this e ect we perform linear stability analysis and establish criteria for instability, then suggesting criteria on how to avoid it. 25 This instability analysis allows us to claim that the manuscript contributes to the growing body of work towards establishing regimes for safe operations of the TCL ensembles, i.e. seeking for operations which allow to mitigate various parasitic e ects. It is important to emphasize, however, that the oscillations re- ported in this manuscript are not related to (but rather imposed on the top of ) other oscillations already discussed in the literature and associated with irreg- ular patterns of consumption and synchronization following demand response signals [43, 44]. We conclude that aggregators and other participants of the energy markets should be aware of this newly reported discretization-caused in- stability as it may be dangerously enhanced, if not mitigated proper, in the case of increasing level of uctuations caused, for example, by increase of renewable penetration. Let us now turn to a brief discussion of the path forward. Even though the paper constitutes a signi cant step towards realistic operation of TCL en- sembles, more work is needed to adapt our results to practical setting of the demand response implementations. We envision relaxing various assumptions made in this study to simplify the analysis, such as accounting for asymmetry in heating and cooling, accounting for variations of parameters on the level of individual devices, etc. More detailed physical modeling at the device scale may and should include in the future modeling and monitoring of the air-quality, i.e. CO concentration, particulate matter concentration, aerosols, humidity, e.g. as discussed in [45]. Finally, we would like to emphasize that demand response is a general energy management tool, which is not restricted to power systems, and is in fact of an even greater utility for integrated energy systems [46, 47]. The mean- eld ap- proach may also be extended to other infrastructure systems such as battery-, water-, waste-, and oil-product systems dependent on exible consumers en- gaged in communications-light demand-response services. 26 Acknowledgments The work at LANL was carried out under the auspices of the National Nu- clear Security Administration of the U.S. Department of Energy under Contract No. DE-AC52-06NA25396, and it was partially supported by DOE/OE/GMLC and LANL/LDRD/CNLS projects. The work at Skoltech was supported by the Skoltech NGP Program (Skoltech-MIT joint project). Appendix A. Continuous limit of discrete master equation Here we discuss the master equation describing an ensemble of loads in continuous space-time limit. We parameterize the state of a device by the tuple, fx; g, where =";# marks the state of a load and x marks the indoor temperature. Let us denote the temperature di erence between neighbouring nodes, x, the number of nodes out of the comfort zone n , the number of out nodes in the comfort zone n , the total number of nodes n = n + n , and in out in the discrete time step t. We assume that n n , i.e. n ! 1. Then, out in out the discrete in space and time master equation (7) takes the following form: (x; t + t) = (1 2) (x + x; t) + (x + 2x; t) + (x; t); x x < x ; > " " " " # " (x; t + t) = (1 2) (x x; t) + (x 2x; t) + (x; t); x < x x ; > # # # # # " (x; t + t) = (1 2 q (t)) (x + x; t) + (x + 2x; t) + (x; t); x < x ; " # " " " # (A.1) > (x; t + t) = (1 2 q (t)) (x x; t) + (x 2x; t) + (x; t); x < x; # " # # # " > (x; t + t) = (1 2) (x + x; t) + q (t) (x; t) + (x + 2x; t) + (x; t); x x; > " " " # " " " (x; t + t) = (1 2) (x x; t) + q (t) (x; t) + (x 2x; t) + (x; t); x x : # # # " # # # One may study this system of equations with discrete derivatives denoted as (n) D , where n is the order of a derivative and x is the target variable. For exam- (1) F (x+x;y) F (x;y) (2) F (x x;y)+F (x+x;y) 2F (x;y) ple D (F (x; y)) = or D (F (x; y)) = . x x 2 x x In these new notations the system of Eqs. (A.1) becomes 27 8 (1) (1) (2) > x x D ( (x; t)) = D ( (x; t)) + D ( (x + x; t)); x x < x ; x x > " " " # " t t > (1) (1) (2) x x D ( (x; t)) = D ( (x x; t)) + D ( (x x; t)); x < x x ; > # x # x # # " t t > (1) (1) (2) q (t) x x # D ( (x; t)) = D ( (x; t)) + D ( (x + x; t)) (x + x; t); x < x ; x x " " " " # t t t (A.2) (1) (1) (2) q (t) x x " >D ( (x; t)) = D ( (x x; t)) + D ( (x x; t)) (x x; t); x < x; # x # x # # " > t t t (1) (1) (2) q (t) x x " >D ( (x; t)) = D ( (x; t)) + D ( (x + x; t)) + (x; t); x x; " x " x " # " > t t t t > 2 (1) (1) (2) q (t) : x x # D ( (x; t)) = D ( (x x; t)) + D ( (x x; t)) + (x; t); x x : # x # x # " # t t t Let us now consider the limit n ! 1; in x ! 0; = const 2 [0; 0:5] t ! 0; r ! 0 n x = const = L; in = const = v; = const = r (the subscript c refers to the continuous case.)(A.3) Notice that di usion related term, O( t), vanishes in the limit. Also, there is no longer a need, in this limit, for the function f , whose role in the quantized model was to ensure non-negativity of the transition matrix elements. This discrete in space and time master equation turns into the following system of 28 continuous Fokker-Planck equations, Eq. (3): >@ (x;t) @ (x;t) " " = v ; x x < x ; > # " @t @x >@ (x;t) @ (x;t) # # = v ; x < x x ; # " > @t @x @ (x;t) @ (x;t) " " = v 2r (N (t)) (x; t); x < x ; c " " # @t @x (A.4) @ (x;t) @ (x;t) # # > = v 2r (1 N (t)) (x; t); x < x; c " # " @t @x @ (x;t) @ (x;t) > " " = v + 2r (1 N (t)) (x; t); x x; c " # " @t @x @ (x;t) @ (x;t) : # # = v + 2r (N (t)) (x; t); x x : c " " # @t @x Appendix B. Consumption in the Steady State (st) Assume that the steady-state of the master equation Eq. (7) is . Then, it satis es (st) (st) (st) p 0 = : (st) (st) We aim to show that N = U = . In order to prove it, let us " 2 0 0 consider the linear transformation T which acts on the state = T and makes the following changes: swaps on and o and reverse order of X . T is also a stochastic matrix Fig. B.13. Other important properties of the matrix T are: T = 1 and TPT = P , where P is the transition matrix in the case without the mean eld control, 1 is the identity matrix. TPT = P directly follows from the symmetry of the P matrix with respect to such a transformation T . This two properties lead to the following commutation relation: TP = PT (B.1) Using this commutation relation one derives P = TP = T PT = T (B.2) 29 In the case when we have only one steady state, T = . The relation is satis ed only if U = 1=2. Therefore, this property is a consequence of the transition matrix symmetry. (Notice that we do not consider here a more complicated case of multiple competing steady states.) x x x x − ↓ ↑ + 0 0 0 0 0 0 0 0 0 X X X X X X X X X 1 2 3 4 5 6 7 8 9 1 1 1 1 1 1 1 1 1 X X X X X X X X X 1 2 3 4 5 6 7 8 9 Figure B.13: Graphical representation of T 0 . Transition probabilities are all equal to 1. Appendix C. Oscillating dynamics Appendix C.1. Variational principle In order to simplify the original master equation one can use a variational principle. to search for approximate solution within a de ned class of functions. Let us start with the master equation p 0 ((t)) 0 (t) = (t + 1): (C.1) and consider the following functional X X L R L R R L R L( ; ) = (t)p 0 ( (t)) (t) (t) (t + 1); (C.2) ; ;t ;t R L where is a probability mass function and is an auxiliary vector (conjugated distribution). Observe that the stationary point of this functional results in the master equation L R X @L( ; ) R R R 0 = = p ( (t)) 0 (t) (t + 1): (C.3) @ (t) k 0 Using this variational principle one can explore di erent class of functions and try to nd the best solution from this family by minimizing of the functional Eq. (C.2). 30 Appendix C.2. Theoretical explanation of instability Let us use the variational formulation to gain a qualitative explanation of the discretization-related instability discussed in the main part of the paper. We derive 2 3 2 3 2 3 N (t) L R " v (t) v (t) " " 6 7 6 7 6 7 . . . 6 7 6 7 6 7 . . . 6 7 6 7 6 7 . . . 6 7 6 7 6 7 6 7 6 7 6 7 N (t) L R " 6 7 6 7 6 7 v (t) v (t) " " n L R 6 7 6 7 6 7 (t) = ; (t) = = (C.4) 6 7 6 7 6 7 (1 N (t)) L R " 6 7 6 7 6 7 v (t) v (t) # # n 6 7 6 7 6 7 6 7 6 7 6 7 . . . . . . 6 7 6 7 6 7 . . . 4 5 4 5 4 5 (1 N (t)) L R v (t) v (t) # # R L R L where v (t), v (t), v (t), v (t) are new variables. This type of variational # # " " ansatz enforces uniform distribution for ON and OFF states along coordinate R R (temperature) form x to x . Changing variables, nv (t) = N (t), nv (t) = + # # " N (t), where n is total number of states, results in the following system of equations h i h i n =2 1 n =2 1 < n 1 out 1 out N (t + t) = f (r (2N (t)) ) N (t) + + f (r (2N (t)) ) N (t); " " " # # n n n n (C.5) h i h i 1 n =2 1 n 1 n =2 1 out out N (t + t) = + f (r (2N (t)) ) N (t) + f (r (2N (t)) ) N (t); # " " # # n n n n where n is the number of states which are outside of the comfort zone. Using out normalization condition, N +N = 1, we reduce the system to a single equation # " n 2 n =2 1 out N (t + t) = (f (r(2N (t)) ) + f (r(2(1 N (t))) )) N (t) " " " " n n 1 n =2 1 out + + f (r(2(1 N (t))) ) : (C.6) n n Consider a small perturbation around the stationary state: N (t) = + N (t). Linearized version of Eq. (C.6) becomes n 2 n =2 1 out N (t + t) = N (t) 2 (1 + )r (C.7) n n Analysis of this relation shows emergence of the 3 distinct regimes in the space- time-quantized model which are interpreted as follows 31 n 2 n =2 1 out • 0 2 (1 + )r 1, mean eld control speeds up relaxation, n n n =2 1 n 2 out • 1 2 (1 + )r 0, mean eld control is too strong and n n it leads to the perturbation alternating its sign at every step, while the absolute value of perturbation is still decaying, n =2 1 n 2 out • 2 (1 + )r 1, mean eld control changes sign and in- n n creases absolute value of perturbation, resulting in the instability. Let us now make a brief comment on the absence of discretization instability in the space-time-continuous model. Consider Eq. (C.7) in the continuous case: out N (t + dt) = N (t) [1 (1 + )r] ; = : (C.8) In this case r ! 0, dt ! 0 and = constant, so that the factor next to N (t) dt in the equation never becomes negative. As a result the continuous dynamics is never unstable. The dynamics never becomes unstable. Notice that the instability occurs in simulations when one uses a too large time step. This e ect is associated with the fact that in the discrete time we may have 1 (1 + ) either negative or positive, and as N is not continuous, the di erence (N (t + dt) N (t)) (slope) is always pointing towards the axis because (N (t +dt) N (t))=N (t) < 0. Note that this observation also applies to the continuous time models, however since N (t) varies continuously it only results in decay (towards zero). The comparison between the numerical analysis and the simple theoretical estimation discussed above is shown in Fig. 10. Appendix D. Some properties of the spectrum of S (i) Let us assume that the spectrum f g and the set of right eigenvectors (i) f g are known. Then we write (i) (i) (i) S = : (D.1) Since, S 0 = 1, one also derives (i) (i) (i) = ; (D.2) 32 P (i) (i) where = . The equality (D.2) implies that there are two types (i) of eigenmodes. The rst type is associated with being arbitrary and P P (i) (i) (i) = 0. The second type occurs when = 1 and as a result, , is not constrained. There is at least one mode of the second type with the fol- lowing property of the corresponding right eigenvector: 6= 0, since it is otherwise impossible to decompose a vector (for which 6= 0) as a linear (i) combination of right eigenvectors f g. Therefore, one reaches the following conclusions about the spectrum of S: 1. there exists at least one mode with = 1 and 6= 0; 2. for those modes that have 6= 1, the corresponding right eigenvectors satisfy: = 0. References [1] N. O'Connell, P. Pinson, H. Madsen, and M. O'Malley, Bene ts and chal- lenges of electrical demand response: A critical review, Renewable and Sus- tainable Energy Reviews 39, 686-699 (2014). doi:10.1016/j.rser.2014.07.098 [2] T. J. Hammons, Integrating renewable energy sources into European grids, International Journal of Electrical Power & Energy Systems 30, 462{475 (2008). doi:10.1016/j.ijepes.2008.04.010 [3] D. K. Critz, S. Busche, and S. Connors, Power systems balancing with high penetration renewables: The potential of demand response in Hawaii, Energy Conversion and Management 76, 609{619 (2013). doi:10.1016/j.enconman.2013.07.056 [4] H. Auer and R. Haas, On integrating large shares of variable re- newables into the electricity system, Energy 115, 1592{1601 (2016). doi:10.1016/j.energy.2016.05.067 [5] M. H. Albadi and E. F. El-Saadany, Demand response in electricity mar- kets: An overview, Proceedings of the 2007 IEEE Power Engineering Soci- ety General Meeting, 24-28 June 2007. doi:10.1109/PES.2007.385728 33 [6] I. Lampropoulos, W. L. Kling, P. F. Ribeiro, and J. van den Berg, History of demand side management and classi cation of demand response control schemes, Proceedings of the 2013 IEEE Power Energy Society General Meeting, pages 1-5. doi:10.1109/PESMG.2013.6672715 [7] M. A. Lynch, S. Nolan, M. T. Devine, M. O'Malley, The impacts of demand response participation in capacity markets, Applied Energy 250, 444{451 (2019). doi:10.1016/j.apenergy.2019.05.063 [8] K. Wohlfarth, M. Klobasa, R. Gutknecht, Demand response in the service sector { Theoretical, technical and practical potentials, Applied Energy 258, 114089 (2020). doi:10.1016/j.apenergy.2019.114089 [9] C. Schiel, P. G. Lind, and P. Maass, Resilience of electricity grids against transmission line overloads under wind power injection at di erent nodes, Scienti c Reports 7, 11562 (2017). doi:10.1038/s41598-017-11465-w [10] A. Gajduk, M. Todorovski, J. Kurths and L. Kocarev, Improving power grid transient stability by plug-in electric vehicles, New Journal of Physics 16, 115011 (2014). [11] M. Yesilbudak and A. Colak, Integration challenges and solutions for renewable energy sources, electric vehicles and demand-side initiatives in smart grids, 7th International IEEE Conference on Renewable En- ergy Research and Applications, ICRERA 2018 8567004, pp. 1407-1412. doi:10.1109/ICRERA.2018.8567004 [12] M. Motalleb, M. Thornton, E. Reihan, and R. Ghorbani, Pro- viding frequency regulation reserve services using demand response scheduling, Energy Conversion and Management 124, 439{452 (2016). doi:10.1016/j.enconman.2016.07.049 [13] C.-Y. Chong and A. S. Debs, Statistical synthesis of power system func- tional load models, Proceedings of the 1979 18th IEEE Conference on Deci- 34 sion and Control including the Symposium on Adaptive Processes 2, 264{ 269 (1979). doi:10.1109/CDC.1979.270177 [14] S. Ihara and F. Schweppe, Physically based modeling of cold load pickup, IEEE Transactions on Power Apparatus and Systems 100, 4142{4150 (1981). doi:10.1109/TPAS.1981.316965 [15] C.-Y. Chong and R. P. Malham e, Statistical synthesis of physically based load models with applications to cold load pickup, IEEE Trans- actions on Power Apparatus and Systems 103, 1621{1628 (1984). doi: 10.1109/TPAS.1984.318643 [16] C. W. Gardiner, Handbook of stochastic methods for physics, chemistry and the natural sciences, 3rd ed. (Springer Series in Synergetics, vol. 13, Berlin: Springer-Verlag, 2004). [17] N. van Kampen, Stochastic Processes in Physics and Chemistry (Third Edition) (Amsterdam: Elsevier, 2007). [18] R. Malham e, R. and C.-Y. Chong, Electric load model synthesis by di usion approximation of a high-order hybrid-state stochastic sys- tem, IEEE Transactions on Automatic Control 30, 854{860 (1985). doi:10.1109/TAC.1985.1104071 [19] R. Malham e and C.-Y. Chong, On the statistical properties of a cyclic di usion process arising in the modeling of thermostat-controlled electric power system loads, SIAM Journal on Applied Mathematics 48, 465{480 (1988). doi:10.1137/0148026 [20] S. El-Ferik and R. P. Malham e, Identi cation of alternating renewal electric load models from energy measurements, IEEE Transactions on Automatic Control 39, 1184{1196 (1994). doi:10.1109/9.293178 [21] N. Lu and D. Chassin, A state-queueing model of thermostatically controlled appliances, IEEE Transactions on Power Systems 19, 1666{1673 (2004). doi:10.1109/TPWRS.2004.831700 35 [22] N. Lu, D. Chassin, and S. Widergren, Modeling uncertainties in aggregated thermostatically controlled loads using a state queueing model, IEEE Transactions on Power Systems 20, 725{733 (2005). doi:10.1109/TPWRS.2005.846072 [23] D. S. Callaway, Tapping the energy storage potential in electric loads to deliver load following and regulation, with application to wind energy Energy Conversion and Management 50, 1389-1400 (2009). doi:10.1016/j.enconman.2008.12.012 [24] D. S. Callaway and I. A. Hiskens, Achieving controllability of electric loads Proceedings of the IEEE 99, 184{199 (2011). doi:10.1109/JPROC.2010.2081652 [25] S. Bashash and H. K. Fathy, Modeling and control insights into demand- side energy management through setpoint control of thermostatic loads, Proceedings of the 2011 American Control Conference, 4546{4553 (2011). doi:10.1109/ACC.2011.5990939 [26] C. Perfumo, E. Kofman, J. H. Braslavsky, and J. K.Ward, Load manage- ment: Model-based control of aggregate power for populations of thermo- statically controlled loads, Energy Conversion and Management 55, 36{48 (2012). doi:10.1016/j.enconman.2011.10.019 [27] D. P. Chassin, J. Stoustrup, P. Agathoklis, and N. Djilali, A new thermostat for real-time price demand response: Cost, comfort and energy impacts of discrete-time control without deadband, Applied Energy 155, 816{825 (2015). doi:10.1016/j.apenergy.2015.06.048 [28] C. Wei, J. Xu, S. Liao, Y. Sun, Y. Jiang, and Z. Zhang, Coordination optimization of multiple thermostatically controlled load groups in distribu- tion network with renewable energy, Applied Energy 231, 456{467 (2018). doi:10.1016/j.apenergy.2018.09.105 36 [29] M. Chertkov and V. Chernyak, Ensemble of thermostatically controlled loads: Statistical physics approach, Scienti c Reports 7, 8673 (2017). doi:10.1038/s41598-017-07462-8 [30] D. M etivier, I. Luchnikov, and M. Chertkov, Power of ensemble diver- sity and randomization for energy aggregation, Scienti c Reports 9, 5910 (2019). doi:10.1038/s41598-019-41515-4 [31] D. M etivier and M. Chertkov, Mean- eld control for ecient mix- ing of energy loads, Physical Review E 101, 022115 (2020). doi:10.1103/PhysRevE.101.022115 [32] A. Rajabi, L. Li, J. Zhang, and J. Zhu, Aggregation of small loads for de- mand response programs { Implementation and challenges: A review, 2017 IEEE International Conference on Environment and Electrical Engineering and 2017 IEEE Industrial and Commercial Power Systems Europe (EEEIC / I&CPS Europe). doi:10.1109/EEEIC.2017.7977631 [33] M. Huang, P.E. Caines, R.P. Malhame, Individual and mass behaviour in large population stochastic wireless power control problems: Centralized and Nash equilibrium solutions, Proceedings of the 42nd IEEE Conference on Decision and Control (CDC), (2003). doi:10.1109/CDC.2003.1272542 [34] G.Y. Weintraub, C.L. Benkard, B. Van Roy, Oblivious Equilibrium: A Mean Field, Advances in Neural Information Processing Systems, MIT Press, (2005). [35] J.-M. Lasry, P.-L. Lions, Jeux a champ Moyens I,II, Comptes Rendus Math ematiques 343, 619{625 (2006). doi:10.1016/j.crma.2006.09.019 [36] J.-M. Lasry, P.-L. Lions, Mean Field Games, Japanese Journal of Mathe- matics, 2, 229{260 (2007). doi:10.1007/s11537-007-0657-8 [37] M. Huang, R.P. Malham e, P.E. Caines, Large Population stochastic dy- namic games: Closed loop McKean-Vlasov and the Nash certainty equiv- 37 alence principle, Communications in Information & Systems 6,221{252 (2006). https://projecteuclid.org/euclid.cis/1183728987 [38] A. Bensoussan, J. Frehse, and P. Yam, Mean Field Games and Mean Field Type Control Theory (Springer-Verlag New York, 2013). [39] ANSI/ASHRAE Standard 55-2013: Thermal environmental conditions forhuman occupancy. [40] P. O. Fanger, Thermal comfort. Analysis and applications in environmental engineering (Danish Technical Press, Copenhagen, 1970). [41] U. Satish, M. J. Mendell, K. Shekhar, T. Hotchi, D. Sullivan, S. Streufert, et al. Is CO an indoor pollutant? Direct e ects of low-to-moderate CO con- 2 2 centrations on human decision-making performance, Environmental Health Perspectives 120, 1671{1677 (2012). [42] J. G. Allen, P. Mac Naughton, U. Satish, S. Santanam, J. Vallarino, J. D. Spengler, Associations of cognitive function scores with carbon diox- ide,ventilation, and volatile organic compound exposures in oce workers: a controlled exposure study of green and conventional oce environments, Environmental Health Perspectives 124, 805{812 (2016). [43] N. A. Sinitsyn, S. Kundu, and S. Backhaus, Safe protocols for generat- ing power pulses with heterogeneous populations of thermostatically con- trolled loads, Energy Conversion and Management 67, 297{308 (2012). doi:10.1016/j.enconman.2012.11.021 [44] N. Mehta, N. A. Sinitsyn, S. Backhaus, B. C. Lesieutre, Safe control of thermostatically controlled loads with installed timers for demand side management, Energy Conversion and Management 86, 784{791 (2014). doi:10.1016/j.enconman.2014.06.049 [45] A. Ryzhov, H. Ouerdane, E. Gryazina, A. Bischi, and K. Turitsyn, Model predictive control of indoor microclimate: existing building stock comfort 38 improvement, Energy Conversion and Management 179, 219{228 (2019). doi:10.1016/j.enconman.2018.10.046 [46] Y. Liu, M. S. Mauter, Assessing the demand response capacity of U.S. drinking water treatment plants, Applied Energy 267, 114899 (2020). doi:10.1016/j.apenergy.2020.114899 [47] A. Hassan, S. Acharya, M. Chertkov, D. Deka and Y. Dvorkin, A hierarchi- cal approach to multi-energy demand response: From electricity to multi- energy applications, Special Issue on Multi-Energy Systems, Proceedings of IEEE 108, 1457{1474 (2020). doi:10.1109/JPROC.2020.2983388
Journal
Computing Research Repository – arXiv (Cornell University)
Published: Aug 3, 2020