Low-complexity control algorithm for decentralised demand response using thermostatic loads

TL;DR

This paper presents a low-complexity decentralized control algorithm for thermostatic loads like refrigerators, enabling efficient aggregate demand management suitable for embedded systems and large-scale simulations.

Contribution

It adapts an existing continuous control approach to a discrete, low-complexity algorithm for decentralized demand response in thermostatic loads.

Findings

Accurately controls aggregate power consumption of refrigerator populations.

Demonstrates high computational efficiency in simulations.

Handles variable-length discrete time steps effectively.

Abstract

Thermostatically controlled loads such as refrigerators are exceptionally suitable as a flexible demand resource. This paper derives a decentralised load control algorithm for refrigerators. It is adapted from an existing continuous time control approach, with the aim to achieve low computational complexity and an ability to handle discrete time steps of variable length -- desirable features for embedding in appliances and high-throughput simulations. Simulation results of large populations of heterogeneous appliances illustrate the accurate aggregate control of power consumption and high computational efficiency.

Equations

$$\frac{\mathrm{d}{T}^a(t)}{\mathrm{d}t}=-{\alpha }^a\left[T^a(t)-T_{\mathrm{off}}^a+c^a(t)\cdot \left(T_{\mathrm{off}}^a-T_{\mathrm{on}}^a\right)\right]$$

$$P^a(t)=P_{\text{on}}^a{c}^a(t).$$

$$E_{{\mathcal{S}}^a}[P^a(t)]=P_0^a.$$

$$P(t)=\sum _{a\in \mathcal{A}}{P}^a(t)$$

$$E_{{\mathcal{S}}^a}[P(t)]=\Pi (t)P_0^a.$$

$$P(t)=\Pi (t)\sum _{a\in \mathcal{A}}{P}_0^a+O(\sqrt{|\mathcal{A}|}),$$

$$f_0(T;{\mathcal{M}}^a)=\frac{k^a}{(T_{\mathrm{off}}^a-T)(T-T_{\mathrm{on}}^a)}$$

$$k^a=\frac{T_{\mathrm{off}}^a-T_{\mathrm{on}}^a}{\log \left(\frac{(T_{\max }^a-T_{\mathrm{on}}^a)(T_{\min }^a-T_{\mathrm{off}}^a)}{(T_{\min }^a-T_{\mathrm{on}}^a)(T_{\max }^a-T_{\mathrm{off}}^a)}\right)}$$

$$\Pi (t)={\Pi }_i,\text{for }t\in (t_{i-1},t_i].$$

$$\overline{T}(t_i)$$

$${\overline{T}}_0=T_{\mathrm{off}}-k\times \log \left(\frac{T_{\max }-T_{\mathrm{on}}}{T_{\min }-T_{\mathrm{on}}}\right)$$

$$z(t)=\frac{{\overline{T}}_0-\overline{T}(t)}{T_{\mathrm{off}}-{\overline{T}}_0},$$

