multistage and stochastic optimization of bidding strategies for a price ...

2 downloads 0 Views 453KB Size Report
where: λ and E are the price-maker bidding price ($/MWh) and quantity (MWh), ...... [30] Read, E.G. George, J.A and MacGregor, A.D. Stochastic dual dynamic ...
Long Term Optimal Allocation of Hydro Generation for a Price-Maker Company in a Competitive Market: Latest Developments and a SDDP approach

Bruno da Costa Flach Pontifical Catholic University of Rio de Janeiro R. Marquês de São Vicente, 225 – Rio de Janeiro, RJ – Brazil [email protected]

Luiz Augusto Barroso PSR Praia de Botafogo, 228 / 1701-A – Rio de Janeiro, RJ – Brazil

[email protected]

Mario Veiga Pereira PSR Praia de Botafogo, 228 / 1701-A – Rio de Janeiro, RJ – Brazil

[email protected]

1

ABSTRACT Since the liberalization of the power industry, there has been a large amount of literature on the determination of optimal bidding decisions for price maker energy producers. The vast majority of the work developed so far has focused on short-term horizons and may be viewed as successful approaches for systems whose operation is generally deterministic. In the case of price maker hydro plants with significant storage capacity, however, the solution of the strategic bidding problem is more subtle. The reason is that hydro reservoirs allow the bidder to postpone energy production if future prices are expected to be higher than the current price. This demands the management of an energy-constrained resource and determines a time-coupling characteristic to the problem, implying that the bidding strategy should ideally take into account the following stages and consider the stochasticity of inflows. These aspects characterize the strategic bidding for price maker hydro agents as a multi-stage stochastic programming problem, with significant computational challenges. The objective of this work is to present a new methodology for the strategic bidding problem of a price-maker hydropower based company, taking into account several hydro plants, time-coupling and stochastic inflow scenarios. The proposed approach considers a deterministic residual demand curve and is based on stochastic dual dynamic programming (SDDP), which has been successfully applied to the least-cost hydrothermal scheduling problem. Since the technique requires the problem to be concave, a piecewise linear approximation of the expected future benefit function is proposed. The application of the methodology is exemplified with a real case study based on the hydrothermal system of El Salvador.

KEYWORDS: Hydroscheduling, Strategic Bidding, Market Power, Stochastic Optimization.

2

1. Introduction Several countries all over the world have implemented reforms in their power sectors, with emphasis on competition and private investment. Although the details of the regulatory frameworks change in each country, the overall organization in most cases follows the same principles, known as the “standard model” [1]. One of its basic features is a short-term electricity market, where energy sales and purchases are settled. In a simplified way, the short-term market works as follows: (i) at the end of each day, generators and loads submit price and quantity bids for the following day; (ii) an equilibrium process is then simulated where the settlement price (or spot price) is adjusted until total dispatched generation is equal to total load; (iii) the energy produced by the dispatched generators (consumed by load agents) is remunerated (charged) based on the energy spot price. The existence of a bid-based dispatch/settlement poses complex technical challenges for energy-producing companies, who are now faced with the problem of having to devise bidding strategies aimed at maximizing their expected net revenue. For example, [2] shows that the optimal strategy for a price-taker producer is to bid the plant’s variable operating cost. However, the same problem gets a lot more complex and interesting for price-maker agents, i.e., those that can affect the equilibrium price with their own offers. The strategic bidding problem for price-maker agents is usually formulated as a bi-level optimization problem, which is in turn transformed into a nonlinear programming model through the use of Karush-Kuhn-Tucker complementarity conditions [3], [4]. Because the resulting problem is non-convex, there has been an intensive search for efficient solution methods with proposed solution schemes ranging from heuristics and iterative procedures to specialized algorithms and specific equilibrium models ([5], [6], [7], [8], [9], [10]). Although these solution schemes have proven to be efficient, it is usually not possible to ensure that the

3

global optimal solution has been found. This has thus motivated the recent development of a new generation of models based on mixed integer linear programming ([11], [12], [13], [14]). Nevertheless, these previous developments may be viewed as successful approaches for systems whose scheduling decisions are generally decoupled in time, which is the case for those that rely predominately on thermal generation. In the case of hydro plants, however, the strategic bidding problem solution is more subtle. The reason is that the hydro reservoirs allow the bidder to postpone energy production if future prices are expected to be higher than the current price. This demands the management of an energy-constrained resource and determines a time-coupling characteristic to the problem, implying that the bidding strategy should ideally take into account the following stages and consider the stochasticity of inflows. These aspects characterize the strategic bidding for hydro agents as a multi-stage stochastic programming problem. 1.1. Strategic bidding for price-taker hydropower companies The objective of a hydro-based generation company in a liberalized electricity market is to develop bidding strategies that maximize their revenues, taking into account the use of its storage capabilities to schedule its energy production over time. If the generation company is a price-taker agent (i.e., it cannot affect the spot price with its bidding strategy) then attention must focus on profit uncertainty caused by uncertainty in spot prices and reservoir inflow. Given uncertainty in these two parameters, the revenue-maximization problem can be formulated as multistage stochastic linear program whose solution approach has been widely studied in different contexts. This problem was first discussed by [33] and [34], where an algorithm based on stochastic dynamic programming was introduced. This approach later became used in generation scheduling models under price uncertainty in Norway ([35], [36]). More recently a large literature on hydroelectric bidding for price-taker companies has become available. They differ from each other essentially on the modeling of system

4

components ([37]), incorporation of risk constraints in the bidding decisions ([38], [39]) and on the modeling of the short-term bidding curve and electricity market ([40]). Reference [41] provides a good overview of this problem. In general, it can be said that there is a consensus on the solution approaches devised for this problem. This situation is very different from the bidding problem for price maker hydroelectric companies, as discussed next. 1.2. Strategic bidding for price-maker hydropower companies Bidding strategies for price-maker hydroelectric companies involves situations where the bidder can affect the market price with its bids. In this case the problem is non-convex, multi stage and stochastic. While it has also been analyzed in the literature during the last few years, the amount of research on this problem is significantly smaller than that in the case of price-taker hydroelectric companies. The classic approach to address it has been through the application of a “direct solution” through the same mathematical-programming solution methods described earlier for one-stage problems in a multi-stage framework. Nevertheless, a deterministic approach is usually taken and future inflows are assumed to be known, such as in the pioneer work done by [8]. This results in a computationally feasible direct solution scheme but it turns out to be an oversimplification in case of hydro companies with large reservoirs that need to be managed some months (or even years) ahead. In these cases, a decision-making process that takes inflow uncertainty into consideration is desirable. Multi-stage and stochastic problems for hydroscheduling have been traditionally solved by stochastic dynamic programming (SDP). It handles the nonlinearity of the bidding problem but the computational effort increases exponentially with the number of state variables – the so-called “curse of dimensionality”. Therefore, its use in strategic bidding problems is restricted to companies with no more than 2 or 3 hydro plants ([15], [16], [17]). On the other hand, decomposition algorithms have proven effective in allowing the solution of large-scale stochastic linear programs (see, for example, [18] and [19]). For

