Optimising cross-border flows in liberalised electricity markets

(1)

2008:110 CIV

M A S T E R ' S T H E S I S

Optimising cross-border flows in liberalised electricity markets

Erik Svensson

Luleå University of Technology MSc Programmes in Engineering

Engineering Physics Department of Mathematics

(2)

Abstract

The electricity markets in Europe are growing together and therefore understanding of cross-border flows becomes more and more important. Due to Market Coupling the method for allocating cross-border transfer capacity is changing, from explicit to implicit auction. This results in strictly price driven cross-border flows. In this thesis a model for forecasting price driven cross- border flows is developed.

The problem of finding the optimal cross-border flows, is modelled as a linear optimisation problem. The flows are found by minimising production cost in all involved areas for 24 hours at a time. The optimisation is subjected to constraints that derive from technical restrictions on the interconnecting cables.

From simulation it is shown that the modelled flows are close to the real flows. The model is also compared with an existing model and the results are at least as good, sometimes even better. This is very good for a fully flexible model.

(3)

1 Introduction

This is a Master Thesis in Applied Mathematics. It is the final part of an MSc programme in Engineering Physics at Lule˚a University of Technology. The thesis deals with the development of a cross-border flow model for European electricity markets. All the work has been done at Vattenfall Trading Services GmbH in Hamburg, Germany.

Cross-border flows are in this thesis assumed to be price driven. This is a direct effect of implicit capacity auctions, explained in section 2.4.4. Due to this, the problem of finding the optimal cross-border flows can be formulated as a linear programming problem.

The thesis is organised in the following way. First an introduction to the electricity market is given. It is, in order to understand what later will be modelled, necessary to gain some basic insights in the mechanisms of the electricity market. In the following chapter, chapter 3, some selected theory on linear programming is presented. The actual modelling is documented in chapter 4 to chapter 6, where description, definition as well as vector and matrix representation of the problem are presented. The last two chapters deal with flow simulations and model conclusions.

(6)

2 Electricity market

Electricity can not be stored and must therefore be produced and consumed at the same time, put in other words, supply must meet demand at all times. This fact, together with liberalisation, leads to an interesting and dynamic market.

Liberalised electricity markets are characterised by free trading on power exchanges and consumers rights to choose their own supplier. Before deregulation of the markets, production, sales, transmission and distribution were managed by one unit. In this monopoly situation the efficiency was questioned and the concept of liberalisation was developed in order to counteract this inefficiency.

The information about the electricity market presented in this chapter comes from [4], [5],[6] and [7].

2.1 Networks

The physical transport of electricity takes place in a power transmission system.

This transmission system consists of different networks and works as the physical link between power plants and end users. The different networks can be divided into two principal levels, the transmission grid and the distribution grid. The transmission grid is used for high voltage long distance transmission, whereas the distribution grid is made for shorter distance transmission and distributing electricity to end users.

Different networks are owned by different network companies. The network company responsible for the transmission grid in a certain area is called the Transmission System Operator, or TSO. All grid operators must ensure that there is electrical balance at all times in their areas, i.e. that production and imports match consumption and exports. Imports and exports exist due to interconnections with other transmission grids, more about that follows in section 2.4.

2.2 Electricity trading

Electricity trading refers to buying and selling electricity on the wholesale market. The trading is carried out at power exchanges or Over-The-Counter (OTC) and can be done both physically and financially. Physical trading is followed by physical delivery, whereas financial trading refers to buying and selling of financial instruments.

2.2.1 Exchange trading

The market places for exchange trading are power exchanges. In Germany, the Leipzig based power exchange EEX consists of three different markets for electricity, two physical spot markets and one future market.

The day-ahead spot market is an electricity auction where contracts for all 24 hours of the following day are traded. Selling and buying bids from the market participants are collected and put together to corresponding hour specific supply and demand curves. The spot price and the traded volumes for each hour is then determined by the intersection between these supply and demand curves.

On the intra-day spot market contracts with delivery on the same or on the following day are traded. Sell and buy orders are updated continuously and

(7)

trades occur when two prices match. The intra-day market offers the participants a last minute opportunity to deal with shortage or excess of capacityand thus contributes to the minimisation of balance costs.

For future market contracts, i.e. contracts for future delivery periods, all trading is done financially which means that no physical delivery takes place.

Future contracts are financial derivatives whose values change in response to changes in underlying variables. The underlying variables in this case are spot market prices. Trading at the future market is done continuously, i.e. a trade occurs when a sell order matches a buy order.

A more detailed introduction to exchange trading can be found in [4]

2.2.2 OTC trading

At the OTC market, trading is done off the exchanges. An OTC trade is a bilateral agreement that is established directly with a counterparty or via an intermediary player, for example a bank or a broker. All different sorts of contracts can be traded on the OTC market, including both physical spot contracts and financial future contracts.

Establishing bilateral agreements between two parties within the same company can sometimes be referred to as trading on an internal market.

2.3 Actors on the electricity market

To increase the understanding of the electricity market the different market actors need to be identified. Figure 2.1 shows the most important actors and the relationship between them.

Transmission Grid Distribution Grid

Power Trading

Company Sales Company

Power Exchange OTC market Electricity

Producer

Electricity User

Physical Flow

Wholesale market Resale market

Monetary Flow

Figure 2.1: Relationships between the market actors

(8)

2.3.1 Electricity producers

The electricity producers are the owners of all the different power plants that generate electricity and put it into the power transmission system. Electricity producers sell their generated electricity on the wholesale market.

2.3.2 Electricity users

The electricity users are the consumers, also referred to as end users. They with- draw electricity from the grid after first establishing resale purchase agreements with electricity suppliers, normally sales companies. End users that consume vast amounts of electricity can also be active on the wholesale market, one example is the German national railway, Deutsche Bahn AG.

