Predictability and Fairness in Load Aggregation and Operations of Virtual Power Plants

In power systems, one wishes to regulate the aggregate demand of an ensemble of distributed energy resources (DERs), such as controllable loads and battery energy storage systems. We suggest a notion of predictability and fairness, which suggests that the long-term averages of prices or incentives offered should be independent of the initial states of the operators of the DER, the aggregator, and the power grid. We show that this notion cannot be guaranteed with many traditional controllers used by the load aggregator, including the usual proportional-integral (PI) controller. We show that even considering the non-linearity of the alternating-current model, this notion of predictability and fairness can be guaranteed for incrementally input-to-state stable (iISS) controllers, under mild assumptions.


I. INTRODUCTION
In response to large-scale integration of intermittent sources of energy [1] in power systems, system operators seek novel approaches to balancing the load and generation. The aggregation of distributed energy resources (DERs) and regulation of their aggregate power output has emerged as a viable balancing mechanism [2, 3, 4, e.g.] for many system operators. 1 For the transmission system operator, the load aggregator provides a single point of contact, which is responsible for the management of risks involved in the coordination of an ensemble of DERs and loads, and the legal and accounting aspects associated with modern grid codes and contracts governing the provision of ancillary services to the grid.
Independent of the form of the ensemble, load aggregation introduces a novel feedback loop [13,14] into regulation J. Marecek and M. Roubalik are at the Czech Technical University, Prague, the Czech Republic.
R. Ghosh and R. Shorten were at University College Dublin, Dublin, Ireland.
R. Shorten is at Imperial College London, South Kensington, UK. F. Wirth is at University of Passau, Germany. Manuscript received October 8, 2021. 1 In Australia, for instance, it is estimated that there are 2.5 million rooftop solar installations with a total capacity of more than 10 GW and 73,000 home energy storage systems with a 1.1 GW capacity. problems in power systems. Recent regulation 2 mandates the use of load aggregation, and hence makes the closed-loop analysis of load aggregation very relevant. A recent pioneering study of the closed-loop aspects of load aggregation by Li et al. [14] leaves three issues open: how to go beyond a linearisation of the physics of the alternating-current (AC) model? How to model the uncertainty inherent in response to the control signal provided by the load aggregator? Finally, what criteria of a good-enough behaviour of such a closed-loop system to consider? The first issue is particularly challenging: While the consideration of losses in the AC model may be necessary [15] for the sake of the accuracy of the model, it introduces a wide range of challenging behaviours [16,17] in the nonlinear dynamics.
In this paper, we address all three issues in parallel. First, we suggest that good-enough behaviour involves "predictability and fairness" properties that are rarely considered in control theory and require new analytical tools. These include: (a) whether each DER and load, on average, receives similar treatment in terms of load reductions, disconnections, or power quality, (b) whether the prices or incentives depend on initial conditions, such as consumption at some point in the past, (c) whether the prices, incentives, or limitations of service are stable quantities, preferably not sensitive to noise entering the system, and predictable.
Building upon the recent work of [18], we show that these criteria can be phrased in terms of properties of a specific stochastic model. In particular, these criteria are achieved whenever one can guarantee the existence of the unique invariant measure of a stochastic model of the closed-loop system. Thus, the design of feedback systems for deployment in the demand-response application should consider both the traditional notions of regulation and optimality and provide guarantees concerning the existence of the unique invariant measure.
Next, we show that many familiar control strategies, in straightforward situations, do not necessarily give rise to feedback systems that possess the unique stationary measure. In theory, we use the model of [18] to show that with any controller C that is linear marginally stable with a pole M M Load aggregator offers ancillary services to the transmission system operator.
A key complication is posed by the nonlinear filter, which models losses, latency in the actuation, etc.
Within the ensemble, DERs, fans of HVAC, pool pumps, etc., may choose to turn on or off in response to the control signal.
Reference signal provided by the transmission system operator (or independent system operator in the USA). s 1 = e qjπ on the unit circle, where q is a rational number, the closed-loop system cannot have a unique invariant measure. In practice, we illustrate this behaviour in simulations using Matpower [19].
As our main result, we show that one can design controllers that allow for predictability and fairness from the point of view of a participant, such as an operator of DERs, even in the presence of losses in the AC model. These results are based on a long work history on incremental input-to-state stability, which we have extended to a stochastic setting.

