An Ising Spin-Based Model to Explore Efficient Flexibility in Distributed Power Systems

This paper analyses customers’ demand flexibility in a local power trading scenario through an Ising spin-based model. We look at quantitative information on the two-way relationships between power exchanges and spin dynamics. To this end, a modified version of the Metropolis-Hastings algorithm was implemented, including a gradient descent through the constraint surface. This allowed us to analyse the system on a large scale (considering the cumulated benefit of all the actors involved) and also from the perspective of total aggregation. In a maximum flexibility scenario, the total aggregation profit increases with the number of aggregators. We also investigate numerically the effect of aggregator boundaries on the spin dynamics.


Introduction
How can customers' collective behaviour affect the efficiency of distributed power systems?Furthermore, could the outcomes of this collective behaviour be exploited in the forthcoming self-organised power systems?In this work, we explore these matters analytically and numerically.An Ising spin-based model allowed us to provide quantitative criteria to couple system's performance and customers' standpoint on being flexible or not in their demand.
Distributed power systems (DPSs) are complex ecosystems encompassing machines, networks, procedures, operators, and customers organised in hierarchical layers [1].After restructuring the power systems, new players have appeared, such as Distribution System Operators (DSOs) which provide electrical demand to local customers [2].Also, customers are increasingly changing their roles from passive (only consumer) to active (generators).Moreover, recent Demand Response (DR) programs require customers to make a timely adjustment of their demand [3].In this interwined subsystems, customers are also exposed to social interactions that can influence their decision as in any other community.In particular, users' flexibility in power consumption has a major influence on the performance of the whole system [4,5].In turn, DPSs powerful restrictions to maintain quality of service and security cause new constraint forces applied to customer behaviour.
On the other hand, collective behaviour has features which are hard to explain by classical statistical methods [6].For example, the so-called herding behaviour or the economic bubbles are phenomena lying outside analyses that neglect large correlations and self-organisation occurring near the critical point [7].In regard to DPSs, a sharp transition in customer's behaviour (i.e., emergence) can trigger dramatic consequences.To the best of our knowledge however, a quantitative analysis of the complexity of customers behaviour in DPSs has not yet been provided.Ising-based models are a promising approach as we demonstrate here.
The Lenz-Ising spin model is known since 1925 [8].In short, the model consists of an arrangement of interacting agents which have two possible states (i.e., up and down).Agents interact with their neighbours locally and are also 2 Complexity exposed to external action and to thermal noise.In the physics metaphor of a ferromagnetic material, the atoms tend to align their quantum spins to minimise energy.However, for high temperatures, the agents flip their spins randomly.There is a critical temperature   (known as Curie temperature) for which the system suffers a sharp change-second-order phase transition-from order to disorder.Near   the system has anomalous behaviour and the correlations among spin states propagate fast to the entire system.As we will show in Section 2, this model is fairly simple but powerful enough to capture most of the features of complex systems (e.g., phase transitions and universality) [9].This has motivated many researchers to apply the Ising model to different fields ranging from ecology [10] to language evolution [11] (see [12] for a discussion of the physical motivation of the model and its limitations).
Given that social systems are both finite and heterogeneous, the Ising model has been exported in different flavours to capture the phenomena under investigation.Perhaps the most known of these Ising-like models is Tom Schelling's model of segregation [13] which has been shown to provide insight into the mechanisms under segregation in U.S. cities.This model has been shown to roughly correspond to an Ising model at  = 0 [14].In this regard, the temperature has been understood as a proxy of tolerance in binary thermodynamic societies; solubility corresponds to integration (i.e., mixing) and the miscibility gap to segregation [15].Ising inspired systems have also been applied to financial markets to explain expectation bubbles and economic crashes [16].For instance, the authors in [17] create a synthetic market where agents can take three actions: buy, sell, or stay inactive.Also, researchers in the power systems and electricity markets domain have made an effort to include part of this complex behaviour in the problem.For instance, in [18] bilateral electricity markets are analysed as complex networks.In [19], the authors use a game theoretic approach to understand the cooperation between small-scale electricity production and consumers.Multiagent systems (MAS) have also been applied to electricity markets to allow decentralized decisionmaking [20].See, for instance, [21], where authors modelled the behaviour of local consumers and producers as active agents in the electricity market based on the distributed control approach.More recently, the entity of aggregators has entered into the modelling scene.These are agents mediating between customers and distribution companies to offer demand bulks at competitive price.In this context, the authors in [22] define a bilevel problem where the upperlevel maximises the profit of the DSO, while the lower-level maximises the profit of each aggregator.Finally, in [23] the authors use a Hopfield neural network to optimise control in power systems and make the whole system self-controlled.Here, there is a super-system composed of several power subsystems which are in turn aggregations of customers.Each customer is simulated as a node in the neural network with two possible states: generation and consumption.Since the Hopfield network is formally equivalent to a spin-glass model, the problem has some formal resemblance with our work.However, our approach is very different.In [23], the aim is to adjust demand and supply to balance the system and let it work autonomously; agent interactions are determined to maintain proper frequency and voltage values only and there is no market.In our work, spins are agents in a decision environment with power exchanges among aggregators, the DSO, and the real-time electricity market (RTEM).Moreover, in this paper, we assume that DSOs act as agents able to participate in the real-time electricity markets to trade realtime power in a two-way fashion.Notice that here we follow the same approach as in [22], where DSO behaves as a proactive market agent able to transact power with the realtime electricity market.This way, besides the RTEM, the DSO is able to trade power with local agents (e.g., aggregators and nodal consumers).Therefore, consumers can play as virtual generation agents according to their behaviour in terms of flexibility.
For the first time, we couple the local power trading problem with an Ising-based model and analyse the interrelationships between them.The Ising model provides quantitative criteria regarding how the diffusion of flexible/nonflexible behaviour impacts the required constraints in the optimisation problem.
This paper is structured as follows.First we describe the local power exchange problem in Section 2 and how the Ising model can be linked to it in Section 3. Then in Section 4 we analyse the problem at the large scale-the social welfare-where we find how the power exchange problem constrains the Ising model in different ways.In Section 5 we decrease the scale and we focus on the problem from the demand aggregator's perspective.Finally we conclude and make our final remarks in Section 6.