End user pay the network owners for the utilisation of the network. This access fee provides the user with the possibility to choose between various suppliers. Electricity users include everything from households to industries.

2.3.3 Power trading companies

Power trading companies are intermediary players, trading electricity on the wholesale market. Power trading companies also participate in cross border trading, discussed further in section 2.4

2.3.4 Sales companies

Sales companies as well as power trading companies are intermediary players.

The major difference is that sales companies in addition to participating in the wholesale market also sell electricity to end users on the resale market. When dealing with households and industries the electricity has to be split up into smaller portions.

2.3.5 TSOs

In Germany the TSOs are obliged to buy all available wind power within their own control areas. A predetermined amount of this electricity must physically be delivered to end users, the rest is sold on the power exchange. This means that german TSOs are active on the wholesale market.

2.4 Cross border flows

When two market areas are connected, electricity can be transported between them, thus opening the possibility for cross-border trading. The cables linking the different networks are called interconnectors. Market participants can use these cables to buy and sell electricity on the foreign side of the border, but in order to do so they need to have utilisation rights for the desired capacity. The capacity is allocated via either explicit or implicit auctions.

2.4.1 Market areas

A market area is a geographical zone where electricity can be transported freely, i.e. without considering network congestion. As a consequence of this, each market area forms a homogenous price zone.

(9)

Normally each country is treated as it’s own market area, but there are of course exceptions. Norway for example consists of three market areas, this is due to congestion in the Norwegian transmission grid. Electricity simply can’t be transported freely all the time and therefore Norway needs to be split, however Norway is controlled by only one TSO.

An interesting different example where four TSOs are joined in one market area is Germany. The four different German transmission grids have enough internal connections so that power is allowed to flow without risk of congestion, hence Germany becomes one market area.

2.4.2 Physical and financial flows

When studying cross border flows it is important to understand the difference between physical and financial flows. Financial flows are commercial schedules that are the result of cross border trading. The physical flows are the corresponding real flows that obey the Kirchhoff’s laws and can be measured in the interconnecting cables.

Physical flows always arise when there are production surpluses and a production deficits in linked networks. The laws of physics apply and therefor electricity will always try to even out the system by flowing along the paths that render least possible resistance. For the border between market area A and market area B this means that even though the scheduled financial flow is fixed for a certain hour the physical flow in the shared interconnector may differ. It could very well be that the electricity also uses other connected paths, e.g. a detour over market area C, in order to get from A to B in the ”cheapest”

possible way. See Figure 2.2.

+ -

A C

B

Scheduled financial flow

Physical cross border flow Physical flow

Figure 2.2: Difference between physical and financial flows

(10)

All physical cables have maximum capacities. The maximum capacities on the interconnecting cables invoke restrictions on the financial flows. It is up to the cable operators, i.e. the TSOs, to determine how the maximum financial flow capacities in either direction will look like. Physical cables also have ramping restrictions. The physical ramping determines the maximum flow change from hour to hour. Financial ramping derives as a consequence from physical ramping and is set by the TSOs.

In this thesis ”flows” will hereafter always refer to financial flows.

2.4.3 Explicit Auction

An explicit auction is when the transmission capacity on an interconnector is auctioned separately from the marketplaces where the electricity is auctioned.

This means that a market participant who wants to buy or sell power in a foreign market also must be sure to have the needed capacity on the corresponding interconnector reserved. The price paid for the capacity is of course a crucial part of cross border trading, the transmission fee reduces the price spread and thereby decreases the possible profit. When capacity is assured, participation at power exchanges and bilateral trades are possible. The allocated capacity can also be left unused.

Explicit auction can be considered an inefficient method of handling the capacity allocations because it sometimes renders ”wrong” directed flows, i.e.

flows from high price areas to low price areas. Explicit auctions are nevertheless very common on the borders across Europe.

2.4.4 Implicit Auction

A more efficient way of dealing with the allocation of interconnector capacity is the implicit auction. With an implicit auction the transmission capacity is included in the electricity auctions at the market place, hence allocation of capacity is dealt with implicitly.

Implicit auction is the common term for two different approaches, market splitting and market coupling.

Market splitting is a cross border trading method where the available capacities on the interconnectors between all involved market areas are utilised by an implicit auction. The electricity and capacity auctions are carried out by one sole power exchange.

The market splitting mechanism at first treats all involved regions as one big market area. Market participants send in their bids and offers on a 24-hour basis and the curves are matched. When there is no congestion the whole area gets the same price and the cross border flows are the flows that lead to the this balanced system. In case of congestion the homogeneous price area is split up. The regions affected by the congestion then become their own market areas with their own prices.

Today market splitting only exists in the nordic area and is carried out by the nordic power exchange Nord Pool.

Market coupling is also a cross border trading method involving an implicit transmission capacitiy auction. It is based on the same general principles as

(11)

market splitting, but with the major difference that it is organised by at least two separate power exchanges.

In a region where market coupling applies, the co-operation between the involved power exchanges ensures that during every hour of the day all available interconnector capacity is utilised in accordance with the markets.

Market coupling is today active in the trilateral coupling of Netherlands, Belgium and France. Other pending coupling projects exist between Denmark and Germany and also between Norway and Netherlands.

Both market splitting and market coupling are auction methods that assure electricity to flow in the ”right” direction. Here, the right direction means from the low price area to the high price area. Entirely price-driven cross border flows will later be used as a model assumption.

(12)

3 Linear Programming

In this chapter a brief introduction to linear programming is given. Two different approaches for solving linear optimisation problems are introduced, the simplex method and interior point methods. The Matlab solver linprog that will be used later on in this thesis primarily uses a variant of an interior point method. This justifies describing interior point methods with more details than the simplex.