5

example, systems with considerable share of hydropower have been using least-cost hydrothermal scheduling tools for at least two decades. In this framework, the objective is to determine an operation strategy of a hydrothermal system that for each stage of the planning period sets generation targets for each plant (hydro releases and thermal production). This strategy should minimize the expected value of the operation cost along the period, composed of fuel cost and penalties for failure of load supply. A DP scheme based on multi-stage stochastic Benders decomposition ([20], [21]) (also known as Stochastic Dual Dynamic Programming or SDDP) has been used to approximate the cost-to-go functions by a set of linear inequalities, known as Benders cuts, avoiding the previously mentioned curse of dimensionality of traditional SDP models. The major advantage is the possibility to represent hydro plants individually. This technique has been applied to the scheduling of large-scale power systems in more than thirty countries, including detailed modeling of system components and transmission networks, in a stochastic and multi-stage framework and with a reasonable computational effort [22]. The SDDP algorithm has led to a number of related methods and applications that are based on the same essential idea, but seek to improve the method by exploiting the structure of particular applications. Therefore, one could immediately wonder about the application of this technique to the strategic bidding problem. However, the SDDP technique requires convexity on the cost-to-go functions, so that it can be approximated by Benders cuts. As mentioned earlier, the strategic bidding problem is nonconvex, thus hindering the straightforward application of this decomposition technique. 1.3. Objective of this work The objective of this work is to review the challenges in devising long term bidding strategies for a price-maker hydroelectric company and to propose a new approach based on stochastic dual dynamic programming. We identify the main difficulties related to the fact that one needs to take into account several hydro plants, time-coupling, inflow uncertainty

6

and the non-convex nature of the relationship between the agent's bids and market outcomes. We start by discussing in detail the application of stochastic dynamic programming to this problem and the computational limitations of this approach. We then proceed to discuss the research direction we pursued, which is focused on the development of a computationally efficient approach to construct the revenue-to-go functions (similar to cost-to-go functions of the least-cost hydroscheduling problem) in a dynamic programming framework. Our approach is based on the application of multi-stage stochastic programming with decomposition techniques. This should allow us to solve larger instances of the problem that would otherwise be computationally infeasible through a direct solution approach but at the cost of some simplifications to the bidding problem. We follow a quantity-only bid approach, which is arguably a coherent model for hydro agents [23][24], and our focus is on the definition of a long-term operating strategy for a single hydroelectric company. Once the long-term uncertainty is represented and an operating policy is calculated, a short-term bidding strategy for day-ahead markets can be devised using MIP algorithms, such as those developed in [11]-[14], that may be adapted to optimize the revenue in the short-term market plus the expected value of the revenue in the future given by the future revenue function calculated by the approaches discussed in this work. This shortterm optimization is not discussed here, our focus is on the modeling of long-term influences. In the same way, we do not discuss the modeling of several price-maker agents, which is left for future work. This means that our model is a particular case of a more general market equilibrium model. This work is organized as follows: sections 2, 3 and 4 review the single-stage and multistage, stochastic strategic bidding problem. These sections introduce the problem in gradual levels of difficulties and discuss existing solution approaches (and their limitations) for the reader not familiar with the subject. Section 5 discusses in detail the proposed methodology,

7

the computational results of which are presented in Section 6. Finally, conclusions and possible future work directions are stated in section 7.

2. Price Maker Strategic Bidding: overview 2.1. Market Clearing Procedure Once the price and quantity bids of the agents are submitted, the market operator (or the independent system operator, depending on the market arrangement) carries out a market clearing procedure (MCP) to determine the spot price and dispatched generation (there is no possibility for an agent to alter his bid after the MCP is carried out). In general, the MCP can be represented by the following linear problem:

Min

λg o + ∑ j =1 c i g i

(dual variable)

subject to

g o + ∑ j =1 g i = d

πd

(1.a)

g i ≤ g j , ∀j

πgj

(1.b)

go ≤ E

πgj

(1.c)

J

J

(1)

where: λ and E are the price-maker bidding price ($/MWh) and quantity (MWh), go is the actual generation of the price-maker agent (MWh), J is the set of price-taker agents, cj and g j

are the bidding price ($/MWh) and quantity of price-taker agent j, gj is the actual generation of the price-taker agent j (MWh) and d is the system demand (MWh). Equation (1.a) represents the load supply and constraints (1.b)-(1.c) are the generation upper bounds. It is implicitly assumed that gj ≥ 0, j = 1, ..., J. Associated with constraints

8

(1.a)-(1.c) of the MCP there is a set of Lagrange multipliers: πd is the spot price, i.e. the marginal cost of increasing load supply ($/MWh), and πgj is the marginal benefit of increasing generator j’s quantity bid . 2.2. Generator net revenue

The price-maker generator receives a gross revenue ($), or operating profit, given by the product of spot price πd and its energy production gj. The net revenue of the price-maker agent, represented by RO, is obtained by deducting the variable operational cost from the total gross revenue: RO = (πd – c) go

(2)

where c is the variable operation cost of its plant ($/MWh). Note that c may be different from λ, the bid price. 2.3. Price-maker bidding strategy

The agent’s objective is to obtain the combination of price and quantity bids which maximizes its total net revenue:

Max

RO = (πd – c) go

subject to

E≤E

(3) (3.a)

where E is the agent’s available capacity (MWh). The strategic bidding problem (3) is a bilevel optimization problem, which besides primal variables (gj) uses dual (πd) results from the economic dispatch – problem (1) – that also depend on the agent’s bid. By using complementarity conditions, one can obtain a one-level non-linear and non-convex

9

formulation of the problem, which is known as an MPEC (mathematical programs with equilibrium constraints), as shown in ([4]-[8]). 2.4. Quantity-only offers

The adoption of a quantity-only bidding model – in which the price bids are set to zero and only quantities are optimized – helps to simplify the problem. The rationale behind this model is that by restricting the quantity offered to the system, the agent may be able to raise the spot price and increase its revenue. In this case, problem (3) can be re-written as (assuming the operational costs are zero):

Max

RO = πd.E

Subject to

E≤E

(4) (4.a)

Similarly, problem (1) can be re-formulated since the quantity offered by the agent will always be automatically dispatched:

Min



J

Subject to



J

j =1

j =1

ci gi gi = d − E

g i ≤ g j , ∀j

(5) (5.a) (5.b)

As discussed in [5,11], the resulting spot price will change discretely depending on the quantity offered by the price maker agent since it will be associated to the variable cost of the most expensive price-taker agent used to close the gap of the residual demand, i.e., the difference between system’s demand and the price-maker supply. Since the agent’s revenue 10

is the result of the spot price multiplied by the dispatched quantity, it will follow a “chainsaw” pattern, as depicted below:

