Power of Ensemble Diversity and Randomization for Energy Aggregation

We study an ensemble of diverse (inhomogeneous) thermostatically controlled loads aggregated to provide the demand response (DR) services in a district-level energy system. Each load in the ensemble is assumed to be equipped with a random number generator switching heating/cooling on or off with a Poisson rate, r, when the load leaves the comfort zone. Ensemble diversity is modeled through inhomogeneity/disorder in the deterministic dynamics of loads. Approached from the standpoint of statistical physics, the ensemble represents a non-equilibrium system driven away from its natural steady state by the DR. The ability of the ensemble to recover by mixing faster to the steady state after its DR’s use is advantageous. The trade-off between the level of the aggregator’s control, commanding the devices to lower the rate r, and the phase-space-oscillatory deterministic dynamics is analyzed. Then, we study the effect of the load diversity, investigating four different disorder probability distributions (DPDs) ranging from the case of the Gaussian DPD to the case of the uniform with finite support DPD. We show that stronger regularity of the DPD results in faster mixing, which is similar to the Landau damping in plasma physics. Our theoretical analysis is supported by extensive numerical validation.


Results
We study the effects of the load inhomogeneity (which we also call disorder, following the statistical physics jargon) on operations of the ensemble, specifically in terms of the ensemble's ability to recover fast from a perturbation after its use by the aggregator for the DR. The ensemble is assumed to be controlled by an aggregator in a communication-minimal way by sending the same signal switching off/on rate to all the consumers simultaneously. We are mainly interested in the regime where both the control and the ensemble variability are weak. The two main messages of the paper (put here in a colloquial format and then quantified formally later) are as follows: (1) There exists an optimal switching rate corresponding to the fastest recovery. Any deviation (increase or decrease) of the rate leads to a slower recovery. (2) Increase of the ensemble variability is advantageous for faster recovery/mixing. Because temporal evolution of the ensemble is at the core of this manuscript analysis, let us define relevant timescales and then restate our main results in a more technical way. We assume that by default (without aggregator), each customer follows a standard bang-bang operation-switching on (off) the cooling device when the temperature exceeds (becomes less than) a preset threshold. We assume that the outside temperature is significantly higher than the switch-on threshold, thus resulting in cycling of the device with its natural timescale τ. The aggregator changes this natural cycling by requesting the consumers to switch on/off with a random delay distributed according to a Poisson distribution with rate r. (Each device is assumed equipped with a random number generator). By default, i.e., without aggregator control, = +∞ r . Weak aggregator control means that τ  r 1. To account for variability within the ensemble, one assumes that devices may have slightly different τ. Formally, one considers τ as the disorder (variability) parameter distributed according the disorder probability distribution (DPD), τ g( ), characterized in terms of its typical value (mean), τ 0 , and the distribution width, ∆. In the following we define and consider four different forms of τ g( )-Gaussian, Lorentzian, Laplace, and finite-support uniformparameterized by τ 0 and ∆, however always assuming (analyzing the disorder case) that the typical control is weak, i.e., τ  r 1 0 , and that the disorder is also weak, i.e., τ ∆  0 . With the timescales and two small dimensionless parameters, τ − r ( ) 0 1 and ∆/τ 0 , defined, we are ready to provide the following more technical, still qualitative but intuitive, explanations for our main results.
(1) When = ∞ r , the system does not decay and temporal evolution of the probability distribution function (PDF) of a device temperature, x, averaged over the ensemble shows a periodic behavior in time, ~λ ±i t exp( ) I , with the period λ τ π = 1/ /(2 ) I . Decrease of r leads to decrease of λ I and simultaneous increase (from zero at = ∞ r ) of the decay rate, λ R . In this oscillatory with a decay regime, temporal behavior of the correction to the stationary probability distribution becomes where ± reflects emergence of two complex-conjugated solutions. At a certain critical value, = r r c , λ I becomes zero, i.e., the two complex-conjugated solutions merge into one (degenerate) solution such that close to the merging point λ λ and λ c is the critical value of λ R achieved at = r r c . The main conclusion of this straightforward qualitative estimate is that the lowest of the two eigenvalues (corresponding to ±1 → −1 and thus to the slowest asymptotic decay) achieves its maximum as a function of r at r c .
(2) In the default regime (no aggregator control), a set of devices with exactly the same τ, i.e., when we set ∆ to zero, would not mix at all, i.e., correction to the stationary probability distribution oscillates and does not decay. Introduction of a small but finite ∆ results in a decay that is largely controlled by τ g( ) in the vicinity of its maximum, i.e., at τ τ ≈ 0 . Specifically, decay of the temperature probability to its stationary value in www.nature.com/scientificreports www.nature.com/scientificreports/ time is controlled by the shifted Fourier transform of the DPD, ∫ τ τ Obviously, details of the decay depend on the shape of the DPD. Of the four model DPDs considered in this manuscript, the Gaussian DPD results in the fastest decay (shifted Gaussian in time), the Lorentzian DPD is a bit slower (exponential in time), and the Laplacian DPD and uniform finite-support DPD are the slowest, with asymptotic 1/t 2 and 1/t decays, respectively. This hierarchy, illustrated in Fig. 1, and the Fourier-transform interpretation suggest that the speed of decay is linked to the regularity of the DPD around its central part. (This phenomenon is reminiscent of the mathematically similar analysis of the Landau damping in plasma physics described by the Vlasov equation [see 12,13 and references therein]. Specifically, we refer here to the fact that the regularity of the initial velocity distribution influences the Landau mixing/damping speed).
These main focal points are detailed and extended in the remainder of the manuscript. Models of the statistical ensemble and of the ensemble inhomogeneity are formulated in "Formulation" section. The basic model of the homogeneous ensemble is analyzed in "Basic Homogeneous Model" section. Effects of the disorder/inhomogeneity are studied in "Disorder" section, where we also compare analytic and numerical results. "Discussions" section is reserved for conclusions and discussion of the path forward.
Formulation of the problem. One characterizes a load by the continuous parameter, x, standing for the temperature, and by the discrete/binary parameter, σ = ↑ ↓ , , indicating whether the air conditioning system/ device of the load (one considers cooling for concreteness) is switched on, σ = ↑, or off, σ = ↓. Conditioned to σ, the dynamics of x follow the deterministic rule where σ | v x ( ) describes the rate of temperature change as a function of the current temperature, x, conditioned to the state of the load's air conditioning device (later in the text referred to simply as "device"). Our basic model is where u is a positive constant. The model is a simplification of a bit richer popular model, e.g., used in 10 , where u is not a constant as in Eq. (2) but a linear function of x. σ in Eqs (1) and (2) is modeled as the following Markovian binary (two-level) stochastic process: www.nature.com/scientificreports www.nature.com/scientificreports/ where dt is the time step (of the properly discretized continuous time limit), r is the rate of exponential (Poisson) switching, and ↑ x , ↓ x marks the size of the temperature band within which no switching occurs, < ↓ ↑ x x . As set above, the basic model has two timescales: one describing deterministic evolution, τ = − ↑ ↓ x x u 2( )/ , which is the time it takes for a device to make a full cycle through the combined σ x ( , ) phase space illustrated in Fig. 2, and 1/r, the typical time of a stochastic jump from σ = ↑ to σ = ↓ or vice versa.
Notice that the selection of the model in Eq. (2) as the basis for this manuscript analysis is dictated not only by its realism but also by considerations of simplicity and our ability to derive analytic results. Specifically, for the case of the asymptotic uniform ensemble consisting of the infinite number of devices with the same characteristics (the same u) and following the same switching protocol (the same r), we are interested in computing and analyzing the evolution in time of the PDFs governed by the system of coupled FP equations following directly from the model definition, given by Eqs (1), (2) and (3): x where θ(y) is unity if > y 0 and zero otherwise. We are seeking a properly normalized solution of Eq. (4): , , , counts proportions of devices that are switched on and off, respectively. As shown in the section "Basic Homogeneous Model", solution of the system of the FP Eq. (6) can be presented explicitly as the spectral expansion in terms of the Lambert-W functions for any initial = t 0 distributions. This analytic expression will allow us to analyze the temporal evolution of the basic homogeneous ensemble in much more detail than 10 for a more complex model, with (1) dependent linearly on x. However, devices contributing realistic ensembles are not necessarily the same in terms of their cooling/heating strength. To model the ensemble diversity, i.e., non-uniform ensemble, one introduces disorder in τ. We assume that τ characterizing a device is drawn independently from one of the following four model DPDs: Gaussian, Lorentzian, Laplace, and uniform (finite support) www.nature.com/scientificreports www.nature.com/scientificreports/ representing different extremes (e.g., in terms of the asymptotics). We parameterize these DPDs via their mean/ max, τ 0 , and variance, ∆, in a similar way to facilitate comparisons. In general, we will assume that τ ∆ ≤ 0 , and in terms of the asymptotic analysis, we will be interested most in the regime of weak disorder, τ ∆  0 . (Notice that the negative values of τ, τ < 0, are not physical. Therefore, when performing asymptotic analysis for the disorder distributions with formally defined infinite support, described by Eqs (7), (8) and (9), one needs to make sure that the fictitious τ < 0 regime does not contribute the asymptotic results). Then, the following average d over the DPDs , 0 , , 0 , will be the focus of our analysis of the inhomogeneous ensembles represented by Eqs (7), (8), (9) and (10). Equation (12) defines the total number of devices observed at the moment of time t in the state ↑ or ↓ (for all τ), Eq. (11) defines similar, however more detailed, object which also differentiate, in addition, between the current temperature, x, of the system. In this manuscript, we pose and answer the following two related questions: • Qualitative Question about the Inhomogeneous Ensemble: Does the disorder accelerate or slow down mixing, i.e., relaxation of the ensemble probability distribution to its steady state? • Quantitative Question about the Inhomogeneous Ensemble: How does the relaxation look depending on the system parameters and the parameters characterizing the PDF of the disorder?
Our choice of the basic model in Eq. (2), resulting in analytic expression for τ | P x t r ( , , ) stated in terms of the explicit spectral series in the section "Basic Homogeneous Model", allows us in the section "Disorder" to answer the quantitative question explicitly and then to use the analytic solution to reach qualitative conclusions.