Also an introduction to a special class of linear programmes, namely network problems, is provided.

3.1 Introduction

Linear programming is a form of mathematical optimisation. The linear objective function is either maximised or minimised within a linearly constrained space. A general linear function can be written

z = c0+ c1x1+ c2x2+ ... + cnxn, which in vector notation yields

z = c^Tx.

The vectors c and x are both of length n.

The constraints, either equality or inequality, are written a^T_ix = b_i, i = 1, 2, ..., m,

˜

a^T_ix ≤ ˜b_i, i = 1, 2, ..., ˜m,

where a_iand ˜a_iare vectors of length n and b_iand ˜b_iare scalars. In vector-matrix notation the constraints become

Ax = b, Ax ≤ ˜˜ b.

Here A and ˜A are matrices with dimensions m × n and ˜m × n respectively. The vectors b and ˜b are of the corresponding lengths m and ˜m.

A linear program can be represented in many different ways. One representation is called the standard form. A linear program on standard form is a minimisation problem subject to equality constraints and with non-negative lower bounds,

min c^Tx subject to

Ax = b x ≥ 0.

All linear programmes can be converted to standard form. A maximisation becomes a minimisation by multiplying the objective function with −1. For free variables and variables with lower bounds other than zero, suitable variable substitutions create equivalent zero lower bounds. If upper bounds exist they can be treated as general constraints, i.e. included in coefficient matrix A.

General constraints involving inequalities can be transformed into equations by adding so called slack variables.

(13)

A point that satisfies all the constraints is called a feasible point. The set of all feasible points is called the feasible region or the feasible set, denoted S.

Each feasible point x ∈ S that cannot be expressed in the form x = αy + (1 − α)z,

with y, z ∈ S, 0 < α < 1, and y, z 6= x, is an extreme point or vertex. Geo- metrically the vertexes can be interpreted as the corner points of the feasible region.

3.2 Simplex Method

The simplex method is one of the most used methods in linear programming. It was developed in the 1940’s by operations research professor George B. Dantzig.

Today a broad variety of altered and enhanced simplex methods exist, although in this section emphasis is laid upon giving the reader a rough idea of how the basic method works. Detailed descriptions on the simplex method can be found in [1] and [2].

The simplex method is an iterative method for solving problems on the standard form. In its search for an optimum it improves the objective function by moving from one basic feasible solution (extreme point) to another. In geometric terms this means that the simplex method moves from vertex to vertex along the edges of the feasible region. The iterative process ends when a vertex is found that has no surrounding vertexes that can improve the objective function further. The optimal solution is hereby said to be found. In linear programming the optimum is always found at a vertex. In special cases there might be an infinite number of optimal solutions that all lie on a certain border of the feasible region. When moving along the edges of the feasible region the objective function is not too difficult to calculate, this contributes to the simplicity of the simplex method. The iterative behaviour of the simplex method can be portrayed by Figure 3.1.

Optimal Solution

Starting Vertex

Feasible Region

Figure 3.1: Moving along the edges of the feasible region

(14)

3.3 Interior Point Methods

Interior-point methods form a group of optimisation algorithms where the points generated at each iteration lie in the interior of the feasible region. This is in contrast to the simplex method where the movement is along the boundary of the feasible region. In this section the primal-dual interior-point method is studied. A modification of this method, called the predictor-corrector method, will also be considered. For more detailed description please consult [1].

3.3.1 Primal-Dual Method

Consider a linear program in standard form (the primal problem) min c^Tx

subject to

Ax = b x ≥ 0.

(3.1)

Here A is assumed to be an m × n matrix of full row rank. The corresponding dual program, with added slacks s, can be written

max b^Ty subject to

A^Ty + s = c s ≥ 0.

(3.2)

A feasible solution is optimal to (3.1) and (3.2) if it also satisfies the complementary slackness condition

x^Ts = 0, which equivalently can be written

XSe = 0.

Here X = diag x is a diagonal matrix whose j-th diagonal term is the j-th element of x, equivalently S = diag s. The vector e = (1 · · · 1)^T, consisting only of ones, is a vector of the same length as x and s.

For each iteration of the primal-dual method the system

Ax = b

A^Ty + s = c

XSe = µe

x, s ≥ 0

(3.3)

is solved for some µ > 0. The value of µ is reduced and the process is repeated until convergence is achieved. A decreasing µ means a decreasing duality gap, c^Tx − b^Ty, and when the duality gap reaches zero the solution is optimal.

From (3.3) it can be shown that the duality gap must equal the complementary slackness condition,

c^Tx − b^Ty = x^Ts = nµ.

Here n represents the length of x and s.

To save computation cost the primal-dual equations (3.3) are not solved exactly at every iteration. The solutions to x, y and s are estimates and need to be improved in each iteration. Suppose that the current solution satisfies

(15)

Ax = b and A^Ty + s = c, but not quite XSe = µe. Then new estimates, x + ∆x, y + ∆y and s + ∆s, that are closer to satisfying XSe = µe can be found.

The requirements for primal and dual feasibility, i.e. satisfying the constraints of (3.1) and (3.2), together with the complementary slackness condition invoke the following restrictions on ∆x, ∆y and ∆s

A∆x = 0

A^T∆y + ∆s = 0

S∆x + X∆s = µe − XSe − ∆X∆Se.

(3.4)

Since the last equation in (3.4) involves the term ∆X∆Se it is not linear. How- ever this term is in comparison small and can therefor be ignored. Solving (3.4) for ∆x, ∆y and ∆s gives

∆y = −(ADA^T)⁻¹AS⁻¹v(µ)