Figure 2.1 – Agent’s revenue vs. amount of energy offered This non-linear and non-convex revenue function is the one to be optimized by a price maker agent. The number of breakpoints of this function is associated to the number of pricetaker agents that might be needed to fulfill the residual demand and, consequently, set the energy spot prices. 2.5. A Mixed Integer Programming (MIP) approach

This problem can be easily modeled as a MIP as follows:

Max



k =1

Subject to



k =1

K

K

π dk E k

Ek ≤ E

(6) (6.a)

Δ k −1 x k ≤ E k ≤ Δ k x k , ∀k

(6.b)



(6.c)

K k =1

xk ≤ 1

11

where πdk is the spot price associated to the kth segment of the agent’s revenue function ($/MWh), Ek is the agent’s quantity bid in the kth interval (MWh), E is the total amount of available energy, Δk is the upper bound of the kth interval into which the revenue function is split (Δ0 = 0) and xk is a binary variable that indicates the interval of the revenue function active at the optimal solution. The basic idea is to associate a binary variable to each interval of the quantity of supplied energy where the spot price remains constant. Hence if the revenue-function has K segments, the problem will have K binary variables. This approach was introduced in [11] and [25] and successfully applied to calculate the optimal bidding strategy of a price-maker company in a single-bus and quantity-only bidding model environment, as illustrated in examples derived from the Spanish system.

3. Deterministic Multi-Stage Bidding 3.1. Time-coupling: the need to look into the future

The problem described in the previous section was only concerned with a single stage. It was implicitly assumed that the bids made at one given stage had no effect on the following stages, i.e., that the problem was decoupled in time. As previously mentioned, a timecoupling characteristic is inherent to the problem of hydro agents due to the existence of water reservoirs, which enable them to transfer available energy from one stage to subsequent ones. The operation of hydro systems with multiple plants requires the introduction of several additional constraints such as the representation of the physical distribution of plants along river basins, storage and turbine capacities, production coefficients and, most importantly, mass balance equations, shown next:

12

vt+1(i) = vt(i) - ut(i) - st(i) + at(i) +



[ut(g) + st(g)]

(7)

g∈G (i )

where vt(i) is the stored volume at the beginning of stage t in plant i, ut(i) and st(i) are, respectively, the turbined and spilled volumes during stage t by plant i, at(i) is the inflow to the reservoir of plant i during stage t and G(i) is the set of plants upstream from plant i. Equation (7) depicts the intrinsic relationship between decisions taken in the present – turbined and spilled volumes – and their consequences in the future – available energy in the form of stored water.

3.2. Energy-constrained bidding problem: a MILP approach

A first solution approach for the multi-stage case is to assume that the agent knows in advance the amount of available energy that may be offered during the horizon, such as done in [8]. In this case, the MILP formulation (6) can be immediately extended to represent this time-coupling and total energy constraint, as presented below (where ρ(i) is the production coefficient of plant i and I indicates the set of hydroelectric plants owned by the agent):

Max

∑ ∑ T

K

t =1

k =1

Subject to

∑ ∑ T

K

t =1

k =1

(8)

π dtk Etk Etk ≤ E

(8.a)

Δ tk −1 xtk ≤ Etk ≤ Δ tk xtk , ∀t , ∀k

(8.b)



k =1

xtk ≤ 1, ∀t

(8.c)



k =1

Etk = ∑ ρ(i ) ⋅ u t (i ), ∀t

(8.d)

K

K

i∈I

[ut(g) + st(g)] vt+1(i) = vt(i) - ut(i) - st(i) + at(i) + g∈G(i) ∑

(8.e)

13

_ vt(i) ≤ v(i), ∀i, ∀t

(8.f)

_ ut(i) ≤ u(i), ∀i, ∀t

(8.g)

However, the computational effort required by integer programming algorithms increases very rapidly with the number of binary variables, which in this problem is equal to K (number of segments of the revenue function) x T (number of time stages). Therefore, for a problem with a moderate number of stages and segments of the revenue function, the computational effort can be significant or even prohibitive. Because of this drawback, the use of dynamic programming emerges as an alternative solution approach and will be discussed next. 3.3. Energy-constrained bidding problem: a DP approach

A usual alternative for problems that deal with sequential decisions is the use of a dynamic programming (DP) procedure. The DP works by decomposing a T-stage problem into T onestage problems (where the term stage is used to refer to time steps of any length). In this case, the amount of energy available at the beginning of a given stage should be defined as the state variable. The algorithm starts by solving the last-stage problem for a set of discretized values of available energy. The resulting pairs {amount of available energy; obtained revenue} are then interpolated providing a future benefit function (FBF) for the previous stage. This function is also known as revenue-to-go function or future revenue function and indicates the expected return the company will have from the stage t until the end of the horizon T. Analogously, the problem of stage T–1 is then solved for the same set of values of the quantity of available energy, but taking into consideration the potential revenue that might be obtained in stage T as a function of the amount of energy that is left to it, as given by the 14

FBF. The algorithm continues backwards to the first stage building the FBF recursively. This function summarizes all the necessary information about the future behavior of the system for the definition of the actually implementable decision in the first stage. A pseudo-code description of the procedure is presented below:

Initialize the FBF of the last stage βT+1(ET+1) = 0 For each t = T, T-1, ..., 1 For each value of available energy Etm = Et1, …, EtM

Solve a one-stage problem similar to (8) considering Etm as the amount of available energy and βt+1(Et+1) as the revenue-to-go function for stage t: βt(Etm) = Max subject to



k =1



k =1

K

K

π dtk Etk + βt+1(Et+1) Etk ≤ E mt

Et+1 = Etm –



K k =1

Etk

[ut(g) + st(g)] vt+1(i) = vt(i) - ut(i) - st(i) + at(i) + g∈G(i) ∑

Constraints (8.b) – (8.d), (8.f), (8.g) End

Create a complete FBF βt(Et) for the previous stage by interpolating the values obtained above. End

15

The representation of the interpolated FBF in the optimization problems presented above requires the use of binary variables, in a way similar to that used to represent the revenue function, and is omitted for the sake of ease of presentation.

4. Stochastic Multi-Stage Bidding 4.1. Challenges for hydro bidding

The assumptions made in the development of the last section allow the calculation of quantity-offers for energy-constrained resources. However, they are still simplistic to allow its practical application to real hydro systems. Water inflows that arrive to hydro plants do so continuously throughout the horizon and not at a single time period. This aspect could still be easily incorporated into the algorithm, but the fact is that inflows to the reservoirs can only be modeled as random variables, which makes the problem inherently stochastic. 4.2. Scenario-based approach

A simple and direct approach for the multi-stage stochastic case is again to formulate it as a large-scale mathematical programming problem with an underlying scenario-tree to represent possible realizations of the random variables [26]. However, an increasing number of plants requires an even larger increase in the number of scenarios to accurately describe the stochastic processes involved. Besides, the number of inflow branches grows exponentially with the length of the time horizon and the resulting stochastic optimization problem quickly becomes computationally intractable. Given that our aim is to solve problems where the planning horizon may be up to some years ahead, the traditional scenario tree approach would soon result impractical. This has motivated the development of solution approaches based on a state-space formulation, described next

