Nonlinear Dynamics in Decentralized Electrical Markets

(1)

MSc Computational Science

Master Thesis

Nonlinear Dynamics in Decentralized

Electrical Markets

by

Nikita Granger

November 22, 2018

Academic Supervisor

and Examiner:

Prof. Dr. B.D.

Kandhai

Second Assessor:

Dr. R. Quax

(2)

Abstract

Due to the similarities between decentralized electrical markets and stock markets, we adapted the Brock Hommes 1998 model to create an agent based model for smart grids that has a lower dimensionality and lower computational complexity than traditional agent based models. Using the wealth of techniques that exist for studying non-linear systems, we investigated drivers of instability in decentralized electrical markets, and managed to find it for certain prices of market

information, and certain belief types. We demonstrated that batteries should not be centralized, but rather distributed evenly over users to reduce volatility in the presence of large supply shocks, and that the same holds true for the number of agents present in the market. Finally we also showed that batteries cannot entirely eliminate seasonality due to opportunity cost, but this

(3)

1 Introduction

In 2007 the European Union declared its goal to achieve "20-20-20" [20]. This means a 20% reduction in carbon emissions, an increase in the production from renewable electricity sources (RES) to 20% of all electricity consumed, and to reduce electricity consumption by 20% by the year 2020 compared to 1990. Most RES will be either wind or solar power [1], and since production from this source of electricity is outside the control of its owners a redesign of electricity markets and transmission grid is required in order to ensure stable and reliable delivery of electricity. For this reason the European Union set up a Smart Grid Task Force [1] in 2009 to help facilitate the transition to the smart grid. According to the goals of this task force the following changes will need to be implemented in order to transition to the smart grid:

• All households and businesses can install windmills and solar panels on their property. This way market participants can be both consumers and producers of electricity.

• Market participants will be able to store electricity. This will accommodate variable produc-tion from RES.

• Market participants will be able to directly buy and sell to other market participants in their neighborhood. In this way markets will be distributed instead of centralized. This will also lead to the creation of micro markets based on neighborhoods that could have its own price level, separate from the wholesale market.

• Nominations and metering will be automated. This, coupled with the ability to store elec-tricity, will allow real time trading of electricity.

• Demand will respond to changes in supply instead of being completely inelastic.

Similar initiatives are currently taking place in other parts of the world, including the United States of America, Brazil, India, China, and South Korea [21]. On top of facilitating the introduc-tion of RES, smart grids will also facilitate a free market for electricity by allowing more flexible trading in real time by streamlining the nomination of transactions. Market efficiency should also theoretically improve with the introduction of more competition and a breaking up of the oligopoly that exists in centralized markets.

These differences will mean that today’s electrical market will have to be completely overhauled, and due to its real time trading, the ability to buy and hold, and large number of heterogeneous market participants, the new electrical market will be much closer in operation to a stock market. For this reason, we believe that employing a traditional agent based model, similar to the ones we have seen used in wholesale markets, and a few models that currently exist for decentralized markets, will encounter many problems, including:

1. High computational requirements due to the large number of agents. Where previously a market would contain 30-40 agents, it would now have to model the interactions of millions of agents to capture the market realistically.

2. More complicated decentralized electrical markets lead to higher level of dimensionality, since there will be a larger quantity and variety of agents. This makes traditional agent based models more difficult to study, since there will be a larger number of parameters and dimensions to consider while running simulations.

For these reasons, and due to the similarities between a decentralized electrical market and a stock market, we believe it is appropriate to borrow from the wealth of techniques that are used to model nonlinear dynamics in stock markets. The main approach we will consider in our paper was developed by Brock and Hommes in 1998 [22], and considers the movement in asset prices to be a consequence of heterogeneous beliefs about the future of the market. The advantage of this approach is that we end up with a system of difference equations that represent the aggregate belief of the market, rather than that of individual agents, reducing the dimensionality of the model significantly compared to traditional agent based models. This reduces the number of possible parameters that influence behavior in our model, which makes it easier to investigate various phenomena in the market. Secondly, aggregating beliefs in this way also reduces computational complexity significantly, which can be a problem in large agent based models. Finally, there is a rich body of literature regarding how to study the stability of such a system.

(5)

Using this approach we will investigate break points and transitions in market behavior occurs, with the goal of giving policy makers advice and a better understanding of what factors cause volatility in decentralized electrical markets. Some more concrete questions include how does the distribution of batteries impact the stability of the market, and how many agents should there be in a single market for effective trading?

We will begin this paper with a literature review that will discuss the research that has already been done in this field and in fields relevant to the construction of our model. Next we will describe the methodology and the construction of our model. We will then finish off with a series of experiments that will explore nonlinear phenomena and break points in decentralized electrical markets.

2 Literature Review

There are many different aspects of the smart grid that could be modeled: demand response (e.g. for electric vehicles (EVs), household appliances), production of virtual power plants and other distributed power generation facilities, neighborhood trading, interaction with the wholesale market and controlling battery storage just to name a few. Our model will assume that there is no wholesale market, and thus all trading is done between neighbors. For this reason in the literature review we will mostly focus on research that addresses neighborhood trading and we will only touch briefly on other aspects of smart grid markets. Additionally, we will also discuss some other papers that served as inspiration in the construction of our model.

2.1 Agent Based Models for Wholesale Markets

As a result of the small number of market participants, agent based modeling has become a very popular approach to modeling contemporary, centralized wholesale markets. A wide variety of different agent based models exist for this, such as PowerACE [1], AMES [2], MASCEM [3], and EMCAS [4] to name some of the more prominent models from Europe and America. Since none of these are decentralized market models, they are generally not of direct interest to us, but they do represent the principle inspiration for the agent based models that do exist now in smart grids, so a general overview on how these models work will be discussed. All these models consist of many different interacting agents. Common examples of agents that are included in a typical model are the transmission system operator, the market operator, retailers, wholesalers, larger consumers, market makers, and regulators. Market participants, such as retailers, wholesalers, large consumers, market makers, all make bids and asks in the futures market, and the market operator sets the price based on the bids and the asks. The transmission system operator and the regulator decide on how the electricity will be transported across the grid. Much of the work done on smart grid market models is simply an extension of these original models. These extensions usually take the form of adding agents that represent an aggregate set of consumers or producers which individually have smart grid characterstics, such as demand responsiveness. Examples include Ensslen et al. 2014 [5], who use an agent in PowerACE to represent the aggregated behavior of all electrical vehicle owners, and Gottwalt et al. 2011 [6] who do the same thing except with households. Alternatively, there exists literature that assumes that smart grid agents interact exclusively with retailers, instead of the wholesale market directly. An example of this includes Zhou et al. 2011 [7] who models the demand responsiveness of commercial buildings, but assumes a simple function for the cost of electricity where the price changes depending on the quantity of electricity consumed.

2.2 Agent Based Models for Smart Grids

One of the earliest papers on agent based modeling for distributed electrical markets [9], the Electricity Market Multi-agent System (EMMAS) model proposed by Gnansounou et al. provides a very rough sketch of what elements should be taken into account while constructing an agent based model for a smart grid market, but it does not delve into the many details that would need to be sorted out in order to actually construct a model. According to the EMMAS model, there should exist many different types of agents, including consumers, generators, transmission system operators, and regulators to name a few. Each agent possesses different assets, such as transmission cables, and generators, as well as possessing many different behaviors, such as the ability to buy and sell. The paper offers an entity relationship-style schema to explain what assets and what