Hierarchical Model and Agent Interactions in the Distributed Power Grid
In our setting, we consider a real-time interplay structure among five kinds of actors: consumers, bus-loads, aggregators, DSOs, and RTEM.In Figure 1, we describe schematically the power exchange problem and how it couples with the spin model.In this way, the power transaction between bus-loads and aggregators is two-way according to the flexibility behaviour from the demand side.Moreover, there are two-way power exchanges between aggregator and DSO.However, the power transaction between DSO and bus-loads is one-way: only from DSO to bus-loads.In other words, bus-loads can only buy real-time power from the DSO, while they are able to buy/sell power from/to aggregators.Also, only DSO can participate in the RTEM.
Each bus-load connects a set of customers whose state in terms of demand corresponds to a binary state of a spin (here represented as an arrow up/down mimicking the spin    metaphor).From the demand side, different types of power balancing or other constraints can force cooperation among customers either at interbus level (blue arrows) or at intrabus level (pink arrow) in Figure 1(b).In this work, we use a 33bus reference system (Figure 1(c)) from [1].This is composed of 32 loads and an entry bus-indicated as "s/s"-playing the role of a slack bus connected to the main grid; power exchanges between the DSO and RTEM are done through this slack bus.The whole demand in the distributed power system is partitioned by aggregators.In our schematic we show the case of three operators as in [22].In this work, we will use the data in [1] (shown in Table 1) as a reference for the maximum scheduled loads at each bus.From this data and assuming that a home might consume around 13 kW [24], we estimate the approximate number of customers per bus (assuming one customer per home) which makes a total of  ≈ 289 customers.Finally in Figure 1(d), we show a 2D grid arrangement where each of the 17 × 17 = 289 cells represents a customer.Furthermore, we assume that in the bus-loads there exists some notion of topological closeness (e.g., spatial or electrical proximity) which is partially mapped onto the square lattice.

Bus
This means that each customer is linked at least to one bus which is in turn connected to at least one aggregator.Below, we describe each agent in the system, specifying both its constraints and its objective function.
where we have introduced consistently    = ∑ ∈     and    = ∑ ∈     .For the scheduled load at each bus , we use the expression: where    represents the power at the source bus and  max  is the normalised expected load at bus  (see Table 1).The perbus flexible component splits itself into the power exchanged with both DSO ( DSO2L

𝑗𝑡
) and aggregator ( L2A  ): According to the schematic shown in Figure 1(a) loads can only buy from the DSO and hence we have On the other hand, the bus-aggregator power exchanges are bidirectional.If  L2A  < 0 demand at bus  is buying from aggregator, the flexibility decreases.However, if  L2A  > 0 then the flexibility of the bus is increased.This can be interpreted as virtual generation injected into the aggregator decreasing the scheduled load at time .
Eventually, the amount of flexibility would be constrained in different ways to ensure self-sustainability and dynamic flexibility of all loads.In this work, we allow two types of flexibility constraint: On one hand, (7) is the definition of the shiftable-loads (i.e., loads that can be shifted over time).On the other hand, (6) increases the self-sustainability of the distributed power system and converges the problem to the optimum social welfare.As we will see both constraints lead to different scenarios.
The optimisation at each bus  can be expressed as the trade-off between load bought from DSO at price  DSO2L and virtual generation sold to its aggregator  at price  L2A  integrated over time: This function must be minimised from the perspective of the loads' profit.However, both aggregators and DSO have their own priorities as we show below.