16

4.3. Stochastic dynamic programming approach

Another possible way to deal with these extensions to the problem is the use of stochastic dynamic programming (SDP). Similarly to the DP approach described in Section 3.3, the algorithm defines a future benefit function that makes it possible to compare the immediate benefit of offering energy at a given stage with the expected benefits of storing it for future use. For the purposes of our work, the operation of the hydro plants at each stage will be affected by the initial level of each reservoir and by the inflow scenario. Since the inflow at a given time period is usually observed to possess serial correlation with past inflows, the latter should also be defined as state variables. A pseudo-code description of the procedure is presented below:

Initialize the FBF of the last stage βT+1(vT+1, aT) = 0 For each t = T, T-1, ..., 1 For each storage level vt = vt1 , ..., vtm , ... , vtM For each past inflow scenario at-1 = at1−1 , ..., atk−1 , ... , atK−1 For each inflow scenario of stage t conditioned on the past inflow scenario at = at1 ,

..., atl , ... , atL Solve the one-stage problem considering an initial storage level equal to vtm and the inflow scenario atl :

β tl ( vtm , atk−1 )=

Max



k =1

subject to



k =1

K

K

π dtk Etk + βt+1( vtl+1 , atl ) Etk ≤ E mt

17

vtl+1 (i ) = vtm (i ) − u t (i ) − st (i ) + atl (i ) +

∑ [u ( g ) + s ( g )]

g∈G ( i )

t

t

Constraints (8.b) – (8.d), (8.f), (8.g) End

Calculate the expected value of the revenues obtained across the conditioned inflow scenarios: βt( vtm , atk−1 ) =



L

l =1

(

p kl × β tl vtm , atk−1

)

End

Create a complete FBF βt(vt, at-1) for the previous stage by interpolating the values obtained above {βt( vtm , atk−1 ), m = 1, ..., M; k = 1, ..., K} End End

where pkl is the probability of inflow scenario l conditioned on the past inflow scenario k and βt(vt, at–1) is the expected future revenue from stage t until the end of the horizon for given values of initial storage vt and past inflow scenario at–1. In the SDP algorithm, for each stage (t = T, T-1, …, 1), storage level (0%, ..., 100%) and conditioned hydrological scenario, the bidding process for the hydro company is optimized. The SDP scheme is straightforward to implement and has been used for several years in some hydro-dominated countries for least-cost scheduling. It has been also applied to hydrobidding – [15] is one of the first works on modeling hydro operation in deregulated markets using DP. Later, references [16] and [17] presented a game-theoretical model for a hydro system with two hydro companies following an SDP algorithm in a multistage and stochastic framework.

18

However, the existence of various plants and the possible relationship of the current inflow scenario with that of several previous time stages create the need for additional “FOR” loops in the algorithm outlined above. As a consequence, the enumeration of all combinations of initial storage values and previous inflows in the SDP recursion leads to an exponential increase in the required computational effort, limiting its application to hydro companies owning one or two plants, as the studies carried out in [15], [16] and [17]. For illustrative purposes, if an agent owns 5 plants and the storage capacity and previous inflow of each plant were to be discretized into 20 levels, one would have to solve 2010 mixed integer programs.

5. Proposed methodology Because of the drawbacks of SDP, it becomes necessary to develop computationally tractable algorithms for practical applications. A possible approach, adopted in some countries, is to reduce system dimensionality by aggregating several reservoirs into one “equivalent reservoir” that should mimic the main characteristics of the cascade. In some cases, this is coupled with the use of partial dynamic programming schemes (typically, calculation of separate cost-to-go functions for each basin). 5.1. Existing multi-stage decomposition techniques for hydroscheduling

In [20], [21] and [18], an approach based on the analytical representation of a convex costto-go function, known as stochastic dual dynamic programming (SDDP) and based on multistage stochastic Benders decomposition, was proposed in the context of least-cost hydroscheduling. This algorithm constructs feasible dynamic programming policies using an outer approximation of a (convex) future cost function computed using Benders cuts. The policies defined by these cuts can be evaluated using simulation, and their performance measured against a lower bound on their expected cost. This provides a convergence criterion that may be applied to terminate the algorithm when the estimated cost of the candidate policy is close enough to its lower bound. This algorithm has been widely applied in several 19

countries ([22], [27], [28], [29]). The SDDP scheme does not require the explicit discretization of the whole state space and, as a consequence, it alleviates the computational requirements of the stochastic DP recursion. It has led to a number of related methods and applications ([30], [31]) that are based on the same essential idea, but seek to improve the method by exploiting the structure of particular applications. The apparent advantages of the algorithm make it natural to wonder about its application to the hydro bidding problem. The non-concavity of the revenue function (and, consequently, of the FBF), however, makes it impossible to approximate it by a set of linear inequalities. The next section suggests two simplifications to the problem that would make it amenable to the application of the methodology. 5.2. Concavization of the revenue function

A common assumption in related works ([4], [8]) that tackle the strategic bidding problem via complementarity models is a linear relationship between residual demand and the spot price. It can be shown that this hypothesis is equivalent to assuming that the operational cost of (price-takers) thermal generators can be represented by a quadratic function. In this case, problem (5) may be re-written as follows:

Min

cg2

Subject to

g=d–E

(9) (9.a)

Problem (9)-(9a) is the equivalent of problem (5) when considering the approximation of the piecewise linear supply curve of the price-taker generators by a quadratic function. By following simple algebraic steps, one arrives at the following expression of the spot price as a function of the amount of energy offered by the hydro price-maker agent:

20

πd = 2c·(d – E)

(10)

The revenue obtained by the agent as a function of its own quantity bid is now represented by the equation below, roughly sketched in Figure 5.1:

R(E) = πd · E = 2cdE – 2cE2

(11)

100

REVENUE [$]

80

60

40

20

0 E0 E1 E2 E3 E4 E5 E6 E7 E8 E9 E10 E11 E12 E13 E14 E15 E16 E17 E18 E19 E20

AMOUNT OF ENERGY OFFERED [MWh]

Figure 5.1 – Revenue function under the assumption of quadratic costs of thermal generators

Observe that the resulting revenue function is obviously concave and, thus, may be approximated by a piecewise linear function, as needed for the SDDP algorithm. The application of the SDDP algorithm might also be possible through an alternative simplification to the original case – where the total thermal operational cost is piecewise linear – by observing that the function that must be approximated by the algorithm is not exactly the one that relates revenue to the amount of energy left for later periods, but the one 21