(6)

behaviors belong to which agent. Trading and price discovery is done as it would be in a regular market, through agents making passive bids and asks, and aggressors buying or selling according to the bids and asks on the market. Agents should also be able to forecast the price, and make bids or asks based on these forecasts. The forecasting mechanism for this is not discussed in detail in the paper. The proposed EMMAS model is very reminiscent of the traditional wholesale market models we discussed earlier, and in many ways is a bridge between traditional electrical market models and smart grid market models.

Kahrobaee et al. 2013 [10] and 2014 [11], are some of the only general agent based models for neighborhood trading, and hence will be discussed in greater depth. Both papers take an entirely agent-based approach to the modeling of smart grid energy markets, where the agents use a decision tree to decide on how much to buy and sell. The decision to buy or sell and to charge or discharge their battery is determined by 1) whether there is a deficit or surplus of electricity produced in a given hour for a particular agent, and 2) whether the deficit/surplus is less than or greater than the average deficit/surplus historically for each agent. These two statistics are then put into the flow chart seen in Figure1 to determine the actions of the agents. In both papers the price for electricity (electricity purchase rate) is determined by a fixed formula that takes the following form

EP R(t) = α1eα2l(t)+ α3eα4l(t), (1)

where the α’s are constants, and l(t) is the load (or quantity demanded) at time t. The excess electricity produced is sold at a different price (electricity selling rate), which is given as

ESR(t) = βEP R(t), (2)

which is the EPR, but adjusted by some constant β, which takes a value between zero and one. This means that the selling rate is typically smaller than the purchasing rate. Neither paper discusses why these formulas are chosen.

The principle difference between the 2013 and 2014 papers is the introduction of different neigh-borhoods that each consist of a variety of different agents. There are three different agent types in the 2014 model: residential, commercial and industrial. First each agent declares its surplus/deficit and interacts with its neighbors, and after distributing electricity across the neighborhood, each neighborhood trades with the other neighborhoods, distributing electricity evenly across all of them. Any excess or shortage of electricity is then addressed by a centralized utility.

The conclusions in both papers are that shifting the electricity consumption results in lower prices for electricity for all consumers. The introduction of neighborhood trading also leads to lower prices for all consumers compared to buying directly from a retailer.

(7)

Figure 1: Decision Tree for Agents. def (tj) is the deficit/surplus for participant j at time t. UMR

and OMR are the percentages of the surplus sold off or the deficit bought versus taken or stored in the battery. The figure is taken from Kahrobaee et al. 2014 [11].

2.3 Models for Physical Smart Grids

GRIDLAB-D [12], is a simulation environment that was commissioned by the United States De-partment of Energy to study smart grids. The creators of this model make very few assumptions of the structure of a smart grid market, and instead just leave it to the user of the model, with only the physics of a grid being implemented. As a result, many of the examples given of projects that use this model focus more on engineering aspects of smart grids. This model can, however, be used as a foundation for an economic model by adding trading households or neighbors on top of it, and is our motivation for mentioning it here.

A more recently proposed model is Kolen et al. 2018 [19], which is a collaborative effort between E.ON and RWTH Aachen University. This model is an extension of the GRIDLAB-D model, but improves on the computational performance of GRIDLAB-D by distributing the calculations over multiple nodes rather than having all calculations done on a single processor using shared memory. This model relies heavily on work already done on RepastHPC in order to function, which is a distributed version of Repast, a common modeling toolkit.

(8)

2.4 Options Pricing for Storage in Energy Markets

The introduction of the smart grid means the introduction of new technologies, such as batteries. This will allow, in contrast to now, the ability to store electricity on a large scale. It is important to understand the consequences of storage in a commodity market, and the closest analogue to the introduction of batteries in electrical markets is the presence of storage in natural gas markets. De Jong 2016 [14], gives an overview on gas storage and the value it brings to a market. In this paper they identified three different ways to measure value for gas storage: intrinsic value, rolling intrinsic, and optional value. Intrinsic value is equivalent to the spread in the futures contracts for a particular time period. Rolling intrinsic is also the spread in futures contracts, but also with the consideration that these spreads might change before maturity. The rolling intrinsic can be thought of as an expectation of spreads in futures contracts. Finally there is optional value. A storage unit, such as a gas reservoir or a battery, can be thought of as a swing option. A swing option is like an American or Bermudan option, but can be exercised multiple times, up until we reach the capacity of the battery. The maturity of this option is usually determined by the storage amounts required by the owner. If the owner of the storage is required to have a certain amount of the commodity in storage for delivery at a particular time, then that particular time can be considered the time of maturity. The actual process of pricing a swing option is given in an earlier paper by Boogert and De Jong 2008 [13], where they determined the price of gas storage by using least squares Monte Carlo, and then dynamically optimizing over functions.

Swing options, which are the basis for gas storage valuation, currently exist in electrical markets, and are exclusively sold as an OTC derivative, mainly for hedging purposes. A popular model for the pricing of electrical swing options in the context of wholesale futures power market is outlined in Haarbrucker and Kuhn 2009 [16], which use stochastic programming and control to price their options. Due to batteries still being expensive and impractical to use in a real world context, we were not able to find any literature on employing swing options to price battery storage in electrical markets. Another issue hampering the study of batteries in the smart grid is the lack of stochastic models used to describe price dynamics that would allow us to price the batteries.

2.5 Break Points and Tipping Points in Electrical Markets

One of the goals of our research is to find tipping points and what we will call break points in smart grid markets. Break points are similar in nature to tipping points, but are more qualitatively defined. The only paper we are aware of that investigated and discovered the same phenomena in a smart grid context is the Dallinger 2011 paper [8], where by adding an aggregated agent that represented all electric vehicles into the PowerACE model, Dallinger managed to demonstrate that there could be what he called ’avalanche effects,’ or break points in the load (demand) for electricity, when enough of the car charging agents assume that the price will increase in the following period. It is important to note that the PowerACE model is a model commissioned by the European Union to study the wholesale, centralized, electric market. This means that the observed break points did not take place in the context of a decentralized smart grid market, but rather a traditional wholesale market that exists today but with some elements of a smart grid market. All agents make a forecast of the market, and place their bids into the futures market according to what they think the price will be. Individual agents use reinforcement learning to determine which forecast procedure leads to the most accurate prediction.

2.6 Brock Hommes 1998

Since our model is largely based on the Brock and Hommes model proposed in 1998, we will describe this model in depth. The general idea behind this model is that dynamics in stock prices are driven by the changing beliefs of agents. We model the percentage of agents with different beliefs over time, and how those beliefs change as one belief becomes a better predictor of the market than the others.

The Brock Hommes model begins by first considering a world with one risky asset, and one risk-free asset. The wealth of an investor is determined by a portfolio of these two assets

Wt+1= RWt+ (pt+1+ yt+1− Rpt)zt, (3)

where the term Wtis the wealth at time t, R is the risk-free rate of return (plus one), pt is the

(9)

risky assets at time t. Variables that are in bold are stochastic in nature, and since our dynamic processes are adapted, everything after time t is unknown, and hence stochastic. We define wealth in this way to model our risky asset. Since investors can make a guaranteed profit investing in a risk-free asset, they will only invest in the risky asset if they are properly compensated for their risk. Using this notation, we clearly see the division between the risk-free return that we are guaranteed, and the excess return above the risk-free rate, which compensates us for taking on risk. There are two ways we are compensated for taking risk: through an increase in the price of the stock, which is represented by pt+1 and through dividends (profits earned by the company)