II. RELATED WORK
There is a rich history of work on dynamics and control in power systems, summarised in several excellent textbooks of [20,21,22,17], and [23]. Without aiming to present a comprehensive overview, one should like to point out that traditionally, load frequency control (LFC) is performed using controllers with poles on the unit circle, such as proportionalintegral (PI) controllers [24, 25, e.g.]. There are many alternative proposals, some of them excellent [26, 27, 28, 29, 30, 31, 32, 33, 34,  Although the idea of demand response has been studied for several decades [35, e.g.], only recently it has been shown to have commercial potential [36,37,38,39]. Notably, with a capacity as low as 20-70 kWh, load aggregators may break even [39] in the spot market in developed economies, although many grid codes still require load aggregators to aggregate a substantially more enormous amount of active power (e.g., 1 MW).
In a recent pioneering study, Li et al. [14] considered the closed-loop system design while quantifying the available flexibility from the aggregator to the system operator, albeit without the non-linearity of the alternating-current model.
While some of the more recent approaches utilize modelpredictive control [54, e.g.], many of the early proposals extend the use of PI. In contrast, we utilise iterated function systems [57,58,59], a class of Markov processes, to design novel controllers following [18].

III. NOTATION, PRELIMINARIES AND OUR MODEL
Let us introduce some definitions and notation: Definition 1. A function γ : R + → R + is is said to be of class K if it is continuous, increasing and γ(0) = 0. It is of class K ∞ if, in addition, it is proper, i.e., unbounded.
Definition 2. A continuous function β : R + × R + → R + is said to be of class K L , if for all fixed t the function β(·, t) is of class K and for all fixed s, the function β(s, ·) is is non-increasing and tends to zero as t → ∞.
In keeping with [18], our model is based on a feedback signal π(k) ∈ Π ⊆ R sent to N agents. Each agent i has a state x i ∈ R ni and an associated output y i ∈ D i , where the latter is a finite set. The non-deterministic response of agent i is within a finite set of actions: A i = {a 1 , . . . , a Li } ⊂ R ni . The finite set of possible resource demands of agent i is D i : We assume there are w i ∈ N state transition maps W ij : R ni → R ni , j = 1, . . . , w i for agent i and h i ∈ N output maps H i : R ni → D i , = 1, . . . , h i for each agent i. The evolution of the states and the corresponding demands then satisfy: where the choice of agent i's response at time k is governed by probability functions p ij : Π → [0, 1], j = 1, . . . , w i , respectively p i : Π → [0, 1], = 1, . . . , h i . Specifically, for each agent i, we have for all k ∈ N that Additionally, for all π ∈ Π, i = 1, . . . , N it holds that The final equality comes from the fact p ij , p i are probability functions. We assume that, conditioned on {x i (k)}, π(k), the random variables {x i (k +1) | i = 1, . . . , N } are stochastically independent. The outputs y i (k) each depend on x i (k) and the signal π(k) only.
Let Σ be a closed subset of R n with the usual Borel σalgebra B(Σ). We call the elements of B(Σ) events. A Markov chain on Σ is a sequence of Σ-valued random vectors {X(k)} k∈N with the Markov property, which is the equality of a probability of an event conditioned on past events and probability of the same event conditioned on the current state, i.e., we always have where G is an event and k ∈ N. We assume the Markov chain is time-homogeneous. The transition operator P of the Markov chain is defined for x ∈ Σ, G ∈ B (Σ) by If the initial condition X(0) is distributed according to an initial distribution λ, we denote by P λ the probability measure induced on the path space, i.e., space of sequences with values in Σ. Conditioned on an initial distribution λ, the random variable X(k) is distributed according to the probability measure λ k which is determined inductively by λ 0 = λ and for G ∈ B.
A probability measure µ is called invariant with respect to the Markov process {X(k)} if it is a fixed point for the iteration described by (6), i.e., if An invariant probability measure µ is called attractive, if for every probability measure ν the sequence {λ k } defined by (6) with initial condition ν converges to µ in distribution.
In an iterated function systems with place dependent probabilities [60,61,62], we are given a set of maps {f j : Σ → Σ | j ∈ J }, where J is a (finite or countably infinite) index set. Associated to these maps, there are probability functions The state X(k+1) at time k+1 is then given by f j (X(k)) with probability p j (X(k)), where X(k) is the state at time k. Sufficient conditions for the existence of a unique attractive invariant measure can be given in terms of "average contractivity" [57,58,59]. Any discrete-time Markov chain can be generated by an iterated function system with probabilities see Kifer [

IV. THE CLOSED-LOOP MODEL AS ITERATED RANDOM FUNCTIONS
Let us now model the system of Figure 1 by a Markov chain on a state space representing all the system components. In general, one could consider: F : C : In addition, we have (e.g., Dini) continuous probability functions p ij , p il : Π → [0, 1] so that the probabilistic laws (4) are satisfied. If we denote by X i , i = 1, . . . , N, X C and X F the state spaces of the agents, the controller and the filter, then the system evolves on the overall state space (13) where each of the maps F m is of the form (14) and the maps F m are indexed by indices m from the set By the independence assumption on the choice of the transition maps and output maps for the agents, for each multi-index m = ((1, j 1 ), . . . , (N, j N ), (1, l 1 ), . . . , (N, l N )) in this set, the probability of choosing the corresponding map F m is given by p ili (π(k)) =: q m (π(k)). (16) We have obtained an iterated function system of the previous subsection III.

A. Our Objective
We aim to achieve regulation, with probability 1 to some reference r > 0, given by the installed capacity. Further, we aim for predictability, in the sense that for each agent i there exists a constant r i such that where the limit is independent of initial conditions. This can be guaranteed, when one has a unique invariant measure for the closed-loop model as iterated random functions. The requirement of fairness is then for the limit r i to coincide for all 1 ≤ i ≤ N .

B. A Motivating Negative Result
The motivation for our work stems from the fact that controllers with integral action [76,77], such as the Proportional-Integral (PI) controller, fail to provide unique ergodicity, and hence the fairness properties. In many applications, controllers with integral action are widely adopted, including leading control systems [78] in power generation. A simple PI control can be implemented as: which means its transfer function from e to π, in terms of the Z transform, is given by Since this transfer function is not asymptotically stable, any associated realisation matrix will not be Schur. Note that this is the case for any controller with any sort of integral action, i.e., a pole at z = 1. Consider the feedback system in Figure 1, where F : y →ŷ is a finite-memory moving-average (FIR) filter. Assume the controller C is a linear marginally stable single-input singleoutput (SISO) system with a pole s 1 = e qjπ on the unit circle where q is a rational number, j is the imaginary unit, and π is Archimedes' constant 3.1416. In addition, let the probability functions p il : R → [0, 1] be continuous for all i = 1, . . . , N, l = 1, . . . , m i , i.e., if π(k) is the output of C at time k, then P(x i (k + 1) = a l ) = p il (π(k)). Then the following holds. (17) is discrete, then the closed-loop system cannot be uniquely ergodic. (iii) For rational A i ⊂ Q for all i = 1, . . . , N , r ∈ Q, and rational coefficients of the FIR filter F , the closed-loop system cannot be uniquely ergodic.
This suggests that it is perfectly possible for the closed loop both to perform its regulation function well and, at the same time, to destroy the ergodic properties of the closed loop. We demonstrate this in a realistic load-aggregation scenario in Section VI.

V. THE MAIN RESULT
For a linear controller and filter, [18] have established conditions that assure the existence of a unique invariant measure for the closed-loop. Let us consider non-linear controllers and filters, as required in power systems, and seek unique ergodicity conditions. Incremental stability is well-established concept to describe the asymptotic property of differences between any two solutions. One can utilise the concept of incremental inputto-state stability, which is defined as follows: Definition 3 (Incremental ISS, [79]). Let U denote the set of all input functions u : Z ≥k0 → R d . Suppose, F : R d × R n → R n is continuous, then the discrete-time non-linear dynamical system is called (globally) incrementally input-to-state-stable (incrementally ISS), if there exist β ∈ K L and γ ∈ K such that for any pair of inputs u 1 , u 1 ∈ U and any pair of initial condition ξ 1 , ξ 2 ∈ R n : Definition 4 (Positively Invariant Set). A set X ⊆ R n is called positively invariant under the system (21) if for any initial state x(k 0 ) = ξ ∈ X we have φ(k) ∈ X ∀k ≥ k 0 .
Definition 5 (Uniform Exponential Incremental Stability). The system (21) is uniformly exponentially incrementally stable in a positively invariant set X if there exists κ ≥ 1 and λ > 1 such that for any pair of initial states ξ 1 , ξ 2 ∈ X , for any pair of u 1 , u 2 ∈ U , and for all k ≥ k 0 the following holds If X = R n , then the system (21) is called uniformly globally exponentially incrementally stable.
With this notation, we can employ contraction arguments for iterated function systems [58] to prove: Theorem 2. Consider the feedback system depicted in Figure  1. Assume that each agent i ∈ {1, · · · , N } has a state governed by the non-linear iterated function system where: (i) we have globally Lipschitz-continuous and continuously differentiable functions W ij and H ij with global Lipschitz constant l ij , resp. l ij (ii) we have Dini continuous probability functions p ij , p il so that the probabilistic laws (4) are satisfied (iii) there are scalars δ, δ > 0 such that p ij (π) ≥ δ > 0, p ij (π) ≥ δ > 0 for all π ∈ Π and all (i, j) (iv) the composition F m (14) of agents' transition maps and the probability functions satisfy the contraction-onaverage condition. Specifically, consider the space X a := N i=1 X i of all agents and the maps F a,m : X a → X a , together with the probabilities p m : Π → [0, 1], m ∈ M, form an average contraction system in the sense that there exists a constant 0 <c < 1 such that for all x,x ∈ X a , x =x, π ∈ Π we have Then, for every incrementally input-to-state stable controller C and every incrementally input-to-state stable filter F compatible with the feedback structure, the feedback loop has a unique, attractive invariant measure. In particular, the system is uniquely ergodic.
Proof. To prove the result, we have to show that the assumptions of Theorem 2.1 in [58] are satisfied. To this end, we note that the maps F m satisfy, by (i), the necessary Lipschitz conditions. Further, the probability maps have the required regularity and positivity properties by (ii) and (iii). Thus, it only remains to show that the maps F m defined in (14) have the average contraction property. The existence of a unique attractive invariant measure then follows from [58] and unique ergodicity follows from [57].
First we note that the filter and the controller are incrementally stable systems that are in cascade. It is shown in [79,Theorem 4.7] that cascades of incrementally ISS systems are incrementally ISS. The proof is given for continuous-time systems but readily extends to the discrete-time case. It will therefore be sufficient to treat filter and controller together as one system with state z.
To show the desired contractivity property, consider two distinct points x = (x a , z) = ((x i ), z),x = (x a ,ẑ) = ((x i ),ẑ) ∈ X. Let L be an upper bound for the Lipschitz constants of the output functions H ij , H f . The assumptions on incremental ISS ensure that for each index m As the maps F a,m satisfy the average contractivity condition, it is an easy exercise to see that also the sets of iterates satisfy the contractivity condition. The result then can be obtained using Theorem 15 of [80], where it is shown that incrementally stable systems are contractions, using the converse Lyapunov theorem of [81]. The proof is concluded by invoking Theorem 2.1 in [58].
Subsequently, one may wish to design a controller that is a priori incrementally input-to-state-stable, e.g. considering extending the work on ISS-stable model-predictive controller [82] or model-predictive controllers [

VI. NUMERICAL RESULTS
Our numerical results are based on simulations utilising Matpower [19], a standard open-source toolkit for power-systems analysis in version 7.1 running on Mathworks Matlab 2019b. We use Matpower's power-flow routine to implement a nonlinear filter, which models the losses in the alternating-current model.
The model consists of an ensemble of N DERs or partiallycontrollable loads divided into two groups that differ in their probabilistic response to the control signal provided by the load aggregator. In particular, following the closed-loop of Figure 1: The response of each system S i to the control signal π(k) provided by the load aggregator at time k takes the form of probability functions, which suggest the probability a DERs is committed at time k as a function of the control signal π(k). In our simulations, the two functions satisfying (8) are: where we have implemented both PI and its lag approximant for the control of an ensemble of DERs or partiallycontrollable loads. The binary-valued vector: captures the output of the probability functions (26) in a particular realization.
The commitment of the DERs within the ensemble is provided to Matpower, which using a standard power-flow algorithm computes the active power output P (k) of the individual DERs: which is aggregated into the total active power output p(k) = P (k) of the ensemble. Then, a filter F is applied: Notice that here we both smooth the total active power output with a moving-average filter and accommodate the losses in the AC model. The error between the reference power output r and the filtered valuê p(k) is then used as the input for the controller. Output of the controller π(k) is a function of error e(k) and an inner state of the controller x c (k). Signal of the controller is given by the PI or its lag approximant: The procedure is demonstrated in the following pseudocode.

A. A Serial Test Case
First, let us consider a simple serial test case depicted in Figure  2, where there are a number of buses connected in series, as they often are in distribution systems. In particular, the buses are: • a slack bus, such as a distribution substation, • 5 buses with 12 DERs of 5 MW capacity connected to each bus, with probability function g i2 of (26) modelling the response of these first 60 DERs to the signal, • 5 buses with 12 DERs of 5 MW capacity connected to each bus, with probability function g i1 of (26) modelling the response of these 60 DERs to the signal, • a single load bus with a demand of 500 MW.
We aim to regulate the system to the reference signal r, which is 300 MW plus losses at time k = 0. Initially, the first 60 (= 5·12) generators are off and the second following 60 (= 5·12) generators are on. Considering the ensemble is composed of 120 DERs, with a total capacity of 600 MW, this should yield a load of less than 200 MW on the slack bus.
We repeated the simulations 300 times, considering the time horizon of 2000 time steps each. The results in Figures 3-5 present the mean over the 300 runs with a solid line and the region of mean ± standard deviation as a shaded area. Notice that we are able to regulate the aggregate power output (cf. Figure 3). With the PI controller, however, the state of the controller, and consequently the signal and the power generated at individual buses over the time horizon, are determined by the initial state of the controller (cf. Figures  4-5).

B. A Transmission-System Test Case
Next, let us demonstrate the analogous results on the standard IEEE 118-bus test case with minor modifications. Notably, the ensemble is connected to two buses of the transmission system. At bus 10, 1 generator with a maximum active power output of 450 MW was replaced by 8 DERs with 110 MW of power each, out of which 4 DERs used probability function g i1 , where x 01 = 660. The other four DERs used function g i2 , where x 02 = 110 and ξ = 100. At bus 25, 1 generator with maximum active power output of 220 MW was replaced by four generators with active power output of 110 MW, out of which two generators used probability function g i1 , and the other two used function g i2 with r = 660.
As in the serial test case, we have repeated the simulations 500 times, considering the time horizon of 2000 time steps each. Figure 6 yields analogical results to Figure 4 in the serial test case, where a different initial state x c (0) of the controller yields a very different control signal produced by the PI controller. This is not the case for the lag approximant. Consequently, some DERs may be activated much more often than others, all else being equal, solely as an artifact of the use of PI control. Many further simulation results illustrating this behaviour are available in the Supplementary material on-line.

VII. CONCLUSIONS AND FURTHER WORK
We have introduced a notion of predictability of fairness in load aggregation and demand-response management, which relies on the concept of unique ergodicity [83], which in turn is based on the existence of a unique invariant measure for a stochastic model of a closed-loop model of the system. Related notions of unique ergodicity have recently been used in socialsensing [83], and two-sided markets [84]; we envision it may have many further applications in power systems.
The notions of predictability and fairness, which we have introduced, may be violated, if controls with poles on the unit circle (e.g., PI) are applied to load aggregation. In particular, the application of the PI controller destroys the ergodicity of the closed-loop model. Considering the prevalence of PI controllers within load frequency control, this is a serious concern. On the other hand, we have demonstrated that certain other controllers guarantee predictability and fairness.  An important direction for further study is to consider unique ergodicity of stability-constrained ACOPF [28,29,30,31,32,33,34], or related controllers based on stability-constrained model-predictive control (MPC). We believe there are a wealth of important results ahead.

APPENDIX
For the convenience of the reader, we attach an overview of the notation and additional numerical results in the following appendices.
A. An Overview of Notation value of y ( k) filtered by filter. p i (π) probability functions. p i (π) probability functions.
internal state of agent i at time k. W ij transition maps. Lipschitz constants. ν an invariant distribution. β a class K L function. γ a class K function. u(k) a binary vector indicating the commitment of DERs in (28). p(k) filtered value of the active power output of DERs in (30).
active power output of DERs in (29). p(k) aggregate active power output of the ensemble at time k. e(k) error defined in (31).

B. Complete Numerical Results
We provide numerical results for three test cases, two of which have been discussed in the main body of the text. First, however, we present results for one additional, intentionally simple test case.

1) A Simple Test Case with Transmission Losses in the Filter:
The first one is a small 3-bus system depicted in Figure 1. At bus 1, there are N = 10 DERs, five of which use probability function g i1 with x 01 = 200. The other five DERs at bus 1 use function g i2 , where x 01 = 40 and ξ = 100. Bus 2 is considered slack; therefore, the generator at bus 2 balances the power deviation of bus 1 and losses in order to fulfil power demand at bus 3. The simulations were run for both PI and its lag approximant for several different initial values of the controller's internal state x c (0). Figures 8-?? demonstrate that both controllers can regulate the ensemble around required power r = 200. However, the PI controller's behaviour is strongly dependent on the initial condition of its internal state, whereas the lag approximant provides predictable behaviour in the long run. These results are in good agreement with simulation results in [18]. Fig. 7. A simple test case. Fig. 8. Results of simulations on the system in Figure 7: Power on buses 1 and 2 as a function of time for the two controllers and two initial states xc of each of the two controllers. Notice that both PI and its lag approximant are able to regulate the aggregate power produced by the ensemble of loads. Fig. 9. Results of simulations on the system in Figure 7: Observations available to the controller, as a function of time, for the two controllers and two initial states xc of each of the two controllers. Notice that both PI and its lag approximant receive similar inputs.

C. A Serial Test Case with Losses in the Filter
In our second example, we consider a serial test-case, with a number of buses connected in series. (A line-graph topology, in the language of graph theory.) The buses are, in this order, • a slack bus, which could represent a distribution substation • 5 buses with 12 DERs of 5 MW capacity each, and the "type 2" probability function. • 5 buses with 12 DERs of 5 MW capacity each, and the "type 1" probability function.
• a single load bus with a 500 MW demand.
Overall, the ensemble is composed of 120 DERs, with a total capacity of 600 MW. Initially, the first 5 · 12 generators are off and the second following 5 · 12 generators are on. We aim to regulate the power output of the ensemble to 60 running DERs (totalling 300 MW), which should yield a load exceeding 200 MW on the slack bus (considering the losses).
In our first set of simulations on the serial test case, we repeat 500 times the simulations considering the time horizon of 2000 time steps each. r = 300, x 0 = 300. The filter is (y(k − 1) + y(k))/2 + losses(k), whereby the losses hence propagate into the signal π(k).
As in the simple test case, we are able to regulate the aggregate power output (cf. Figure 12). With PI controller, however, the state of the controller, and consequently the signal and the power at individual buses are determined by the initial state of the controller (cf. Figures 13-14).
In our second set of simulations on the serial test case, we repeat 500 times the simulations considering the time horizon of 2000 time steps each. Instead of r = 300, x 0 = 300 of the first set of simulations, we consider r to be 300 plus losses at time Fig. 11. A serial test case. k = 0, and similarly x 0 to be 300 plus losses at time k = 0. Otherwise, the system's behaviour is similar.
As above, we are able to regulate the aggregate power output (cf. Figure 15). With PI controller, however, the state of the controller, and consequently the signal and the power at individual buses are determined by the initial state of the controller (cf. Figures 16-17).

D. A Transmission-System Test Case with Losses in the Filter
The third example demonstrates the procedure's applicability to the standard IEEE 118-bus test case with minor modifications. Notably, the ensemble is connected to two buses of the transmission system. At bus 10, 1 generator with a maximum active power output of 450 MW was replaced by 8 DERs with 110 MW of power each, out of which 4 DERs used probability function g i1 , where x 01 = 660. The other four DERs used function g i2 , where x 02 = 110 and ξ = 100. At bus 25, 1 generator with maximum active power output of 220 MW was replaced by four generators with active power output of 110 MW, out of which two generators used probability function f i1 , and the other two used function f i2 with r = 660. Bus 69 is the slack bus, similarly to the generator at bus 3 in the previous example. Figures 18-20 yield analogical results to the first example, for both controllers and different initial conditions x c (0).