Decentral Smart Grid Control

DR may have huge benefits in future energy systems see, e.g., [7, 8, 25, 29, 30]. Here, we briefly summarize the most important economic aspects, following [7], which hold regardless of the technical implementation: (1) consumers may reduce their electricity bill by shifting their demand to periods of low prices. (2) In addition, DR may reduce the global costs of the power system as it allows for a more efficient use of the existing infrastructure and avoids costs for additional infrastructure. (3) DR may improve system stability by avoiding dangerous peaks of the demand and thus reduce the probability of power outages. (4) Finally, DR may improve market performance and reduce the price volatility. In addition to these points, DR is particularly important for future power grids because its implementation can potentially allow a higher penetration with renewable energies [31].

In the current proposal of a Decentral Smart Grid Control there is no central computer which controls the demand of the consumers and no central exchange to determine the electricity price. The control is realized in a decentralized way using local frequency measurements, thus requiring no long-distance communication. Is this sufficient to provide an efficient market, i.e., to reach an economic equilibrium?

To answer this question we first note that the stable stationary operation of a power grid requires that all machines rotate in synchrony, i.e., the frequency must be the same everywhere

$$\begin{eqnarray}&&{{\omega }_{j}}(t)\mathop{=}\limits^{!}\langle \omega \rangle \quad {\rm for}\;{\rm all}\;j\in \{1,\ldots ,N\}.\end{eqnarray} \tag{ 8 }$$

Otherwise the power flow between two nodes j and k

$$\begin{eqnarray}&&{{P}_{jk}}(t)={{K}_{jk}}{\rm sin} ({{\theta }_{k}}(t)-{{\theta }_{j}}(t)),\end{eqnarray} \tag{ 9 }$$

would be oscillating and average out over time. The synchronous state must be dynamically stable for small disturbances to be damped out [15, 32] and the grid may self-organize to a synchronous state with steady power flows [18, 22, 33]. Substituting the condition ${{\omega }_{j}}(t)=\langle \omega \rangle$ into the equations of motion (1) shows that the synchronous state is determined by the equation

$$\begin{eqnarray}&&{{\kappa }_{j}}\langle \omega \rangle ={{S}_{j}}({{p}_{j}})-{{D}_{j}}({{p}_{j}})-\mathop{\sum }\limits_{k=1}^{N}{{K}_{jk}}{\rm sin} ({{\theta }_{k}}-{{\theta }_{j}})\quad {\rm for}\;{\rm all}\;j\in \{1,\ldots ,N\}.\end{eqnarray} \tag{ 10 }$$

Summing up the equations for all j and using ${{K}_{ij}}={{K}_{ji}}$ yields

$$\begin{eqnarray}&&\mathop{\sum }\limits_{j}{{S}_{j}}({{p}_{j}})=\mathop{\sum }\limits_{j}{{D}_{j}}({{p}_{j}})+\mathop{\sum }\limits_{j}{{\kappa }_{j}}\langle \omega \rangle .\end{eqnarray} \tag{ 11 }$$

This shows that a dynamical equilibrium of the grid also implies the economic equilibrium of the market, i.e., the supply equals the demand including transmission losses. Hence, we have to analyze the dynamic stability of the combined techno-economical system to evaluate its stability properties. This will be done in detail in section 4.

For now, we assume that the grid is in equilibrium with $\langle \omega \rangle =0$ at a price ${{p}_{\Omega }}$. To analyze the stability of this equilibrium, we linearize the supply and demand curves close to the equilibrium:

