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

Authors: Francisco Prieto-Castrillo; Amin Shokri Gazafroudi; Javier Prieto; Juan Manuel Corchado
Journal: Complexity
 (2018)

---

## Abstract

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.


---

## Body

## 1. 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-calledherding 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 anddown). Agents interact with their neighbours locally and are also exposed to external action and to thermal noise. In the physics metaphor of a ferromagnetic material, the atoms tend to align their quantumspins to minimise energy. However, for high temperatures, the agents flip their spins randomly. There is a critical temperature Tc (known as Curie temperature) for which the system suffers a sharp change—second-order phase transition—from order to disorder. Near Tc 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 anduniversality) [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 T=0 [14]. In this regard, the temperature has been understood as a proxy of tolerance in binarythermodynamic 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, orstay 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 decision-making [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 ofaggregators 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 upper-level 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 real-time 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 real-time 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 Section2 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.


## 2. 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 Figure1, we describe schematically the power exchange problem and how it couples with the spin model. The power system at distribution level involves a hierarchical structure as shown in (a). Here, only DSO is able to exchange power PtRT with the RTEM at price λtRT. At time t each bus j interchanges power with both DSO (PjtDSO2L) and aggregator k (PjtL2A) at prices λDSO2L and λktL2A, respectively. Then, a generic aggregator k exchanges an amount of PktA2DSO power with DSO at price λktA2DSO. 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.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.
(a)
Hierarchical structure in the DPS (b)
Bus-loads and spin interactions (c)
33-bus reference system and aggregators (d)
Lattice arrangement of customersEach 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 Figure1(b). In this work, we use a 33-bus 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 N≈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.Table 1
Bus maximum scheduled load from [1] and number of homes assuming that each home demands approximately 13 kW as in [24].
Bus L j m a x [kW] Homes 1 100 8 2 90 7 3 120 9 4 60 5 5 60 5 6 200 15 7 200 15 8 60 5 9 60 5 10 45 3 11 60 5 12 60 5 13 120 9 14 60 5 15 60 5 16 60 5 17 90 7 18 90 7 19 90 7 20 90 7 21 90 7 22 90 7 23 420 32 24 420 32 25 60 5 26 60 5 27 60 5 28 120 9 29 200 15 30 150 12 31 210 16 32 60 5More rigorously, given a setC of N consumers C={1,…,N}, a set B of nJ electrical buses B={1,…,nJ}, and a set A of nA potential aggregators A={1,…,nA}, we arrange two partitions: in-bus-communities C=⨄j∈BBj and bus-aggregations B=⨄k∈AAk. This way Bj represents the subset of consumers with loads electrically connected to bus j and Ak stands for the subset of buses whose load is aggregated by aggregator k. These quantities enable us to define the following maps: (1) customer-bus map: Δcb:C×B→{0,1};(i,j)↦χ(Bj)(i) and (2) bus-aggregator map: Δba:B×A→{0,1};(j,k)↦χ(Ak)(j), where χA(x) represents the indicator function. To lighten notation, we will use Δij or Δjk to represent Δcb(i,j) or Δba(j,k), respectively, depending on the context. Notice that because Bj and Ak partition the sets C and B, respectively, it holds that(1)∑j∈BΔij=1,∑k∈AΔjk=1,∀i,j∈C,B.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.


### 2.1. Bus-Loads and Consumers

Each consumeri∈C has an associated load lit at time t. Since each bus j∈B connects Bj customers, the load at each bus j at time t is Ljt=∑i∈Bj‍lit. Additionally, each consumer load can be split into a scheduled amount litc and a flexible portion litf, which represents how customers can act as either upward or downward flexible loads: lit=litc-litf,∀i∈C,t. The sign convention here is that a positive flexibility litf>0 tends to reduce the scheduled load, whereas a negative flexibility litf<0 would increase customer’s expected demand. In other words, if litf>0 the corresponding customer decreases his day-ahead scheduled electrical demand in the real time. However, if litf<0 his real-time electrical demand is more than his day-ahead scheduled demand. From these expressions, we find that(2)Ljt=Ljtc-Ljtf,∀j∈B,t,where we have introduced consistently Ljtc=∑i∈Bj‍litc and Ljtf=∑i∈Bj‍litf. For the scheduled load at each bus j, we use the expression:(3)Ljtc=PtsLjmax∑j∈BLjmax,∀j∈B,t,where Pts represents the power at the source bus and Ljmax is the normalised expected load at bus j(see Table 1). The per-bus flexible component splits itself into the power exchanged with both DSO (PjtDSO2L) and aggregator (PjtL2A):(4)Ljtf=PjtL2A-PjtDSO2L,∀j∈B,t.According to the schematic shown in Figure 1(a) loads can only buy from the DSO and hence we have(5)PjtDSO2L≥0,∀j∈B,t.On the other hand, the bus-aggregator power exchanges are bidirectional. If PjtL2A<0 demand at bus j is buying from aggregator, the flexibility decreases. However, if PjtL2A>0 then the flexibility of the bus is increased. This can be interpreted asvirtual generation injected into the aggregator decreasing the scheduled load at time t.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:(6)∑j∈BLjtf=0,∀t,(7)∑tLjtf=0,∀j∈B.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 busj can be expressed as the trade-off between load bought from DSO at price λDSO2L and virtual generation sold to its aggregator k at price λktL2A integrated over time:(8)OFj∈Ak=λDSO2L∑tPjtDSO2L-∑tλktL2APjtL2A.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.


