Ensemble of Thermostatically Controlled Loads: Statistical Physics Approach

Thermostatically controlled loads, e.g., air conditioners and heaters, are by far the most widespread consumers of electricity. Normally the devices are calibrated to provide the so-called bang-bang control – changing from on to off, and vice versa, depending on temperature. We considered aggregation of a large group of similar devices into a statistical ensemble, where the devices operate following the same dynamics, subject to stochastic perturbations and randomized, Poisson on/off switching policy. Using theoretical and computational tools of statistical physics, we analyzed how the ensemble relaxes to a stationary distribution and established a relationship between the relaxation and the statistics of the probability flux associated with devices’ cycling in the mixed (discrete, switch on/off, and continuous temperature) phase space. This allowed us to derive the spectrum of the non-equilibrium (detailed balance broken) statistical system and uncover how switching policy affects oscillatory trends and the speed of the relaxation. Relaxation of the ensemble is of practical interest because it describes how the ensemble recovers from significant perturbations, e.g., forced temporary switching off aimed at utilizing the flexibility of the ensemble to provide “demand response” services to change consumption temporarily to balance a larger power grid. We discuss how the statistical analysis can guide further development of the emerging demand response technology.


A. Fokker-Planck Equations
In the hard model limit, i.e., at r → ∞, the FP Eqs. (8,9) of the main text turn into supplemented by the so-called absorbing boundary conditions i.e., conditions re-enforcing termination of a "particle"/air conditioner after it reaches its respective switching boundary. The "switching-induced" fluxes in Eqs. (8,9) of the main text are defined according to The "instantaneous" switch terms on the right-hand sides (rhs) of Eqs. (1,2,4) reflect the absorbing nature of the jumping process, i.e., the immediate jump/switch once the hard boundary for the respective ↑, ↓ state is achieved, and then Eq. (4) re-enforces relationships between switching rates and respective probability fluxes at the hard boundaries. Note that the current fluxes on the rhs of Eq. (4) do not include the "advection"-induced (τ -dependent) terms, because these vanish as a result of the boundary condition (4). Notice that a system of equations for the hard model equivalent to Eqs. (1,2,3,4) was derived for the first time in [1].

B. Stationary Probability Distributions
In the stationary state of the hard model, the fluxes of probability in the ↑ and ↓ states are piecewise-constant: where J is a positive constant, while P ↑ (x) at x ≥ x ↑ and P ↓ (x) at x ≤ x ↓ admit, in accordance with Eq. (11) of the main text, the following equilibrium (Boltzmann) forms: ∼ exp(−k(x − x − ) 2 /(2κ)). Then, solving Eq. (6) and accounting for the normalization condition, Eq. (12) of the main text, one arrives at the following explicit expressions where J is a positive constant, while P ↑ (x) at x ≥ x ↑ and P ↓ (x) at x ≤ x ↓ admit, in accordance with Eq. (11) of the main text, the following equilibrium (Boltzmann) forms: ∼ exp(−k(x − x − ) 2 /(2κ)). Then, solving Eq. (6) and accounting for the normalization condition, Eq. (12) of the main text, one arrives at the following explicit solution of the Fokker-Planck equations for the hard model . (10)
Substituting Eq. (12) into Eqs. (15,16), we arrive at four linear homogeneous relationships on four c constants: where M (λ) is thus an explicitly defined 4 × 4 matrix with coefficients that are parametrically dependent on λ. The allowed λ n are solutions of the zero-determinant condition for the respective matrix Even though the determinant is an explicit function of λ, resolving the zero-determinant equation explicitly does not seem feasible and we rely here on a numerical exploration of the zero-determinant Eq. (19), which is illustrated in Fig. (2) of the main text.
For the sake of completeness, notice also that the eigenvalue problem was analyzed in [2], where the author guessed that λ n = 2n/τ , making the first argument in the Kummer's confluent hypergeometric function a negative integer (in which case the hypergeometric functions are reduced to polynomials). However, a direct check shows that Eq. (19) is not satisfied for this guess at a finite κ.