$$\begin{eqnarray}\begin{array}{rcl} {{S}_{j}}({{p}_{j}}) & = & {{S}_{j}}({{p}_{\Omega }})+{\underbrace{{{\left. \frac{{\rm d}{{S}_{j}}}{{\rm d}p} \right|}_{{{p}_{\Omega }}}}}_{=:{{s}_{j}}}}({{p}_{j}}-{{p}_{\Omega }}) \\ {{D}_{j}}({{p}_{j}}) & = & {{D}_{j}}({{p}_{\Omega }})+{\underbrace{{{\left. \frac{{\rm d}{{D}_{j}}}{{\rm d}p} \right|}_{{{p}_{\Omega }}}}}_{=:{{d}_{j}}}}({{p}_{j}}-{{p}_{\Omega }}). \\ \end{array}\end{eqnarray} \tag{ 12 }$$

The modeling or measurement of supply and demand curves then reduces to the measurement of the elasticity of supply and demand. Generally, the supply increases with the price while the demand decreases such that ${{s}_{j}}\geqslant 0$ and ${{d}_{j}}\leqslant 0$ and thus in particular

$$\begin{eqnarray}&&{{s}_{j}}-{{d}_{j}}\gt 0\ {\rm for}\;{\rm all}\;j.\end{eqnarray} \tag{ 13 }$$

Here, we used course-graining, i.e., not every consumer is represented by one demand function but several consumers are aggregated to form one node in the network and hence supply and demand curves are assumed to be smooth.

The new aspect of our proposal of a Decentral Smart Grid Control is the direct encoding of the electricity price in the grid frequency. Thus, the grid frequency must be allowed to vary in certain boundaries so that fluctuations of the price are directly related to fluctuations of the frequency. Currently, frequency variations are limited due to technical reasons [34]. In the European grid $\pm 200\;{\rm mHz}$ are acceptable in normal operation and up to $\pm 800\;{\rm mHz}$ can occur in extreme cases for short times before emergency measures such as load shedding are initiated [35]. This sets the order of magnitude at which the frequency may vary. Accordingly, we consider a frequency range of $(50\pm 0.5)\;{\rm Hz}$.

In figure 1 we analyze the possible variations of the price and the frequency in more detail. Panel (a) shows a histogram of plausible values of the consumer price, if this price is strictly coupled to the variable spot market price for Germany in 2012 [36]. To obtain plausible consumer prices, we add 9 ct kWh⁻¹ for distribution and service, 7 ct kWh⁻¹ fees plus 19% VAT on the total. The variations of the electricity price directly relate to variations of the grid frequency as described above. We consider a linear price curve for all nodes

$$\begin{eqnarray}&&p({{\bar{\omega }}_{j}})={{p}_{\Omega }}-\epsilon \times {{\bar{\omega }}_{j}},\end{eqnarray} \tag{ 14 }$$

with ${{p}_{\Omega }}=24.1$ ct kWh⁻¹, as shown in panel (b1). For a slope $\epsilon =10$ (ct kWh⁻¹)/2π Hz, the price curve maps the operational range $({{\bar{\omega }}_{j}}+\Omega )/2\pi \in [49.5,50.5]$ Hz to a price interval $p\in [19.1,29.1]\;{\rm ct}\;{\rm kW}{{{\rm h}}^{-1}}$, which covers 98% of the observed fictitious consumer prices. A histogram of the resulting frequency variations is shown in panel (c1). In 2% of all time slots prices outside of this interval were recorded which can even become negative.

Zoom In Zoom Out Reset image size

To treat such events accordingly, a nonlinear function must be chosen which maps a fixed finite frequency interval to all possible prices, i.e., to the real line. Still, the slope of this function should be bounded around the operational point $\Omega =50$ Hz. These requirements are satisfied by an inverse sigmoidal function. As a particular example, we here consider the function

$$\begin{eqnarray}&&p({{\bar{\omega }}_{j}})={{p}_{\Omega }}-\frac{\alpha }{2}{{{\rm tanh} }^{-1}}(\frac{{{{\bar{\omega }}}_{j}}}{\pi \beta }),\end{eqnarray} \tag{ 15 }$$