### 2.2. 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:(9)PktA2DSO=∑j∈AkPjtL2A,∀k∈A,t.Also there is a price model for these exchanges which is constrained in the following way:(10)δktλktL2A≤λktA2DSO≤λtRT,where δkt≥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 aggregator-to-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 δkt=1.1, λktL2A and λtRT related to the NordPool market. We summarise the respective values in Tables 2 and 3.Table 2
Day-ahead market power transactionsPts and real-time market price λtRT obtained from [22].
Time (h) P t s [kW] λ t R T [€/kW] 1 1114,50 0,13 2 1114,50 0,12 3 1300,25 0,15 4 1114,50 0,11 5 2972,00 0,30 6 2972,00 0,32 7 3343,50 0,35 8 3715,00 0,40 9 3715,00 0,42 10 6315,50 0,66 11 6687,00 0,71 12 6687,00 0,74 13 6315,50 0,69 14 3715,00 0,50 15 3715,00 0,41 16 3715,00 0,40 17 3715,00 0,42 18 5572,50 0,60 19 5944,00 0,65 20 6315,50 0,67 21 6501,25 0,70 22 2972,00 0,35 23 1857,50 0,28 24 1486,00 0,15Table 3
Aggregator prices table obtained from [22].
Time (h) λ k = 1 , t L 2 A [€/kW] λ k = 2 , t L 2 A [€/kW] λ k = 3 , t L 2 A [€/kW] 1 0,05 0,08 0,06 2 0,05 0,08 0,07 3 0,05 0,09 0,07 4 0,04 0,07 0,05 5 0,11 0,18 0,15 6 0,12 0,20 0,16 7 0,13 0,22 0,17 8 0,15 0,24 0,19 9 0,16 0,25 0,20 10 0,24 0,41 0,33 11 0,26 0,42 0,36 12 0,28 0,43 0,37 13 0,25 0,40 0,32 14 0,18 0,26 0,21 15 0,15 0,24 0,20 16 0,14 0,22 0,18 17 0,15 0,25 0,19 18 0,20 0,36 0,30 19 0,21 0,36 0,29 20 0,22 0,41 0,30 21 0,24 0,42 0,33 22 0,12 0,22 0,16 23 0,11 0,19 0,15 24 0,06 0,09 0,07The optimisation function for each aggregatork can be expressed as the accumulated balance over time between power PktA2DSO sold to DSO at price λktA2DSO and power PjtL2A bought from demand j, at price λktL2A:(11)OFk=∑t∑j∈AkλktL2APjtL2A-∑tλktA2DSOPktA2DSOwhich as in the case of demand must be minimised to reach a profitable situation from the aggregators’ perspective.


### 2.3. 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 power-flow 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:(12)PtRT+∑k∈APktA2DSO=∑j∈BPjtDSO2L.Therefore, the DSOs optimisation function can be expressed as the accumulated trade-off among these quantities, times the respective prices over time:(13)OFDSO=∑t∑k∈AλktA2DSOPktA2DSO+∑tλtRTPtRT-λDSO2L∑t∑j∈BPjtDSO2L.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.


## 2.1. Bus-Loads and Consumers