that connects revenue to the state variables of the DP recursion. Therefore, we need to map the function “revenue vs. energy offered” into a “revenue vs. storage” function. This is straightforward: if we consider a simple example with one reservoir, this function may be obtained by observing that the amount of energy that may be produced at a certain stage depends on the initially available volume of water and on the inflows at the current stage. Thus, the relation between the revenue obtained by the agent and the storage level of the reservoir at the beginning of a time period will be given by a function such as the one in the figure below:

Figure 5.2 – Agent’s revenue vs. initial storage level

The function takes this shape because, even though it might be possible to produce more energy for certain levels of initial stored volume, it is more advantageous to restrict the amount of energy offered and store this extra energy in order to prevent spot prices from getting lower and negatively affecting its profits. These intervals correspond to the horizontal segments (in Figure 5.2) where the revenues are constant. From the point where a higher production – even with a lower spot price – becomes more attractive, the revenue function starts to increase again. Once the sum of the initial storage level and the inflow scenario are

22

enough to produce the amount of energy that corresponds to the global maximum revenue, the function becomes flat. This results on a much smoother function, that may be approximated by its piecewise linear concave hull. Again, that would allow the application of the SDDP methodology. 5.3. Application of Stochastic Dual Dynamic Programming

Independently of the adopted approach, one ends up with a piecewise linear concave hull of the revenue function, such as the one depicted below:

Figure 5.3 – Length and slope of the linear segments that approximate the revenue function Each one of the P segments of this approximation can be characterized by its slope r(p) and by the length of its projection on the E-axis, emax(p). The one-stage problem may now be formulated as: r ( p ) × et ( p )

Max



Subject to

(operating constraints)



P

(12)

p =1 t

J j =1

ρ(i ) ⋅ uT (i ) − ∑ p =1 et ( p ) = 0 P

(12.a)

23

et(p) ≤ etmax(p)

Vp

(12.b)

where constraint (12.a) links the physical production of the plants with the variables that are used to approximate the revenue function. For the sake of simplicity, the production coefficient in this work is assumed to be constant. It is well-known this is a non-linear factor. It can be, however, approximated by a piecewise linear function [42] without any implications to the proposed approach. However, we believe that even with the constant head assumption, our methodology still may be valuable in capturing the long term influences – in the form of the future benefit functions – of present decisions. We surely introduce an error into the problem by approximating the revenue function by its concave hull. The quality of the approximation depends, basically, on the difference between the marginal costs of consecutive price-taker agents (when sorted in ascending merit order of variable cost offer 1 ) – the smaller the difference between the marginal costs of consecutive price-taker agents, the smaller is the approximation error. Anyhow, we believe this is a valid trade-off since it allows us to apply a computationally-efficient methodology to determine the market-based operation policy. We also conjecture that in most real electricity markets there is not a big difference between the marginal costs of “consecutive” price-taker agents (usually thermal plants in hydro-based systems). This is the case of several Latin American countries and Turkey, as per the personal experience of the authors. An important limitation is the assumption that costs/bids for the price-taker agents are known throughout the horizon – although they may have different values at each time stage. This means that we “manually” pre-calculate the bidding curve for each stage based on demand and price-taker costs/bids and capacities. The main difference between the SDP procedure explained in Section 4 and the proposed

1

Price takers are represented by their capacities and marginal production costs. In addition, there is a generator whose capacity is infinite and whose marginal cost is equal to the deficit cost that represents the costs of energy shortages. 24

SDDP approach is related to the way each algorithm constructs the future benefit function when decomposing the multi-stage problem into a series of one-stage problems. While SDP requires the future revenue function to be evaluated at numerous points so that it can be reasonably interpolated (discrete approach), SDDP represents it as a piecewise linear function which can be constructed around a much smaller number of states. The approximation to the FBF can be constructed by observing that its slope around a certain point can be obtained from the solution of the single-stage problem. The hydro bidding problem of the last stage – in which the FBF is zero – is presented below for given values of initial storage of each reservoir and of the inflow scenario (the inflow scenarios and discretized values for the stored volumes in reservoirs used in the problem will become clearer when we present the whole structure of the algorithm in Section 5.4). βlT (vm,T, ak,T-1) = Max

Max



P

r ( p ) × et ( p ) (13)

p =1 t

Subject to

v l ,T +1 (i ) + uT (i ) + sT (i ) − ∑ j∈G ( i ) (uT ( j ) + sT ( j ) ) = v m,T (i ) + a l ,T (i )

(13.a)

v l ,T +1 (i ) ≤ v(i )

(13.b)

uT (i ) ≤ u (i )

(13.c)



J j =1

ρ(i ) ⋅ uT (i ) − ∑ p =1 et ( p ) = 0

et(p) ≤ etmax(p)

P

Vp

(13.d) (13.e)

where: vm,T(i)

mth discretized value for the initial stored volume in reservoir i at stage T

ak,T–1

kth past inflow scenario

ρ(i)

production coefficient of reservoir i 25

al,T(i)

lth inflow scenario at stage T in reservoir i

v( i )

maximum stored volume of reservoir i

u( i )

maximum turbined volume of reservoir i

uT(i)

turbined volume in reservoir i at stage T

sT(i)

spilled volume in reservoir i at stage T

vl,T+1

initial (final) volume at the beginning (end) of stage T+1 (T)

From the theory of linear programming, there is a set of dual variables associated with the constraints of problem (13) at its optimal solution which represent the partial derivatives of the objective function with respect to the right-hand sides of these constraints. In particular, the dual variables associated with the water balance equations (13.a) can be utilized to calculate the partial derivatives with respect to the state variables, as shown below for the case of the volume stored in reservoir i at the beginning of stage T: ∂β Tl = π vlmT (i ) m ∂vT (i )

(9)

Likewise, the rate of variation of the objective function with respect to the current stage inflow is given by the same simplex multiplier: ∂β Tl

∂a (i ) l T

= π vlmT (i )

( 10 )

The partial derivatives of the objective function with respect to the inflows of the previous stages – which are actually the state variables – can be obtained from by the application of the chain rule, as in the expression below: ∂β Tl ∂a

l T −1

(i )

=

∂β Tl

∂a (i ) l T

×

∂a T (i ) ∂aT (i ) = π vlmT (i )× ∂aT −1 (i ) ∂aT −1 (i )

(11)

26

Note that the derivative of the current stage inflow with respect to the ones in previous time periods may be calculated from the parameters of any time-series model applied to the process. Thus, when solving a one-stage problem, one obtains a point of the FBF and the partial derivatives with respect to the state variables. In that way, a hyperplane that is tangent to the FBF of the previous stage can be generated. This hyperplane corresponds to a linear approximation of the expected value of the benefits obtained in the neighborhood of the initial state considered (vmT, akT-1) for a given inflow scenario and is mathematically represented by: ∂β Tl ∂β Tl m β = β l ,t (v m,T , a m ,T −1 ) + m × vT − vT + l × aT −1 − aTl −1 ∂vT ∂aT −1

(

)

(

)

(12)

By repeating this process for different initial states, a set of hyperplanes is generated. This set approximates the FBF and can be incorporated into the problems of the previous stage as a series of linear constraints. 5.4. The SDDP Algorithm