which maps all angular frequency deviations in the interval ${{\bar{\omega }}_{j}}\in (-\pi \beta ,+\pi \beta )$ to a price $p\in (-\infty ,+\infty )$. The operational range can thus be fixed beforehand, and emergency measures can be specified, if frequencies outside this range are measured. Choosing $\alpha =\beta \epsilon$ yields the same slope of the price curve around the reference frequency as the linear price function (14). Indeed, using $\beta =1$ Hz the statistics of the observed frequencies in figure 1 hardly change in comparison to the linear price function, see panels (b2) and (c2). The corresponding frequencies change significantly only for extreme price events, which now map to the operational range $({{\bar{\omega }}_{j}}+\Omega )/2\pi \in [49.5,50.5]$ Hz as desired. A similar statistic is found for other sigmoidal functions. The precise form of such a nonlinear price function can be designed by the grid operators on the basis of actual market and consumption data.

We note that figure 1 is based on the German spot market prices, which show huge fluctuations compared to other energy markets [39]. Furthermore, one major effect of a comprehensive DR would be to suppress extreme price fluctuations anyway [7, 14]. Hence, we expect that the fluctuations shown in the figure represent extremes, and that they would be significantly smaller in energy systems with DR and (virtual) storage.

We first investigate a DR that is fast compared to the grid dynamics. Assuming an instantaneous adaption of the demand, i.e., $T=\tau =0$ in equations (6) and (7), the effective power $P_{j}^{{\rm mech}}(t)$ in equation (1) becomes a function of the current angular frequency deviation ${{\omega }_{j}}$. In particular, we consider a linear relation of price and angular frequency deviation

$$\begin{eqnarray}&&{{p}_{j}}:={{p}_{\Omega }}-\epsilon {{\omega }_{j}}(t).\end{eqnarray} \tag{ 16 }$$

Linearizing the supply and demand curves around ${{p}_{\Omega }}$ as in equation (12) yields

$$\begin{eqnarray}&&{{M}_{j}}\frac{{{{\rm d}}^{2}}{{\theta }_{j}}}{{\rm d}{{t}^{2}}}+{{\kappa }_{j}}\frac{{\rm d}{{\theta }_{j}}}{{\rm d}t}={{S}_{j}}({{p}_{\Omega }})-{{D}_{j}}({{p}_{\Omega }})-\epsilon ({{s}_{j}}-{{d}_{j}})\frac{{\rm d}{{\theta }_{j}}}{{\rm d}t}-\mathop{\sum }\limits_{k=1}^{N}{{K}_{jk}}{\rm sin} ({{\theta }_{k}}-{{\theta }_{j}}).\end{eqnarray} \tag{ 17 }$$

By equation (13), an instantaneous economic response thus increases the effective damping constant to

$$\begin{eqnarray}&&\kappa _{j}^{{\rm eff}}={{\kappa }_{j}}+\epsilon ({{s}_{j}}-{{d}_{j}})\gt {{\kappa }_{j}}.\end{eqnarray} \tag{ 18 }$$

Therefore, it always lowers the return times after perturbations, see equation (3).

A second, more realistic scenario is that the Decentral Smart Grid Control is much slower than the grid dynamics. Here, we consider a discrete time control system, where the price function is the same for all nodes and given by (6). Both supply and demand are updated periodically with periodicity $T\gg {{\sum }_{j}}{{M}_{j}}/{{\sum }_{j}}{{\kappa }_{j}}$. Assuming that the grid is dynamically stable for the given parameters, the angular frequency deviation will relax to

$$\begin{eqnarray}&&\langle \omega \rangle =\frac{\Delta P}{{{\sum }_{j}}{{\kappa }_{j}}}\quad {\rm with}\quad \Delta P=\mathop{\sum }\limits_{j}{{S}_{j}}(p)-{{D}_{j}}(p).\end{eqnarray} \tag{ 19 }$$

On this basis, a new market price ${{p}^{\prime }}$ is calculated. Assuming an affine-linear price function as in equation (14) and using the linearized supply and demand curves (12), we find