Each consumeri∈C has an associated load lit at time t. Since each bus j∈B connects Bj customers, the load at each bus j at time t is Ljt=∑i∈Bj‍lit. Additionally, each consumer load can be split into a scheduled amount litc and a flexible portion litf, which represents how customers can act as either upward or downward flexible loads: lit=litc-litf,∀i∈C,t. The sign convention here is that a positive flexibility litf>0 tends to reduce the scheduled load, whereas a negative flexibility litf<0 would increase customer’s expected demand. In other words, if litf>0 the corresponding customer decreases his day-ahead scheduled electrical demand in the real time. However, if litf<0 his real-time electrical demand is more than his day-ahead scheduled demand. From these expressions, we find that(2)Ljt=Ljtc-Ljtf,∀j∈B,t,where we have introduced consistently Ljtc=∑i∈Bj‍litc and Ljtf=∑i∈Bj‍litf. For the scheduled load at each bus j, we use the expression:(3)Ljtc=PtsLjmax∑j∈BLjmax,∀j∈B,t,where Pts represents the power at the source bus and Ljmax is the normalised expected load at bus j(see Table 1). The per-bus flexible component splits itself into the power exchanged with both DSO (PjtDSO2L) and aggregator (PjtL2A):(4)Ljtf=PjtL2A-PjtDSO2L,∀j∈B,t.According to the schematic shown in Figure 1(a) loads can only buy from the DSO and hence we have(5)PjtDSO2L≥0,∀j∈B,t.On the other hand, the bus-aggregator power exchanges are bidirectional. If PjtL2A<0 demand at bus j is buying from aggregator, the flexibility decreases. However, if PjtL2A>0 then the flexibility of the bus is increased. This can be interpreted asvirtual generation injected into the aggregator decreasing the scheduled load at time t.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:(6)∑j∈BLjtf=0,∀t,(7)∑tLjtf=0,∀j∈B.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 busj can be expressed as the trade-off between load bought from DSO at price λDSO2L and virtual generation sold to its aggregator k at price λktL2A integrated over time:(8)OFj∈Ak=λDSO2L∑tPjtDSO2L-∑tλktL2APjtL2A.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.


## 2.2. 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:(9)PktA2DSO=∑j∈AkPjtL2A,∀k∈A,t.Also there is a price model for these exchanges which is constrained in the following way:(10)δktλktL2A≤λktA2DSO≤λtRT,where δkt≥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 aggregator-to-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 δkt=1.1, λktL2A and λtRT related to the NordPool market. We summarise the respective values in Tables 2 and 3.Table 2
Day-ahead market power transactionsPts and real-time market price λtRT obtained from [22].
Time (h) P t s [kW] λ t R T [€/kW] 1 1114,50 0,13 2 1114,50 0,12 3 1300,25 0,15 4 1114,50 0,11 5 2972,00 0,30 6 2972,00 0,32 7 3343,50 0,35 8 3715,00 0,40 9 3715,00 0,42 10 6315,50 0,66 11 6687,00 0,71 12 6687,00 0,74 13 6315,50 0,69 14 3715,00 0,50 15 3715,00 0,41 16 3715,00 0,40 17 3715,00 0,42 18 5572,50 0,60 19 5944,00 0,65 20 6315,50 0,67 21 6501,25 0,70 22 2972,00 0,35 23 1857,50 0,28 24 1486,00 0,15Table 3
Aggregator prices table obtained from [22].
Time (h) λ k = 1 , t L 2 A [€/kW] λ k = 2 , t L 2 A [€/kW] λ k = 3 , t L 2 A [€/kW] 1 0,05 0,08 0,06 2 0,05 0,08 0,07 3 0,05 0,09 0,07 4 0,04 0,07 0,05 5 0,11 0,18 0,15 6 0,12 0,20 0,16 7 0,13 0,22 0,17 8 0,15 0,24 0,19 9 0,16 0,25 0,20 10 0,24 0,41 0,33 11 0,26 0,42 0,36 12 0,28 0,43 0,37 13 0,25 0,40 0,32 14 0,18 0,26 0,21 15 0,15 0,24 0,20 16 0,14 0,22 0,18 17 0,15 0,25 0,19 18 0,20 0,36 0,30 19 0,21 0,36 0,29 20 0,22 0,41 0,30 21 0,24 0,42 0,33 22 0,12 0,22 0,16 23 0,11 0,19 0,15 24 0,06 0,09 0,07The optimisation function for each aggregatork can be expressed as the accumulated balance over time between power PktA2DSO sold to DSO at price λktA2DSO and power PjtL2A bought from demand j, at price λktL2A:(11)OFk=∑t∑j∈AkλktL2APjtL2A-∑tλktA2DSOPktA2DSOwhich as in the case of demand must be minimised to reach a profitable situation from the aggregators’ perspective.


## 2.3. 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 power-flow 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:(12)PtRT+∑k∈APktA2DSO=∑j∈BPjtDSO2L.Therefore, the DSOs optimisation function can be expressed as the accumulated trade-off among these quantities, times the respective prices over time:(13)OFDSO=∑t∑k∈AλktA2DSOPktA2DSO+∑tλtRTPtRT-λDSO2L∑t∑j∈BPjtDSO2L.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.


## 3. An Ising Spin Model for Customers’ Flexibility

