Research Article  Open Access
Francisco PrietoCastrillo, Amin Shokri Gazafroudi, Javier Prieto, Juan Manuel Corchado, "An Ising SpinBased Model to Explore Efficient Flexibility in Distributed Power Systems", Complexity, vol. 2018, Article ID 5905932, 16 pages, 2018. https://doi.org/10.1155/2018/5905932
An Ising SpinBased Model to Explore Efficient Flexibility in Distributed Power Systems
Abstract
This paper analyses customers’ demand flexibility in a local power trading scenario through an Ising spinbased model. We look at quantitative information on the twoway relationships between power exchanges and spin dynamics. To this end, a modified version of the MetropolisHastings 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.
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 selforganised power systems? In this work, we explore these matters analytically and numerically. An Ising spinbased 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 socalled herding behaviour or the economic bubbles are phenomena lying outside analyses that neglect large correlations and selforganisation 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. Isingbased models are a promising approach as we demonstrate here.
The LenzIsing 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 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—secondorder 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 Isinglike 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 [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 smallscale 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 lowerlevel 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 selfcontrolled. Here, there is a supersystem 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 spinglass 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 realtime electricity market (RTEM). Moreover, in this paper, we assume that DSOs act as agents able to participate in the realtime electricity markets to trade realtime power in a twoway 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 Isingbased 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.
2. Hierarchical Model and Agent Interactions in the Distributed Power Grid
In our setting, we consider a realtime interplay structure among five kinds of actors: consumers, busloads, aggregators, DSOs, and RTEM. In Figure 1, 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 with the RTEM at price . At time each bus interchanges power with both DSO () and aggregator () at prices and , respectively. Then, a generic aggregator exchanges an amount of power with DSO at price . In this way, the power transaction between busloads and aggregators is twoway according to the flexibility behaviour from the demand side. Moreover, there are twoway power exchanges between aggregator and DSO. However, the power transaction between DSO and busloads is oneway: only from DSO to busloads. In other words, busloads can only buy realtime power from the DSO, while they are able to buy/sell power from/to aggregators. Also, only DSO can participate in the RTEM.
(a) Hierarchical structure in the DPS
(b) Busloads and spin interactions
(c) 33bus reference system and aggregators
(d) Lattice arrangement of customers
Each busload 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 customers. Finally in Figure 1(d), we show a 2D grid arrangement where each of the cells represents a customer. Furthermore, we assume that in the busloads there exists some notion of topological closeness (e.g., spatial or electrical proximity) which is partially mapped onto the square lattice.

More rigorously, given a set of consumers , a set of electrical buses , and a set of potential aggregators , we arrange two partitions: inbuscommunities and busaggregations . This way represents the subset of consumers with loads electrically connected to bus and stands for the subset of buses whose load is aggregated by aggregator . These quantities enable us to define the following maps: customerbus map: and busaggregator map: , where represents the indicator function. To lighten notation, we will use or to represent or , respectively, depending on the context. Notice that because and partition the sets and , respectively, it holds thatThis 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. BusLoads 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 , which represents how customers can act as either upward or downward flexible loads: . The sign convention here is that a positive flexibility tends to reduce the scheduled load, whereas a negative flexibility would increase customer’s expected demand. In other words, if the corresponding customer decreases his dayahead scheduled electrical demand in the real time. However, if his realtime electrical demand is more than his dayahead scheduled demand. From these expressions, we find thatwhere 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 is the normalised expected load at bus (see Table 1). The perbus flexible component splits itself into the power exchanged with both DSO () and aggregator ():According to the schematic shown in Figure 1(a) loads can only buy from the DSO and hence we haveOn the other hand, the busaggregator power exchanges are bidirectional. If demand at bus is buying from aggregator, the flexibility decreases. However, if 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 selfsustainability 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 shiftableloads (i.e., loads that can be shifted over time). On the other hand, (6) increases the selfsustainability 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 tradeoff between load bought from DSO at price and virtual generation sold to its aggregator at price 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.
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:Also there is a price model for these exchanges which is constrained in the following way:where represents a lower bound threshold for the aggregatortoDSO price compared to loadtoaggregator price. This ensures aggregators’ profit, which makes it reasonable for them to be part of the market. Also, the aggregatortoDSO price is limited by the DSOtomarket price; this upper bound acts as a price control, limiting the bidding price of aggregators below the realtime price. Here we will use reference values from [22] for , and 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 sold to DSO at price and power bought from demand , at price :which 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 powerflow balance, the power sold by the DSO to the whole demand—busloads—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 tradeoff 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.
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 and between each spin and a global magnetic field . The coupling constant measures the strength of spintospin interactions.The notation refers to pairs of spins belonging to the same radius of action or neighbourhood. When (ferromagnetism) spins tend to align in the same direction and if (antiferromagnetism) the spins tend to align in opposite directions. For , there is no spintospin interaction. The external action of a positive field will also foster positive spin alignments (and the other way around for ). For a given temperature the probability for finding a spin configuration is proportional to the Boltzmann factor: , where stands for the Boltzmann constant. It is usual to take units such that and . An important magnitude is the magnetisation which measures the macroscopic effect of the spin states. In the socalled Mean Field approximation its value is .
One way to implement the Ising model numerically is through the MetropolisHastings 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 spin system shown in Figure 1(d). For 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 although in a smoothed and rounded way.
When using this model one has to fix the interpretation for the following elements: spin states, spintospin interactions, external field, magnetisation, and temperature. In our setting, we can parametrize the eagerness of consumer to follow the flexibility program at time with two state variables: flexible and not flexible: . A possible and simple model for is thenThe spintospin interaction can be a proxy for communication among customers and the external field can simulate a topdown 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 and 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: (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 perbus load ). Therefore, the flexible component of the load at bus at time renderswhere represents the number of customers in bus and stands for the bus magnetisation . Let us also define dimensionfree variables as follows: and . Further, we assume that the realtime power exchanged from DSO and RTEM is a factor of the total load the DSO has to supply to buses. . With these definitions and (4), we findNotice that the power flow from aggregator to DSO in (9) can be expressed asAlso, from (1), (12), and (18) one can express the DSO power constraints asSince and , inequality 5 is expressed asNotice that constraints (17) and (19) can be combined intoFinally, flexibility constraints in (6), (7) can be expressed through (16) asBy 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 :(where we have used ). The optimisation functions for each aggregator from (9) and (11) renderwhere we have used . Finally, from (9) and (13), we findBy defining and , 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 to ensure the maximised profits of all players (DSO, aggregator, and busloads) in the distributed power system. The global optimisation function can be expressed asNotice how from this global scale the specifics about power exchanges between DSOtodemand, demandstoaggregators, and aggregatorstoDSO 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 ( 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 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 . This can be set into a linear programming problem for the variables . On the other hand, in scenario the objective function vanishes regardless of consumer’s flexibility. However, inbuscommunities (i.e., spin communities) are forced to constrain their spin state every hour so that the whole system renders . This requires a level of coordination among the bus communities at every hour .
4.1. ShiftableLoads
The constraint in (29) can be implemented as follows. First we solve the optimisation problem in (29). The output is the perbus magnetisations . Notice that the constraint can be factorised in . This way, minimising is equivalent to solvingFor 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 busflexibility 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 spincommunity to follow the constraint:This must hold for each buscommunity of size for all times. By redefining the spin states as with , the constraint in (32) is equivalent toSince the spin system is discrete, the left hand side of (33) is a positive natural number: . Therefore, there is no solution if is not a rational number: . If (i.e., ), there is a solution when (1) and is even; the solution consists in having half of the spins up and half down;(2) with solution ;(3) with solution ;(4) and is a multiple of .
Notice that in Figure 3 all solutions are for all times except for , where we found . In the first case the only possible solutions are if and if . However, in the latter case of customers will not be able to flip their states in order to attain this magnetisation; although we are in case , none of the available buses connects a number of homes which is a multiple of (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 in the Ising Hamiltonian of (14). Then we can measure how customer’s flexibility increases by setting , evolve the system for different field intensities, and see if we reach . In Figure 4 we show the final magnetisation for a pair of reached from a starting configuration with all spins down () 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 that makes the system reach at the theoretical critical temperature (blue dotted line in the inset). This is a way of forcing the socalled herding behaviour by imposing a topdown 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. SelfSustainability
Now let us explore a different situation where loads are forced to adjust so that 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 shiftableloads 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 . 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 , 132 with , and 32 with . On the other hand, the minimisation results in 132 locations with , 125 with , and also 32 with . Therefore, the total magnetisation in both cases is 7 and −7 for the maximisation or minimisation problem, respectively.
(a) Maximisation: Mj
(b) Minimisation: Mj
(c) Maximisation: spin solutions
(d) Minimisation: spin solutions
As in Scenario , spin communities with constraints or can be obtained by switching the spin states regardless of the size of the community. However for the communities (shown in white color in Figure 5), a solution is only found when is even. Therefore, buses in the maximisation and 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 buscommunity (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).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 and odd by rewiring a random customer from another bus with . 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 spins.(3)Evolve the spin system with both the unconstrained Metropolis and with the CMH algorithm for a range of temperatures .(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 constraintsfree 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 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.
(a) Unconstrained evolution:
(b) Unconstrained evolution:
(c) Constrained evolution:
(d) Constrained evolution:
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.
5. The Aggregators’ Perspective
From (17) and (25) we find that the total aggregation optimisation function renderswhich must attain a minimum value to maximise total aggregation profit. Each aggregatortoDSO transaction has a price bounded by constraint (10) and as stressed before: . Since we also have , 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 busloads. Hence, optimising the problem from the aggregators’ point of view propels a bottomup power transaction from busloads to aggregators, from aggregators to DSO, and from DSO to the realtime market, hierarchically. Moreover, the exchanged power benefit between the DSO and busloads 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 busloads to aggregators and with the power sent from DSO to busloads. 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 busloads to join this energy management approach.
One way to motivate busloads 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 . Also notice that, in optimality conditions, the variables will attain their maximum values in (10): . The modified aggregation optimisation function is then expressed asBelow 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 busloads, (i.e., there is no heterogeneity in the distributed power system), the optimisation function does not depend on the loadtobus mapping . Also, if the selfsustainability flexibility constraint (22) is added to the problem, the function will be zeroed. This way, selfsustainability 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 realtime 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 , , and the optimisation function in (37) renders . In this case the function is zeroed and again independent of the loadtobus mapping if we impose the shiftableloads 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 LoadtoAggregator Mappings
If we assume that customers are maximally flexible we can set in the optimisation function and conjecture which bustoaggregator 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 always holds the minimum price, the optimal solution will map all buses to . 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 , where is a slack threshold for the uniform bustoaggregator 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 (maximum homogeneity of load split) we find that for combinations of aggregators there are no optimal solutions (black dots). The reason for this is that for and, say, , 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 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 , , corresponds to our optimal solution from (38) for and :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.
5.2. Aggregator Boundaries and Their Effects
Since aggregators operate on busload 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 (neighbourhoodfree) with another evolution where spins interact only with neighbours belonging to the same aggregator. We start with a feasible solution for maximum flexibility and then we measure how thermal noise worsens the solution.
In Figure 8 we show these results by comparing the boundaryfree 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 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 for (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 (spintospin 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 phaselike transitions. However, that is left for a future work.
6. Summary and Discussion
In this work we have analysed the customers’ demand flexibility in a local power trading problem through an Ising spinbased 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. Spintospin 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 twoway 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 MetropolisHastings 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 (shiftableloads and selfsustainability). Each leads to two different scenarios which in turn motivate different types of analyses in our spin system. In the shiftableloads scenario the maximum welfare requires each buscommunity 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 selfsustainability 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: shiftableloads in the fixed tariff case and selfsustainability 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 indepth 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.
Nomenclature
:  Number of customers 
:  Set of customers 
:  Set of busloads 
:  Customers connected to bus 
:  Set of aggregators 
:  Buses managed by aggregator 
:  Number of buses 
:  Number of aggregators 
:  Number of hours 
:  Number of customers connected at bus 
:  = 1(0) if customer is connected (not connected) to bus 
:  = 1(0) if bus is managed (not managed) by aggregator 
:  Ising Hamiltonian 
:  Flexibility for customer 
:  State of spin 
:  Energy of the spin system 
:  Temperature of the spin system 
:  Critical temperature 
:  Spintospin interaction strength 
:  Magnetic external field 
:  Magnetisation of the spin system 
:  One configuration in the spin system 
:  Set of possible spin configurations 
:  Load of customer at time 
:  Expected load of customer at time 
:  Flexible proportion of load of customer at time 
:  Expected load in bus at time 
:  Flexible proportion of load in bus at time 
:  Power at the source bus at time 
:  Normalised expected load at bus 
:  Power exchange for buy/sell between bus and aggregator 
:  Realtime power exchange between DSO and RTM at time 
:  Power exchange or buy/sell between aggregator and aggregator DSO 
:  Power bought to DSO by demand at bus (customers connected to bus ) 
:  Normalised power bought to DSO by 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 
:  Price for buy/sell between demand at bus and aggregator at time 
:  Realtime market price for DSO and RTM power exchanges 
:  Price for DSO and load power exchanges 
:  Realtime market price for DSO and aggregator power exchanges 
:  Profit guarantee factor of aggregator at time 
:  Additional constraint imposed to the MetropolisHastings algorithm 
:  Optimisation function for busload 
:  Optimisation function for all buses 
:  Optimisation function for aggregator 
:  Optimisation function for all aggregators 
:  Optimisation function in the social welfare scenario 
:  Slack threshold to relax constraints in (38). 
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was partially supported by the Regional Ministry of Education from Castilla y León (Spain) and the European Social Fund under the MOVIURBAN project with Ref. SA070U16. The authors also acknowledge the European Commission H2020 MSCARISE2014: Marie SklodowskaCurie project DREAMGO; Enabling Demand Response for Short and RealTime Efficient and Market Based Smart Operation; An Intelligent and RealTime Simulation Approach with Ref. 641794. Amin Shokri Gazafroudi also acknowledges the support by the Ministry of Education of Junta de Castilla y León and the European Social Fund through a grant from predoctoral recruitment or research personnel associated with the research project “Arquitectura Multiagente para la Gestión Eficaz de Redes de Energía a Través del Uso de Técnicas de Inteligencia Artificial” of the University of Salamanca.
References
 N. Mithulananthan, Q. Duong Hung, and Y. Kwang, Intelligent Network Integration of Distributed Renewable Generation, Springer International Publishing, 2017.
 A. S. Gazafroudi, F. PrietoCastrillo, T. Pinto, and J. M. Corchado, “Organizationbased multiagent system of local electricity market: Bottomup approach,” Advances in Intelligent Systems and Computing, vol. 619, pp. 281–283, 2017. View at: Publisher Site  Google Scholar
 M. A. Fotouhi Ghazvini, J. Soares, O. Abrishambaf, R. Castro, and Z. Vale, “Demand response implementation in smart households,” Energy and Buildings, vol. 143, pp. 129–148, 2017. View at: Publisher Site  Google Scholar
 S. S. Torbaghan, N. Blaauwbroek, P. Nguyen, and M. Gibescu, “Local market framework for exploiting flexibility from the end users,” in Proceedings of the 13th International Conference on the European Energy Market, EEM 2016, pp. 1–6, Portugal, June 2016. View at: Publisher Site  Google Scholar
 A. Gazafroudi, F. PrietoCastrillo, T. Pinto, J. Prieto, J. Corchado, and J. Bajo, “Energy flexibility management based on predictive dispatch model of domestic energy management system,” Energies, vol. 10, no. 9, p. 1397, 2017. View at: Publisher Site  Google Scholar
 G. Ajmone Marsan, N. Bellomo, and A. Tosin, Complex Systems and Society Modeling and Simulation, Springer International Publishing, New York, NY, USA, 2013. View at: Publisher Site  MathSciNet
 A. Carro, R. Toral, and M. S. Miguel, “Markets, herding and response to external information,” PLoS ONE, vol. 10, no. 7, Article ID e0133287, 2015. View at: Publisher Site  Google Scholar
 S. G. Brush, “History of the LenzIsing model,” Reviews of Modern Physics, vol. 39, no. 4, pp. 883–893, 1967. View at: Publisher Site  Google Scholar
 R. V. Sole, Phase Transitions, Springer International Publishing, 2011.
 R. Schlicht and Y. Iwasa, “Forest gap dynamics and the Ising model,” Journal of Theoretical Biology, vol. 230, no. 1, pp. 65–75, 2004. View at: Publisher Site  Google Scholar  MathSciNet
 D. Nettle, “Is the rate of linguistic change constant?” Lingua, vol. 108, no. 23, pp. 119–136, 1999. View at: Publisher Site  Google Scholar
 A. M. Guenault, Statistical Physics, Springer International Publishing, Dordrecht, Netherlands, 2007. View at: Publisher Site
 T. C. Schelling, “Dynamic models of segregation,” Journal of Mathematical Sociology, vol. 1, pp. 143–186, 1971. View at: Publisher Site  Google Scholar
 D. Stauffer, “Social applications of twodimensional Ising models,” American Journal of Physics, vol. 76, no. 45, pp. 470–472, 2008. View at: Publisher Site  Google Scholar
 J. Mimkes, “Binary alloys as a model for the multicultural society,” Journal of Thermal Analysis and Calorimetry, vol. 43, no. 2, pp. 521–537, 1995. View at: Publisher Site  Google Scholar
 S. Bornholdt, “Expectation bubbles in a spin model of markets: Intermittency from frustration across scales,” International Journal of Modern Physics C, vol. 12, no. 5, pp. 667–674, 2001. View at: Publisher Site  Google Scholar
 P. Sieczka and J. A. Holyst, “A threshold model of financial markets,” Acta Physica Polonica A, vol. 114, no. 3, pp. 525–530, 2008. View at: Publisher Site  Google Scholar
 E. Bompard and Y. Ma, “Modeling bilateral electricity markets: A complex network approach,” IEEE Transactions on Power Systems, vol. 23, no. 4, pp. 1590–1600, 2008. View at: Publisher Site  Google Scholar
 W. Lee, L. Xiang, R. Schober, and V. W. S. Wong, “Direct electricity trading in smart grid: A coalitional game analysis,” IEEE Journal on Selected Areas in Communications, vol. 32, no. 7, pp. 1398–1411, 2014. View at: Publisher Site  Google Scholar
 S. M. Sajjadi, P. Mandal, T.L. Tseng, and M. VelezReyes, “Transactive energy market in distribution systems: A case study of energy trading between transactive nodes,” in Proceedings of the 48th North American Power Symposium, NAPS 2016, pp. 1–6, Springer International Publishing, USA, September 2016. View at: Publisher Site  Google Scholar
 A. Jokic, P. P. J. Van Den Bosch, and R. M. Hermans, “Distributed, pricebased control approach to marketbased operation of future power systems,” in Proceedings of the 2009 6th International Conference on the European Energy Market, EEM 2009, pp. 1–6, Leuven, Belgium, May 2009. View at: Publisher Site  Google Scholar
 C. Zhang, Q. Wang, J. Wang, P. Pinson, J. M. Morales, and J. Ostergaard, “Realtime procurement strategies of a proactive distribution company with aggregatorbased demand response,” IEEE Transactions on Smart Grid, vol. 3053, no. c, p. 1, 2016. View at: Publisher Site  Google Scholar
 T. Ishii and K. Yasuda, “Hierarchical decentralized autonomous control in superdistributed energy systems,” IEEJ Transactions on Electrical and Electronic Engineering, vol. 2, no. 1, pp. 63–71, 2007. View at: Publisher Site  Google Scholar
 B. Canizes, T. Pinto, J. Soares, Z. Vale, P. Chamoso, and D. Santos, “Smart city: A GECADBISITE energy management case study,” Advances in Intelligent Systems and Computing, vol. 619, pp. 92–100, 2017. View at: Publisher Site  Google Scholar
 L. David and B. Kurt, A Guide to Monte Carlo Simulations in Statistical Physics, Springer International Publishing, New York, NY, USA, 2005. View at: MathSciNet
 V. E. Brkzin, “The Onset of Phase Transitions in Finite Systems,” Physikertagung Heidelberg, vol. 42, no. 7, pp. 6–8, 1986. View at: Publisher Site  Google Scholar
 S. Bornholdt and F. Wagner, “Stability of money: Phase transitions in an Ising economy,” Physica A: Statistical Mechanics and its Applications, vol. 316, no. 14, pp. 453–468, 2002. View at: Publisher Site  Google Scholar
 D. O. Cajueiro, “Enforcing social behavior in an Ising model with complex neighborhoods,” Physica A: Statistical Mechanics and its Applications, vol. 390, no. 9, pp. 1695–1703, 2011. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2018 Francisco PrietoCastrillo et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.