$$\begin{eqnarray}&&{{p}^{\prime }}-{{p}_{\Omega }}=-\epsilon \;\frac{{{\sum }_{j}}{{s}_{j}}-{{\sum }_{j}}{{d}_{j}}}{{{\sum }_{j}}{{\kappa }_{j}}}\;(p-{{p}_{\Omega }}).\end{eqnarray} \tag{ 20 }$$

This yields an oscillating dynamics of the market price, which is stable if and only if

$$\begin{eqnarray}&&\epsilon \lt {{\epsilon }_{{\rm cr}}}=\frac{{{\sum }_{j}}{{\kappa }_{j}}}{{{\sum }_{j}}{{s}_{j}}-{{\sum }_{j}}{{d}_{j}}},\end{eqnarray} \tag{ 21 }$$

setting a strict upper limit for the slope of the price function. The potential instability for $\epsilon \gt {{\epsilon }_{{\rm cr}}}$ is caused by an overreaction of the suppliers and consumers to market incentives. Similar rebound effects can generally occur in DR systems [8], such that this problem is not specific to the current proposal.

An example for the possible dynamics of the techno-economical system is shown in figure 2 for a model grid with three generators and three consumers with values inspired by the IEEE 9-bus test grid (see panel (e)). We assume that the price elasticities

$$\begin{eqnarray}&&{{E}_{S}}={{\left. \frac{{\rm d}{{S}_{j}}/{\rm d}p}{{{S}_{j}}/p} \right|}_{{{p}_{\Omega }}}}\quad {\rm and}\quad {{E}_{D}}={{\left. \frac{{\rm d}{{D}_{j}}/{\rm d}p}{{{D}_{j}}/p} \right|}_{{{p}_{\Omega }}}}\end{eqnarray} \tag{ 22 }$$

are the same for all nodes and given by ${{E}_{S}}\;=\;+0.3$ and ${{E}_{D}}=-0.3$ in agreement with empirical results [38–40]. The price and subsequently the demand and supply are adapted after a time step of T = 60 s. We assume an equilibrium price ${{p}_{\Omega }}=24$ ct kWh⁻¹ and a damping constant ${{\kappa }_{j}}=0.2/s\times {{M}_{j}}$ with ${{M}_{j}}={{10}^{4}}\;{\rm kg}{{{\rm m}}^{2}}\times \Omega$. All transmission lines have the same transmission capacity K = 200 MW. For these parameters, the system is stable if and only if $\epsilon \lt {{\epsilon }_{{\rm cr}}}\approx 3\;({\rm ct}/{\rm kWh})/(2\pi \;{\rm Hz})$. Otherwise the prices and the grid frequency will diverge after a small perturbation as shown in figure 2.

Zoom In Zoom Out Reset image size

We note that this result can be interpreted as an application of the famous cobweb theorem from microeconomics [41] to Decentral Smart Grid Control. The generalized equilibrium condition (11) can be seen as the intersection of the loss curve ${{\sum }_{j}}\kappa \langle \omega \rangle$ and the net supply curve $\Delta P(\langle \omega \rangle )={{\sum }_{j}}{{S}_{j}}(p(\langle \omega \rangle ))-{{D}_{j}}(p(\langle \omega \rangle ))$. The loss takes the role of an effective demand function with fast relaxation, while the net supply is adapted much slower. The cobweb model then states that the economic system will relax to an equilibrium, if the slope of the loss curve is larger than the slope of the net supply curve, which yields the stability condition (21).