Aggregators.
The first thing to notice is that each aggregator is able to sell to the DSO all the virtual generation collected among the set of buses he operates on: Also there is a price model for these exchanges which is constrained in the following way: where   ≥ 1 represents a lower bound threshold for the aggregator-to-DSO price compared to load-to-aggregator price.This ensures aggregators' profit, which makes it reasonable for them to be part of the market.Also, the aggregatorto-DSO price is limited by the DSO-to-market price; this upper bound acts as a price control, limiting the bidding price of aggregators below the real-time price.Here we will use reference values from [22] for   = 1.1,  L2A  and  RT  related to the NordPool market.We summarise the respective values in Tables 2 and 3.
The optimisation function for each aggregator  can be expressed as the accumulated balance over time between power  A2DSO  sold to DSO at price  A2DSO  and power  L2A  bought from demand , at price  L2  : which as in the case of demand must be minimised to reach a profitable situation from the aggregators' perspective.

Distribution System
Operator.In our model, DSO is an agent able to exchange power directly with all other agents in the system.As stressed, the DSO is the only player who can exchange power with the RTEM.In the powerflow balance, the power sold by the DSO to the whole demand-bus-loads-arrives from the power transacted with both aggregators and market.We express this in the following equation: Therefore, the DSOs optimisation function can be expressed as the accumulated trade-off among these quantities, times the respective prices over time: This function will be minimised from the DSO's perspective to maximise its profit.So far we have described the main actors in the power exchange scenario with power balance and price constraints along with the respective optimisation functions for each agent.Now we describe the Ising model and our interpretation of its constituents in this context.

An Ising Spin Model for Customers' Flexibility
The Ising Hamiltonian ( 14) describes the interaction among entities (i.e., agents or spins) given their state   ∈ {−1, 1} and between each spin and a global magnetic field .The coupling constant  measures the strength of spin-to-spin interactions.The notation ⟨, ⟩ refers to pairs of spins belonging to the same radius of action or neighbourhood.When  > 1 (ferromagnetism) spins tend to align in the same direction and if  < 1 (antiferromagnetism) the spins tend to align in opposite directions.For  = 0, there is no spin-to-spin interaction.The external action of a positive field  > 0 will also foster positive spin alignments (and the other way around for  < 0).For a given temperature  the probability for finding a spin configuration Γ = {  } is proportional to the Boltzmann factor: exp(−/  ), where   stands for the Boltzmann constant.It is usual to take units such that  = 1 and   = 1.An important magnitude is the magnetisation  = (1/) ∑  =1   which measures the macroscopic effect of the spin states.In the so-called Mean Field approximation its value is One way to implement the Ising model numerically is through the Metropolis-Hastings algorithm [25], which belongs to the family of the Markov Chain Monte Carlo (MCMC) methods.Applied to our case it can be understood as a random walk over the configuration space Ω = {Γ} that converges to the Boltzmann distribution.In Algorithm 1 we show the pseudocode of a slight variant of the classical algorithm where we have included the possibility for implementing a constraint  at each step.Here  can represent the constraints in ( 6) and (7).The idea is that the spin shift is performed also when the constraint is minimised in absolute value.The number of iterations is chosen so that the final configuration reaches equilibrium.In operative terms this means that the correlations among the spin states are negligible at this point.
The size of the lattice is another important factor when using the model.On one hand, a small lattice with free boundary conditions-the one used here-shows border effects which are not present in large systems or in systems with other boundary conditions (e.g., periodic).In particular, this can affect the number of iterations to decorrelate the system and the potential transitions of state.There is an interesting debate on whether finite systems can undergo phase transitions or not [26].Clearly, a finite system cannot reproduce a singularity in a purely mathematical sense using a finite number of sums.As it is shown in [26] phase transitions in finite systems tend to be rounded and smoothed near the critical points.In Figure 2 we show the autocorrelation function for the spin states between the initial and evolved configuration for increasing number of iterations.We start with a random initial configuration for the 17 × 17 spin system shown in Figure 1(d).For > 5 iterations it is possible to achieve a reasonable value of decorrelation.In the inset, we show the phase transition in  beyond   .Notice how the numerical implementation is still able to reproduce qualitatively the transition  = ±1 → 0 although in a smoothed and rounded way.
When using this model one has to fix the interpretation for the following elements: The spin-to-spin interaction can be a proxy for communication among customers and the external field  can simulate a top-down directive forcing customers to follow some policies or Demand Response (DR) programs that are adopted by players in the top layers of the system (e.g., DSO and aggregators).On the other hand, in the context of financial markets [27] interprets  as a measure for the deviation from the fundamental value; if  is large agents are afraid of trading unless they have strong support from their private opinions or from their neighbours.In our case  = 1 and  = −1 represent that the system has either maximum or minimum flexibility.A value of  close to 0 is interpreted as a perfect balance where half of the customers are flexible and the other half is not.Finally, the temperature can represent customers' uncertainty about the effectiveness of the flexibility program that can impact on the reliability and security of the system.We can also associate  with noise in customers' information channels; too much noise destroys the spontaneous magnetisation (max/min flexibility previously achieved) [15].Also, when the "social temperature" [28] is high it destroys large domains (subpopulations) and hence it favours the mixing of options.
Let us also assume that all consumers connected to a bus have the same amount of scheduled load:    =    > 0, ∀ ∈    , ,  (notice that this does not reduce the generality of the problem; we can always adjust the number of customers per bus to represent the scheduled per-bus load    ).Therefore, the flexible component of the load at bus  at time  renders where Notice that the power flow from aggregator  to DSO in ( 9) can be expressed as Also, from (1), (12), and (18) one can express the DSO power constraints as Since    ≥ 0, ∀ ∈ ,  and  ≥ 0, ∀ ∈ , inequality 5 is expressed as Notice that constraints ( 17) and ( 19) can be combined into Finally, flexibility constraints in ( 6), ( 7) can be expressed through (16) as By combining constraints (10), ( 17), (19), (20), with either (22) or (23) and different objective functions, we can build scenarios from different perspectives and scales as we show below.