which is represented by yt+1. Using this definition of wealth we can determine the conditional

expectation of wealth at time t + 1 given that our filtrationFt(i.e. the time that has already been

observed by agents) is up until time t to be

E(Wt+1|Ft) = E(RWt+ (pt+1+ yt+1− Rpt)zt|Ft) = RWt+ E(pt+1+ yt+1|Ft)zt− Rptzt, (4)

we can also deduce the conditional variance in the same manner

V (Wt+1|Ft) = V (RWt+ (pt+1+ yt+1− Rpt)zt|Ft) = V (pt+1+ yt+1− Rptzt|Ft)z2t, (5)

and in the Brock Hommes model they assume for simplicity that the conditional variance of the excess returns is set always a constant for all agents, so that V (p_t+1+ y_t+1− Rptzt|Ft) ≡ σ2. From

here we can then construct a utility function that is essentially a trade-off between the conditional risk and conditional return of our total wealth. This is given as

Ut= E(Wt+1|Ft) −

a

2V (Wt+1|Ft), (6)

where a is a constant we can set to determine the risk aversion of a particular agent, and we will assume it is the same across all agents. When we determine Equation6using the results we got from Equations4 and5, we get

Ut= RWt+ E(pt+1+ yt+1− Rpt|Ft)zt−

a 2σ

2_z2

t. (7)

This utility curve now gives us a trade-off between the expected return of the asset, and the riskiness of the asset. When we optimize over this function we end up with the optimal amount of the stock we should purchase to maximize our risk adjusted returns. We can also control the role risk plays by adjusting the aversion parameter a. If a = 0, then we say our agent is risk neutral, and risk plays no role in his actions. If a > 0, then we say the agent is risk averse, and actively avoids risk. Finally if a < 0, then our agent is risk loving, which means he wants greater risk exposure. Since this utility function is a upside down facing parabola when our agents are risk averse (a > 0), we can find the optimal trade-off between risk and return by simply taking a derivative with respect to zt, equating our expression to zero and isolating the ztterm

dUt

dzt

= E(p_t+1+ y_t+1− Rpt|Ft) − aσ2zt= 0 ⇔ zt=

E(p_t+1+ y_t+1− Rpt|Ft)

aσ2 . (8)

Equation8is now our demand curve at time t as it is a downward sloping function that depends on price. Since different agents have different opinions on the conditional expectation of the market, we can add an additional condition to find the demand curve to be

zht=

E(p_t+1+ y_t+1− Rpt|Ft, h)

aσ2 , (9)

where h the type of the agent (i.e. the type of beliefs he holds). We can then determine the aggregate demand of the market by adding together all demand curves for all types of agents multiplied by the percentage of each type of agent. This is the same as taking the weighted average of all demand curves. The price can be determined by setting aggregate supply equal to aggregate demand

X

h∈types

(10)

where nhtis the proportion of people with demand curve zht, and zst is the aggregate supply. The

logic behind this derivation of the equilibrium price is that when we find the price that leads to supply equaling demand, there is neither a surplus nor a deficit of stocks on the market, hence the price is optimal. Any other price will lead to a shortage, or a surplus of stocks on the market. If expectations are rational, then this equilibrium price is equivalent to the fundamental price. This is the first of three approaches we have for finding the equilibrium price, but they will only yield the same result if the expectations of the future are rational. We will use different approaches depending on what is most convenient. The second way to derive the equilibrium price is by first assuming that there is only one belief in the market

E(p_t+1+ y_t+1− Rpt|Ft)

aσ2 − zst= 0 ⇔ Rpt= E(pt+1+ yt+1) − aσ 2_z

st. (11)

The logic behind this assumption is that we assume that on average people are correct. Next we simplify this expression further by assuming there are no new shares issued into this market, so that zst= 0, which gives us

Rpt= E(pt+1+ yt+1|Ft). (12)

Next we can assume that the dividend process yields E(yt+1|Ft) = ¯y, which is a constant. We can

also set our expected price equal to the fundamental price pt+1= ¯p, and with this we can derive

the following

R ¯p = ¯p + ¯y, (13)

and the price can be determined to be ¯p = _R−1¯y , which is the fundamental market price. This pricing approach can be thought of as a no arbitrage argument, or in other words that the returns from both the risk-free and risky assets should, on average, be the same, or else you can make a profit going long in one and short selling the other. This definition will be used to determine what the deviation of the market price is from the fundamental price for the risky asset, since it is the simplest formulation of the fundamental price.

Going back to Equation 12, we can also reintroduce the multiple different agents here to get the third derivation for the equilibrium price, which is as follows

Rpt=

X

h∈types

nhtE(pt+1+ yt+1|Ft, h). (14)

This is the same no arbitrage argument as before, but now we take the average expected price across all beliefs. Again, if all expectations are rational this becomes equal to the fundamental price.

2.6.1 Deviation from Fundamental Price

Using Equation13we now know the fundamental market price. Combining this with the result in Equation14, we can now express the market as a deviation from the fundamental price. There are two advantages of doing this, firstly we know when the market enters a bull or a bear market with respect to the fundamental price, and secondly we can use this to enforce stability in our market. We will now continue to work with the deviation from the market price

xt= pt− ¯pt. (15)

This deviation from the fundamental price means that the expected price can be decomposed into the following

E(p_t+1+ y_t+1− Rpt|Ft, h) = E(¯pt+1+ yt+1|Ft) + f (xt−1, xt−2, ..., xt−L|h), (16)

where E(¯p_t+1+y_t+1|Ft) is the conditional expected fundamental price filtered to time t, and is the

same for all traders. We have seen a particular case of this conditional expectation in Equation13. Function f is a deterministic function of deviations from the fundamental price, and this differs by trader depending on their beliefs. Using the same approach as in Equation14, the market clearing deviation from the fundamental price can be given as

(11)

Rxt=

X

h∈types

nhtf (xt−1, xt−2, ..., xt−L|h), (17)

where the fundamental price is ignored because it persists regardless of the beliefs of traders, so the clearing only applies to the deviation from the fundamental price.

The mean-variance optimization that was done in Equation8, can also be done here in exactly the same manner as before, with the only difference being that the deviation from the fundamental price is used rather than the fundamental price

max z E(pt+1+ yt+1− Rpt|Ft, h)z − a 2z 2_{V (p} t+1+ yt+1− Rpt|Ft, h) = max z ρhtz − a 2z 2_σ2_, ₍₁₈₎

where we set E(pt+1+ yt+1− Rpt|Ft, h) =zt and V (pt+1+ yt+1− Rpt|Ft, h) = σ2. Since the

optimization in Equation18 is the same as the one we had for Equation 8 except for a constant which drops out anyway when we take the derivative, the demand curve remains the same when we are working with excess returns alone. In the Brock Hommes model they denote the demand curve for Equation18 as z(ρht), where ρht is the expectation for the excess returns by an agent

with type h belief at time t.

2.6.2 Changing Beliefs in the Brock Hommes Model

Agents change their beliefs based on which forecast type has the best performance. Performance is defined as a function of excess returns, which will be expressed in the following way

pt+1+ yt+1− Rpt= xt+1+ p∗t+1+ yt+1− Rxt− Rp∗t =

= xt+1− Rxt+ p∗t+1+ yt+1− E(p∗t+1+ yt+1|Ft) + E(p∗t+1+ yt+1|Ft) − Rp∗t, (19)