New risks may emerge when the Decentral Smart Grid Control acts on similar time scales as the dynamics of the grid such that the two system become truly interdependent. We first consider the case of a delayed feedback, i.e., we consider the price function (7) with $\tau \geqslant 0$ but without averaging (T = 0), i.e., consumers measure their local frequency and try to adapt as fast as possible but need their intrinsic time τ to react. To obtain analytic solutions for the interdependent techno-economical system, we use a very simple system: we linearize the supply and demand curves (12) and consider only two nodes with equal technical and economical parameters, i.e., ${{M}_{1}}={{M}_{2}}$, ${{\kappa }_{1}}={{\kappa }_{2}}$, ${{s}_{1}}={{s}_{2}}$ and ${{d}_{1}}={{d}_{2}}$. Defining our new variables as the phase difference $\Delta \theta ={{\theta }_{1}}-{{\theta }_{2}}$ and the angular frequency difference of the two nodes $\Delta \omega ={{\omega }_{1}}-{{\omega }_{2}}$, the equations of motion read

$$\begin{eqnarray}\begin{array}{rcl} \frac{{\rm d}}{{\rm d}t}\Delta \theta & = & \Delta \omega \\ \frac{{\rm d}}{{\rm d}t}\Delta \omega & = & 2P-\alpha \Delta \omega -2K\cdot {\rm sin} \left( \Delta \theta \right)-\gamma \Delta {{\omega }_{\tau }}, \\ \end{array}\end{eqnarray} \tag{ 23 }$$

where we introduced the abbreviations $\alpha =\kappa /M$, $K={{K}_{12}}/M$, $\gamma =\epsilon ({{s}_{1}}-{{d}_{1}})/M$, $2P=({{S}_{1}}({{p}_{\Omega }})-{{D}_{1}}({{p}_{\Omega }})-{{S}_{2}}({{p}_{\Omega }})+{{D}_{2}}({{p}_{\Omega }}))/M$ and $\Delta {{\omega }_{\tau }}=\Delta \omega (t-\tau )$. In the following, time is measured in seconds, α and γ in s⁻¹ and P and K in s⁻². This notation is similar to the one used in [19]. Delayed differential equations need a history function as initial condition, which we chose to be $\Delta \omega (t\lt 0)=\Delta {{\omega }_{0}}\cdot (1+0.1{\rm tanh} (t/2))$ for our dynamical simulations. Furthermore, we used standard mathematica routines [42] to solve the ordinary and delayed differential equations.

The delayed adaption of supply and demand can both stabilize and destabilize the grid dynamics as shown in figure 3. The physical reason of this effect can be easily understood. The frequency-adaptive 'effective' power ${{P}_{{\rm eff}}}\left( t \right)=2P-\gamma \Delta {{\omega }_{\tau }}$ in equation (23) can be seen as a resonant driving acting on an oscillating system. Such a driving term will either damp or amplify the oscillations depending on whether the driving is in-phase or out-of-phase. The phase shift of this driving term is directly given by the delay τ, such that the stability of the system depends crucially on the value of τ. To illustrate this result, we compare the dynamics for two different values of τ in figure 3. For ${{\tau }_{1}}=0.75\;{\rm s}$ the driving is in-phase, the oscillations are amplified and the grid becomes unstable with potentially fatal results, whereas ${{\tau }_{2}}=1.5\;{\rm s}$ stabilizes the system.

Zoom In Zoom Out Reset image size

To obtain a global view of the stability and the role of the system parameters we analyze the dynamical stability around the steady-state operation of the grid given by the fixed point

$$\begin{eqnarray}&&\left( \Delta {{\theta }^{*}},\Delta {{\omega }^{*}} \right)=\left( {\rm arcsin} \left( \frac{P}{K} \right),0 \right).\end{eqnarray} \tag{ 24 }$$

A solution exists only if the transmission capacity is larger than the power which has to be transmitted, $K\gt P$ [32]. The linear stability of a dynamical system is determined by the eigenvalues of the Jacobian matrix. For a delayed system [43–45], we have to calculate the Jacobian of both the non-delayed system, based on equation (23),