Analytic solution for the Basic Homogeneous
; ; ; ; x τ ; where †  , the adjoint of , and the standard L 2 scalar product between two vectors P and G are defined according to It is straightforward to check that the eigenvectors, defined by Eqs (16) and (18) . Now closing the loop in Eq. (13) and linking the a coefficients there to the initial condition, ). We discuss consequences of the analysis on special features of the spectral problem, long time analysis (of the gap), and sensitivity of the asymptotic solution to the parameters in the following three subsections.
Features of the Spectral Problem. We find it useful to identify a number of significant features of the spectral problem defined by Eqs (13-21): (1) λ = + 0 0; and all other eigenvalues have a positive real part that grows with |k|, i.e., the spectrum is discrete, positive, and ordered. 1; . All other eigenvalues (with a nonzero imaginary part) have a real part that is larger than λ −1;− . Therefore, in the "low switching rate" regime, the solution decays with time (no oscillations). (4) When one of the two parameters, τ or r, is fixed, one finds that the largest value of λ 0;− is achieved at the bifurcation point, where β τ =  r /4 reaches β = C/4 c , i.e., given fixed r and changing τ, or fixed τ and changing r, mixing is the fastest at τ = C/r. (5) Moreover, given r is fixed, dRe(λ 0;−1 )/dτ, i.e., the rate of change with τ of the real part of the leading eigenvalue is positive/negative when rτ is smaller/larger than C. (6) When r is fixed and β τ = r /4 is sent to zero, one finds that λ → − r 0; . Indeed, in this regime, all devices that are in the allowed range, ∈ ↓ ↑ x x x [ , ], move fast to their respective boundaries; thus, at small τ, PDF decay is controlled primarily by the Poisson jumps/switchings. (7) When r is fixed and β τ = r /4 is sent to ∞ (or alternatively when τ is fixed and β is sent to ∞), one arrives at the following asymptotic: . One concludes that in the asymptotic regime of the "highest switching rate", the temporal evolution of the PDF becomes oscillatory and relaxation to the steady state slows down asymptotically to zero.
Evolution of the (four) leading non-zero eigenvalues (containing the smallest real part) with the dimensionless parameter β τ = r /4 and fixed τ is illustrated in Fig. 3. It is worth noting that the behavior of our system described , respectively, evaluated at τ τ = 0 . The coefficients of interest show the following asymptotics at small e τ =  r 1/( ) 0 : plane. Markers indicate β = .
.  Basic Model with Disorder. Averaging over the disorder according to Eq. (11) with only the leading = k 0, "−" term in Eq. (13) is justified when the spectral gap condition (see the section "Long Time Asymptotic Analysis: Gap Condition") is verified. Furthermore, we assume τ ∆  / 1 0 so that the integral Eq. (11) is concentrated for τ located around τ 0 so that τ τ τ | − |  0 0 . Then, taking into account the large time asymptotic (Eq. (22)) and assuming that the Taylor series expansion (Eqs (26) and (27) is legitimate (when τ > r C 0 ), one arrives at Versions of Eq. (32) for the four example probability distributions of the disorder (Eqs (7-10)) computed for Particle Simulations and Comparison with the Theory. To test our analytic results, we performed particle simulations of the dynamics of Eqs (1), (2) and (3). One first associates with each of N devices its own relaxation time, τ, drawn i.i.d. from one of the DPDs defined by Eqs (7), (8), (9) and (10). (Negative values of τ are rejected). Initially, at = t 0, all devices are set to = ↓ x x and σ = +1, corresponding to the "worst case", i.e., least mixed, initial distribution. Then the dynamics, advanced discretely and independently for each device, are implemented according to the following rules. At the beginning of each time interval, t, the state of each of the N devices, characterized by σ and x, is advanced in time according to the first-order (in time) version of Eqs (1) and (3). At each t, we monitor ↑ N t ( ), which is the total number of the devices in the state +1, also corresponding to the instanta-www.nature.com/scientificreports www.nature.com/scientificreports/ neous energy consumption of the ensemble (under the model assumption that each device, when switched on, consumes the same amount of energy).
The results of the straightforward particle simulations are illustrated in Fig. 4 for four different DPDs. To facilitate comparison, we juxtapose the results of the simulations with the corresponding analytic predictions given by Eqs (33), (35), (37) and (39). We observe very good agreement between the theory and the simulations at short and intermediate times. The conclusion is reached based on comparison of the amplitude and the frequency of the oscillations and the relaxation rate of ↑ N t ( ). Note that the theory results are derived in the asymptotic, weak disorder regime described by Eq. (32) and its complex-conjugated expression corresponding to the same (worst case) initial condition as in the simulations; hence, there are no fitting parameters. We also observe that at sufficiently large t, controlled by the finite (not infinite) size of the ensemble, the theory and the simulations start to deviate. Indeed, when | − | ↑ N t ( ) 1/2 becomes of the order of 1/ N , fluctuations associated with the finiteness of the ensemble start to dominate results of the simulations. In the simulations with = N 10 5 , this threshold is reached 3 . Comparing the four subfigures in Fig. 4 with each other is useful because it illustrates dependence of the ensemble mixing on different types of disorder.

Discussions
The main conclusion of the manuscript is that both types of randomizations, smoothing out the bang-bang control via Poisson-delayed switching and introducing diversity of loads in the ensemble, result in acceleration of the mixing/recovery following a heavy DR use of the ensemble. Specifically, we have shown via rigorous analysis and numerical simulations that (a) increasing the level of control (decreasing the switching rate) is advantageous only at sufficiently large rates, > r r c ; and (b) diversity of the devices' natural timescale (speed of cooling/heating), which is more "regular" (e.g., distributed according to the Gaussian DPD), is advantageous in leading to a faster mixing (more efficient recovery). In this paper we study quantitatively effect of the load diversity on the speed of the ensemble restoration to normal (steady distribution) following a significant DR perturbation. To the best of our knowledge, this is the first systematic analysis showing and explaining that disorder does provide a significant acceleration. It is also worth mentioning that similar analysis and effects (acceleration of relaxation to the steady state due to disorder) is expected in other applications where ensemble of particles follow the same stochastic signal. This expectation is consistent with numerical observations made in the context of signal tracking in 16 , where it was reported that diversity of loads results in a faster tracking. Moreover, many models of population www.nature.com/scientificreports www.nature.com/scientificreports/ dynamics 17,18 and more generally of non-equilibrium statistical mechanics 19 , has mathematical structure resembling one analyzed in the manuscript thus suggesting that there should be a wide range of possible applications where disorder in the model parameters results in faster mixing.
Encouraged by the reported results, we plan to extend the study in the following directions: • Complex Modeling. We envision considering more complex models of both the individual device dynamics and the ensemble compilation. For the former, different switching rates (for switching on and off) and more general dependence of the relaxation speed u on x are two practical complications that can be included in the analysis. For the latter (richer disorder), most significant generalization corresponds to adding disorder/ inhomogeneity in other model parameters, such as switching on/off temperatures. Our working hypothesis is that these modifications/generalizations will lead to (possibly significant) quantitative but not qualitative changes in the predictions. • Mean-Field, Nonlinear Control. Switching rate, r, communicated by the aggregator to consumers, was constant in the model discussed above. It is interesting to experiment with changing the rate, in particular allowing it to depend on the current state of the ensemble, i.e., on the instantaneous probability distribution in the (x, σ) space. This Mean-Field control improves greatly the relaxation time, as shown by the team in 20 . This intricate scenario is related to developing and extending the study to the so-called mean-field games and control 21 . • Optimal Control. This manuscript has focused primarily on analysis of the stochastic ensemble with a control.
However, the control in this setting was not optimal but rather preset. The natural evolution of this analysis (which would also complicate it) consists of a two-level formulation where solution of the problem analyzed here is also optimized. For example, one minimizes a cumulative cost including DR tasks (such as tracking time-evolving consumption signal from the system operator) and the mixing/recovery characteristics of the ensemble investigated above. • Discrete Phase Space. Given practical constraints in the device resolution, it is natural to reduce the hybrid (continuous-discrete) state space of the analyzed model to a purely discrete space simply by binning the temperature. Moreover, following the logic of 22,23 , it is practically appropriate to also consider the resulting Markov Process (MP) model in discrete time. In fact, this MP formulation is also practically advantageous for analysis of the aforementioned optimal control, where the problem becomes of the Markov decision process (MDP) type, as in 9,24,25 . We would argue that the MP and MDP approaches are naturally appropriate and algorithmically attractive to account for the randomization effects analyzed in the manuscript. • Data-Driven Control. Individual devices included in the aggregation may change their behavior, which then should be accounted for through data-driven identification of a device and ensemble parameters 26 . To track changes in real time and then account for them in the control, one would naturally resort to the data-driven approaches of the reinforcement learning type 27 , combining learning and control and aimed at developing on-line algorithms for optimal control.