∆s = −A^T∆y

∆x = S⁻¹v(µ) − D∆s,

(3.5)

where D ≡ S⁻¹X and v(µ) = µe − XSe. These formulas are not defined if there are elements of x and s that are zero.

The algorithm now looks as follows: Given strictly feasible initial estimates on the primal, dual and dual slack variables, the directions ∆x, ∆y and ∆s are computed in accordance to (3.5). New estimates, x + ∆x, y + ∆y and s + ∆s are defined, µ is reduced and (3.3) is solved. This is repeated until µ has been reduced to a sufficiently low value.

How the parameter µ is updated is of course of great importance. Principally it could be said that the bigger the decrease the better, meaning fewer iterations and faster convergence. However, with large decreases the new solutions run the risk of not being strictly feasible, i.e. failing to sastisfy x > 0 and s > 0. This is also the reason why µ = 0 isn’t chosen directly from the start. In order to keep x and s positive at all times a step length approach can be applied when defining the new estimates to the soultion.

3.3.2 Predictor-Corrector Method

A modification of the primal-dual method is the already mentioned predictor- corrector method. Instead of ignoring the nonlinear term ∆X∆Se in (3.4), this method incorporates it in a 2-step procedure for finding the directions.

In the first step, the predictor step, (3.4) is solved in a similar way as before, i.e. with the nonlinear term and the term µe ignored. The intermediate solutions to ∆x, ∆y and ∆s found here are used only to calculate an approximate value on ∆X∆Se.

In the second step, the corrector step, the approximate value of ∆X∆Se computed in the predictor step is used to substitute the nonlinear term in equation system (3.4). The system is now solved and new directions, ∆x, ∆y and

∆s, are found.

This approach might seem to double the computational cost of one iteration of the primal-dual method. Luckily this is not the case. The corrector step only slightly increases the cost because it uses the same factorisation of ADA^T as in the predictor step. The positive effect of the predictor-corrector method is that it often decreases the number iterations needed by the primal-dual method.

(16)

This method is also sometimes called Mehrotra’s predictor-corrector method, simply because Sanjay Mehrotra was the first to apply this technique in the context of interior point methods.

3.4 Network Problems

Network problems are linear optimisation problems defined on networks. They form an important special class of linear programming due to their many special properties. By exploiting these special properties the simplex method can be implemented more efficiently. A network consists of two types of objects: nodes and arcs. The nodes are numbered from 1 to m, where m is the total number of nodes in the network. The arcs connect the nodes and the number of arcs in the network is n. An arc connecting two nodes, say i and j, is denoted (i, j).

The arcs are assumed to be directed, this means that the arc (i, j) is not the same as (j, i). A simple network is shown in Figure 3.2

1 5

4 2

3 6

Figure 3.2: Sample network

There are some different types of network problems, e.g. the transportation problem, the assignment problem and the shortest path problem, more about these can be found in [1]. The problem that will be treated here is the so called minimum cost network flow problem. It is a linear program on the form

min c^Tx subject to

Ax = b l ≤ x ≤ u.

The vectors x and c both consist of n elements. For each arc (i, j) there is a flow variable x_i,j in x and a cost coefficient c_i,j in c. The flow variable x_i,j records the flow in arc (i, j) and the cost coefficient c_i,j determines the cost for shipping one unit of flow over (i, j). The vectors l and u represent lower and upper bounds on x.

Each row of the constraint matrix A corresponds to the flow balance at each of the network’s m nodes. The matrix A will therefor be an m × n matrix, where

(17)

the net flow at each node must be equal to either the supply or the demand at the very same node. Whether supply or demand is to be considered is a question of definition.

Let the net flow at each node be defined as outflow minus inflow. With this convention a positive net flow describes supply and a negative demand. For the i-th node the flow balance then becomes

X

j

xi,j−X

k

xk,i = bi,

where bi> 0 describes supply and bi< 0 describes demand. Demand is simply treated as negative supply. The vector b then collects all supply (or demand) bi for the nodes in the network.

The matrix A and the vectors x and b for the sample network in Figure 3.2 would be written

x = x12 x13 x15 x24 x34 x36 x45 x64 ^T

A =







1 1 1 0 0 0 0 0

−1 0 0 1 0 0 0 0

0 −1 0 0 1 1 0 0

0 0 0 −1 −1 0 1 −1

0 0 −1 0 0 0 −1 0

0 0 0 0 0 −1 0 1





 , b =





 b1

b2

b3

b4

b5

b6





 .

Each column in A represents an arc (i, j). Because the arc (i, j) carries flow from i to j, +1 appears in row i and -1 in row j. All the other entries of the column must be zero. The matrix A can be said to describe the topology of a network. Sometimes in network flow problems A is called the node-arc incidence matrix, [2].

More information on network problems can be found in [1] and [2].

(18)

4 Problem description

The purpose of this thesis is to solve the problem of how to model future cross- border flows between an arbitrary number of market areas, I. The flows are assumed to be entirely price driven, which is a direct consequence of implicit capacity auctions, see section 2.4.4. Implicit capacity auction on all modelled borders therefore becomes a prerequisite.

The flows should be modelled over 24 hour periods, as cross-border flows are determined at day-ahead auctions, see chapter 2. The model should be able to produce hourly flow values for any desired time interval.

4.1 Merit order and price setting

In order to model price driven flows a price setting mechanism is needed. As mentioned in section 2.2.1 the real pricing is done by matching supply and demand curves composed of sell and buy bids at power exchanges. In this section a price setting model using merit orders and inelastic demand is introduced.