Large Scale Optimisation: Flexibility and Social Welfare
With the definitions above and from (8), we find the optimisation for each bus : (where we have used OF  = ∑ ∈ Δ  OF ∈  ).The optimisation functions for each aggregator from ( 9) and (11) render where we have used Finally, from ( 9) and ( 13), we find By defining OF  ≡ ∑ ∈ OF  and OF  ≡ ∑ ∈ OF  , the total load and total aggregation objective functions are found straightforward.Finally our objective is to maximise the social welfare of the distributed power system.Hence, the optimisation function is defined as OF  = OF  + OF  + OF DSO to ensure the maximised profits of all players (DSO, aggregator, and bus-loads) in the distributed power system.The global optimisation function OF  can be expressed as Notice how from this global scale the specifics about power exchanges between DSO-to-demand, demands-toaggregators, and aggregators-to-DSO cancel out.The only relevant quantity which remains is the aggregated power exchanges between the DSO and RTEM integrated over time.
From constraint (21) it holds which links the objective function to the spin flexibility.What this equation is telling us is that social welfare increases (OF  is minimised) when flexibility increases.Now depending on which additional constraint we use for the flexible amount (( 22) or ( 23)), there are two main scenarios to analyse: These two forms for constraining customers' flexibility lead to very different strategies in terms of consumers' interaction and cooperation (see arrows in Figure 1(b)).The constraint in 1 can be achieved by asking spin communities to adjust their flexibility over time independently of other communities; the only requirement is that at the end of the 24 h cycle each bus  renders ∑       = 0.This can be set into a linear programming problem for the variables   .On the other hand, in scenario 2 the objective function vanishes regardless of consumer's flexibility.However, in-buscommunities (i.e., spin communities) are forced to constrain their spin state every hour so that the whole system renders ∑ ∈        = 0.This requires a level of coordination among the bus communities at every hour .