$$\begin{eqnarray}{{J}_{0}}=\left( \begin{array}{ccccccccccccccc} \frac{\partial }{\partial \Delta \theta }\left( \frac{{\rm d}}{{\rm d}t}\Delta \theta \right) & \frac{\partial }{\partial \Delta \omega }\left( \frac{{\rm d}}{{\rm d}t}\Delta \theta \right) \\ \frac{\partial }{\partial \Delta \theta }\left( \frac{{\rm d}}{{\rm d}t}\Delta \omega \right) & \frac{\partial }{\partial \Delta \omega }\left( \frac{{\rm d}}{{\rm d}t}\Delta \omega \right) \\ \end{array} \right)=\left( \begin{array}{ccccccccccccccc} 0 & 1 \\ -2\;K{\rm cos} \left( \Delta {{\theta }^{*}} \right) & -\alpha \\ \end{array} \right),\end{eqnarray} \tag{ 25 }$$

and the derivatives for the delayed term

$$\begin{eqnarray}{{J}_{\tau }}=\left( \begin{array}{ccccccccccccccc} \frac{\partial }{\partial \Delta {{\theta }_{\tau }}}\left( \frac{{\rm d}}{{\rm d}t}\Delta \theta \right) & \frac{\partial }{\partial \Delta {{\omega }_{\tau }}}\left( \frac{{\rm d}}{{\rm d}t}\Delta \theta \right) \\ \frac{\partial }{\partial \Delta {{\theta }_{\tau }}}\left( \frac{{\rm d}}{{\rm d}t}\Delta \omega \right) & \frac{\partial }{\partial \Delta {{\omega }_{\tau }}}\left( \frac{{\rm d}}{{\rm d}t}\Delta \omega \right) \\ \end{array} \right)=\left( \begin{array}{ccccccccccccccc} 0 & 0 \\ 0 & -\gamma \\ \end{array} \right),\end{eqnarray} \tag{ 26 }$$

where we consider exponential solutions, see [45]. The stability eigenvalues λ are then determined by the solution of the characteristic equation

$$\begin{eqnarray}&&{\rm det} ({{J}_{0}}+{{{\rm e}}^{-\lambda \tau }}{{J}_{\tau }}-\lambda {\mathbb{1}})={{\lambda }^{2}}+\alpha \lambda +2\;K{\rm cos} (\Delta {{\theta }^{*}})+\lambda \gamma {{{\rm e}}^{-\lambda \tau }}=0.\end{eqnarray} \tag{ 27 }$$

Small perturbations induce an oscillatory motion with eigenfrequency ${\rm Im}(\lambda )$ and period $2\pi /{\rm Im}(\lambda )$. The amplitude grows or decreases exponentially as ${{{\rm e}}^{\operatorname{Re}(\lambda )t}}$ such that $\operatorname{Re}(\lambda )\gt 0$ is the condition for a dynamic instability. This is possible only if the frequency adaption is strong enough compared to the damping of the system,

$$\begin{eqnarray}&&\operatorname{Re}(\lambda )\geqslant 0\quad {\rm is}\;{\rm possible}\;{\rm only}\;{\rm if}\quad \gamma \geqslant \alpha ,\end{eqnarray} \tag{ 28 }$$

as shown below. When the delay τ changes, we observe a periodic pattern of stable and unstable parameter values as shown in figures 4 and 5. As explained above, destabilization occurs in the case of an in-phase driving which happens when τ is an integer multiple of the period of the eigenoscillations of the system.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To proof these statements we first note that $\operatorname{Re}(\lambda )\lt 0$ for $\gamma =0$ or $\tau =0$ as long as $\alpha \gt 0$. We now consider the parameter values where $\operatorname{Re}(\lambda )$ changes its sign such that the system becomes unstable. Decomposing the characteristic equation (27) into real and imaginary parts and setting $\operatorname{Re}(\lambda )=0$ yields the conditions

$$\begin{eqnarray}&&-{\rm Im}{{(\lambda )}^{2}}+2K{\rm cos} (\Delta {{\theta }^{*}})+\gamma {\rm Im}(\lambda ){\rm sin} \left( \tau {\rm Im}(\lambda ) \right)=0,\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&\alpha {\rm Im}(\lambda )+\gamma {\rm Im}(\lambda ){\rm cos} \left( \tau {\rm Im}(\lambda ) \right)=0.\end{eqnarray} \tag{ 30 }$$