The total electricity supply in a market area can be described with a so called merit order curve. For each hour of the day a merit order is put together ranking all available plants based on production cost. This results in monotonically increasing step functions where the cheaper plants are found to the left and the more expensive plants to the right, see Figure 4.1. Because many plants of the same energy type often share the same production cost they can be put together into bigger blocks, these blocks will in this thesis be referred to as plant clusters.

Figure 4.1: German merit order

As seen in the figure a grey shaded margin is added to the plants’ production

(19)

costs. This margin is individual for each plant and represents the spread between the price at which a plant sells its power at the market and the production cost.

The modelling of this margin does not lie within the framework of this thesis and therefore margin-complete merit orders are taken for granted as input to the model. This means that the price associated with each plant can be understood as production cost plus margin.

To the very far left in the German merit order in Figure 4.1 a quite con- siderable amount of zero priced capacity can be found. These are the so called must-run plants that comprises energy types like wind, CHP, hydro and waste.

The demand in a market area can be approximated with a vertical line, see Figure 4.1. This is due to the fact that short-term demand is very inelastic.

Here the demand represents the actual consumption in a market area and not the buy bids seen on the wholesale market auctions.

The electricity price is now set by the intersection between the supply and the demand curves, in this case the intersection between merit order curve and inelastic demand. The plants to the left of the demand are in use and the plants to the right are not. The price is therefor set by the last, most expensive, plant needed to cover the demand.

4.2 Network of market areas

The modelling problem can be viewed as a network with I nodes and n arcs.

The nodes represent the market areas and the arcs the interconnectors. When electricity flows between the market areas the supply and demand situations change and as a consequence of that prices change. The network is simply a system of market areas that interact with each other via cross-border flows, see Figure 4.2.

MA I

MA 2 MA 1

Figure 4.2: Interacting system of market areas

The flows between the market areas are for every hour limited by maximum

(20)

transfer capacities and ramping restrictions, determined by the TSOs. The restrictions exist due to technical interconnector constraints, see section 2.4.2.

In order to make decent future forecasts of the cross-border flows, an accurate model for future maximum capacities will be required, although this will not be dealt with in this thesis.

To sum this chapter up a model with the following input and output for all desired hours has to be constructed:

Input

• Merit order for all market areas

• Demand for all market areas

• Maximum Capacity and Ramping values for all cables

• Topology for the network Output

• Hourly cross-border flows for all cables

(21)

5 Problem definition

The problem described in chapter 4 can be modelled as a linear optimisation problem. A mathematical definition of this problem will be derived step by step in this chapter.

The objective of the optimisation is to minimise the system-wide cost of electricity production over 24 hours. By minimising the total cost, the price differences will also be minimised, which means that the optimal price driven flows will be found.

The question is now how the system-wide cost can be described. One sug- gestion is to integrate the merit order and treat that as a cost function.

Fcost= Z

fmerit order

The integral of a function in 1 dimension specified in a certain interval can be interpreted as the area under the graph in the same interval. See A in Figure 5.1.

Mw

€/Mwh Merit order curve

Demand

∆x P∆x

A

∆A

x'

Figure 5.1:

Consider the small capacity range, ∆x, in Figure 5.1. Assume for simplicity that ∆x is chosen where no steps in the merit order occur. The cost for producing ∆x MWh, i.e. ∆x MW during 1 hour, is calculated by multiplying the amount of capacity, ∆x, with the price, P_∆x, and the time, 1 hour. The product, P_∆x· ∆x, can be represented by the rectangular area ∆A and describes the cost. This means that the area A, which can be seen as the sum of all ∆A:s between 0 and x⁰, describes the cost to produce x⁰ MWh. The area A is also the same as the integral of the merit order from 0 to x⁰.

(22)

5.1 The objective function

In order to minimise the system-wide cost of production, an appropriate objective function has to be defined.

The cost of production in a market area is equal to the integral of, or the area under, the merit order. Since the merit order is a step function, the integral from 0 to a certain capacity can be expressed as a finite sum. Before this sum can be formulated some merit order related variables need to be defined, see Figure 5.2

Mw

€/Mwh Merit order curve

Demand

C t,i,2 C t,i,3 P t,i,1

C t,i,j C t,i,1

P t,i,j

P t,i,2 P t,i,3

Figure 5.2: Defining variables

Ct,i,j = Capacity of the j:th plant cluster

Pt,i,j = Price corresponding to the j:th plant cluster

The indices t, i and j are used to indicate hour, market area and plant cluster.

Now a sum that describes the cost in one market area for one hour can be created as,

J_t,i⁰

X

j=1

P_t,i,j· Ct,i,j, (5.1)

where J_t,i⁰ is the number of active plant clusters in the current market area.

By introducing the weight variable w_t,j,i∈ [0, 1] which can ”turn on” and ”turn off” any chosen plant cluster, J_t,i⁰ is set to J_t,i. Here J_t,i denotes the number of available plant clusters in the current market area. Adding wt,i,j to (5.1) gives

J_t,i

X

j=1

P_t,i,j· C_t,i,j· w_t,i,j. (5.2)

(23)

The next step is to expand expression (5.2) so that the cost for the whole system is covered. Since the cost for each single market area already is explained, all there is to do is to sum them up, which gives

I

X

i=1





Jt,i

X

j=1

Pt,i,j· Ct,i,j· wt,i,j



, (5.3)

where I is the number of market areas.

The last step in establishing a proper objective function is to take the time perspective into consideration. As mentioned before the optimisation needs to contain all 24 hours. This means expanding (5.3) with yet another sum,

24

X

t=1





I

X

i=1





Jt,i

X

j=1

Pt,i,j· Ct,i,j· wt,i,j







. (5.4)

Now an expression for the system-wide cost for all 24 hours has been derived.

The objective will hereby be defined as finding the weights, wt,i,j, minimising (5.4). The objective function is written,