Shiftable-Loads.
The constraint in (29) can be implemented as follows.First we solve the optimisation problem in (29).The output is the per-bus magnetisations   .Notice that the constraint can be factorised in .This way, minimising OF  is equivalent to solving For all  ∈ .Therefore the solutions   do not depend on  and we can write   ≡   , ∀ ∈ .In Figure 3 we show the solution along with the RTEM prices scaled in the following way: . This scaling makes both quantities comparable.Notice how the bus-flexibility solution follows the RTEM prices over time.This means that at demand peak times (10-13 h) and from (18-21 h) it is necessary to increase customers' flexibility.
Then we might ask how the spin system can comply with the total magnetisation imposed by the electrical constraints at every time  (see Figure 3).We need each spin-community   , ∀ ∈   to follow the constraint: This must hold for each bus-community of size   for all times.By redefining the spin states as   ≡ (1 +   )/2 with   ∈ {0, 1}, the constraint in (32) is equivalent to Since the spin system is discrete, the left hand side of (33) is a positive natural number: N + ∋  ≡ ∑ ∈    .Therefore, there is no solution if   is not a rational number: ,   = /; ,  ∈ Z, ,  ̸ = 0), there is a solution when (1)  = 0 and   is even; the solution consists in having half of the spins up and half down; (2)  =  with solution   = 1, ∀; (3)  = − with solution   = 0, ∀; (4) / ∈ (−1, 0) and   is a multiple of 2.
Notice that in Figure 3 all solutions are   = ±1 for all times except for  = 18, where we found   = −0.4.In the first case the only possible solutions are   = 1, ∀ if   = 1 and   = 0, ∀ if   = −1.However, in the latter case of   = −0.4= −2/5 customers will not be able to flip their states in order to attain this magnetisation; although we are in case (4), none of the available buses connects a number of homes which a multiple of 5 (Table 1).This evidences a potential limitation of a discrete model when it is coupled to the local power exchange problem.In the presence of electrical constraints, a discrete system will in general not be able to follow the continuum limit every time.A possible strategy to cope with this situation is to include an external effect (field) in the customer interactions to force their flexibility.In the Ising parlance, this is equivalent to setting  > 0 in the Ising Hamiltonian of ( 14).Then we can measure how customer's flexibility increases by setting  =18 = −1, evolve the system for different field intensities, and see if we reach  =19 = 1.In Figure 4 we show the final magnetisation for a pair of (, ) reached from a starting configuration with all spins down ( = −1) after 10k iterations using the Metropolis algorithm without any constraint.Magnetisations are averaged from 5000 Monte Carlo runs for each (, ) combination.As a reference we find the field strength  = 1.5 that makes the system reach   = 1 at the theoretical critical temperature   (blue dotted line in the inset).This is a way of forcing the so-called herding behaviour by imposing a top-down approach.Notice how the external action ramps up customers' flexibility.However, if temperature is high, the noise will destroy these dynamics and the intended constraint can not be reached.4.2.Self-Sustainability.Now let us explore a different situation where loads are forced to adjust so that ∑ ∈        = 0 at every time.A simple way to tackle this problem is by seeking solutions where the total magnetisation is either maximised or minimised.Hence, the problem can be set in terms of an Integer Linear Programming (ILP) problem of the form: Analogous for  in the shiftable-loads case, now the optimal solutions   do not depend on time and we can set   ≡   , ∀.For each case (max or min) we have the range of values for   = {−1, 0, 1}.In Figures 5(a) and 5(b) we show both solutions with our spin arrangement (Figure 1(c)).The maximisation solutions from the system in (34) result in 125 locations with   = −1, 132 with   = 1, and 32 with   = 0.
On the other hand, the minimisation results in 132 locations with   = −1, 125 with   = 1, and also 32 with   = 0. Therefore, the total magnetisation in both cases is 7 and −7 for the maximisation or minimisation problem, respectively.As in Scenario 1, spin communities with constraints   = 1 or   = −1 can be obtained by switching the spin states regardless of the size of the community.However for the   = 0 communities (shown in white color in Figure 5), a solution is only found when   is even.Therefore, buses  = 10, 28 in the maximisation and  = 3, 10 for the minimisation, respectively, will lack a solution as indicated by a red label in the figure.We have also verified this numerically by solving the following simple Binary Linear Programming (BLP)  problem for each bus-community (notice that minimising instead of maximising the problem will provide the same result; the only difference among solutions is the switching of states in the bus communities when   is an even number). max We show the spin solutions from (35) in Figures 5(c) and 5(d).
Since these constraints must hold for all times, it is interesting to check how robust the system is for complying with such magnetisations.To this end we test the robustness as follows: (1) Regularise buses with   = 0 and   odd by rewiring a random customer from another bus with   ̸ = 0.This results in a feasible spin configuration compatible with the constraints in (30).
(2) Perturb the feasible spin solution with strength  by switching the state of round(/ ⋅ 100) spins.
(3) Evolve the spin system with both the unconstrained Metropolis and with the CMH algorithm for a range of temperatures /  = 0.5, 0.6, . . ., 1.5.
(4) Measure the value of the constraint  = ∑         in both cases.
(5) In the resulting ensemble find the realisation with minimum  and monitor how || deviates from 0.
We found that the constraints-free evolution renders  values of 3 orders of magnitude larger than those obtained with the constrained version (CMH).The results of this experiment are shown in Figure 6.Here we compare on a normalised scale how the constraints deviate from 0 as we perturb the feasible spin solution with increasing perturbation strength.We measure the perturbation by the Hamming distance between the original and the perturbed configuration.Every T/Tc 0.5 0,2 0 ,4 0 ,6 0 ,8 1,0 0,0 distance from solution (norm.Hamming)  point represents the average value of 20 Monte Carlo replicas for the same parameters.As expected, the unconstrained evolution fails to provide a systematic trend in these dynamics since once the constraint is broken, there are no mechanisms to bend the spin dynamics towards regions with increasing feasibility.However, we find a pattern for the constrained evolution: as we deviate from the solution the value for  increases monotonically (worsening the feasibility).Notice that if perturbations are large, the minimum  is only found for high temperatures, since as we increase thermal fluctuations, the system will be able to explore different configurations and discover lower  values.The effect of the external field is different in both cases.In the unconstrained case the field only tilts up the values and reduces the noise.This is useless in our case since the field action reinforces the perturbations worsening the constraint.In the constrained case, however, the field linearises the point pattern on our scaling (we have added a dotted grey line for reference) and this improves feasibility.Therefore, an external action in the community would not be useful for recovering flexibility unless additional mechanisms are implemented to penalise deviations from the feasible region.These actions can be implemented smoothly because the drift from optimality increases monotonically with the perturbation.