The Ising Hamiltonian (14) describes the interaction among entities (i.e., agents or spins) given their state si∈-1,1 and between each spin and a global magnetic field b. The coupling constant J measures the strength of spin-to-spin interactions.(14)H=-J∑i,jsisj-b∑isi.The notation i,j refers to pairs of spins belonging to the same radius of action or neighbourhood. When J>1 (ferromagnetism) spins tend to align in the same direction and if J<1 (antiferromagnetism) the spins tend to align in opposite directions. For J=0, there is no spin-to-spin interaction. The external action of a positive field b>0 will also foster positive spin alignments (and the other way around for b<0). For a given temperature T the probability for finding a spin configuration Γ={si} is proportional to the Boltzmann factor: exp⁡(-H/kBT), where kB stands for the Boltzmann constant. It is usual to take units such that J=1 and kB=1. An important magnitude is the magnetisation M=1/N∑i=1N‍si which measures the macroscopic effect of the spin states. In the so-called Mean Field approximation its value is Tc=2J/(kblog⁡(1+(2))).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 f at each step. Here f 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 reachesequilibrium. In operative terms this means that the correlations among the spin states are negligible at this point.Algorithm 1:Constrained Metropolis-Hastings (CMH).
INPUT:Δcb, Γ, f, T, num.Iters OUTPUT:Γt′ ( 1 )   i t e r ← 0 ( 2 ) while i t e r ≤ n u m . I t e r s do ( 3 ) n← pick a random spin from Γ ( 4 ) ΔE←2γnt∑i,n‍γit{energy change if spin flipped} ( 5 ) P(ΔE)←exp(-ΔE/T){Boltzmann’s probability for ΔE} ( 6 ) γit′←γit(1-2δin){new configuration if spin flipped} ( 7 ) Mjt′←1/nj∑Bj‍γit′{bus magnetization in new configuration} ( 8 ) iffMjt′≤fMjtthen ( 9 )  x← random number ∈[0,1] from uniform distribution ( 10 )    if(ΔE≤0orx≤P(ΔE)) then ( 11 )   Γt′←Γt∖{γnt}∪{-γnt}{flipspin} ( 12 )   end if ( 13 ) end if ( 14 ) iter←iter+1 ( 15 ) end whileThe 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 >5k iterations it is possible to achieve a reasonable value of decorrelation. In the inset, we show the phase transition in M beyond Tc. Notice how the numerical implementation is still able to reproduce qualitatively the transition M=±1→0 although in a smoothed and rounded way.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 inM beyond Tc (red line).When using this model one has to fix the interpretation for the following elements:(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 i to follow the flexibility program at time t with two state variables: (1) flexible γit=1 and (2) not flexible: γit=-1. A possible and simple model for litf is then(15)litf=litcγit.The spin-to-spin interaction can be a proxy for communication among customers and the external field b 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 M as a measure for the deviation from the fundamental value; if M is large agents are afraid of trading unless they have strong support from their private opinions or from their neighbours. In our case M=1 and M=-1 represent that the system has either maximum or minimum flexibility. A value of M 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 T 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:litc=ljtc>0,∀i∈Bj,j,t(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 Ljtc). Therefore, the flexible component of the load at bus j at time t renders(16)Ljtf=ljtcMjtnj,where nj=Bj represents the number of customers in bus j and Mjt stands for thebus magnetisationMjt≡(1/nj)∑i∈Bj‍γit. Let us also define dimension-free variables ∀j∈B,t as follows: ηjt≡PjtL2A/Ljtc and ψjt≡PjtDSO2L/Ljtc. Further, we assume that the real-time power exchanged from DSO and RTEM is a factor of the total load the DSO has to supply to buses. PtRT≡∑j∈B‍θjtLjtc. With these definitions and (4), we find(17)ηjt-ψjt=Mjt,∀j∈B,t.Notice that the power flow from aggregator k to DSO in (9) can be expressed as(18)PktA2DSO=∑j∈BΔjkljtcηjtnj.Also, from (1), (12), and (18) one can express the DSO power constraints as(19)∑j∈Bnjljtcθjt+ηjt-ψjt=0,∀t.Since ljtc≥0,∀j∈B,t and nj≥0,∀j∈B, inequality 5 is expressed as(20)ψjt≥0,∀j∈B,t.Notice that constraints (17) and (19) can be combined into(21)∑j∈Bnjljtcθjt+Mjt=0,∀t.Finally, flexibility constraints in (6), (7) can be expressed through (16) as(22)∑j∈BljtcnjMjt=0,∀t,(23)∑tljtcMjt=0,∀j∈B.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.


## 4. Large Scale Optimisation: Flexibility and Social Welfare

With the definitions above and from (8), we find the optimisation for each bus j:(24)OFj=∑t,k∈AnjljtcΔkjλDSO2Lψjt-λktL2Aηjt,∀j∈B,t,(where we have used OFj=∑k∈A‍ΔkjOFj∈Ak). The optimisation functions for each aggregator from (9) and (11) render(25)OFk=∑t,j∈BnjljtcΔkjλktL2A-λktA2DSOηjt,∀k∈A,t,where we have used ∑j∈Ak[·]=∑j∈B‍Δkj[·]. Finally, from (9) and (13), we find(26)OFDSO=∑t,k∈A,j∈BnjljtcΔkjλktA2DSOηjt+λtRTnAθjt-λDSO2LnAψjt,∀t.By defining OFL≡∑j∈B‍OFj and OFA≡∑k∈A‍OFk, the total load and total aggregation objective functions are found straightforward. Finally our objective is to maximise thesocial welfare of the distributed power system. Hence, the optimisation function is defined as OFW=OFL+OFA+OFDSO to ensure the maximised profits of all players (DSO, aggregator, and bus-loads) in the distributed power system. The global optimisation function OFW can be expressed as(27)OFW=∑tλtRT∑j∈Bnjljtcθjt.Notice how from this global scale the specifics about power exchanges between DSO-to-demand, demands-to-aggregators, 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(28)OFW=-∑tλtRT∑j∈BnjljtcMjt which links the objective function to the spin flexibility.What this equation is telling us is that social welfare increases (OFW 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: (29)
S
1
(Shiftable-Loads):
O
F
W
=
-
​
∑
j
∈
B
n
j
∑
t
λ
t
R
T
l
j
t
c
M
j
t
​
∑
t
l
j
t
c
M
j
t
=
0
,
∀
j
∈
B
,
(30)
S
2
(Self-Sustainability):
O
F
W
=
0
​
∑
j
∈
B
n
j
l
j
t
c
M
j
t
=
0
,
∀
t
. 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 S1 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 j renders ∑t‍ljtcMjt=0. This can be set into a linear programming problem for the variables Mjt. On the other hand, in scenario S2 the objective function vanishes regardless of consumer’s flexibility. However, in-bus-communities (i.e., spin communities) are forced to constrain their spin state every hour so that the whole system renders ∑j∈B‍njljtcMjt=0. This requires a level of coordination among the bus communities at every hour t.


### 4.1. 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 Mjt. Notice that the constraint can be factorised in j. This way, minimising OFW is equivalent to solving(31)min⁡-∑tλtRTPtsMjts.t.∑tPtsMjt=0.For all j∈B. Therefore the solutions Mjt do not depend on j and we can write Mt≡Mjt,∀j∈B. In Figure 3 we show the solution along with the RTEM prices scaled in the following way: λ¯tRT=2(λtRT-min⁡λtRT)/(max⁡λtRT-min⁡(λtRT))-1∈[-1,1]. 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.Figure 3
Resulting per-bus magnetisations. Blue dotted line shows the values for the normalised RTEM pricesλ¯tRT. The red dot represents a situation where the discrete spin system is unable to comply with the electrical constraint by flipping spin states.Then we might ask how the spin system can comply with the total magnetisation imposed by the electrical constraints at every timet (see Figure 3). We need each spin-community γit,∀i∈Bj to follow the constraint:(32)∑i∈Bjγit=njMt,∀j,t.This must hold for each bus-community of size nj for all times. By redefining the spin states as sit≡(1+γit)/2 with sit∈{0,1}, the constraint in (32) is equivalent to(33)∑i∈Bjsit=nj21+Mt,∀j,t.Since the spin system is discrete, the left hand side of (33) is a positive natural number: N+∋x≡∑i∈Bjsit. Therefore, there is no solution if Mt is not a rational number: Mt∉Q. If Mt∈Q (i.e., Mt=a/b;a,b∈Z,a,b≠0), there is a solution when (1)
a=0 and nj is even; the solution consists in having half of the spins up and half down;(2)
a=b with solution sit=1,∀i;(3)
a=-b with solution sit=0,∀i;(4)
a/b∈(-1,0) and nj is a multiple of 2b.Notice that in Figure3 all solutions are Mt=±1 for all times except for t=18, where we found Mt=-0.4. In the first case the only possible solutions are sit=1,∀i if Mt=1 and sit=0,∀i if Mt=-1. However, in the latter case of Mt=-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 is 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 settingb>0 in the Ising Hamiltonian of (14). Then we can measure how customer’s flexibility increases by setting Mt=18=-1, evolve the system for different field intensities, and see if we reach Mt=19=1. In Figure 4 we show the final magnetisation for a pair of (T,b) reached from a starting configuration with all spins down (M=-1) after 10k iterations using the Metropolis algorithm without any constraint. Magnetisations are averaged from 5000 Monte Carlo runs for each (T,b) combination. As a reference we find the field strength b=1.5 that makes the system reach Mt=1 at the theoretical critical temperature Tc (blue dotted line in the inset). This is a way of forcing the so-calledherding 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.Figure 4
Ramping up to meet maximum flexibility. Each curve represents the final magnetisation for a pair of(T,b) reached from a starting configuration with all spins down (M=-1). The inset shows the curve for b=1.5 (grey dots are numerical realisations and red line represents the averaged value) and the critical temperature (blue vertical line).


### 4.2. Self-Sustainability

Now let us explore a different situation where loads are forced to adjust so that∑j∈B‍njljtcMjt=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:(34)maxmin∑jMjts.t.∑jnjljtcMjt=0-1≤Mjt≤1Mjt∈Z.Analogous for j in the shiftable-loads case, now the optimal solutions Mjt do not depend on time and we can set Mj≡Mjt,∀t. For each case (max or min) we have the range of values for Mj=-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 Mj=-1, 132 with Mj=1, and 32 with Mj=0. On the other hand, the minimisation results in 132 locations with Mj=-1, 125 with Mj=1, and also 32 with Mj=0. Therefore, the total magnetisation in both cases is 7 and −7 for the maximisation or minimisation problem, respectively.Figure 5
Electrical constraints and spin solutions in a self-sustainability scenario. (a, b) Max/min solutions withMj=-1 (red), Mj=0 (white), and Mj=1 (blue); numbers indicate the bus number the customers belong to. (c, d) Spin solutions with γjt=-1 (red), γjt=1 (green), and grey color if there is no solution.
(a)
Maximisation:Mj (b)
Minimisation:Mj (c)
Maximisation: spin solutions (d)
Minimisation: spin solutionsAs in ScenarioS1, spin communities with constraints Mj=1 or Mj=-1 can be obtained by switching the spin states regardless of the size of the community. However for the Mj=0 communities (shown in white color in Figure 5), a solution is only found when nj is even. Therefore, buses j=10,28 in the maximisation and j=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 nj is an even number).(35)max⁡∑i=1njsis.t.∑i=1njsi=nj21+Mjsi∈0,1.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 howrobust the system is for complying with such magnetisations. To this end we test the robustness as follows:(1)
Regularise buses withMj=0 and nj odd by rewiring a random customer from another bus with Mj≠0. This results in a feasible spin configuration compatible with the constraints in (30).(2)
Perturb the feasible spin solution with strengthq by switching the state of round(q/N·100) spins.(3)
Evolve the spin system with both the unconstrained Metropolis and with the CMH algorithm for a range of temperaturesT/Tc=0.5,0.6,…,1.5.(4)
Measure the value of the constraintF=∑j‍njljtcMj in both cases.(5)
In the resulting ensemble find the realisation with minimumF and monitor how F deviates from 0.We found that the constraints-free evolution rendersF 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 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 F increases monotonically (worsening the feasibility). Notice that if perturbations are large, the minimum F is only found for high temperatures, since as we increase thermal fluctuations, the system will be able to explore different configurations and discover lower F values.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 ofF is found.
(a)
Unconstrained evolution:b=0 (b)
Unconstrained evolution:b=5 (c)
Constrained evolution:b=0 (d)
Constrained evolution:b=5The 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.


## 4.1. 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 Mjt. Notice that the constraint can be factorised in j. This way, minimising OFW is equivalent to solving(31)min⁡-∑tλtRTPtsMjts.t.∑tPtsMjt=0.For all j∈B. Therefore the solutions Mjt do not depend on j and we can write Mt≡Mjt,∀j∈B. In Figure 3 we show the solution along with the RTEM prices scaled in the following way: λ¯tRT=2(λtRT-min⁡λtRT)/(max⁡λtRT-min⁡(λtRT))-1∈[-1,1]. 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.Figure 3
Resulting per-bus magnetisations. Blue dotted line shows the values for the normalised RTEM pricesλ¯tRT. The red dot represents a situation where the discrete spin system is unable to comply with the electrical constraint by flipping spin states.Then we might ask how the spin system can comply with the total magnetisation imposed by the electrical constraints at every timet (see Figure 3). We need each spin-community γit,∀i∈Bj to follow the constraint:(32)∑i∈Bjγit=njMt,∀j,t.This must hold for each bus-community of size nj for all times. By redefining the spin states as sit≡(1+γit)/2 with sit∈{0,1}, the constraint in (32) is equivalent to(33)∑i∈Bjsit=nj21+Mt,∀j,t.Since the spin system is discrete, the left hand side of (33) is a positive natural number: N+∋x≡∑i∈Bjsit. Therefore, there is no solution if Mt is not a rational number: Mt∉Q. If Mt∈Q (i.e., Mt=a/b;a,b∈Z,a,b≠0), there is a solution when (1)
a=0 and nj is even; the solution consists in having half of the spins up and half down;(2)
a=b with solution sit=1,∀i;(3)
a=-b with solution sit=0,∀i;(4)
a/b∈(-1,0) and nj is a multiple of 2b.Notice that in Figure3 all solutions are Mt=±1 for all times except for t=18, where we found Mt=-0.4. In the first case the only possible solutions are sit=1,∀i if Mt=1 and sit=0,∀i if Mt=-1. However, in the latter case of Mt=-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 is 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 settingb>0 in the Ising Hamiltonian of (14). Then we can measure how customer’s flexibility increases by setting Mt=18=-1, evolve the system for different field intensities, and see if we reach Mt=19=1. In Figure 4 we show the final magnetisation for a pair of (T,b) reached from a starting configuration with all spins down (M=-1) after 10k iterations using the Metropolis algorithm without any constraint. Magnetisations are averaged from 5000 Monte Carlo runs for each (T,b) combination. As a reference we find the field strength b=1.5 that makes the system reach Mt=1 at the theoretical critical temperature Tc (blue dotted line in the inset). This is a way of forcing the so-calledherding 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.Figure 4
Ramping up to meet maximum flexibility. Each curve represents the final magnetisation for a pair of(T,b) reached from a starting configuration with all spins down (M=-1). The inset shows the curve for b=1.5 (grey dots are numerical realisations and red line represents the averaged value) and the critical temperature (blue vertical line).


## 4.2. Self-Sustainability

Now let us explore a different situation where loads are forced to adjust so that∑j∈B‍njljtcMjt=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:(34)maxmin∑jMjts.t.∑jnjljtcMjt=0-1≤Mjt≤1Mjt∈Z.Analogous for j in the shiftable-loads case, now the optimal solutions Mjt do not depend on time and we can set Mj≡Mjt,∀t. For each case (max or min) we have the range of values for Mj=-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 Mj=-1, 132 with Mj=1, and 32 with Mj=0. On the other hand, the minimisation results in 132 locations with Mj=-1, 125 with Mj=1, and also 32 with Mj=0. Therefore, the total magnetisation in both cases is 7 and −7 for the maximisation or minimisation problem, respectively.Figure 5
Electrical constraints and spin solutions in a self-sustainability scenario. (a, b) Max/min solutions withMj=-1 (red), Mj=0 (white), and Mj=1 (blue); numbers indicate the bus number the customers belong to. (c, d) Spin solutions with γjt=-1 (red), γjt=1 (green), and grey color if there is no solution.
(a)
Maximisation:Mj (b)
Minimisation:Mj (c)
Maximisation: spin solutions (d)
Minimisation: spin solutionsAs in ScenarioS1, spin communities with constraints Mj=1 or Mj=-1 can be obtained by switching the spin states regardless of the size of the community. However for the Mj=0 communities (shown in white color in Figure 5), a solution is only found when nj is even. Therefore, buses j=10,28 in the maximisation and j=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 nj is an even number).(35)max⁡∑i=1njsis.t.∑i=1njsi=nj21+Mjsi∈0,1.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 howrobust the system is for complying with such magnetisations. To this end we test the robustness as follows:(1)
Regularise buses withMj=0 and nj odd by rewiring a random customer from another bus with Mj≠0. This results in a feasible spin configuration compatible with the constraints in (30).(2)
Perturb the feasible spin solution with strengthq by switching the state of round(q/N·100) spins.(3)
Evolve the spin system with both the unconstrained Metropolis and with the CMH algorithm for a range of temperaturesT/Tc=0.5,0.6,…,1.5.(4)
Measure the value of the constraintF=∑j‍njljtcMj in both cases.(5)
In the resulting ensemble find the realisation with minimumF and monitor how F deviates from 0.We found that the constraints-free evolution rendersF 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 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 F increases monotonically (worsening the feasibility). Notice that if perturbations are large, the minimum F is only found for high temperatures, since as we increase thermal fluctuations, the system will be able to explore different configurations and discover lower F values.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 ofF is found.
(a)
Unconstrained evolution:b=0 (b)
Unconstrained evolution:b=5 (c)
Constrained evolution:b=0 (d)
Constrained evolution:b=5The 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.