The second equation can be solved for τ with the result

$$\begin{eqnarray}&&\tau ={\rm arccos} (-\alpha /\gamma )+\frac{2\pi }{{\rm Im}(\lambda )}n,\quad n\in \mathbb{N}.\end{eqnarray} \tag{ 31 }$$

One directly sees the periodicity in the delay τ, where the period $2\pi /{\rm Im}(\lambda )$ is equal to the period of the eigenoscillations of the system. Furthermore, the ${\rm arccos} (-\alpha /\gamma )$ is real only if $\gamma \geqslant \alpha$, which yields a necessary condition for the destabilization by delay. Note that this statement is equivalent to the one from section 4.2, namely the damping of the system has to be larger than the price influence to always guarantee stability.

In addition to the local stability analysis, we also consider the grid dynamics after a large perturbation. We integrate the equations of motion for a long time period until ${{t}_{{\rm max} }}=600\;{\rm s}$ with initial conditions randomly drawn from a subset of the phase space $Q=[-\pi ,\pi ]\times [-30,30]\;{\rm Hz}$, see also [46], and evaluate whether the system relaxes to a stable operation, i.e., whether it converges to the fixed point (24), or not. As a criterion for stable operation we assume an angular frequency deviation of $\Delta \omega \left( t={{t}_{{\rm max} }} \right)\leqslant 0.1\;{\rm Hz}$. The results are visualized in figure 5, panels (a)–(c), where we plot the stability in a color code (light green: stable, dark blue: unstable) as a function of the initial location in phase space for different values of τ. The results confirm the finding of the linear stability analysis. Depending on the value of τ the global stability is altered dramatically. For $\tau =0.8\;{\rm s}$ the fixed point is linearly unstable and the system does not relax for almost all initial conditions. On the contrary, a delay of $\tau =5\;{\rm s}$ leads to an almost perfect stability; the grid converges for almost all initial states. Notably, the regions of initial conditions leading to stable or unstable behavior are not clearly separated because the actual boundary of these regions has a rather complex geometric structure already in the non-delayed case.

The global stability of a fixed point of a dynamical system can be quantified by the volume of its basin of attraction. The 'basin size' ${{V}_{{\rm Basin}}}$ can be determined numerically using a Monte Carlo method as the ratio of initial conditions converging to a stable operation to the total number of initial conditions. Figure 5 panel (d) shows how the basin volume depends on the delay time τ in comparison to the linear stability exponent $\operatorname{Re}(\lambda )$. As expected, ${{V}_{{\rm Basin}}}$ tends to zero in the case of linear instability, $\operatorname{Re}(\lambda )\geqslant 0$. Maxima of ${{V}_{{\rm Basin}}}$ of different height are observed in the stable parameter regions, including an almost perfect stabilization for $\tau \approx 5\;{\rm s}$. However, these maxima do not coincide with the minima of $\operatorname{Re}(\lambda ).$ As both the linear stability and the basin volume predict similar delays to be problematic for the system, we focus on the computationally easier to handle linear stability in the following.

Measuring the local grid frequency will generally take some time in a real-world system. We thus consider the dynamics of the elementary model grid for the price function (7) including both delay $\tau \geqslant 0$ and averaging over a period $T\geqslant 0$. The equation of motion for the angular frequency difference $\Delta \omega ={\rm d}\Delta \theta /{\rm d}t$ then reads

$$\begin{eqnarray}&&\frac{{\rm d}}{{\rm d}t}\Delta \omega =2P-\alpha \Delta \omega -2K{\rm sin} \left( \Delta \theta \right)-\frac{\gamma }{T}\int _{t-T}^{t}\Delta \omega (t^{\prime} -\tau ){\rm d}t^{\prime} .\end{eqnarray} \tag{ 32 }$$