The Aggregators' Perspective
From (17) and (25) we find that the total aggregation optimisation function renders which must attain a minimum value to maximise total aggregation profit.Each aggregator-to-DSO transaction has a price  A2DSO  bounded by constraint (10) and as stressed before: .Since we also have   ≥ 0, this problem will be in general unbounded in   (the power that customers buy from DSO).From the aggregator's profit perspective, the larger the quantity (  +  ), the higher their benefit (notice that in this case the flexibility constraints in ( 22) and ( 23) do not bound the   values).This way, maximum flexibility and maximum   will render an optimal benefit to aggregators without considering the optimum profit of the bus-loads.Hence, optimising the problem from the aggregators' point of view propels a bottom-up power transaction from bus-loads to aggregators, from aggregators to DSO, and from DSO to the real-time market, hierarchically.Moreover, the exchanged power benefit between the DSO and bus-loads is missing because this is irrelevant to aggregators.Consequently, as (36) shows, aggregators ask for the maximum loads' flexibility.In other words, the power flexibility is balanced with the power exchanged between the bus-loads to aggregators and with the power sent from DSO to bus-loads.Hence, from the aggregators' perspective, this will push the system to increase the flexibility from demand side and transacted power from the DSO to the loads.Notice that although this case is profitable for the DSO, it is not a profitable way for bus-loads to join this energy management approach.
One way to motivate bus-loads to join this setting is by letting   be a parameter instead of a variable in this optimisation problem.In particular-and without loss of generality-we can set   = 0. Also notice that, in optimality conditions, the variables  A2DSO Below we explore two analytical limits of (37) and how flexibility constraints in ( 22) and ( 23) affect the respective solutions.
(1) Aggregator Homogeneous Prices Limit.If all aggregators offer the same price for their transactions with the bus-loads,  L2A  →  L2A  (i.e., there is no heterogeneity in the distributed power system), the optimisation function does not depend on the load-to-bus mapping OF  → ∑  ( L2A  −  RT  ) ∑         .Also, if the self-sustainability flexibility constraint ( 22) is added to the problem, the function will be zeroed.This way, self-sustainability constraints do not have any effect on the total aggregation benefit because it makes the distributed power system an independent energy system.Hence, the system does not depend on the real-time electricity market to provide its local demand.
(2) Stationary Limit in Aggregator Prices.On the other hand, consider a situation where there is no price evolution over time (i.e., a fixed tariff scenario).In this case we get  L2A  →  L2A   ,  RT  →  RT , and the optimisation function in (37) renders OF  → ∑  Δ    ( L2S  −  RT ) ∑       .In this case the function is zeroed and again independent of the loadto-bus mapping if we impose the shiftable-loads flexibility constraint (23).This is because in a fixed tariff scenario shifting demand over time is not rational, and it does not render any extra profit for the distributed power system.

Optimal Load-to-Aggregator Mappings.
If we assume that customers are maximally flexible we can set   = 1 in the optimisation function and conjecture which bus-toaggregator mappings are more effective in different scenarios.In [22] the authors analysed the problem with 3 aggregators and other different configurations built by merging these 3 operators with aggregator prices in Table 3.Since for all times aggregator  1 always holds the minimum price, the optimal solution will map all buses to  1 .However, this is unrealistic since it is unlikely that all demand can be monopolised by a single aggregator.One way to cope with this is to force aggregators to split load as ∀ ∑  Δ  ≤ round(  /  ) + , where  ∈ N is a slack threshold for the uniform bus-to-aggregator mapping.By adding this constraint, the optimisation can be tackled by solving the BLP problem in the Δ  variables In Figure 7 we show the normalised value of the objective function in (38) as we increase the number of aggregators (in %) and for different  values.First we notice a global trend: the optimisation increases (lower objective function values) with the number of aggregators.Also  helps in this reduction by pulling the solutions downwards and by smoothing the peaks as it increases.For low  values the constraints are more stiff and the feasible space becomes fragmented.In the limit  = 0 (maximum homogeneity of load split) we find that for combinations of {5, 6,10,13,14,15,22,23,24,25,26,27,28,29, 30, 31} aggregators there are no optimal solutions (black dots).The reason for this is that for  = 0 and, say,   = 5, we are forcing to evenly distribute 32 buses among 5 aggregators so that at least one bus belongs to one aggregator.But each aggregator can only group round(32/5) = 6 buses at maximum.On the other hand, as  increases the peaks are smoothed, since the feasible space increases too.To test our model we compare our results (39) with the mapping in [22] displayed in Figure 1(d).The cardinality of  1 ,  2 ,  3 = (11, 11, 10) corresponds to our optimal solution from (38) for  = 0 and   = 3: 3,6,7,13,23,24,28,29,30  From these results we conclude that, in a maximum flexibility scenario, the total aggregation profit increases with the number of aggregators.In general, depending on the bus scheduled loads and the number of aggregators, demand cannot be split 100% evenly among aggregators; there must be some slack mechanisms to allow for flexibility in this mapping.Since customers group into buses and buses are mapped to aggregators, this discussion is relevant to the customer's flexibility as we show below with our Ising model.