## 5. The Aggregators’ Perspective

From (17) and (25) we find that the total aggregation optimisation function renders(36)OFA=∑tjnjljtcψjt+Mjt∑kΔkjλktL2A-λktA2DSOwhich must attain a minimum value to maximise total aggregation profit. Each aggregator-to-DSO transaction has a price λktA2DSO bounded by constraint (10) and as stressed before: λktA2DSO≥δktλktL2A. Since we also have ψjt≥0, this problem will be in general unbounded in ψjt (the power that customers buy from DSO). From the aggregator’s profit perspective, the larger the quantity (ψjt+Mjt), the higher their benefit (notice that in this case the flexibility constraints in (22) and (23) do not bound the ψjt values). This way, maximum flexibility and maximum ψjt 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ψjt be a parameter instead of a variable in this optimisation problem. In particular—and without loss of generality—we can set ψjt=0. Also notice that, in optimality conditions, the variables λktA2DSO will attain their maximum values in (10): λktA2DSO=λtRT. The modified aggregation optimisation function is then expressed as(37)OF¯A=∑tjnjljtcMjt∑kΔkjλktL2A-λtRT.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, λktL2A→λtL2A (i.e., there is no heterogeneity in the distributed power system), the optimisation function does not depend on the load-to-bus mapping OF¯A→∑t‍(λtL2A-λtRT)∑j‍njljtcMjt. 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 λktL2A→λkL2A, λtRT→λRT, and the optimisation function in (37) renders OF¯A→∑kj‍Δjknj(λkL2S-λRT)∑t‍ljtcMjt. In this case the function is zeroed and again independent of the load-to-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.


