Complexity / 2018 / Article
Special Issue

Complex Optimization and Simulation in Power Systems

View this Special Issue

Research Article | Open Access

Volume 2018 |Article ID 5905932 |

Francisco Prieto-Castrillo, Amin Shokri Gazafroudi, Javier Prieto, Juan Manuel Corchado, "An Ising Spin-Based Model to Explore Efficient Flexibility in Distributed Power Systems", Complexity, vol. 2018, Article ID 5905932, 16 pages, 2018.

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

Academic Editor: João Soares
Received25 Jan 2018
Accepted19 Mar 2018
Published08 May 2018


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.

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-called herding behaviour or the economic bubbles are phenomena lying outside analyses that neglect large correlations and self-organisation occurring near the critical point [7]. In regard to DPSs, a sharp transition in customer’s behaviour (i.e., emergence) can trigger dramatic consequences. To the best of our knowledge however, a quantitative analysis of the complexity of customers behaviour in DPSs has not yet been provided. Ising-based models are a promising approach as we demonstrate here.

The Lenz-Ising spin model is known since 1925 [8]. In short, the model consists of an arrangement of interacting agents which have two possible states (i.e., up and down). Agents interact with their neighbours locally and are also exposed to external action and to thermal noise. In the physics metaphor of a ferromagnetic material, the atoms tend to align their quantum spins to minimise energy. However, for high temperatures, the agents flip their spins randomly. There is a critical temperature (known as Curie temperature) for which the system suffers a sharp change—second-order phase transition—from order to disorder. Near the system has anomalous behaviour and the correlations among spin states propagate fast to the entire system. As we will show in Section 2, this model is fairly simple but powerful enough to capture most of the features of complex systems (e.g., phase transitions and universality) [9]. This has motivated many researchers to apply the Ising model to different fields ranging from ecology [10] to language evolution [11] (see [12] for a discussion of the physical motivation of the model and its limitations).

Given that social systems are both finite and heterogeneous, the Ising model has been exported in different flavours to capture the phenomena under investigation. Perhaps the most known of these Ising-like models is Tom Schelling’s model of segregation [13] which has been shown to provide insight into the mechanisms under segregation in U.S. cities. This model has been shown to roughly correspond to an Ising model at [14]. In this regard, the temperature has been understood as a proxy of tolerance in binary thermodynamic societies; solubility corresponds to integration (i.e., mixing) and the miscibility gap to segregation [15]. Ising inspired systems have also been applied to financial markets to explain expectation bubbles and economic crashes [16]. For instance, the authors in [17] create a synthetic market where agents can take three actions: buy, sell, or stay inactive. Also, researchers in the power systems and electricity markets domain have made an effort to include part of this complex behaviour in the problem. For instance, in [18] bilateral electricity markets are analysed as complex networks. In [19], the authors use a game theoretic approach to understand the cooperation between small-scale electricity production and consumers. Multiagent systems (MAS) have also been applied to electricity markets to allow decentralized 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 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 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 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 real-time interplay structure among five kinds of actors: consumers, bus-loads, aggregators, DSOs, and RTEM. In Figure 1, we describe schematically the power exchange problem and how it couples with the spin model. 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 bus-loads and aggregators is two-way according to the flexibility behaviour from the demand side. Moreover, there are two-way power exchanges between aggregator and DSO. However, the power transaction between DSO and bus-loads is one-way: only from DSO to bus-loads. In other words, bus-loads can only buy real-time power from the DSO, while they are able to buy/sell power from/to aggregators. Also, only DSO can participate in the RTEM.

Each bus-load connects a set of customers whose state in terms of demand corresponds to a binary state of a spin (here represented as an arrow up/down mimicking the spin metaphor). From the demand side, different types of power balancing or other constraints can force cooperation among customers either at interbus level (blue arrows) or at intrabus level (pink arrow) in Figure 1(b). In this work, we use a 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 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 bus-loads there exists some notion of topological closeness (e.g., spatial or electrical proximity) which is partially mapped onto the square lattice.