min

24

X

t=1





I

X

i=1





J_t,i

X

j=1

P_t,i,j· Ct,i,j· wt,i,j







. (5.5)

5.2 Flow balance constraint

In all market areas supply must meet demand at all times. This fact is presented and defined as the first constraint.

Consider a market area, i, with a demand, Demandt,i, at a certain hour, t.

By using the weights defined in section 5.1, the supply in i can be written

Supplyt,i=

Jt,i

X

j=1

Ct,i,j· wt,i,j. (5.6)

Now (5.6) is used to write the electricity balance equation, supply equals demand, for i,

Jt,i

X

j=1

Ct,i,j· wt,i,j= Demandt,i. (5.7) If the market area i is interconnected with other market areas, import and export flows can exist. These flows are defined as

F lowt,i→k = Export-flow, from i to k

F lowt,k→i = Import-flow, from k to i, (5.8) where k indicates the market areas which are connected with i. This definition implies that for each physical cable a distinction between the two different flow directions is made. One real cable can simply be regarded as two cables in the model. The advantage is that negative flows don’t have to be considered.

(24)

By summing up all the export and import flows, they can be incorporated in the balance equation. Exports are treated as additional demand and imports as additional supply, this gives

J_t,i

X

j=1

C_t,i,j· wt,i,j− X

k∈Ki

F low_t,i→k+ X

k∈Ki

F low_t,k→i= Demand_t,i, (5.9)

∀ t, i,

where K_i is the set of indices representing market areas that are connected with i. Equation (5.9) is the flow balance constraint.

5.3 Cable capacity constraint

Each interconnecting cable has two maximum transfer capacities, one for each direction. The two maximum capacities limit the maximum values of the flows in (5.8), thus

max (F low_t,i→k) = M axCap_t,i→k

max (F low_t,k→i) = M axCap_t,k→i. (5.10) The time dependence refers to the fact that the maximum transfer capacities do not have to be the same from hour to hour, see section 2.4.2. Maintenance and cable outages are factors that may invoke further reductions on the maximum capacities. The cable capacities will also differ from summer to winter.

The two equations in (5.10) can in fact be reduced to only one,

M axCapt,i→k = max (F lowt,i→k) . (5.11) It is not necessary with two equations because the flow in a cable is an export flow and import flow at the same time.

From (5.11) the following inequality derives,

F lowt,i→k≤ M axCapt,i→k. (5.12) Combining (5.12) with the fact that all flows are non-negative, see section 5.2, gives the cable capacity constraint,

0 ≤ F low_t,i→k≤ M axCapt,i→k, (5.13)

∀ t, i, k ∈ K_i.

Due to similar form another constraint will be presented here. The weight variables that were introduced in section 5.1 are bounded between 0 and 1, hence

0 ≤ wt,i,j≤ 1, (5.14)

∀ t, i, j.

(5.14) is called the weight constraint.

(25)

5.4 Ramping constraint

The maximal change in flow from hour to hour in one interconnector is limited to a certain value, the ramping. This means that the flow for hour t depends on the flow from the previous hour, t − 1. By merging the previously defined flows, F lowt,i→kand F lowt,k→i, a net flow for each cable can be expressed as

N etF lowt,i→k= F lowt,i→k− F lowt,k→i. (5.15) For positive values on N etF low_i→k,t the flow goes from i to k, for negative it goes from k to i.

Let the non-negative variable R_t,i→k denote the ramping. Then for the net flow of hour t the following inequalities must hold,

N etF lowt,i→k≤ N etF lowt−1,i→k+ Rt,i→k, (5.16) N etF lowt,i→k≥ N etF lowt−1,i→k− Rt,k→i. (5.17) Inequality (5.16) suggests that the net flow of t can not be bigger than the net flow of t − 1 plus the ramping Rt,i→k. In the same manner (5.17) states that the net flow of t also can not be smaller than the net flow of t − 1 minus the ramping Rt,k→i. Since the flow change can have two directions, two different ramping variables are used. Rt,i→kdoes not necessarily have to be the same as Rt,k→i,t.

When (5.15) is inserted in (5.16) and (5.17) the ramping constraint can be formulated as

−F lowt−1,i→k+ F low_t−1,k→i+ F low_t,i→k− F lowt,k→i ≤ R_t,i→k, F low_t−1,i→k− F low_t−1,k→i− F low_t,i→k+ F low_t,k→i ≤ R_t,k→i,

(5.18)

∀ t, i, k ∈ Ki.

(26)

5.5 Summary

By gathering and arranging the objective function, (5.5), and all the constraints, (5.9), (5.13), (5.14) and (5.18), a final definition of the problem can be established. The objective is finding the weights, wt,i,j, and the flows, F lowt,i→kand F lowt,k→i, minimising

min P24 t=1

PI i=1

P^Jt,i

j=1P_t,i,j· C_t,i,j· w_t,i,j

subject to:









 P^Jt,i

j=1Ct,i,j· wt,i,j−P

k∈KiF lowt,i→k

+P

k∈K_iF lowt,k→i= Demandt,i, ∀ t, i

−F lowt−1,i→k+ F lowt−1,k→i

+F lowt,i→k− F lowt,k→i≤ Rt,i→k, ∀ t, i, k ∈ Ki

F lowt−1,i→k− F lowt−1,k→i

−F lowt,i→k+ F lowt,k→i≤ Rt,k→i, ∀ t, i, k ∈ Ki

0 ≤ F low_t,i→k≤ M axCapt,i→k, ∀ t, i, k ∈ Ki

0 ≤ w_t,i,j≤ 1, ∀ t, i, j

(5.19)

where

t = index for hour (1-24) i = index for market area j = index for plant cluster