### 5.1. Optimal Load-to-Aggregator Mappings

If we assume that customers are maximally flexible we can setMjt=1 in the optimisation function and conjecture which bus-to-aggregator 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 A1 always holds the minimum price, the optimal solution will map all buses to A1. 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 ∀k∑j‍Δjk≤round(nJ/nA)+q, where q∈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 Δjk variables(38)min⁡z=∑tjnjljtc∑kΔkjλktL2A-λtRTs.t.∑kΔjk=1,∀j∑jΔjk≤roundnJnA+q,∀kΔjk∈0,1.In Figure7 we show the normalised value of the objective function in (38) as we increase the number of aggregators (in %) and for different q values. First we notice a global trend: the optimisation increases (lower objective function values) with the number of aggregators. Also q helps in this reduction by pulling the solutions downwards and by smoothing the peaks as it increases. For low q values the constraints are more stiff and the feasible space becomes fragmented. In the limit q=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 q=0 and, say, nA=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 q 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 A1, A2, A3=(11,11,10) corresponds to our optimal solution from (38) for q=0 and nA=3:(39)A1=1,3,6,7,13,23,24,28,29,30,31,A2=2,4,5,8,9,17,18,19,20,21,22,A3=10,11,12,14,15,16,25,26,27,32.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.Figure 7
Bus-to-aggregator map efficiency. The normalised value ofz for the optimal solution is shown for different number of aggregators and q. Black circles represent feasible solutions for q=0.