$$R(t)=\{\begin{array}{ll}{T}_{\max },& \text{if }\overline{T}(t)\ge {\overline{T}}_0\\ T_{\min },& \text{if }\overline{T}(t)<{\overline{T}}_0\end{array}$$

$$\beta (t;R)$$

$$=\frac{(\Pi (t)-1)-z(t)}{z(t)-\zeta (R(t))}$$

$$\zeta (R)=\frac{{\overline{T}}_0-R}{T_{\mathrm{off}}-{\overline{T}}_0}.$$

$$s(t)$$

$$=1-z(t)/\zeta (R(t)).$$

$$z_i=\sum _{j=0}^i({\Pi }_j-1)\left({\mathrm{e}}^{-\alpha (t_i-t_j)}-{\mathrm{e}}^{-\alpha (t_i-t_{j-1})}\right).$$

$$z_i=z_{i-1}{\mathrm{e}}^{-\alpha \Delta t_i}+({\Pi }_i-1)(1-{\mathrm{e}}^{-\alpha \Delta t_i}).$$

$$\hat{R}(t)=R_i,\text{for }t\in (t_{i-1},t_i].$$

$$R_{i+1}=\{\begin{array}{ll}{T}_{\max },& \text{if }{z}_i\le 0,\\ T_{\min },& \text{if }{z}_i>0.\end{array}$$

$$R_i^{-}$$

$$R_i^{+}$$

$${\zeta }_i^{\pm }$$

$$s_i^{\pm }$$

$$(1-w){\overline{T}}_0+w{T}_{\min }\le \overline{T}(t)\le (1-w){\overline{T}}_0+w{T}_{\max }.$$

$$w\zeta (T_{\max })\le z(t)\le w\zeta (T_{\min }).$$

$${\Pi }_{i+1}:=\{\begin{array}{ll}\max ({\Pi }_{i+1},1+w\zeta (T_{\max })),& \text{if }{z}_i\le w\zeta (T_{\max })\\ \min ({\Pi }_{i+1},1+w\zeta (T_{\min })),& \text{if }{z}_i\ge w\zeta (T_{\min })\\ {\Pi }_{i+1},& \text{otherwise}\end{array}$$

$$\begin{gathered} \left(\frac{{\overline{T}}_0-T_{\min }}{T_{\max }-T_{\min }}\right)\left(\frac{T_{\mathrm{off}}-T_{\max }}{T_{\mathrm{off}}-{\overline{T}}_0}\right)\le {\Pi }_{i+1}\le \\ \left(\frac{T_{\mathrm{off}}-T_{\max }}{T_{\mathrm{off}}-{\overline{T}}_0}\right)+\frac{(T_{\max }-{\bar{T}}_0)(T_{\max }-T_{\mathrm{on}})}{(T_{\max }-T_{\min })(T_{\mathrm{off}}-{\bar{T}}_0)}. \end{gathered}$$

$$\begin{gathered} \left(\frac{T_{\max }-{\overline{T}}_0}{T_{\max }-T_{\min }}\right)\left(\frac{T_{\mathrm{off}}-T_{\min }}{T_{\mathrm{off}}-{\overline{T}}_0}\right)\le {\Pi }_{i+1}\le \\ \left(\frac{T_{\mathrm{off}}-T_{\min }}{T_{\mathrm{off}}-{\overline{T}}_0}\right)+\frac{({\bar{T}}_0-T_{\min })(T_{\min }-T_{\mathrm{on}})}{(T_{\max }-T_{\min })(T_{\mathrm{off}}-{\bar{T}}_0)}. \end{gathered}$$

$$T_{\mathrm{low}}(t)=R(t)-(R(t)-T_{\min })s(t),$$

Full text

I Introduction

The physical characteristics of refrigerators and other thermostatically controlled loads (TCLs) make them exceptionally suitable as a low-cost provider of flexibility to the grid: their power consumption can be shifted by tens of minutes without noticeable effects on cooling performance. This flexibility can then be used for the provision of response and reserve services, to reduce extreme load levels and to alleviate ramping constraints [1]. Given the large number of devices involved, large populations can be effectively controlled in a decentralised manner using randomised control schemes, as proposed e.g. in [2, 3, 4]. Practical implementation must also consider constraints on implementation, computation and communication, as discussed e.g. in [5, 6].

A robust decentralised control scheme for heterogeneous TCLs was introduced in [7], and extended in [8] to allow for scenarios in which TCLs collectively absorb energy from the grid. The control strategy has the desirable feature that it requires only one-way broadcast information, yet achieves tracking of a reference signal that is exact in expectation (i.e. exact for large numbers) without violating temperature limits in individual devices. However, its continuous time formulation in integral form was not conducive to implementation on embedded controllers or for rapid simulation of many devices at once.

This paper addresses that shortcoming by deriving a discrete-time control algorithm that implements the control strategy that was presented in [7, 8]. The algorithm is particularly suitable for implementation on devices with computational constraints. Specifically, it avoids numerical integration and uses greedy time steps of variable size, so that real-time performance requirements are relaxed. First, a discretisation procedure is described, for the case of a piecewise-constant control signal, and the controller is reformulated in natural coordinates. Then, expressions are derived for each of the on/off-switching processes governing the behaviour of the device. Finally, an explicit control algorithm is provided, and its efficacy is illustrated with Python-based simulations for a heterogeneous population of refrigerators.

II Preliminaries

II-A Appliance model

Throughout this paper, we consider the first order TCL model, expressed by the following differential equation for the temperature $T^{a}$ of the compartment of appliance $a$:

[TABLE]

For ease of exposition, we shall refer to refrigerating applicances throughout, although the same model and control strategy can be used for other TCLs, e.g. space heaters. Here, $c^{a}(t)\in\{0,1\}$ is the state of the compressor, and $T^{a}_{\mathrm{off}}$ and $T^{a}_{\mathrm{on}}$ are the asymptotic temperatures in the off and on states, respectively. The power consumption $P^{a}(t)$ of the appliance is assumed to be dominated by the compressor power consumption $P^{a}_{\textrm{on}}$, so that

[TABLE]

In the steady state (no control actions), the appliance is subject to a hysteresis controller that switches to the on state whenever an upper temperature bound $T_{\mathrm{max}}^{a}$ is reached, and to the off state when $T_{\mathrm{min}}^{a}$ is reached. This results in a periodic cycling of the power consumption, with an average power level $P^{a}_{0}$. Let us consider that each appliance $a$ has a model $\mathcal{M}^{a}=\{\alpha^{a},P^{a}_{\textrm{on}},,T_{\mathrm{off}}^{a},T_{\mathrm{on}}^{a},T_{\mathrm{min}}^{a},T_{\mathrm{max}}^{a}\}$ that is known to us, but an internal state $\mathcal{S}^{a}=\{T^{a}(t),c^{a}(t)\}$ that is not. We assume the latter is distributed according to a steady state distribution. Then, the power consumption is in effect a random process, with the expectation (at each time $t$)

[TABLE]

II-B Aggregate power modulation

The objective of the TCL demand response controller is to control the aggregate power consumption

[TABLE]

of a collection of appliances $\mathcal{A}$. In this paper we consider the control approach introduced in [7], which modulates the power consumption using a broadcast reference signal $\Pi(t)$ in the following way. Each device individually adapts its power consumption to $\Pi(t)$ in order to satisfy

[TABLE]

Clearly, $\Pi(t)=1$ represents the steady state (3) and changes in $\Pi(t)$ are immediately reflected in the expected power consumption. Moreover, the controller maintains independence between appliances (conditional on the control signal $\Pi(\cdot)$) so that the central limit theorem can be applied to the total power consumption (4), resulting in

[TABLE]

where the last term is a random process that decreases in relative importance to the first term as the set of appliances increases. We note that this is the case even for heterogeneous appliances.

The ability to closely track a reference signal was first demonstrated in [7]. In [8], the control signal was generated using a mixture of off-line scheduling and real-time control, and in [9] various frequency-sensitive controllers to locally compute $\Pi(t)$ (e.g. a simple droop controller) were implemented.

II-C Distribution-referred control

The control approach introduced in [7] and extended in [8] isdistribution-referred approach, because it implements control via the probability distribution of temperatures for appliances with a known model $\mathcal{M}^{a}$ and unknown state $\mathcal{S}^{a}$. It defines a family of alternative distributions $f_{z}(T;\mathcal{M}^{a})$ that varies continuously in the parameter $z$, containing as a special case the steady state temperature distribution [7, Eqs. (31)-(32)]

[TABLE]

with

[TABLE]

The controller consists of two major elements, which are evaluated by each appliance in order:

Choose the distribution parameter $z(t)$ such that the power consumption tracks the reference signal $\Pi(t)$ according to (5). Determine the collective device switching actions required to keep device temperatures aligned with $f_{z(t)}(T)$, and identify the temperature limits $T_{\mathrm{min}}(t)\geq T_{\mathrm{min}}$ and $T_{\mathrm{max}}(t)\leq T_{\mathrm{max}}$. 2. 2.

Based on the actual appliance state $\mathcal{S}^{a}$, compute stochastic control actions, in the form of on/off switching. Switching events can be initiated in three distinct ways:

• •

Deterministic switching when temperature limits $T_{\mathrm{min}}(t)$ or $T_{\mathrm{max}}(t)$ are exceeded.

• •

Continuous-time stochastic switching at intermediate temperatures in order to shape the temperature distribution.

• •

Instantaneous stochastic switching on discrete changes of power setpoints, or when the controller switches between energy-provision and energy-absorption modes [8].

The temperature distribution and appliance switching phases for the discrete time control strategy are addressed in Sections III and IV, respectively. We will henceforth drop the appliance superscript $a$, because the control steps are executed locally within each appliance (or independently for each appliance in a simulation). Note that this implies a single model $\mathcal{M}$ is used in the derivations, but the results remain valid for portfolios of heterogeneous devices (each with their own model).

II-D Discretisation procedure

Moving from a continuous time formulation to a discrete time formulation, we consider a partitioning of the timeline by the ordered sequence of times $\{t_{i}\}$, indexed by the integer $i$, at which the controller is invoked. They define time intervals $(t_{i+1},t_{i}]$ with durations $\Delta t_{i}=t_{i}-t_{i-1}$. Note that the duration $\Delta t_{i}$ refers to the interval prior to $t_{i}$, and the intervals may have variable size. The reference signal $\Pi(t)$ is assumed to be piecewise constant, defined by

[TABLE]

The controller thus receives at $t_{i}$ a new reference power level $\Pi_{i+1}$ that must be applied for the upcoming interval $(t_{i},t_{i+1}]$.

Although the discontinuous changes of reference power will trigger switching events at $t_{i}$, the other switching events may occur at any time $t$. In the discretised approximation of the continuous time controller, they will be synchronised with the control execution times $t_{i}$ as follows. It is assumed that switching is immediate (at $t_{i}$).

• •

A violation of the temperature limits will trigger corrective switching as soon as it is detected.

• •

Switching due to toggling between energy-provision and energy-absorption modes is implemented as soon as a change in regime is detected.

• •

Continuous time stochastic switching is implemented by approximating the integrated switching rate (i.e. the switching probability) over $\Delta t_{i}$ using the trapezoidal method, and executing any switching events at $t_{i}$ (the end of the interval).

The algorithm is thus implemented in a ‘backward’ fashion, meaning that at time $t_{i}$, the algorithm implements switching actions resulting from reference changes at $t_{i}$, and those accumulated over the preceding interval $(t_{i-1},t_{i}]$. The advantage of this approach is that the interval $\Delta t_{i}$ can be chosen opportunistically: the controller does not need to know in advance when the next time step will take place. This is convenient, for example when computational limitations cause a delay in intended invocation time, or when the time step adapts to sudden changes in grid frequency. It should be pointed out that this ‘backward’ integration does not delay the response to changes in reference power, which is implemented immediately at $t_{i}$.

III Distribution behaviour

This section focuses on the first part of the controller. It computes the desired evolution of the probability distribution of temperatures of fridges with model $\mathcal{M}$, when tracking a piecewise constant reference $\Pi_{i}$. The derivation is initially performed in continuous time. The results are subsequently expressed in natural coordinates and restated in a form that is suitable for discrete-time evaluation.

III-A Aggregate physics

The average temperature of a TCL population is affected by the desired power consumption $\Pi(t)$ according to [7, Eq. (26)]. With the convention that $t_{-1}=-\infty$ and $\Pi_{0}=1$ (assuming an initial steady state), it follows that

[TABLE]

where the steady state average temperature $\overline{T}_{0}$ is computed using [7, Eqs. (23) and (32)-(33)] as

[TABLE]

with $k$ defined in (8). We define the dimensionless variable

[TABLE]

to parameterise the distributions $f_{z(t)}(T)$, and simplify the notation in what follows. Note that it is related to the $\sigma$ variable used in [8] as $z=\sigma-1$.

III-B Controller modes

The algorithm in [7] implicitly generates the family of temperature distributions $f_{z(t)}(T)$ by the net heating rate $v(T,t)$, which is determined by averaging over devices in the off (heating) and on (cooling) states at time $t$ and temperature $T$. The heating rate is controlled by a parameter $\beta(t)$ through $v(T,t)=\alpha\beta(t)(T-T_{\mathrm{max}})$. The effect of this heating rate profile is a temperature distribution that contracts to the pivot temperature $T_{\mathrm{max}}$ in order to provide energy to the grid - and reverses this process to recover the energy supplied. In [8] it was coupled to a ‘mirrored’ controller that is capable of absorbing energy from the grid by contracting to the pivot temperature $T_{\mathrm{min}}$. Switching between the two controller modes takes place whenever $\overline{T}(t)$ crosses $\overline{T}_{0}$ (when $z(t)$ crosses 1).

A generalised formulation covering both regimes is obtained by defining a heating rate of the form $v(T,t;R)=\alpha\beta(t;R)(T-R(t))$, where $R\in\{T_{\mathrm{min}},T_{\mathrm{max}}\}$ is a reference temperature, which acts as a pivot temperature for the controller, with the property $v(R,t;R)=0$. The reference temperature is defined as follows:

[TABLE]

III-C Control parameter

The control parameter $\beta(t;R)$ is determined by the desired reference power $\Pi(t)$ according to [7, Eq. (36)]:

[TABLE]

where

[TABLE]

The denominator in the definition of $\beta$ reflects, in dimensionless form, the energy limits of the TCL aggregate. Note also that $\beta$ switches sign depending on the value of $R(t)$.

III-D Distribution scaling

The controller has the effect of scaling the steady state temperature distribution $f_{0}(T)$ around the pivot temperature $R(t)$, such that the distribution does not exceed the temperature bounds $T_{\mathrm{min}}$ and $T_{\mathrm{max}}$ [7]. The extent of this scaling at time $t_{i}$ is compactly represented by the scale parameter

[TABLE]

III-E Discretisation

We now consider the restriction of the continuous time controller to the set of discrete times $t_{i}$. We replace the coordinate $z(t)$ by its discretisation $z_{i}=z(t_{i})$, which is computed from (10) as

[TABLE]

Updates to $z_{i}$ are efficiently implemented using $z_{0}=0$ (for a steady state initialisation) and the recursive relation

[TABLE]

The discretised controller switches modes only at instants $t_{i}$, so $R(t)$ is approximated by the delayed function

[TABLE]

with

[TABLE]

Because our analysis focuses on the control time $t_{i}$, where $\hat{R}(t)$ and $\Pi(t)$ are potentially discontinuous, we introduce $\pm$-notation for the left and right limits at $t_{i}$:

[TABLE]

Similar definitions using left and right limits naturally apply to $\zeta(R)$, $s(t)$ and $\beta(t;R)$:

[TABLE]

III-F Energy and power constraints

The ability of the aggregate appliances to sustain a low or high power level is determined by operating temperature bounds of the appliance, applied to the distribution-averaged temperature: $\overline{T}(t)\in(T_{\mathrm{min}},T_{\mathrm{max}})$ (no feasible solutions for the distribution $f_{z}(T)$ exists outside of this domain). However, because operation near the limits is infeasible in practice, due to diverging switching rates, we shall use a restricted range of operating temperatures that is scaled with a fraction $w<1$ around the steady state operating temperature $\overline{T}_{0}$:

[TABLE]

Rewriting this in terms of $z(t)$ and $\zeta(\cdot)$, we get

[TABLE]

Small excursions out of this temperature band will be permitted, but if this happens, the requested power level $\Pi_{i+1}$ will be restricted to not exacerbate the excursion, using the relation (17). This leads to the update rule for $\Pi_{i+1}$:

[TABLE]

In addition to energy constraints related to the distribution-averaged temperature, the controller is subject to instantaneous power constraints that result from the maximum rate of change of the distribution. If the controller is in energy provision mode ($z_{i}\leq 0$), these power constraints are given by [8]

[TABLE]

If the controller is in energy absorption mode ($z_{i}>0$), the power constraints are given by

[TABLE]

IV Device switching

The desired evolution of the temperature distribution can be used to compute the necessary control actions of individual appliances. This section identifies such control actions using the three types of switching events identified in Section II-C. These are computed as a function of the time of evaluation $t_{i}$, the compressor state $c_{i}\in\{0,1\}$ during the preceding interval $(t_{i-1},t_{i}]$, and the current device temperature $T_{i}$ (assumed to be measured in the appliance at time $t_{i}$).

IV-A Forced switching

TCLs are forced to switch on or off when their temperatures exceed the permitted interval $[T_{\mathrm{low}}(t),T_{\mathrm{high}}(t)]$. From the linear scaling of the temperature distributions around the pivot temperature $R(t)$ with a factor $s(t)$, it follows that

[TABLE]

At $t_{i}$, the refrigerator must act if these bounds are violated at the start of the next time interval:

[TABLE]

IV-B Continuous-time switching

We now consider the continuous-time stochastic switching rates from on to off states ($r^{1\rightarrow 0}(t)$) and vice versa ($r^{0\rightarrow 1}(t)$), required to maintain the desired shape of the temperature distribution. The switching rates for the energy provision mode are defined in [7, Eqs. (48)-(52)]. Here, we generalise these expressions to cover both energy provision and absorption modes ($R(t)\in\{T_{\mathrm{min}},T_{\mathrm{max}}\}$) and simplify them using the $z$-coordinate transformation. Finally, we specialise the expressions for trapezoidal integration with piecewise constant power references.

The derivative of $\beta$ can be simplified by substitution using (12) and (15), resulting in:

[TABLE]

This substitution can be used in [7, Eq. (51)] (adjusted for general $R$). Further simplification follows from setting $\mathrm{d}\Pi(t)/\mathrm{d}t=0$ (because we consider piecewise constant sections between $t_{i}$). We compute the intermediate quantity $\Xi(t)$, using the identity found in [7, Eq. (38)], again taking left and right limits due to discontinuity at $t_{i}$.

[TABLE]

The stochastic transition rates at the left and right limits to $t_{i}$ are computed from $\Xi^{\pm}_{i}$ using [7, Eq. (49)] (adjusted for general $R$), resulting in

[TABLE]

Midpoint integration between adjacent time instants $t_{i-1}$ and $t_{i}$ is used to determine the resulting switching probabilities, where switching is implemented at $t=t_{i}$:

[TABLE]

Note that the rates at both ‘inner’ edges of the interval $\Delta t_{i}$ are used: the ‘$+$’ side at $t_{i-1}$ and the ‘$-$’ side at $t_{i}$.

IV-C Instantaneous switching

Finally, consider the instantaneous stochastic switching at time $t_{i}$ due to mode changes (energy absorption, energy delivery) or changes in $\Pi(t)$. This results in a discontinuous change in the net heating rate $v(T,t)$, which can only be achieved by a fraction of devices switching on or off at $t_{i}$. Following [8], we compute the probability of switching from the on to off state at time $t_{i}$, for a refrigerator that is currently on, as

[TABLE]

IV-D Combined stochastic switching

Formally, the continuous-time (28) and instantaneous (29) switching probabilities should be evaluated in sequence, because the former occurs during the interval $(t_{i-1},t_{i}]$ and the latter at time $t_{i}$. This would account for the possibility that an appliance switches off and on again within a single interval, or vice versa. Here, we assume that the switching probability associated with the continuous-time process is small to allow us to evaluate both probabilities in a single step.

[TABLE]

V Algorithm and results

The discrete time algorithm for updating the compressor state derived in sections III and IV, is summarised in pseudocode in Algorithm 1. The algorithm was implemented in Python 3.7 using the numba package to benefit from just-in-time compilation for considerable speedups.

For simulations, thermal model parameters were taken from [8, domestic refrigerator class]: $\alpha=1/7200s$; $T_{\mathrm{max}}=7\degree C$; $T_{\mathrm{min}}=2\degree C$; $T_{\mathrm{on}}=-44\degree C$; $T_{\mathrm{off}}=20\degree C$. Heterogeneous appliances were generated from these parameters by individually multiplying them with a random factor that was uniformly distributed between $0.8$ and $1.2$. All appliances had a maximum power consumption $P^{a}_{\mathrm{on}}=70W$ and operating range $w=0.9$ (not binding for the parameters used). Each appliance was randomly initialised as follows. The compressor was set to the on state with a probability equal to the steady state duty cycle $P^{a}_{0}/P^{a}_{\mathrm{on}}$ and the temperature was initialised according to the steady state probability distributions $f_{0}(T|c^{a}_{0}=1)\propto 1/(T-T_{\mathrm{on}})$ and $f_{0}(T|c^{a}_{0}=0)\propto 1/(T_{\mathrm{off}}-T)$.

Figure 1 shows a reference signal (top), with a length of 5 hours, that demonstrates a variety of features. The middle and bottom panels show the compressor state $c_{i}$ and temperature $T_{i}$, respectively, of a single appliance that tracks the reference signal. These results illustrate the apparently weak relation between the reference signal and single device dynamics (middle), and the ability of the controller to strictly respect the temperature bounds (bottom).

Next, heterogeneous populations of $1{,}000$ and $100{,}000$ appliances were simulated, tracking the same reference signal. Figure 2 illustrates the convergence of the aggregate response to the reference signal as the number of independent appliances increases. The top panel shows absolute power consumption per appliance; the bottom panel the deviation from the reference.

VI Conclusions and future work

This paper has derived a discrete time TCL controller for decentralised demand response. The results illustrate the ability to accurately track a reference signal with a large population of heterogeneous appliances. Moreover, Algorithm 1 has low computational complexity, which permits implementation on embedded hardware with severe computational constraints, or it can be used to achieve efficient simulations. The simulation of 100,000 devices for 5 hours using 10s time steps took only 36 seconds (using an Intel i5-7360U CPU under macOS 10.14.2). Moreover, the ability to use variable time steps can further alleviate real time constraints.

Lab testing of the algorithm in a modified refrigerator is currently ongoing. Both in the lab and using simulations, it is of interest to investigate the robustness of the controller against perturbations from the idealised setting [10]. Initial simulation experiments suggest that the performance of the controller is quite robust to misspecification of the thermal model. Nevertheless, it is interesting to consider how the controller could be enhanced with a means for an appliance to learn and test its own thermal model.