The integration can be carried out in a straightforward way by using (23) such that we obtain the modified delayed dynamical system

$$\begin{eqnarray}&&\frac{{\rm d}}{{\rm d}t}\Delta \omega =2P-\alpha \Delta \omega -2K{\rm sin} \left( \Delta \theta \right)-\frac{\gamma }{T}\left[ \Delta \theta \left( t-\tau \right)-\Delta \theta \left( t-T-\tau \right) \right].\end{eqnarray} \tag{ 33 }$$

To evaluate the stability of the steady state (24), we have to calculate the eigenvalues λ for a system with both the delay τ and the delay $\tilde{T}=T+\tau$. The Jacobian of the non-delayed system is given by (25) as above and for the delay terms only $\frac{\partial }{\partial \Delta {{\theta }_{\tau }}}\left( \frac{{\rm d}}{{\rm d}t}\Delta \omega \right)=-\frac{\gamma }{T}$ and $\frac{\partial }{\partial \Delta {{\theta }_{\tilde{T}}}}\left( \frac{{\rm d}}{{\rm d}t}\Delta \omega \right)=\;+\frac{\gamma }{T}$ are non-zero. The stability eigenvalues λ are then given by the roots of the characteristic equation [43]

$$\begin{eqnarray}&&p\left( \lambda \right)=-{{\lambda }^{2}}+\alpha \lambda +2K\sqrt{1-\frac{{{P}^{2}}}{{{K}^{2}}}}+\frac{\gamma }{T}\left( {{{\rm e}}^{-\lambda \tau }}-{{{\rm e}}^{-\lambda \left( T+\tau \right)}} \right)=0.\end{eqnarray} \tag{ 34 }$$

Figure 6 shows the real part of the stability eigenvalue $\operatorname{Re}(\lambda )$ as a function of the delay τ and the averaging time T. Instabilities are observed for certain values of τ, if T is small as discussed above. But for a sufficiently high T, the system stays stable regardless of the time τ. The actual values of the stability exponent $\operatorname{Re}(\lambda )$ for large $T\gg \gamma$ are comparable to the one of the system without any price adaptation.

Zoom In Zoom Out Reset image size

The results shown here are very interesting: while a delay in adaptation poses a stabilization risk to the grid, averaging the measured signal for a certain time removes the short time fluctuations that could resonantly drive the system and thereby guarantee a stable operation state. Still, the nodes can adapt to changes of the generation on all time scales slower then T, which provides an effective DR management system.

The larger the grid becomes, the more complex behavior it is able to show. Here, we consider a network with four nodes to test our results. Two consumers ($P=(S-D)/M\lt 0$) are supplied by two generators ($P=(S-D)/M\gt 0$). Each generator is coupled to a consumer to balance the power production/consumption. In addition, the generators are coupled to each other. As above we assume that all technical and economic parameters are the same for all nodes of the network. The equations of motion can be read in appendix B.

Risks from resonant driving emerge in a grid with delayed response and T = 0 as discussed above. The grid becomes unstable in certain regions of parameter space which are periodically spaced on the τ-axis. In a complex network, there are generally many different eigenmodes with different frequencies. The parameter regions where these modes can be excited generally overlap, such that the grid becomes unstable for most values of the delay τ as shown in the dark blue curves of panels (c), (d) in figure 7. Only for a very fast response $\tau \to 0$ the dynamic is stabilized.

Zoom In Zoom Out Reset image size

We conclude that either a delayed adaption must be avoided, or alternatively, an averaging over a sufficiently long period T must be introduced as well. Figure 7 panel (b) shows the stability exponents $\operatorname{Re}(\lambda )$ as a function of τ and T for the four-node network. Increasing T to values well above τ restores stability also for 'dangerous' values of τ as shown in the light green curves of panels (c), (d) in the figure.