### 5.2. Aggregator Boundaries and Their Effects

Since aggregators operate on bus-load collections they might impose certainboundaries on the customers (and hence, in the associated spin system). These boundaries can be, for instance, by-contract 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 M=1 and then we measure how thermal noise worsens the solution.In Figure8 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 M=0 equilibrium as shown in Figure 2 (inset). However, the interesting dynamics occur precisely near criticality. Here we can monitor how sharply the flexibility (M) drops as T increases. The free-neighbourhood solution drops fast to M=0 for T/Tc>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 T 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 asemicritical exponent to these phase-like transitions. However, that is left for a future work.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 eachT and lines show averaged values (Tc shown in blue vertical line). The inset zooms in the solutions where there is a phase-like transition at T/Tc (grey dotted line).


## 5.1. Optimal Load-to-Aggregator Mappings

If we assume that customers are maximally flexible we can setMjt=1 in the optimisation function and conjecture which bus-to-aggregator 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 A1 always holds the minimum price, the optimal solution will map all buses to A1. 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 ∀k∑j‍Δjk≤round(nJ/nA)+q, where q∈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 Δjk variables(38)min⁡z=∑tjnjljtc∑kΔkjλktL2A-λtRTs.t.∑kΔjk=1,∀j∑jΔjk≤roundnJnA+q,∀kΔjk∈0,1.In Figure7 we show the normalised value of the objective function in (38) as we increase the number of aggregators (in %) and for different q values. First we notice a global trend: the optimisation increases (lower objective function values) with the number of aggregators. Also q helps in this reduction by pulling the solutions downwards and by smoothing the peaks as it increases. For low q values the constraints are more stiff and the feasible space becomes fragmented. In the limit q=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 q=0 and, say, nA=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 q 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 A1, A2, A3=(11,11,10) corresponds to our optimal solution from (38) for q=0 and nA=3:(39)A1=1,3,6,7,13,23,24,28,29,30,31,A2=2,4,5,8,9,17,18,19,20,21,22,A3=10,11,12,14,15,16,25,26,27,32.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.Figure 7
Bus-to-aggregator map efficiency. The normalised value ofz for the optimal solution is shown for different number of aggregators and q. Black circles represent feasible solutions for q=0.


## 5.2. Aggregator Boundaries and Their Effects

Since aggregators operate on bus-load collections they might impose certainboundaries on the customers (and hence, in the associated spin system). These boundaries can be, for instance, by-contract 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 M=1 and then we measure how thermal noise worsens the solution.In Figure8 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 M=0 equilibrium as shown in Figure 2 (inset). However, the interesting dynamics occur precisely near criticality. Here we can monitor how sharply the flexibility (M) drops as T increases. The free-neighbourhood solution drops fast to M=0 for T/Tc>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 T 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 asemicritical exponent to these phase-like transitions. However, that is left for a future work.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 eachT and lines show averaged values (Tc shown in blue vertical line). The inset zooms in the solutions where there is a phase-like transition at T/Tc (grey dotted line).


## 6. 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 top-down 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 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); asT 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.



---

*Source: 5905932-2018-05-08.xml*