k = index for interconnected market area I = number of market areas in the model

Jt,i = number of plants clusters in each market area for each hour Ki = Set containing i:s interconnected market areas, k

Pt,i,j = Specific price for j:th plant cluster of each supply function Ct,i,j = Specific capacity for j:th plant cluste of each supply function wt,i,j = Weights determining which plant clusters that are active Demandt,i = Fixed demand for each market area and hour

F lowi→k,t = Flow from i to k at hour t F lowk→i,t = Flow from k to i at hour t

(27)

6 Vector and matrix representation

The problem definition, (5.19), from the previous chapter will in this chapter be translated into vectors and matrices.

The objective function and the constraints of (5.19) have to be constructed so that they fit the following definition,

min c^Tx

subject to







Aeqx = beq

Ax ≤ b lb ≤ x ≤ ub.

(6.1)

MatLab’s solver linprog accepts linear problems written on the form given by (6.1). Since linprog is the primary choice of solver for the optimisation problem in this thesis it makes sense to express (5.19) in terms of (6.1).

First the vector x is constructed. The vector x must contain all the optimisation variables, i.e. all the weights, w_t,i,j, and all the flows, F low_t,i→k, F low_t,k→i. We start with collecting the weights and the flows for a given hour, t, in two vectors,

wt=





 w_t,1,1 wt,1,2

... wt,1,J_t,1

wt,2,1

wt,2,2

... w_t,2,J_t,2

... ...

wt,I,1

wt,I,2

... wt,I,J_t,I







, (6.2)

Flow_t=







F lowt,1→2

F lowt,2→1

F low_t,1→3 F low_t,3→1

... F low_t,i→k F lowt,k→i

...







. (6.3)

(28)

The vector x can now be defined for all 24 hours by using (6.2) and (6.3),

x =





 w1

w2

... w24

Flow1

Flow2

... Flow₂₄







. (6.4)

6.1 Objective function

The aim is to find the cost vector, c, that makes the scalar product c^Tx equivalent with the sum definition of the objective function, (5.4).

Analogous with the approach for x, an hour-specific vector can be assembled,

c_t=







Pt,1,1· Ct,1,1

Pt,1,2· Ct,1,2

...

Pt,1,Jt,1· Ct,1,Jt,1

P_t,2,1· Ct,2,1

P_t,2,2· C_t,2,2 ...

Pt,2,J_t,2· Ct,2,J_t,2

... ...

Pt,I,1· Ct,I,1

Pt,I,2· Ct,I,2

... P_t,I,J_t,I· Ct,I,Jt,I







. (6.5)

The vector ctcorresponds to the hour-specific weight vector, wt. Since no flows are involved in the objective function the final cost vector has to be filled out

(29)

with zeros, i.e.,

c =





 c₁ c₂ c3

... c24

0 0 0 ... 0







. (6.6)

The length of the additional zero-vector must be the same as the length of





 Flow1

Flow₂ ... Flow₂₄





 .

The objective function c^Tx is now described by the scalar product of the two vectors (6.6) and (6.4).

(In the actual modelling, the zeros will be replaced with small positive num- bers. Small means small in relation to the costs, P_t,i,j· Ct,i,j. This is done to prevent spontaneous flows from arising.)

6.2 Equality constraint

For all hours and market areas, the flow balance constraint, (5.9), must be satisfied. This constraint ends up in a system of equations that can be expressed as the matrix equality Aeqx = beq. Here Aeqx contains the supply and the flows of (5.9) whereas b_eq contains the demand.

Let’s start to look at the supply. The supply is given by the capacities of the plant clusters multiplied with their weights. Consider the vector,

C_t,i=





 C_t,i,1 C_t,i,2

... Ct,i,J_t,i







. (6.7)

This vector collects all the available plant cluster capacities for market area i at hour t. By using (6.7) a matrix that comprises all the I market areas can be written

Ct=





 C^T_t,1

C^T_t,2 . ..

C^T_t,I







. (6.8)

(30)

Each row of C_tnow contains available plant capacities for one market area. The blank entries in (6.8) are zeros that have been excluded in order to get a better graphical overview. All blank matrix entries will hereafter represent zeros.

The supply defined defined by equation (5.6) is given by Ctwt, hence







Supplyt,1

Supplyt,2

... Supply_t,I







=





 C^T_t,1

C^T_t,2 . ..

C^T_t,I







· wt. (6.9)

Now a matrix for the flows needs to be constructed. The node-arc incidence matrix from section 3.4 describes the topology of a network, i.e. how the nodes are connected to each other. The topology of this problem is represented by a node-arc incidence matrix denoted N AIM . An example of N AIM , for a network with 6 market areas (nodes) and 8 non-negative cables (arcs), looks like this,







−1 1 −1 1 0 0 0 0 0 0

1 −1 0 0 −1 1 0 0 0 0

0 0 1 −1 0 0 −1 1 −1 1

0 0 0 0 0 0 1 −1 0 0

0 0 0 0 1 −1 0 0 0 0

0 0 0 0 0 0 0 0 1 −1





 .

In this node-arc incidence matrix −1 denotes outflows and +1 denotes inflows.

This is fully in line with (5.9), where outflows are negative and inflows positive.

Create a new matrix by horizontal concatenation of Ct and N AIM ,





 C^T_t,1

C^T_t,2 . ..

C^T_t,I

−1 1 0 0 · · · 0

1 −1 −1 1 · · · 0 ... ... ... ... . .. ... 0 0 1 −1 · · · −1







| {z }

N AIM

. (6.10)