Bus [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 5

More rigorously, given a set of consumers , a set of electrical buses , and a set of potential aggregators , we arrange two partitions: in-bus-communities and bus-aggregations . 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: customer-bus map: and bus-aggregator 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. Bus-Loads and Consumers

Each consumer has an associated load at time . Since each bus connects customers, the load at each bus at time is . Additionally, each consumer load can be split into a scheduled amount and a flexible portion , 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 day-ahead scheduled electrical demand in the real time. However, if his real-time electrical demand is more than his day-ahead 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 per-bus 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 bus-aggregator 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 self-sustainability and dynamic flexibility of all loads. In this work, we allow two types of flexibility constraint:

On one hand, (7) is the definition of the shiftable-loads (i.e., loads that can be shifted over time). On the other hand, (6) increases the self-sustainability of the distributed power system and converges the problem to the optimum social welfare. As we will see both constraints lead to different scenarios.

The optimisation at each bus can be expressed as the trade-off between load bought from DSO at price 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 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 , and related to the NordPool market. We summarise the respective values in Tables 2 and 3.

Time (h) [kW] [€/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,15

Time (h) [€/kW] [€/kW] [€/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,07

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 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:Therefore, the DSOs optimisation function can be expressed as the accumulated trade-off among these quantities, times the respective prices over time:This function will be minimised from the DSO’s perspective to maximise its profit.

So far we have described the main actors in the power exchange scenario with power balance and price constraints along with the respective optimisation functions for each agent. Now we describe the Ising model and our interpretation of its constituents in this context.

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 spin-to-spin 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 spin-to-spin 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 so-called Mean Field approximation its value is .

One way to implement the Ising model numerically is through the Metropolis-Hastings algorithm [25], which belongs to the family of the Markov Chain Monte Carlo (MCMC) methods. Applied to our case it can be understood as a random walk over the configuration space that converges to the Boltzmann distribution. In Algorithm 1 we show the pseudocode of a slight variant of the classical algorithm where we have included the possibility for implementing a constraint at each step. Here can represent the constraints in (6) and (7). The idea is that the spin shift is performed also when the constraint is minimised in absolute value. The number of iterations is chosen so that the final configuration reaches equilibrium. In operative terms this means that the correlations among the spin states are negligible at this point.

INPUT: , , , ,
while do
pick a random spin from
change if spin
ann’s probability for
configuration if spin
magnetization in new
if then
   random number from uniform distribution
    if ( or ) then
   end if
end if
end while

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, spin-to-spin 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 spin-to-spin interaction can be a proxy for communication among customers and the external field can simulate a top-down directive forcing customers to follow some policies or Demand Response (DR) programs that are adopted by players in the top layers of the system (e.g., DSO and aggregators). On the other hand, in the context of financial markets [27] interprets as a measure for the deviation from the fundamental value; if is large agents are afraid of trading unless they have strong support from their private opinions or from their neighbours. In our case 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 per-bus 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 dimension-free variables as follows: and . 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. . 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 bus-loads) in the distributed power system. The global optimisation function can be expressed asNotice 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 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, in-bus-communities (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. Shiftable-Loads

The constraint in (29) can be implemented as follows. First we solve the optimisation problem in (29). The output is the per-bus magnetisations . Notice that the constraint can be factorised in . This way, minimising 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 bus-flexibility solution follows the RTEM prices over time. This means that at demand peak times (10–13 h) and from (18–21 h) it is necessary to increase customers’ flexibility.

Then we might ask how the spin system can comply with the total magnetisation imposed by the electrical constraints at every time (see Figure 3). We need each spin-community to follow the constraint:This must hold for each bus-community of size for all times. By redefining the spin states as 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 so-called herding behaviour by imposing a top-down approach. Notice how the external action ramps up customers’ flexibility. However, if temperature is high, the noise will destroy these dynamics and the intended constraint can not be reached.

4.2. Self-Sustainability

Now let us explore a different situation where loads are forced to adjust so that at every time. A simple way to tackle this problem is by seeking solutions where the total magnetisation is either maximised or minimised. Hence, the problem can be set in terms of an Integer Linear Programming (ILP) problem of the form:Analogous for in the shiftable-loads case, now the optimal solutions do not depend on time and we can set . For each case (max or min) we have the range of values for . 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.

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 bus-community (notice that minimising instead of maximising the problem will provide the same result; the only difference among solutions is the switching of states in the bus communities when is an even number).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 constraints-free evolution renders values of 3 orders of magnitude larger than those obtained with the constrained version (CMH). The results of this experiment are shown in Figure 6. Here we compare on a normalised scale how the constraints deviate from 0 as we perturb the feasible spin solution with increasing perturbation strength. We measure the perturbation by the Hamming distance between the original and the perturbed configuration. Every point represents the average value of 20 Monte Carlo replicas for the same parameters. As expected, the unconstrained evolution fails to provide a systematic trend in these dynamics since once the constraint is broken, there are no mechanisms to bend the spin dynamics towards regions with increasing feasibility. However, we find a pattern for the constrained evolution: as we deviate from the solution the value for increases monotonically (worsening the feasibility). Notice that if perturbations are large, the minimum is only found for high temperatures, since as we increase thermal fluctuations, the system will be able to explore different configurations and discover lower values.

The effect of the external field is different in both cases. In the unconstrained case the field only tilts up the values and reduces the noise. This is useless in our case since the field action reinforces the perturbations worsening the constraint. In the constrained case, however, the field linearises the point pattern on our scaling (we have added a dotted grey line for reference) and this improves feasibility. Therefore, an external action in the community would not be useful for recovering flexibility unless additional mechanisms are implemented to penalise deviations from the feasible region. These actions can be implemented smoothly because the drift from optimality increases monotonically with the perturbation.

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 aggregator-to-DSO 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 bus-loads. Hence, optimising the problem from the aggregators’ point of view propels a bottom-up power transaction from bus-loads to aggregators, from aggregators to DSO, and from DSO to the real-time market, hierarchically. Moreover, the exchanged power benefit between the DSO and bus-loads is missing because this is irrelevant to aggregators. Consequently, as (36) shows, aggregators ask for the maximum loads’ flexibility. In other words, the power flexibility is balanced with the power exchanged between the bus-loads to aggregators and with the power sent from DSO to bus-loads. Hence, from the aggregators’ perspective, this will push the system to increase the flexibility from demand side and transacted power from the DSO to the loads. Notice that although this case is profitable for the DSO, it is not a profitable way for bus-loads to join this energy management approach.

One way to motivate bus-loads to join this setting is by letting be a parameter instead of a variable in this optimisation problem. In particular—and without loss of generality—we can set . 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 bus-loads, (i.e., there is no heterogeneity in the distributed power system), the optimisation function does not depend on the load-to-bus mapping . 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 , , and the optimisation function in (37) renders . 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 set 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 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 bus-to-aggregator mapping. By adding this constraint, the optimisation can be tackled by solving the BLP problem in the variables

In Figure 7 we show the normalised value of the objective function in (38) as we increase the number of aggregators (in %) and for different values. First we notice a global trend: the optimisation increases (lower objective function values) with the number of aggregators. Also helps in this reduction by pulling the solutions downwards and by smoothing the peaks as it increases. For low values the constraints are more stiff and the feasible space becomes fragmented. In the limit (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 bus-load collections they might impose certain boundaries 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 and then we measure how thermal noise worsens the solution.

In Figure 8 we show these results by comparing the boundary-free evolution with the aggregators’ boundary dynamics with the reference boundaries shown in Figure 1(d). As expected, with very high temperature, noise will dominate interactions and the system will reach the 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 free-neighbourhood 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 (spin-to-spin interactions) will not be strong enough to keep customers’ confidence level in the flexibility program. Consequently, they will randomly decide whether to increase or decrease their demand and the constraint for maximum flexibility will start to deteriorate. However, when the information flow among customers is bounded within regions in the configuration space, this drop is smoothed, since the information does not propagate to the entire system in the same way. From our previous discussion about the limitations of a finite model for reproducing a phase transition, one might attempt to provide a semicritical exponent to these phase-like transitions. However, that is left for a future work.

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); as grows shared opinions as to whether demand should be increased will not be strong enough to keep customers’ confidence level in the flexibility program. Consequently, they will randomly decide whether demand should be increased or decreased and the constraint for maximum flexibility will start to deteriorate. However, when the information flow among customers is bounded within regions in the configuration space, it is smoothed, since the information does not propagate to the entire system in the same way.

In a future work we will make a more in-depth analysis of the relationships between more general spin system topologies and the flexibility constraints. Also we will develop a more thorough study of the phase transitions in presence of aggregator boundaries and their implications.


:Number of customers
:Set of customers
:Set of bus-loads
: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
:Spin-to-spin 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
:Real-time 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
:Real-time market price for DSO and RTM power exchanges
:Price for DSO and load power exchanges
:Real-time market price for DSO and aggregator power exchanges
:Profit guarantee factor of aggregator at time
:Additional constraint imposed to the Metropolis-Hastings algorithm
:Optimisation function for bus-load
: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.


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 MSCA-RISE-2014: Marie Sklodowska-Curie project DREAM-GO; Enabling Demand Response for Short and Real-Time Efficient and Market Based Smart Operation; An Intelligent and Real-Time 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.


  1. N. Mithulananthan, Q. Duong Hung, and Y. Kwang, Intelligent Network Integration of Distributed Renewable Generation, Springer International Publishing, 2017.
  2. A. S. Gazafroudi, F. Prieto-Castrillo, T. Pinto, and J. M. Corchado, “Organization-based multi-agent system of local electricity market: Bottom-up approach,” Advances in Intelligent Systems and Computing, vol. 619, pp. 281–283, 2017. View at: Publisher Site | Google Scholar
  3. 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
  4. 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
  5. A. Gazafroudi, F. Prieto-Castrillo, 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
  6. 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
  7. 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
  8. S. G. Brush, “History of the Lenz-Ising model,” Reviews of Modern Physics, vol. 39, no. 4, pp. 883–893, 1967. View at: Publisher Site | Google Scholar
  9. R. V. Sole, Phase Transitions, Springer International Publishing, 2011.
  10. 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
  11. D. Nettle, “Is the rate of linguistic change constant?” Lingua, vol. 108, no. 2-3, pp. 119–136, 1999. View at: Publisher Site | Google Scholar
  12. A. M. Guenault, Statistical Physics, Springer International Publishing, Dordrecht, Netherlands, 2007. View at: Publisher Site
  13. T. C. Schelling, “Dynamic models of segregation,” Journal of Mathematical Sociology, vol. 1, pp. 143–186, 1971. View at: Publisher Site | Google Scholar
  14. D. Stauffer, “Social applications of two-dimensional Ising models,” American Journal of Physics, vol. 76, no. 4-5, pp. 470–472, 2008. View at: Publisher Site | Google Scholar
  15. 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
  16. 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
  17. 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
  18. 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
  19. 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
  20. S. M. Sajjadi, P. Mandal, T.-L. Tseng, and M. Velez-Reyes, “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
  21. A. Jokic, P. P. J. Van Den Bosch, and R. M. Hermans, “Distributed, price-based control approach to market-based 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
  22. C. Zhang, Q. Wang, J. Wang, P. Pinson, J. M. Morales, and J. Ostergaard, “Real-time procurement strategies of a proactive distribution company with aggregator-based demand response,” IEEE Transactions on Smart Grid, vol. 3053, no. c, p. 1, 2016. View at: Publisher Site | Google Scholar
  23. T. Ishii and K. Yasuda, “Hierarchical decentralized autonomous control in super-distributed energy systems,” IEEJ Transactions on Electrical and Electronic Engineering, vol. 2, no. 1, pp. 63–71, 2007. View at: Publisher Site | Google Scholar
  24. B. Canizes, T. Pinto, J. Soares, Z. Vale, P. Chamoso, and D. Santos, “Smart city: A GECAD-BISITE energy management case study,” Advances in Intelligent Systems and Computing, vol. 619, pp. 92–100, 2017. View at: Publisher Site | Google Scholar
  25. 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
  26. 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
  27. S. Bornholdt and F. Wagner, “Stability of money: Phase transitions in an Ising economy,” Physica A: Statistical Mechanics and its Applications, vol. 316, no. 1-4, pp. 453–468, 2002. View at: Publisher Site | Google Scholar
  28. 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 © 2018 Francisco Prieto-Castrillo 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.

Related articles

No related content is available yet for this article.
 PDF Download Citation Citation
 Download other formatsMore
 Order printed copiesOrder

Related articles

No related content is available yet for this article.

Article of the Year Award: Outstanding research contributions of 2021, as selected by our Chief Editors. Read the winning articles.