The SDDP algorithm works iteratively and is basically composed of two procedures: the (i) forward and (ii) backward steps. From the (known) initial state, the forward simulation produces a sample of “interesting” states for each stage t = 1,…T by simulating the hydro bidding for a set of different hydrological scenarios. In the backward recursion the algorithm goes back in time (repeat for t = T, ..., 1) applying stochastic DP to the states selected in step (1) and generating Benders cuts for the revenue-to-go function (or “future benefit function”) in the previous stage. These two procedures can be schematically depicted as: i)

Forward Simulation:

define a set of inflow scenarios at = {atq, q = 1, ..., Q} for all time periods

27

for each inflow scenario at = {atq, q = 1, ..., Q} initialize storage value for stage 1 as vtq = v1 (known value) for t = 1, 2, …, T solve the one-stage scheduling problem for initial storage vtq and inflow atq: Max Σp=1...P rt(p).et(p) + βt+1 subject to vt+1(i) = vtq(i) – ut(i) – st(i) + atq(i) vt+1(i) ≤ vmax(i)

Vi

ut+1(i) ≤ umax(i)

Vi

et(p) ≤ emax(p)

Vp

βt+1 ≥ ϕnt+1vt+1 + δnt+1

n = 1, ...,Nt

next calculate total expected revenue zq for scenario q as the sum of the revenues obtained at each stage along the horizon: zq = Σt=1…T Σp=1...P rt(p).et(p)

next where Nt is the number of linear inequalities generated for the FBF of the given time period up to the current iteration and ϕnt+1 and δnt+1 are the vectors of coefficients that characterize the hyperplanes. ii)

Backward recursion:

In the first iteration, initialize FBF of the last stage as zero for t = T, T-1, ..., 1 for each storage value vt = {vtm, m = 1, ..., M}

28

for each inflow scenario at = {atk, k = 1, ..., K} solve the one-stage scheduling problem for initial storage vtm and inflow atk: βtk(vtm) = Max Σp=1...P rt(p).et(p) + βt+1 subject to

(simplex multiplier)

vt+1(i) = vtm(i) – ut(i) – st(i) + atk(i) vt+1(i) ≤ vmax(i)

Vi

ut+1(i) ≤ umax(i)

Vi

et (p) ≤ emax(p)

Vp

βt+1 ≥ ϕnt+1vt+1 + δnt+1

n = 1, ...,Nt

πtk

next calculate the coefficients and constant terms for the mth linear segment of the FBF of the previous stage and add them to the set Nt – 1: ϕtm = Σk=1…K pk×πtk and δtm =Σk=1…K pk×αtk(vtm) – ϕtm×vtm

next next

Once the backward recursion pass is complete, the solution of the first stage problem with the approximated FBF provides an upper bound on the optimal solution to the problem. On the other hand, the expected total revenue across all scenarios used in the simulation gives a lower bound on the optimal solution as well. Since we use a Monte Carlo simulation procedure, we can also calculate the standard deviation of this estimate to construct a confidence interval against which we can compare the previously calculated upper bound and devise a convergence criterion: if the upper bound lies inside the confidence interval, the algorithm terminates; else, another iteration is run starting with a new set of storage values. The structure of the complete algorithm is presented below – where η indicates the number of

29

standard deviations used in the calculation of the confidence interval: Do Until zUPP ∈ (zLOW – η.σ, zLOW – η.σ) Run Simulation Store initial storage values of each stage to be used in the Backward Recursion phase Calculate zLOW = (1/Q) Σt=1…Q zq Calculate σ = ((1/(Q – 1)) Σt=1…Q (zq – zLOW)2)1/2 Run Backward Recursion Calculate zUPP Next iteration

As shown above, at each iteration the approximation to the FBF is refined and upper and lower bounds to the problem solution are calculated until the convergence criteria is met. Detailed analytics can be found in [20] and [21].

6. Case Study This section applies the proposed methodology to a case study from the El Salvador system.

6.1. Characteristics of the power system

El Salvador is a country in Central America whose power system has a total installed capacity of 1425MW, of which 38% (554 MW) are represented by 5 hydro plants: Guajoyo, Cerrón Grande, 5 de Noviembre, 15 de Septiembre and Chaparral. The following table presents some of their basic operational data: TABLE I BASIC OPERATIONAL DETAILS OF THE HYDRO PLANTS IN EL SALVADOR Capacity Maximum Production Turbine Hydro plant [MW] storage coefficient capacity 30

Guajoyo Cerron Grande 5 de Noviembre 15 de Septiembre El Chaparral

90 135 95 170 64

[hm3] 560 2180 72 346 0

[MW/m3/s] 2.09 0.53 0.41 0.29 0.53

[m3/s] 43 255 231 593 121

The figure below displays the schematic representation of the hydro basin in which the plants are located:

Guajoyo

Cerron Gr ande

5 de Noviembre

El Chaparral

15 de Septiembre

Figure 6.1– Hydro basin in which the hydro plants are located

The remaining generation capacity is composed by a set of 17 thermal plants with operating costs ranging from 20 US$/MWh to 200 US$/MWh and capacities ranging from 20 to 150 MW. The system used to have a bidding-scheme for system dispatch and spot price formation, which in 2008 was substituted by a centralized least-cost dispatch. 6.2. Description of the case study

We have simulated the system behavior assuming a competitive market where the 5 hydro plants belong to the same company and is the price maker agent. A study horizon of 5 years in monthly steps was considered. Demand was considered to be constant at 700 average MW per year. In this simulation, we applied our proposed methodology, which consists of two steps: (i) approximate the revenue functions at each month by its concave hull; (ii)

31

approximate the concave hulls by the SDDP algorithm described in Section 5.4. The convergence tolerance was defined by setting η = 1.96, i.e. the lower bound must fall within the 95% confidence interval for the upper bound. 6.3. Step 1: concavization of the revenue functions

Having the information regarding the thermal price-taker agents (marginal costs and capacities) and the demand, it is possible to build the residual demand curves that relate the amount of energy offered at zero-price by the price-maker hydro company to the spot price and to the obtained revenue, as presented in Figure 6.2 and Figure 6.3, respectively. Observe that, in our case, this function is the same for all months because demand is assumed to be constant and there is one price-maker. If demand was not constant, this function would vary per time stage. In Figure 6.3, the approximation of the curve – up to the point where it achieves its maximum value – by its concave hull is also shown. 160 140

Spot price [$/MWh]

120 100 80 60 40 20 0 0

100

200

300

400

500

600

700

E: Amount of energy offered [average MW]

Figure 6.2 – Spot price vs. energy offered by the hydro agent

32

35,000

30,000

R(E): Revenue [$]

25,000

20,000

15,000

10,000

5,000

0 0

100

200

300

400

500

600

700

E: Amount of energy offered [average MW]