Aggregator Boundaries and Their Effects.
Since aggregators operate on bus-load collections they might impose certain boundaries on the customers (and hence, in the associated spin system).These boundaries can be, for instance, bycontract forcing actions or other rules which have a greater or lesser influence on the customers' behaviour.The effect of this on the spin model is that agents will prefer to interact with neighbours who share the same aggregator.In our final experiment we measure the effect of aggregator boundaries by comparing a standard evolution (neighbourhood-free) with another evolution where spins interact only with neighbours belonging to the same aggregator.We start with a feasible solution for maximum flexibility  = 1 and then we measure how thermal noise worsens the solution.
In Figure 8 we show these results by comparing the boundary-free evolution with the aggregators' boundary dynamics with the reference boundaries shown in Figure 1(d).As expected, with very high temperature, noise will dominate interactions and the system will reach the  = 0 equilibrium as shown in Figure 2 (inset).However, the interesting dynamics occur precisely near criticality.Here we can monitor how sharply the flexibility () drops as  increases.The freeneighbourhood solution drops fast to  = 0 for /  > 1.5 (shown as a grey dotted line in the inset of the figure).However, notice how the neighbourhood constrained evolution slows down the worsening of the flexibility.Following our interpretation of temperature as uncertainty (lack of information), as  grows, shared opinions as to whether demand should be increased (spin-to-spin interactions) will not be strong enough to keep customers' confidence level in the flexibility program.Consequently, they will randomly decide whether to increase or decrease their demand and the constraint for maximum flexibility will start to deteriorate.However, when the information flow among customers is bounded within regions in the configuration space, this drop is smoothed, since the information does not propagate to the entire system in the same way.From our previous discussion about the limitations of a finite model for reproducing a phase transition, one might attempt to provide a semicritical exponent to these phase-like transitions.However, that is left for a future work.

Summary and Discussion
In this work we have analysed the customers' demand flexibility in a local power trading problem through an Ising spin-based model.We interpret spin states as customers' standpoint on following a flexibility program or not, which translates into an increase or decrease of their scheduled load.Spin-to-spin interactions simulate how these customer attitudes can change when information is retrieved from neighbouring customers.An external field is a proxy of topdown directives to enforce flexibility criteria and we associate the temperature in the spin system with uncertainty about the flexibility program.Further, we have addressed quantitative information about the two-way relationships between power Complexity exchanges and spin dynamics.These are formalised in terms of constraints, which force the spin system to evolve in different ways.To this end we provide a modified version of the Metropolis-Hastings algorithm including a gradient descent through the constraint surface.This implementation allowed us to analyse the system on a large scale (considering the cumulated benefit of all the actors involved) and also from the perspective of total aggregation.
At the global scale-welfare-we made two reasonable assumptions to limit customers' flexibility (shiftable-loads and self-sustainability).Each leads to two different scenarios which in turn motivate different types of analyses in our spin system.In the shiftable-loads scenario the maximum welfare requires each bus-community to meet a given value of flexibility every time.However, in general this is not possible in a discrete system.We provide conditions for this and also measure how an external field can force customers' flexibility.
Maximum welfare in the self-sustainability scenario requires all customers to meet the flexibility criteria every hour.We monitor the robustness of a feasible solution by perturbing the spin matrix and then measuring the feasibility of the perturbed magnetisation.When spins evolve by decreasing the constraint, there is an improvement of 3 orders of magnitude in the solutions with respect to the classical Metropolis evolution.An external action in the community would not be useful for recovering flexibility unless additional mechanisms are implemented to penalise deviations from the feasible region.
On the aggregation scale we analyse two limits: homogeneous and fixed electricity tariff of aggregators.These are interesting cases since flexibility constraints set the total aggregation optimisation function to zero: shiftable-loads in the fixed tariff case and self-sustainability in the case of homogeneous prices.Here we also address quantitative results of how aggregators can effectively split the total demand in the system and which are the implications of aggregator subcommunities in the spreading of flexible behaviour.We find that, in a maximum flexibility scenario, the total aggregation profit increases with the number of aggregators.In general, depending on the bus scheduled loads and the number of aggregators, demand cannot be split 100% evenly among aggregators; there must be some slack mechanisms to allow flexibility in this mapping.Finally, we check the effect of aggregator boundaries on spin dynamics following our interpretation of temperature as uncertainty (lack of information); as  grows shared opinions as to whether demand should be increased will not be strong enough to keep customers' confidence level in the flexibility program.Consequently, they will randomly decide whether demand should be increased or decreased and the constraint for maximum flexibility will start to deteriorate.However, when the information flow among customers is bounded within regions in the configuration space, it is smoothed, since the information does not propagate to the entire system in the same way.
In a future work we will make a more in-depth analysis of the relationships between more general spin system topologies and the flexibility constraints.Also we will develop a more thorough study of the phase transitions in presence of aggregator boundaries and their implications.Energy of the spin system :