D. Flux Generation in the Hard Model: Dynamical Viewpoint and First-Passage Time
Even though the relaxation dynamics of the hard (r = ∞) model and the diffusionless model are quite different, the LD function of the hard model and its relationship to the spectrum of the FP operator can be identified similarly to how they were identified for the zero diffusivity model in the main text. An appropriate dynamic object in the case of the hard model is the PDF of staying at the level σ for time t, P σ (t), which we can also refer to as the PDF of the first-passage time at the level σ. The PDF of making one full cycle in time t becomes are related to each other according to In the regime of weak diffusion, κ (x ↑ − x ↓ ) 2 /τ , the first-passage time distributions can be computed explicitly using spectral decomposition and utilizing the Wentzel-Kramers-Brillouin approximation of quantum mechanics to resolve the spectrum problem.

A. Relating Eulerian and Lagrangian Views
Explicit expressions discussed in the main text for the probability of advancing along the cycle in time t can also be used to compute an Eulerian object-the probability of observing a device in the state (x, σ) at a time t. This is achieved by relating the probability P (x, σ|x 0 , σ 0 ; t) of the system being in the state (x, σ) at time t, conditioned to be in state (x 0 , σ 0 ) initially, to the probabilityP (t|x 0 , σ 0 → x, σ) of reaching the state (x, σ), starting at (x 0 , σ 0 ), in time t. The former object is whereẋ(x, ↑) = −(x − x − )/τ andẋ(x, ↓) = −(x − x + )/τ are the velocities of the deterministic motion at point x on the discrete level σ. The latter object iŝ where t dc , defined in Eq. (5) of the main text, is the time needed for the device to complete one deterministic cycle within the comfort zone. P (t |x 0 , σ 0 → x, σ), entering Eq. (24), is the probability of the device transitioning from the state (x 0 , σ 0 ) to the state (x, σ) in time t without completing a single cycle. It is assumed in Eq. (24) that all the devices are in the (x ↑ , ↑) state initially, in order to simplify the resulting expressions. Notice that at t → ∞ the initial state will be forgotten, assuming sufficient mixing. It is important to emphasize that Eq. (23) plays a crucial technical role in our approach by relating the distribution of x on the left-hand side (lhs), which is the object of our interest, to the distribution of time on the rhs, which is easy to compute using the Laplace transform technique. Such a simple relationship via just a Jacobian is due to the deterministic nature of dynamics for a given discrete level; the jumps between the discrete levels are the only source of stochasticity.
Obviously, Eq. (24) simplifies when (x, σ) is chosen to coincide with the the initial state, (x ↑ , ↑). We will focus on this case, where P (t |x ↑ , ↑→ x ↑ , ↑) = δ(t ) and thus Eq. (24) becomes with u = |ẋ(x ↑ , ↑)|. Extending P out;n (t) with zero to negative times, t < 0, we can substitute m in Eq. (25) with ∞. Then, the Laplace transform (over time) version of Eq. (25) becomes given that F s (α)F s (β) < 1, where F s (α) is defined in Eq. (34) of the Method Section of the main text. Inverse Laplace transform applied to the rhs of Eq. (26) results in a sum over poles (in the domain of complex s), solving Following the same strategy and with a bit of additional effort, an explicit expression for P (x, σ|x 0 , σ 0 ; t) for a full range of arguments can be derived. Indeed, introducing the notation we arrive at the following extended version of Eq. (27): withP (x, σ|x 0 , σ 0 ; t) being the PDF of being in state (x, σ) at time t provided the starting state was (x 0 , σ 0 ) and a stochastic trajectory never went through state (x ↑ , ↑). A direct computation yields the following natural relationship for the eigenmodes: ξ σ,n (x) = ξ −λn (x, σ) and η σ,n (x) = η −λn (x, σ). Three follow-up remarks are necessary. First, we note that Eq. (28) is fully consistent with the spectral expression given by Eq. (18) of the main text when s → −λ n . To check the consistency, one just needs to change variables in the first and second integrals on the lhs of Eq. (18) of the main text according to x = x+α+x− exp(t/τ ) α+exp(t/τ ) and x = x−β+x+ exp(t/τ ) β+exp(t/τ ) , respectively. In fact, the dynamical derivation of Eq. (28) just presented should be considered as a rigorous proof of the spectral assumptions made in the main text.
Second, we observe that according to Eqs. (37,38) of the Methods Section of the main text, Eq. (28) and preceding remarks on the LD function, S(ω) is directly related to G s , fully defining the spectrum of the FP operator.
Finally, the sum over poles in our description is finite. We do not prove this fact here but instead illustrate it in Appendix III A on a simpler toy model assuming an "instantaneous escape" (from the uncomfortable zone).  1: Graphical solution of the spectral problem for the toy "instantaneous escape" (from the uncomfortable zone) model discussed in Appendix III A. We choose α and β parameters as described before. We also choose τ = 1, r = 2. Plotting conventions are the same as in Fig. 2 of the main text.