where pt+1 = p ∗t+1+xt+1 means that the price is equal to the fundamental price p∗ plus the

deviation from the fundamental price x. Since we know that E(p∗t+1+ yt+1|Ft) = Rp∗t, the last

two terms drop out. The remaining p∗t+1+ yt+1− E(p∗t+1+ yt+1|Ft) is simply the excess returns.

Brock and Hommes assumed that the excess returns process has a martingale term and represent it by δt+1. So this gives us the following simplified expression for excess returns

xt+1− Rxt+ δt+1. (20)

Equation20can be thought of as the summation of the excess returns due to belief driven prices plus the efficient market excess return (which is stochastic, with an average of zero). Now this expression can be used to determine the performance of a specific type of belief h at time t

πht= (xt+1− Rxt+ δt+1)z(ρht), (21)

which can be thought of as the excess returns multiplied by the quantity demanded. Next a memory can be introduced to performance so that agents consider more than just one period’s performance when changing their belief type. It can take the form

Uh,t= πh,t+ ηUh,t−1, (22)

where Uh,t is the utility of agent with belief type h at time t. The η represents the memory

strength, and can be adjusted to determine whether agents are more or less backward looking. Finally Brock and Hommes construct the percentage of agents with a certain belief by using the following logistic function to construct a probability

nh,t=

exp[βUh,t−1]

P

h∈typesexp[βUh,t−1]

, (23)

where nh,tis the percentage of agents with belief type h at time t. β is a parameter that determines

how likely agents are to change their beliefs. A logistic function is a common way to create a probability out of a function. In this case we make the utility associated with one belief type an exponent of e in the numerator, and divide it by the summation of all utilities that are made exponents of e in the same way. In this way we will always have a function between zero and one, and hence we will always have a probability.

(12)

2.6.3 Bifurcations

After trying this model with several different belief types, Brock and Hommes manage to discover that trend chasers (trend followers) can lead to the creation of pitchfork bifurcations, contrarians (people who go against the trend) can lead to doubling bifurcations and opposite biases (people who always predict a constant deviation from the fundamental price) can lead to Hopf bifurcations. Examples of what these price movements look like can be seen in Figures2and3. All these models assume that the fundamentalist belief type is always present in the market.

Figure 2: Price and belief movements for trend chasers and contrarians. X is the deviation from the fundamental price, and m is the percentage of agents that are fundamentalists. Movements take place over time. These figures are taken from Brock Hommes 1998 [22].

(a) Trend chaser movements in price and belief. (b) Contrarian movements in price and belief.

Figure 3: Price and belief movements for opposite biases. X is the price deviation from the fundamental price, n1 is the proportion of fundamentalists, n2 is the proportion of agents with positive biases, and n3 is the proportion with negative biases. Movements take place over time. This figure is taken from Brock Hommes 1998 [22].

2.6.4 Stability Analysis

Besides running the simulations over time there are also other approaches we have of analyzing our dynamic system and its stability. A common approach to studying dynamic systems is to determine which values for parameters result in a system becoming stable, and which lead to divergence. This is done by taking an existing collection of functions and taking the Jacobian with respect to every variable. We can then take the eigenvalues of the resulting square Jacobian matrix. The real part of these eigenvalues is what explains the dynamics of a system, with the real part of the eigenvalues having to be inside the unit circle in order for the system to be stable. How the eigenvalues change

(13)

inside the unit circle determine how the qualitative behavior of the dynamic system changes. A qualtitative change in the behavior of a dynamic system is called a bifurcation. Some examples of bifurcations include

1. If the eigenvalue is +1, then we have a saddle node bifurcation. This is when we have a pair of steady states where one is stable and the other is a saddle. It is possible for there to be more than one steady state in this case. An example of this can be seen in Figure4.

2. If the eigenvalue is -1, then we have a period doubling or flip bifurcation. This results in a cycle being created. An example of this can be seen in Figure5.

3. A pair of complex eigenvalues on the unit circle. This results in Hopf or Neimark-Sacker bifurcation. This leads to potentially complicated invariant circle with quasi-period cycles. An example of a Hopf bifurcation can be seen in Figure6.

Figure 4: This is a saddle bifurcation for the system. Taken from [23].

Figure 5: This is a flip or doubling bifurcation. Taken from [23].

Figure 6: This is a Hopf bifurcation. Taken from [23].

Outside of the three different bifurcations we can also study other bifurcations with the help of a bifurcation diagram. The idea behind a bifurcation diagram is that we can create a graphical

(14)

representation of all equilibria that a dynamic system has. An example of a bifurcation diagram can be seen in Figure7, where the darker areas represent equilibria for different sets of parameters x and r.

Figure 7: An example of a bifurcation diagram from a Lotka-Volterra model. x can be thought of as the population size, while the r parameter is the rate of reproduction. This figure is taken from [23].

Deeper investigation into bifurcations can be analytically difficult, but numerically simple by plotting the resulting dynamic system with the help of a phase portrait. A phase portrait shows the derivatives of a certain dynamical system for different values of that dynamical system. This gives an idea about which direction the general system is moving in. An example of this can be seen in Figure8.

(15)

Figure 8: Example of a phase portrait for a pendulum. The red arrows represent the derivative for that point in space, while the blue lines represent the path taken by the pendulum in space. This is taken from [23].

Another numerical technique that can be used to study the average rate of divergence for a particular state of our dynamic system is by looking at the largest Lyapunov characteristic exponent (LCE), sometimes called the Lyapunov principle exponent (LPE), which is defined as follows

λmax= lim n→∞

1

nlog(||(DxF

n_)δx||), ₍₂₄₎

where Fnis a n dimensional function, and DxFnis the Jacobian of the nth variable of the function

at point x, and δx is a small perturbation for the vector, and || • || denotes the length or the norm of the vector. We can define the stability of our dynamic system by looking at the n different Lyapunov exponents where λ1> λ2> ... > λn. Some possible things that can be measured using

LCEs include

1. If all LCEs are negative, then we have a steady state or a stable cycle.

2. If the largest LCE or LPE is equal to zero, and the remaining are negative, then we have a circular strange attractor.

3. If the largest LCE or LPE is positive, then we have a situation where our attractor is depen-dent on the initial conditions. Under certain conditions our dynamic system can diverge.

Finally in situations where we have strange attractors, we can end up with fractal-type shapes. There are ways of studying these as well using techniques such as box counting. As the name implies, the idea behind box counting is that we see how many n-dimensional boxes with side length we need to draw around a shape in order to cover it entirely. The number of boxes we need to achieve this is called the box counting dimension, and it can be expressed with the following equation

Db= lim →0

log(N ())

log(1) , (25)

where N () is the number of n-dimensional boxes with side needed to cover the entire area of the shape. An example of box counting in fractal geometry can be seen in Figure9. In practice what we will do is first generate a large number of points to draw our strange attractor. Next we will

(16)

calculate N () for different values of , and regress the two against each other in a log(N ())−log(1) plot to find a relationship between the two. From here we calculate relationship for a very small value of to get the box counting dimension Db.

Figure 9: Example of box counting to count the fractal dimension of the von Koch curve. Taken from [24].

3 Methodology

3.1 Differences between the Smart Grid Market and Stock Markets

Before we get into how we adapted the Brock Hommes 1998 model, we will first dicuss the differ-ences between electrical markets and stock marktes. Traditionally electrical markets were modeled using a completely separate set of techniques from other financial markets. Some reason for this include:

1. Electricity is immediately perishable. The biggest consequence of this is that electricity prices fluctuate relatively predictably throughout the day, and hence modeling electrical markets using statistical or stochastic models can be problematic. For this reason most electricity market models are agent based.

2. Due to the fact that transmission system operators (TSOs) need to keep the electrical grid balanced, it is important for them to know exactly all the electricity bought and sold ahead of time in order to be able to balance the grid. This means that all trading is done on the futures market, and not in real time.

3. Electricity markets are generally driven by market fundamentals, and not speculation. Due to restrictions on market participation, there are few speculators and poorly informed market participants.

4. Electricity is consumed and produced continuously, whereas stocks are not.

Moving to the smart grid, however, will change many of these characteristics. The introduction of batteries will mean that electricity will no longer be immediately perishable. This means that the battery can behave as sort of a ’portfolio’ for electricity ownership, and that it will be possible to buy and hold as one would a stock or a bond. In the smart grid metering and balancing will be handled automatically, meaning that trading could feasibly be done in real time, or close to real time, again, making the electrical market more like a stock market. Finally since all households and businesses will be able to participate directly instead of acting through retailers. Anyone could be a market participant, much like a stock or a bond market. These similarities were the inspiration to try to adapt some agent based models that were initially developed for the stock market into the smart grid market. As mentioned before, the model we chose to adapt for use in the smart grid market is the Brock Hommes model proposed in 1998. The reason we chose this model in particular is due to its relative simplicity, and because of the ease with which we could study stability.

Finally the last change that will be made will be the new value created by the addition of batteries. Traditionally the only value to electricity at any given period of time is the associated

(17)

supply and demand at that point in time. The change over to smart grids will mean the introduction of batteries to electrical markets. This creates an interesting new dynamic to electricity markets, which is the intrinsic value and the real option value of storage. The intrinsic value is the spread in the futures contracts of different maturities. For example in Europe electricity is used for heating in countries like France, and as a result tends to be more expensive during the winter than in the summer. That means the value of a battery in France is the spread in prices between the summer and winter months. The owner of the battery can lock in this spread using futures contracts, and reap the profit immediately. The second form of value is the real optionality of the storage. Going back to our example of a battery in France, between the time we go long in the summer and short the winter, we can still trade between those two time periods, as long as at the end of the time period we have enough electricity to honor our promise to deliver the amount stipulated in the futures contract. In summary, the intrinsic value of our storage stems from predictable seasonal trends in electricity markets, while the real option value stems from less predictable short term shocks to the market. We will try to capture both dynamics in our model as they relate to the price of electricity.

3.2 Adapting the Brock Hommes Model to Smart Grid Markets

As mentioned before, one big difference between stock markets and electricity markets is that demand for electricity is generally inelastic. This is due to the fact that market participants have a certain minimum amount of electricity they need to consume regardless of the price. In addition to this minimum agents need to consume, they can also speculate on the future price of electricity using their battery as storage. Borrowing concepts from gas storage markets, agents can buy electricity either on the futures market or on the spot market. The futures market is used to lock in predictable seasonal spreads in the electricity price, while the spot market is used to sort out any sudden, unforeseen changes in the electricity supply or demand in the short term. The demand curve for each agent is given as the following

Z_it= Z_icons+ Z_i,tf utures+ Z_i,tspot, (26) where Zit is the quantity demanded by agent i at time t. Zicons is completely constant for each

agent i, Z_i,tf utures is the quantity demanded of electricity bought on the futures market for agent i at time t, and finally Z_i,tspot is the demand for electricity in the spot market. In this paper we assume that the quantity demanded in the futures market is known beforehand, and therefore is constant for each period t. For the demand in the spot market, however, we will introduce a function. Using the same technique as in Equation8in the Brock Hommes model, we can find the demand function that optimizes utility. The utility curve we will use in for the spot demand will be as follows

U_i,tspot= Z_i,tspot pe_i,t+1− Rpi,t − a(Zi,tspot)

2_V, ₍₂₇₎

where Z_i,tspotis the demanded by agent i at time t, U_i,tspotis the utility of agent i at time t, V is the perceived risk which we assume is the same across all agents, and R is the risk free rate (plus one). Equation27follows the same logic as the utility curve in the Brock Hommes model, where utility is expressed as a trade-off between risk and reward. pe

i,t+1 and pt are the expected price at time

t + 1 for agent i and the price at time t respectively, where ptis found through our market clearing

procedure. We assume that agents are myopic, meaning that they will only consider forecasts for the next period only when considering whether to buy or sell on the spot market now. In the future an extension could be adding more than one step ahead, and then dynamically optimizing over the forecasted prices for all those time periods instead of just one. We assume that storing electricity is a form of investment, that on average should earn the risk free rate, so pe_i,t+1−(1+r)pi,tis basically

the excess returns over the risk free rate. Since, as before, our function U_i,tstorage is parabolic and continuous, and hence is convex with a single maximum, we can find the optimal value by taking the first derivative and setting it equal to zero, which gives us the following result

∂U_i,tspot

∂Z_i,tspot = 0 ⇔ p

e

i,t+1− Rpi,t− 2a(Zi,tspot)V = 0 ⇔ Z spot i,t =

pe

i,t+1− Rpt

2aV , (28)

and that gives us the optimal quantity demanded with respect to any price at time t for agent i. By plugging this back into our full demand curve in Equation26, we get the following

(18)

Zit= Zicons+ Z f utures i,t + pe i,t+1− Rpt 2aV . (29)

Next we will incorporate the constraints introduced by batteries. Batteries allow agents to trade using electricity from their storage, but only until while the capacity of the battery permits. Batteries cannot surpass their maximum capacity or store less than a zero volume of electricity. We can introduce these restrictions to the battery by using a piecewise function

¯ Zit=     

Z_icons+ Z_if utures+ bcap_i − Si,t−1 Zit− Z cons i + Si,t−1≥ bcapi Zt i 0 ≤ Zit− Zicons+ Si,t−1≤ b cap i

Z_icons+ Z_if utures− Si,t−1 Zit− Z cons

i + Si,t−1≤ 0.

So we have a downward sloping linear demand curve, and a limit when we reach the lower and upper capacities of the batteries. The ¯Z_itterm is the demand with the battery constraints imposed, Z_itis the demand from Equation 29, and the other Z terms are the quantity demanded either for consumption, or as a speculative futures or spot market purchase by agent i at time t. bcap_i is the battery capacity of agent i, and Si,t is the amount of electricity stored by agent i at time t.

We can also rewrite the above piecewise function in an aggregated form to benefit computational efficiency by removing the need to sum across all agents individually. We can find the aggregate demand by adding together all the individual demand curves

X i Zit= X i Zicons+ Z f utures i,t + pe i,t+1− Rpt 2aV , (30) Zt= Zcons+ Ztf utures+ P ip e i,t+1− Rpt 2aV , (31)

and then we can take the replace the individual expected excess return with an average across all agents for simplicity to get pe

i,t+1− Rpt= pet+1− Rpt, where pet is the weighted average forecast

across all agents. This gives us the following aggregate demand

Zt= Zcons+ Ztf utures+

n(pe

t+1− Rpt)

2aV . (32)

We can then reintroduce the constraints the same way we did before, but this time without the subscript i denoting the number of the agent

¯ Zt=     

Zcons_{+ Z}f utures_{+ b}cap_{− S}

t−1 Zt− Zcons+ St−1≥ bcap

Zt _{0 ≤ Z}t_{− Z}cons_{+ S}