Nomenclature
T emperatureofthespinsystem   : Criticaltemperature : Spin-to-spin interaction strength : Magnetic external field : Magnetisation of the spin system Γ: One configuration in the spin system Ω: Set 33-bus reference system and aggregators

Figure 1 :
Figure 1: Power exchange and the spin model.Each bus in (a), (c) groups customers interacting through the Ising Hamiltonian with binary states represented by up/down arrows (b).Aggregators share the demand in the bus-loads (c) and customers are arranged into a 2D lattice (d) in a way that preserves the bus distance partially.

( 1 )
spin states, (2) spin-to-spin interactions, (3) external field, (4) magnetisation, and (5) temperature.In our setting, we can parametrize the eagerness of consumer  to follow the flexibility program at time  with two state variables: (1) flexible   = 1 and (2) not flexible:   = −1.A possible and simple model for    is then    =      .

Figure 2 :
Figure 2: Autocorrelation function for the spin states between initial and evolved configurations.Grey dots represent numerical realisations and the blue line represents the average.The inset shows the phase transition in  beyond   (red line).

Figure 3 :
Figure 3: Resulting per-bus magnetisations.Blue dotted line shows the values for the normalised RTEM prices  RT  .The red dot represents a situation where the discrete spin system is unable to comply with the electrical constraint by flipping spin states.

Figure 4 :
Figure 4: Ramping up to meet maximum flexibility.Each curve represents the final magnetisation for a pair of (, ) reached from a starting configuration with all spins down ( = −1).The inset shows the curve for  = 1.5 (grey dots are numerical realisations and red line represents the averaged value) and the critical temperature (blue vertical line).

5 Figure 6 :
Figure 6: Constraint robustness in a self-sustainability scenario.(a, b) Unconstrained evolution.(c, d) Evolution with the CMH version.The size of the points is proportional to the temperature where the minimum value of || is found.

Figure 7 :
Figure 7: Bus-to-aggregator map efficiency.The normalised value of  for the optimal solution is shown for different number of aggregators and .Black circles represent feasible solutions for  = 0.

2 MFigure 8 :
Figure 8: Aggregator boundaries and total magnetisation.Standard evolution (red) compared with another evolution where spins interact only with neighbours belonging to the same aggregator (blue).Dots represent 100 realisations for each  and lines show averaged values (  shown in blue vertical line).The inset zooms in the solutions where there is a phase-like transition at /  (grey dotted line).
Notice that because   and   partition the sets  and , respectively, it holds that 2.1.Bus-Loads and Consumers.Each consumer  ∈  has an associated load   at time .Since each bus  ∈  connects |  | customers, the load at each bus  at time  is   = ∑ ∈    .Additionally, each consumer load can be split into a scheduled amount    and a flexible portion

Table 2 :
[22]ahead market power transactions    and real-time market price  RT  obtained from[22].
of possible spin configurations   : Loadofcustomer at time     : Expected load of customer  at time   Power bought to DSO by demand at bus  (customers connected to bus )   : NormalisedpowerboughttoDSOby demand at bus  at time    : Normalised power exchange for buy/sell between bus  and aggregator at time    : Fraction of power exchange between DSO and RTM reaching bus  in time   L2A  : Price for buy/sell between demand at bus and aggregator  at time   RT  : Real-time market price for DSO and RTM power exchanges  DSO2L : Price for DSO and load power exchanges  A2DSO : Profit guarantee factor of aggregator  at time  : Additional constraint imposed to the Metropolis-Hastings algorithm OF  : Optimisation function for bus-load  OF  : Optimisation function for all buses OF  : Optimisation function for aggregator  OF  : Optimisation function for all aggregators OF  : Optimisation function in the social welfare scenario : Slack threshold to relax constraints in (38).
:  : Real-time market price for DSO and aggregator  power exchanges Complexity 15