III. SOFT MODEL: THE CASE OF SMALL FINITE DIFFUSIVITY
The limit κ → 0 of small (but finite) diffusivity is important because noise/diffusion is always present in the system. The formal side of this statement is linked to the fact that the FP operator changes in transition from the diffusionless case, when it is of the first order (hyperbolic), to the small diffusion case, when it becomes of the second order (elliptic). This change in the operator order has principal consequences for the shapes of the stationary solution, the spectrum, and the dynamic behavior.
The goal of this Section is to discuss these principal changes briefly and informally, leaving more detailed mathematical analysis for future publications.
We start the discussion by mentioning that it is clear from analysis of the steady distribution conducted for the diffusionless regime above that accounting for small but finite diffusion will be important mainly at rτ < 1. Indeed, in this regime, the PDFs of the on/off states are chosen in the vicinity of x ± . On the other hand, when the diffusion is strictly zero the x ± limits are not achievable. Therefore, immediate vicinities of x ± are controlled solely by the diffusion.
Two temporal scales are significant in this regime: first, 1/r, which is the typical time for switching from on to off and vice versa, and second, the time for the system to stay within the on or off state. Notice that the latter time is not just τ , as one might naively suggest. Indeed, given that the majority of devices in the rτ < 1 regime get rather close to x ± while the dynamics in the vicinity of the equilibrium points slow down, the time for a device to traverse the range and then get inside the diffusion-controlled zone on the other end of the range is estimated to be τ log(|x + − x | / √ κτ ). This suggests that the dynamics, spectrum, and structure of eigenvalues, discussed in the main text for the diffusionless regime, are in fact realized within the following intermediate time/transient asymptotic: . We now discuss the long time asymptotic regime, t τ log(|x + − x − |/ √ κτ ), dominated by the diffusion. First, we notice that the stationary distribution, achieved as a result of a balance of the ballistic motion within a state and jumps between the states, is renormalized by small diffusion in the |x − x ± | ∼ √ κr vicinity of x ± . Similar renormalization of the x-shape by diffusion also applies to the rest of the spectrum (beyond the stationary solution). In fact, the spectrum itself is also modified in this long time regime. In particular, if r is sufficiently small, i.e., τ r log(|x + − x | / √ κτ ) 1, one observes a perfect equidistant spectrum λ n = n/τ , where n = 0, 1, · · · , corresponding to separate equilibrations within the σ =↑ and σ =↓ ensembles.

A. Spectral Properties of a Toy "Instantaneous Escape" Model
In this Section of the Supplementary Information, we discuss the spectrum in the simple case of a toy "instantaneous escape" (from the uncomfortable zone) model, assuming that the device returns to the control zone instantaneously. See Fig. (1) for an illustration. In this case, we have F s (α) = F s (β) = F s = − r r + s , and the spectral Eq. (28) simplifies to (s + r) 2 = r 2 e −st dc .
One observes that at r > r cr Eq. (32) has no real solutions, whereas at r ≤ r cr there is exactly one real solution, where r cr is given implicitly by (r cr t dc ) 2 = 4e −(rcrt dc +2) .
This implies, in particular, that r cr < t −1 dc . Next we analyze the convergence of the spectral decomposition of the PDF of finding the device in the state (x, σ) at time t given that the device was in the state (x 0 , ↓) at the moment 0: P (x, ↓ |x 0 , ↓; t) = n e −εnt−iωn ξ ↓,n (x)η ↓,n (x 0 ), where ξ and η are the respective left and right eigenvalues of the FP operator of the toy model. Substituting s with −λ n in Eq. (32) and analyzing the large n part of the spectrum, one derives e −λnt ∼ e ±iπ(2n+1)(t/t dc ) rt dc π(2n + 1)
It follows from Eq. (35) that for given values of x and x 0 the spectral decomposition of P (x, ↓ |x 0 , ↓; t) converges at sufficiently large t. However, if the initial distribution P ↓,0 (x) is represented by a smooth function with compact support, the spectral decomposition converges at all times.