Note that the matrix representing N AIM displays an example of a topology configuration. Now multiply (6.10) with the vector (w^T_t Flow^T_t)^T. The resulting rows of (6.10) represents, for each market area, the left-hand-side of the flow balance constraint, (5.9). For hour t the flow balance constraint then yields,





 C^T_t,1

C^T_t,2 . ..

C^T_t,I

−1 1 0 0 · · · 0

1 −1 −1 1 · · · 0 ... ... ... ... . .. ... 0 0 1 −1 · · · −1







·

wt

Flow_t

=







Demandt,1

Demandt,2

... Demandt,I







, (6.11)

(31)

where the demand vector is the right-hand-side of (5.9). Let this vector be denoted Demand_t. The flow balance constraint is in (6.11) formulated for all market areas, but only for one hour. Hence, the next step is to include all 24 hours,







C1 N AIM

C2 N AIM

. .. . ..

C24 N AIM







| {z }

A_eq

·





 w1

w2

... w24

Flow1

Flow2

... Flow₂₄







| {z }

x

=







Demand₁ Demand₂

... Demand₂₄





 .

| {z }

b_eq

(6.12)

Hereby A_eq and b_eq are identified.

6.3 Inequality constraint

The ramping constraint needs to be expressed as the matrix inequality Ax ≤ b.

The matrix A can be considered as two concatenated matrices, (AweightsAf low), one for weights and one for flows. The part of A corresponding to the weights will be an all zero matrix, this is because the weights aren’t affected by ramping.

The matrix for the flows needs more investigation. Starting with rewriting the ramping constraint (5.18) as:

−1 1 1 −1

1 −1 −1 1







F low_t−1,i→k F lowt−1,k→i

F lowt,i→k

F lowt,k→i







≤Rt,i→k

Rt,k→i

. (6.13)

This only applies for one connection and one hour. For all connections (6.13) becomes

−B BFlowt−1

Flowt

≤ Rt, (6.14)

(32)

where

−B B =







−1 1 1 −1

1 −1 −1 1

−1 1 1 −1

1 −1 −1 1

. .. . ..

−1 1 1 −1

1 −1 −1 1

. .. . ..







and

Rt=





 Rt,1→2

Rt,2→1

Rt,1→3

Rt,3→1

... R_t,i→k R_t,k→i

...





 .

In order to expand (6.14) to a full 24-hour period a new vector needs to be introduced,

Flow0=







F low0,1→2− F low0,2→1

F low0,2→1− F low0,1→2

F low0,1→3− F low0,3→1

F low0,3→1− F low0,1→3

...

F low0,i→k− F low0,k→i

F low_0,k→i− F low0,i→k

...





 ,

which contains flows for t = 0, i.e. the hour before the first hour of the optimisation. Flow₀ is needed to express the ramping for hour 1 and will therefore contain the flows of hour 24 from the optimisation of the day before. Expression (6.14) expanded to hold for all 24 hours can now be written





 B

−B B

. .. . ..

−B B











 Flow1

Flow2

Flow3

... Flow24







≤







R1+ Flow0

R2

R3

... R24







. (6.15)

By adding the weights and the previously mentioned zero matrix the inequality (6.15) becomes

(33)







0 B

0 −B B

. .. . .. . ..

0 −B B







| {z }

A





 w₁ w₂ ... w24

Flow1

Flow2

... Flow24







| {z }

x

≤







R1− Flow0

R2

R3

... R24







| {z }

b

,

(6.16) where A and b can be identified.

6.4 Upper and lower bound

The upper and lower bounds are rather uncomplicated, they are both described by two column vectors of the same length as x. The cable capacity constraint (5.13) and the weight constraint (5.14) must be combined and written on the form, lb ≤ x ≤ ub.

First consider the weights. Inequality (5.14) can at hour t for all market areas be written,

0 ≤ wt≤ 1, (6.17)

where 0 and 1 are two vectors containing respectively only zeros and only ones.

For the flows the inequality (5.13) becomes

0 ≤ Flow_t≤ MaxCap_t, (6.18)

where 0 once again is the zero vector and

MaxCap_t=







M axCap_t,1→2 M axCapt,2→1

M axCapt,1→3

M axCapt,3→1

... M axCapt,i→k

M axCapt,k→i

...







. (6.19)

When (6.17) and (6.18) are expanded to 24 hours and combined, the upper and lower bound can be identified as,

(34)





 0 0 ... 0 0 0 ... 0







| {z }

lb

≤





 w₁ w₂ ... w24

Flow1

Flow2

... Flow24







| {z }

x

≤





 1 1 ... 1 MaxCap₁ MaxCap₂

... MaxCap₂₄







| {z }

ub

. (6.20)

(35)

7 Simulation

To test the model, a simulation of cross-border flows for one day in the nordic areas is preformed. The nordic area, involving Denmark, Finland, Norway and Sweden, is chosen because implicit auctions are applied on all borders and because input data is available. The day which will be considered is 5 April, 2008.

In the nordic area, 7 different market areas with a total of 9 cross-border connections exist, see Figure 7.1. In order to simulate the flows, for 24 hours,

FI NO2

NO1

SE

NO3

DK1

DK2 Interconnections

Figure 7.1: Market areas and interconnections in the Nordic area on all these interconnectors,the optimisation problem described and defined in chapters 4 - 6 needs to be solved. Remember (6.1),

min c^Tx

subject to







Aeqx = beq

Ax ≤ b lb ≤ x ≤ ub.

For this linear program Matlab’s linear solver linprog is used. The linprog solver can use different algorithms depending on the users wishes.

7.1 Special features

The output from the model will be compared with output from an already existing model for the nordic region. The nordic model doesn’t take the ramping constraint into account and in order to make the comparison fair, the ramping

Optimising cross-border flows in liberalised electricity markets

M A S T E R ' S T H E S I S