t−1≤ bcap

Zcons_{+ Z}f utures_{− S}

t−1 Zt− Zcons+ St−1≤ 0,

where Zt _{comes from Equation} ₃₂_{, and ¯}_Zt _{is the piecewise aggergate demand.} _{The other Z}

terms are the aggregate amount consumed and purchased on the futures market, while bcap _{is the}

aggregate battery capacity, and St−1is the aggregate amount stored at time t − 1. An important

consequence of our aggregation is that our aggregate demand is now entirely linear in nature. When we sum across each demand curve individually, it is possible to end up with a non-linear demand curve if agents have different slopes to their demand, as well as different battery capacities. In this sense, the aggregation represents a linear approximation of the non-aggregated model.

So now we have our demand curve, both in individual and aggregate form. In order to find the market price we also need a supply curve. It is unreasonable to assume that supply will always equal zero as was done in Brock Hommes, so we need to consider a more elaborate supply function. Since we assume all our production is generated by renewable sources our supply is outside the control of the agents, so it is completely price inelastic. We also assume supply will always be sufficient to satisfy the quantity demanded by agents. We will vary the supply function we use in different experiments, but the general form they will take will appear as follows

Gt=

X

i

Z_icons+ t, (33)

where tis an error and E(t) = 0. With this supply function we can now find the market clearing

price by setting each part of our piece wise demand function equal to supply, starting with the middle part of the demand function, which is when the storage is not a constraint for demand

(19)

X

i

Z_icons+ t=

X

i

Z_icons+ Z_i,tf utures+p

e i,t+1− Rpt 2aV ⇔ X i t= X i Z_i,tf utures+p e i,t+1− Rpt 2aV , (34)

and then we set Z_i,tf uturesequal to ctand we setP_i to get n, and we can then isolate the price

ptto get the following

n − ct= X i pe_i,t+1− Rpt 2aV , (35) (n − ct)(2aV ) = X i pe_i,t+1− Rpt, (36) (n − ct)(2aV ) = X i pe_i,t+1− nRpt, (37) 1 R (ct− n)(2aV ) n + X i pe i,t+1 n = pt, (38)

then assume pt = p∗, or in other words that the market price is correct. We also take the

expectations E(pe i,t+1) = p∗ or P ip e i,t n = p

∗ _{so that on average expectations are correct, then we}

get p∗= 1 R (ct− n)(2aV ) n + p ∗_, ₍₃₉₎ p∗− 1 Rp ∗₌ 1 R (ct− n)(2aV ) n , (40) r Rp ∗₌(ct− n)(2aV ) nR , (41) p∗= (ct− n)(2aV ) nr , (42)

as the analytical solution to the equilibrium price. Since we assume that on average a household’s production of electricity is always enough to satisfy supply, the price of electricity is on average always zero. In Equation42 we see that the more agents we have in our market, and the higher the risk free rate, the less volatile the price becomes. As we would expect, Volatility and errors made in the futures market (i.e. when the futures market does not correctly predict seasonality) increase price volatility.

When the constraints are active we assume that one of four things can happen:

1. The batteries are full and the production is just enough to satisfy consumption.

2. The batteries are empty and the production is just enough to satisfy consumption.

3. The batteries are full and the production is more than is required to satisfy consumption.

4. The batteries are empty and the production is less than is required to satisfy consumption.

In the first two scenarios the price of electricity is equal to zero. In this situation every household can satisfy its own demand using their own store electricity, and hence no trading is required. For the third and fourth situations we assume that there is a market failure. Situation three results in prices diverging to negative infinity, while in situation four the prices diverge to positive infinity. With these options the fundamental prices become

(20)

p∗=                ∞ Z_it− Zcons i + Si,t−1≥ bcapi 0 Zt i− Zicons+ Si,t−1= b cap i 2aV (ct− ¯t) nr b cap i > Z t i− Z cons i + Si,t−1> 0, 0 Zt i− Zicons+ Si,t−1= 0 −∞ Z_it− Zcons i + Si,t−1≤ 0.

Or the subscripts i can be dropped off when you are working with aggregates.

3.3 Adding the Opportunity Cost of the Battery

In the first model the battery capacity does not serve any purpose outside of defining the boundaries where the market can exist. Storage facilities in other commodity markets, such as oil and gas [17], have a considerable impact on the price and this dynamic should be introduced into our model. Additionally, since our agents manage their own batteries, they have to take into account the opportunity cost of their actions. For example, if they sell off everything they have stored, they will not have any electricity in the event of an even more acute shortage later on. Likewise, if they fill their battery to full capacity, and prices become negative, they lose the opportunity to make money from charging their battery. This opportunity cost should play a role in the agent’s calculations. We will begin with the same demand function we had in the previous model, which is given by Zit= Zi,tcons+ Z f utures i,t + Z spot i,t , (43)

but has a different utility function for the spot market, which will now become

U_ispot= −γ 2(Si,t+1− bcap_i 2 − ∂V (t + 1, pe i,t+1, Si,t+1) ∂Si,t+1

)2+ Z_i,tspot(pe_i,t+1− Rpt) − a(Zi,tspot)

2_{V, (44)}

where Si,t is the storage at time t for agent i, γ is a constant that determines how important the

opportunity cost of the battery is relative to the excess returns, and bcap_i is the battery capacity for agent i. This new −γ₂(Si,t+1−

bcap_i 2 −

∂V (t+1,pe_i,t+1,Si,t+1)

∂Si,t+1 )

2 _{term is introduced to influence the}

price of electricity when the battery deviates too much from half way full. The derivative term

∂V (t+1,pe_i,t+1,Si,t+1)

∂Si,t+1 is present in order to move the optimal fullness of the battery from half full

when the situation calls for it. For example, when the price is negative, it would make sense to have the batteries fuller than usual, which is captured in this partial derivative.

Next we can simplify Equation44 by first removing the St+1 term from Equation 44, which we

can write as St+1= St+ Z spot i,t + Z f utures i,t − Z cons

i and get the following

U_ispot = −γ 2(St+Z spot i,t +Z f utures i,t −Z cons i − bcapacity 2 − ∂V (t + 1, Pt+1, St+1) ∂St+1

)2+Z_i,tspot(pe_i,t+1−Rpt)−a(Z spot i,t )

2_V.

(45) Next we will discuss how to approximate the ∂V (t+1,p

e

i,t+1,Si,t+1)

∂Si,t+1 term by pricing the battery as

an option, and agents will strive to optimize the value of their battery as well.

3.4 Deriving the Derivative of the Value of Storage

In the previous derivation we stated the need to find a function to represent ∂V_∂S, where V is the value of the storage, and S is the amount stored. We will approximate the function numerically, using the technique described in Boogert 2008 [13], which was already discussed in the literature review, but will be discussed in greater depth here.

3.4.1 Stochastic Difference Equation

Initially we require a stochastic difference equation (since we are in discrete time) from which we will perform Monte Carlo simulations. The stochastic difference equation will model the spot market and will have the following characteristics:

(21)

1. We will assume that on average supply and demand will be correctly resolved in the futures market, and hence there will be no trading required in the spot market. That means that the price should be, on average, zero. If the market deviates from zero, it will eventually revert back to it.

2. There will be a small amount of noise in the market.

3. There will occasionally be large jumps in prices.

4. Since the generation of electricity is entirely generated by renewable sources of electricity like solar and wind power, we will assume that the stochastic difference equation will have auto correlation.