Figure 6.3 – Revenue vs. energy offered by the hydro agent and its concave approximation This concave hull will be then approximated by the SDDP algorithm. As we have mentioned earlier, the quality of the approximation depends, basically, on the difference between the marginal costs of consecutive price-taker agents (when sorted in ascending merit order of marginal cost). As seen in Figure 6.3, for the El Salvador case, this difference is quite small, which ensures a good approximation. 6.4. Application of SDDP

We have applied the algorithm presented in Section 5.4 considering 50 inflow scenarios. For comparison purposes, the results of the model based on the maximization of the hydro agent’s revenue are plotted along those of a model that aims at minimizing total operating costs. Seven iterations were run and CPU time was 10 minutes for the revenue maximization model in a Pentium IV 3 Ghz 1 Mb RAM memory. 6.5. Results

The following graph shows the average total hydro generation along with the inflow energy (water inflow transformed into energy):

33

Figure 6.4 – Average total hydro generation and inflow energy It can be seen that both models have similar behaviors regarding the transfer of energy from wet seasons to dry seasons. However, the revenue-maximization model determines that the amount of energy produced by the hydro agent is never above 220GWh in any period – which is equivalent to approximately 300 average MW during a month and corresponds to the value where revenue function reaches its maximum (see Figure 6.3). The effects of the different hydro generation profiles on spot prices can be observed in Figure 6.5. While the cost-minimization model is able to significantly reduce spot prices during wet seasons, its revenue-maximization counterpart keeps them high by lowering the amount of energy produced and thus requiring the dispatch of more expensive thermal plants.

Figure 6.5 – Average spot price Finally, we observe the evolution of the average revenue obtained by the hydro agent in 34

each period. It can be seen that by exercising its market power, the company obtains an obviously larger income stream.

Figure 6-1 – Average revenue obtained by the hydro agent

Finally, as it can also be observed (Fig 6.4), the blue (cost minimization) and red (revenue maximization) lines or do not have the same area. The “excess” hydro energy (that which would cause the spot price and revenue to decrease) is spilled because we set no penalty for spillage. This could obviously be changed with no impacts on the proposed solution methodology. This same behavior was observed in other papers such as [16].

7. Conclusions and future work While there have been numerous articles on the literature regarding the determination of optimal bidding strategies for price-maker thermal and small hydro plants, substantially less effort was devoted to the case of large hydro portfolios. The need to account for stochastic inflows and the combinatorial explosion of traditional solution methods turn this into a very hard problem. Our work is intended to partially contribute to these developments by presenting a detailed overview of current methodologies for the strategic bidding problem of a price-maker hydropower based company, taking into account several hydro plants, timecoupling and multiple inflow scenarios. The limitations involved with the existing solution

35

procedures motivated us to propose an approach based on stochastic dual dynamic programming (SDDP), which has been previously applied to the least-cost hydrothermal scheduling problem. Since the technique requires the problem to be concave, a piecewise linear approximation of the revenue function is proposed. This allowed us to adapt the SDDP algorithm and solve larger instances of the problem that would otherwise be computationally infeasible – as presented in the case study of the El Salvador system. Of course, this benefit is obtained at the cost of an approximation of the future benefit function. Since our model deals with medium to long term horizons and the bidding process actually occurs on a daily basis, we envision its application as part of a sequence of optimization runs, as done in the traditional hydroscheduling procedures. First, the proposed model may be run in monthly time steps with a horizon of a few years so that long term effects are captured by the generated approximations of the future benefit functions. After that, the model could be run again with a horizon comprising a few months divided into weekly time steps and the final-stage future benefit function taken from the appropriate time period of the previous run. Finally, these results could then be refined by a short time model, which represents the exact non-convex revenue function in a MIP framework and represent the long-term effects in two possible ways: by setting generation targets for each plant for the end of the horizon or by incorporating the generated linear constraints that represent the future benefit function. Both results are obtained from the long-term model. In terms of future work, we highlight three possible areas for research: the first relates to the representation of a market equilibrium resulting from the interaction of several pricemakers agents. The model discussed in this work considers a single price-maker company. A natural extensions is then the modeling of a dynamic interaction between all price-makers in a game-theoretic framework when more than one price-maker is considered. We suggest to tackle this by calculating a Nash equilibrium by means of an iterative algorithm following a

36

successive approximation approach. The second research line consists in overcoming the assumption that bids of price-takers are deterministic. In this case, we suggest representing these bids by means of scenarios. Finally, a third extension is related to the fact that some price-takers might adjust their bidding strategies as a function of the price-maker’s behavior. This is more important to price-taker hydro companies, whose bidding strategy (even if costbased) is also coupled in time and handled by means of SDP algorithms.

8. References [1]

Hunt, S. “Making Competition Work in Electricity”, Published by John Wiley & Sons, 2003.

[2]

Gross, G. and Finlay, D., “Generation supply bidding in perfectly competitive electricity markets”, Computational & Mathematical Organizations Theory,Vol. 6, pp. 83-98, 2000.

[3]

Hobbs B. H., Metzler C. B. and Pang J. S., “Strategy Gaming Analysis for Electric Power Systems: An MPEC Approach”. IEEE Trans. Power Syst., v. 15, pp. 638-645, May 2000.

[4]

Hobbs, B. F. and Helman, U., “Complementarity-based equilibrium modeling for electric power markets”, in D. Bunn, ed., Modeling Prices in Competitive Electricity Markets, J. Wiley, 2005.

[5]

Conejo A. J., Contreras J., Arroyo J. M. and Torre S., “Optimal response of an oligopolistic generating company to a competitive pool-based electric power market”, IEEE Trans. Power Syst, Vol. 17, No.2, 2002.

37

[6]

Ramos, A., Ventosa, M. and Rivier, M., “Modeling Competition in Electric Energy Markets by Equilibrium Constraints” Utilities Policy, Vol. 7 Issue 4 pp 233-242, Apr. 1999.

[7]

Weber, J. and Overbye T., “An individual welfare maximization algorithm for electricity markets”, IEEE Trans. Power Syst., Vol. 17, No.3, 2002

[8]

Bushnell, J., “A mixed complementarity model of hydrothermal electricity competition in the western United States”, Operations research, Vol. 51, No. 1, pp. 80-93, 2003.

[9]

I. Otero-Novas, C. Meseguer, and J. J. Alba, “An iterative procedure for modeling strategic behavior in competitive generation markets,” in Proc. 13th Power Syst. Conf., Trondheim, Norway, July 1999.

[10] Villar, J.; Rudnick, H.; “Hydrothermal market simulator using game theory: assessment of market power”, IEEE Transactions on Power Systems. Volume 18, Issue 1, Feb. 2003 Page(s):91 – 98. [11] Torre, S., Arroyo, J. M., Conejo, A. J. and Contreras, J., “Price Maker Self Scheduling in a Pool-Based Electricity Market: A MILP Approach”, IEEE Trans. Power Syst, Vol. 17, No.4, 2002. [12] Baillo, A., Ventosa, M., Rivier M. and Ramos, A., “Optimal Offering Strategies for Generation Companies Operating in Electricity Spot Markets”, IEEE Trans. Power Syst, Vol. 19, No.2, 2004. [13] M. V. Pereira, S. Granville, M. Fampa, R. Dix, L. A. Barroso. “Strategic Bidding Under Uncertainty: A Binary Expansion Approach”. IEEE Trans. on Power Syst., Vol. 20, No 1, 2004.