These four characteristics can then be captured in the following difference equation

∆Pt= −aPt+ ψ(σZ) + (1 − ψ)Pt−1+ IpJ, (46)

where Ptis the price at time t, ψ is a constant to determine the strength of the auto correlation,

a is the rate of mean reversion, σ is the variance, Z is a standard normal distribution, I is an indicator function that takes one with a probability p, and J is a jump, that we will assume is uniformally distributed.

Using Equation46we can then run Monte Carlo simulations that we will assume are different possible paths the market could take. We can then optimize over these paths.

3.4.2 Objective Function

Since pricing storage is fundamentally a dynamic optimization problem, we will need an objective function. The objective function will be taken directly from Boogert and de Jong 2008 [13], and will be as follows

π(t, Pt, vt) = arg min ∆v∈D(t,vt)

h(Pt, ∆vt) + C(t, Pt, vt+ ∆vt), (47)

where ∆vt is the change in the volume stored between time periods t and t − 1, h(Pt, ∆vt) is a

function that represents the payoff in period t, and has the following specification

h(Pt, ∆v) =      −Pt∆vt Charge on period t 0 Do nothing Pt∆vt Discharge on period t.

C(t, Pt, vt+ ∆v) is the discounted continuation function, and will be approximated by a linear

regression for every period, but the final period, similar to what is done by Longstaff-Schwarz for American options [15]. We will discuss the specification of the regression a little later. The final period will have its continuation function set to zero.

Finally we have the function D(t, v(t)) which is the constraints imposed on the charging and discharging. In this case the only constraints we have on charging and discharging is the capacity of the battery, which is given by the following

D(t, vt) = {∆vt|vmin(t + 1) ≤ vt+ ∆vt≤ vmax(t + 1)}. (48)

Additional constraints can be imposed as well for the rate of charging and discharging, but we will assume that the time frame we are considering is sufficient for this not to be a constraint.

This means that our final payoff will be as follows

Y (t + 1, Pt+1, vt+1) = h(Pt+1, π(t + 1, Pt+1, vt+1)) + e−δY (t + 2, Pt+2, vt+2+ π(t + 1, Pt+1, vt+1)),

(49) and we recursively walk backwards in time until we find the payoff Y (0, P0, v0), and we take the

average across all the different simulated paths for the prices to find the average value of our storage.

(22)

3.4.3 Regression Specification

Since we are already using simple market beliefs for our agents, we will also use simple forecasts for our regressions. These will be as follows

Y (t, Pt, vt) = β0+ β1Pt+ β2vt, (50)

or in other words, just an ordinary least squares regression on the price and the volume. The βs are our parameters that will be estimated during the regression. This gives us the following pay off function