38

[14] Bakirtzis, A.G.; Ziogos, N.P.; Tellidou, A.C.; Bakirtzis, G.A.; “Electricity Producer Offering Strategies in Day-Ahead Energy Market With Step-Wise Offers”,

IEEE

Transactions on Power Systems, Volume 22, Issue 4, Nov. 2007 Page(s):1804 - 1818 . [15] Read, E.G. and Scott, T.J. Modelling Hydro Reservoir Operation in a Deregulated Electricity Sector International Transactions in Operations Research 3, 3-4, 243253,1996. [16] Kelman, R., Barroso, L. A. and Pereira M. V., “Market power assessment in hydrothermal systems”, IEEE Trans. Power Syst, Vol. 16, No.3, 2001. [17] Barquín, J., García-González J., Ubeda, J.R. “Water value in competitive markets: Dynamic programming and game theory”, 6th International Conference on Probabilistic Methods Applied to Power Systems., 2000, Funchal, Madeira (Portugal). [18] Birge, J. R. . “Decomposition and Partitioning Methods for Multi-Stage Stochastic Linear Programs.” Operations Research 33 (1985), 989–1007. [19] Gassmann, H. I.. “MSLiP: A Computer Code for the Multistage Stochastic Linear Programming Problem.” Mathematical Programming 47 (1990), 407–423. [20] Pereira, M. V., Pinto and L. M. V. G., “Multi stage stochastic optimization applied to energy planning”, Mathem.Programming, 52, 359-375, 1991. [21] Pereira, M.V.; Pinto, L. M. V. G. “Stochastic Optimization of a Multireservoir Hydroelectric System: A Decomposition Approach.” Water Resources Research, v.21, n.6, pp.779-792, June 1985. [22] Granville, S., Oliveira, G. C., Thomé, L. M., Campodónico, N., Latorre, M., Pereira, M. V. and Barroso L. A., “Stochastic optimization of transmission constrained and large scale Hydrothermal Systems in a Competitive Framework”, Proceedings of the IEEE General Meeting, Toronto, 2003.

39

[23] J Bushnell et al, “An International Comparison of Models for Measuring Market Power in Electricity”, Energy Modeling Forum, Stanford University, March, 1999 [24] Borenstein, S., J. Bushnell and C. Knittel “Market Power in Electricity Markets: Beyond Concentration Measures” The Energy Journal 20 (4), 1999. [25] J. García-González, J.Barquín., "Self-unit Commitment of thermal units in a competitive electricity market", Proceedings of IEEE PES Summer Meeting, Vol 4, pp: 2278-2283, Seattle (USA), July 2000 [26] Birge, J.R., and Louveaux, F. Introduction to Stochastic Programming, Springer Verlag, New York, 1997. [27] Gorenstin, B.G.

Campodonico, N.M.

da Costa, J.P.

Pereira, M.V.F. , “Stochastic

optimization of a hydro-thermal system including network constraints”

IEEE

Transactions on Ppwer Syst., Vol 7, No.; 2, page(s): 791-797, 1992. [28] Rotting, T.A.

Gjelsvik, A., “Stochastic dual dynamic programming for seasonal

scheduling in the Norwegian power system”, IEEE Transactions on Power Syst., Vol 7, No 2 page(s): 273-279, 1992. [29] Halliburton, T.S., “An Optimal Hydrothermal Planning Model for the New Zealand Power system”, Australian Journal of Electrical & Electronics Engineering, Volume 1 Issue 3 (2004). [30] Read, E.G. George, J.A and MacGregor, A.D. Stochastic dual dynamic programming with lagged variables. Proceedings of the 30th Annual Conference of the O.R. Society of N.Z. August 1994 145-53 [31] Donohue, C. J., and Birge, J.R. “The Abridged nested Decomposition Method for Multistage Stochastic Linear Programs with Relatively Complete Recourse”, Algorithmic Operations Research, 1, 20-30, 2006.

40

[32] Wallace, S.; Fleten, S-E; “Stochastic programming in energy”, In "Stochastic Programming", A. Ruszczynksi and A. Shapiro (eds), Vol. 10 in the series Handbooks in Operations Research and Management Science, North-Holland, 2003. [33] Gjelsvik A., Belsnes M., Haugstad A. (1998). An algorithm for stochastic medium-term hydrothermal scheduling under spot price uncertainty. Proceedings of 13th Power Systems Computation Conference, Trondheim, Norway. [34] Gjelsvik A., Wallace S.W. (1996). Methods for stochastic medium-term scheduling in hydro-dominated power systems. EFI TR A4438, Norwegian Electric Power Research Institute, Trondheim. [35] Flatabo, N., Haugstad A, Mo B., Fosso O. (1998). Short and Medium-term Gerneration Scheduling in the Norwegian Hydro System under a Competitive Power Market Structure, Proceedings of EPSOM Conference, Zurich. [36] Mo, B.; Gjelsvik, A.; Grundt, A.; Karesen, K. (2001): Optimisation of hydropower operation in a liberalised market with focus on price modelling. Power Tech Proceedings, Porto, Portugal. [37] Cabero J, Baillo A, Cerisola S, Ventosa M, Garcıa-Alcalde A, Peran F, et al. A medium-term integrated risk management model for hydrothermal generation company. IEEE Trans Power Syst 2005;20(3):1379–88. [38] Fleten S.-E, Wallace S.W., Ziemba, W.T. (2002): Hedging electricty portfolios via stochastic programming. In: Decision Making Under Uncertainty: Energy and Power, Greengard, C., and A. Ruszczynski (ed.), 71-94. [39] J. García-González, E. Parrilla, A. Mateo, "Risk-averse profit-based optimal scheduling of a hydro-chain in the day-ahead electricity market," European Journal of Operational Research. vol. 181, no. 3, pp. 1354-1369, September 2007

41

[40] Pritchard, G., Zakeri, G. (2003): Market offering strategies for hydroelectric generators. Operations Research 51, 602-612. [41] L. Hongling, J.Chuanwen, Z.Yan, A review on risk-constrained hydropower scheduling in deregulated power market, Renewable and Sustainable Energy Reviews 12 (2008) 1465–1475. [42] Diniz, A. L.; Maceira, M. E. P., “A Four-Dimensional Model of Hydro Generation for the Short-Term Hydrothermal Dispatch Problem Considering Head and Spillage Effects”, IEEE Transactions on Power Systems, Volume 23, Issue 3, Aug. 2008 Page(s):1298 – 1308

42