Y (t + 1, Pt+1, vt+1) = h(Pt+1, π(t + 1, Pt+1, vt+1)) + e−δ(β0+ β1Pt+ β2(vt+2− π(t + 1, Pt+1, vt+1).

(51)

3.5 Approximating the Pay-off and Greeks

The pay-off and derivative functions depend on the parameters, and for the purposes of illustrating our technique for approximating the pay-off and Greeks we will use parameters given in Table1.

Table 1: Parameters used to calculate the value of the option.

Parameter Value

Time Steps 10

Risk-free Rate (annualized) 0.05

Battery Capacity 100

Volatility 3

ψ (from SDE) 0.3

α (from SDE) 10

Jump Size (from SDE) ±100

We will approximate the pay-off and the derivative of the pay-off with respect to the initial storage using a regression. Although there are many different parameters we could include in our regression, we will only alter the starting price and the initial storage because these two are the most relevant for our model. The impact on the value of storage for different values of initial price and initial storage can be seen in Figure10. All values found using this approach have a confidence interval of less than 0.01.

(23)

We will then perform the following regression to get an approximation of the pay-off function

V = β0+ β1S + β2P + β3P S, (52)

where V is the value of the storage, S is the initial stored amount, and P is the initial price. The coefficients are -198.705 for β0, -50.616 for β2 and 0.999 for β3. β1 is insignificant at a 95%

confidence interval, so we will remove that variable in future calculations. The R2_{is 0.99998, which}

is a very tight fit suggesting that this is the data generating process.

Using a finite difference we can approximate the value of ∂V_∂S, and we can see the results for this derivative in Figure11for different values of starting storage and staring price. The approximation of our derivative becomes less accurate for larger values of the start price, especially once we exceed a starting price of 1500. Prior to the starting price of 1500, our confidence interval remains below one.

Figure 11: The derivative of the value of the battery with respect to the initial amount stored for different prices.

The derivative of the value of the storage with respect to the initial storage amount is given as

∂V

∂S = β0+ β1P, (53)

where the parameters are -0.013 for β0 and 4.9997 for β1, and β0 is not significantly different

from zero on a 95% confidence interval, while β1 is not significantly different from five on a 95%

confidence interval. For this reason we can drop off the β0term, and just keep the β1P term. The

R2 _{is 0.99968, a very tight fit, suggesting that this is likely the data generating process. Another}

observation is that the β1term is always the same, regardless of the battery capacity of the agent in

question. Having a larger battery just leads to a proportional increase in the value of the battery, and not a change in the sensitivty with respect to the initial storage.

3.5.1 Summary of Options Pricing Methodology

In summary we will do the following steps to calculate the pay-off and derivative of the pay-off with respect to different parameters (known in finance as Greek functions) for our options:

1. Simulate the stochastic difference equation given in Equation46.

2. Optimize the consumption in the final period. We assume that the continuation function here is zero, so this is just a static optimization.

3. Run a regression on the final period using the specification given by Equation50. This will be the continuation function. Then we perform the optimization using Equation51 to find the optimal value π.

(24)

4. We then perform this recursively until we get to time zero. We then take the average of all pay-offs at time zero to get the value of the battery.

5. We then repeat this process for different values of the starting price, and different values of the initial amount stored in our battery. We use finite differencing to approximate the Greeks. We then perform regressions52 and53to get an approximation of our pay-off and Greeks.

3.6 Continuing the Derivation of the Optimal Utility

Using β1,iPt+1 or β1,ipei,t+1 as an approximation for our derivative

∂V (t+1,pe

i,t+1,Si,t+1)

∂Si,t+1 . This gives

us the following expression

U_ispot = −γ 2(Si,t+ Z spot i,t + Z f utures i,t − Z cons i − bcap_i 2 − β1p e i,t+1) 2_{+ Z}spot i,t (p e

i,t+1− Rpi,t) − a(Z spot i,t )

2_V.

(54) Next we will find the value of Z_i,tspot which optimizes our utility by taking the derivative of U_ispot with respect to Z_i,tspot, which gives us the following result

− γSi,t− γZ spot i,t − γZ f utures i,t + γZ cons i + γ bcap_i 2 + γβ1p e i,t+1+ p e i,t+1− Rpt− 2aV Z spot i,t = 0, (55)

and when we isolate Z_i,tspot we get

Z_i,tspot= γ( bcap_i 2 − Si,t− Z f utures i,t + Z cons i ) + pei,t+1(γβ1+ 1) − Rpt 2aV + γ . (56)

The entire demand is then given by

Z_it= Z_icons+ Z_if utures+γ(

bcap_i

2 − Si,t− Z f utures

i,t + Zicons) + pei,t+1(γβ1+ 1) − Rpt

2aV + γ . (57)

Notice that if we set γ = 0 then we end up with exactly the same demand function we derived for the first model. Finally we introduce the battery constraints

¯ Z_it=      Zcons i + Z f utures i + b cap

i − Si,t−1 Zit− Zicons+ Si,t−1≥ b cap i Z_it 0 ≤ Z_it− Zcons i + Si,t−1≤ bcapi Zcons i + Z f utures

i − Si,t−1 Zit− Zicons+ Si,t−1≤ 0,

where ¯Zt

i is the piecewise individual demand, and Zitcomes from Equation57.

Again, we can rewrite the individual demand curve into an aggregated one to reduce computational complexity. We can calculate this by summing across all n agents

X i Z_it=X i Z_icons+ Z_if utures+γ( bcap_i 2 − Si,t− Z f utures

2aV + γ . (58)

By summing we drop the subscript i from the demand for consumption and the demand for the futures. We will also add together all battery capacities and storage and get the aggregate of both

Zt= Zcons+ Zf utures+γ( bcap 2 − St− Z f utures t + Zcons) + P ip e i,t+1(γβ1+ 1) − Rpt 2aV + γ . (59)

Finally we average of all values for pe

i,t+1(γβ1+1)−Rpt, which will be written as pet+1(γβ1+1)−Rpt,

where pe

t+1 is the average expected price across all agents at time t. This gives us the following

Zt= Zcons+ Zf utures+γ( bcap 2 − St− Z f utures t + Zcons) + n(pet+1(γβ1+ 1) − Rpt) 2aV + γ . (60)

(25)

3.6.1 Deriving the Market Price with Opportunity Cost of Battery

Next we can calculate the market clearing price by setting the demand equal to supply, the same method we used previously for our first model. We assume that the supply follows a similar function to before, that being consumption plus some disturbance term . The gives us the following equation X i Z_icons+ t= X i

Z_icons+ Z_i,tf utures+γ(

bcap_i 2 − Si,t− Z f utures i,t + Z cons i ) + p e i,t+1(γβ1+ 1) − Rpt 2aV + γ , (61) and as before we can remove the Z_icons terms by subtracting them from both sides to get

X i t= X i Z_i,tf utures+γ( bcap_i 2 − Si,t− Z f utures i,t + Z cons i ) + p e i,t+1(γβ1,i+ 1) − Rpt 2aV + γ . (62) Again, we setP iZ f utures

i,t equal to a ct, and we set Pit equal to n, where n is the number of

individual agents. This gives us

n = ct+ X i γ(b cap i 2 − Si,t− Z f utures

2aV + γ , (63) (n − ct)(2aV + γ) = X i γ(b cap i 2 − Si,t− Z f utures i,t + Z cons i ) + p e i,t+1(γβ1+ 1) − Rpt, (64) (n − ct)(2aV + γ) = X i γ(b cap i 2 − Si,t− Z f utures i,t + Z cons i ) + X i pe_i,t+1(γβ1+ 1) − nRpt, (65) (n − ct)(2aV + γ) − X i γ(b cap i 2 − Si,t− Z f utures i,t + Z cons i ) = X i pe_i,t+1(γβ1+ 1) − nRpt, (66) (n − ct)(2aV + γ) − X i γ(b cap i 2 − Si,t− Z f utures i,t + Z cons i ) − X i pe_i,t+1(γβ1+ 1) = −nRpt, (67) (ct− n)(2aV + γ) + X i γ(b cap i 2 − Si,t− Z f utures i,t + Z cons i ) + X i pei,t+1(γβ1+ 1) = nRpt, (68) (ct− n)(2aV + γ) +P_iγ( bcap_i 2 − Si,t− Z f utures i,t + Zicons) + P ip e i,t+1(γβ1+ 1) nR = pt. (69)

This is market price. Next we assume that the fraction

P

ip e i,t

n is on average correct, so it is on

average equal to the market clearing price:

P

ip e i,t

n = p

∗_{. We also assume that the market price p} t

is correct, and can be set to p∗ = pe. This gives us

(ct− n)(2aV + γ) +Piγ( bcap_i 2 − Si,t− Z f utures i,t + Zicons) nR = pt− P ip e i,t+1(γβ1+ 1) nR , (70) (ct− n)(2aV + γ) +Piγ( bcap_i 2 − Si,t− Z f utures i,t + Z cons i ) nR = pt− P ipei,t+1 n γβ1+ 1 R , (71)

(26)

and we set pt= p∗ and P ip e i,t n = p ∗ _{to get} (ct− n)(2aV + γ) +Piγ( bcap_i 2 − Si,t− Z f utures i,t + Zicons) nR = p ∗_{− p}∗γβ1+ 1 R , (72) (ct− n)(2aV + γ) +Piγ( bcap_i 2 − Si,t− Z f utures i,t + Z cons i ) nR = p ∗_{(1 −} γβ1+ 1 R ), (73)

we then isolate p∗ to find the fundamental market price

(ct− n)(2aV + γ) +Piγ( bcap_i 2 − Si,t− Z f utures i,t + Z cons i ) nr − γβ1 = p∗. (74)

When we introduce the constraints from the battery capacities, we end up with the following piece wise function

p∗=                  ∞ Z_it− Zcons i + Si,t−1≥ bcapi , 0 Zt i − Zicons+ Si,t−1= b cap i , (ct−n)(2aV +γ)+Piγ( bcap_i

2 −Si,t−Zi,tf utures+Z cons i ) nr−γβ1 b cap i > Z t i − Z cons i + Si,t−1> 0, 0 Zt i − Zicons+ Si,t−1= 0, −∞ Zt i − Zicons+ Si,t−1≤ 0.

This is the same result that we had in the first model if we set γ = 0.

For the sake of computational efficiency we can also use an aggregated version of this model, where we do not sum over all agents and instead just take the aggregated storage, consumption, battery capacity, and electricity bought in the futures market. If we do this, we get the following derivation for the market price

p∗= (ct− n)(2aV + γ) + γ( bcap 2 − St− Z f utures t + Zcons) nr − γβ1 , (75)

where we aggregate battery capacities and amounts stored for all market participants. When we reintroduce the battery constraints we get the following price

p∗=                ∞ Zt_{− Z}cons_{+ S} t−1≥ bcap, 0 Zt_{− Z}cons_{+ S} t−1= bcap (ct−n)(2aV +γ)+γ(bcap2 −St−Zf uturest +Z

cons₎ nr−γβ1 , b cap_{> Z}t_{− Z}cons_{+ S} t−1> 0, 0 Zt_{− Z}cons_{+ S} t−1= 0, −∞ Zt− Zcons_{+ S} t−1≤ 0.

3.6.2 Aggregate Supply for a Smartgrid Market

As mentioned earlier, we cannot make the same assumption as in Brock Hommes that supply will equal zero. For this reason we need to specify a supply curve for the production of electricity. Since we assume that supply is entirely generated by renewables, such as wind and solar power, we assume that it is completely inelastic in nature, since agents have no (or very little) control of their production. This means that production is either a constant, or a function of something that is not price. There are two main features we want to capture in our generation model, one being seasonality of production, the other being sudden shocks (jumps) in the production of electricity.

Seasonality is inherent in both wind and solar energy [18], as can be seen in Figure 12 taken from the Energy Information Administration [25]. We can approximate the annual seasonality of production by using a simple sinusoidal function, such as a sine or cosine function, as is given in Equation76 Gt= X i Z_icons+ γ sin(t n), (76)

Nonlinear Dynamics in Decentralized Electrical Markets

MSc Computational Science

Master Thesis