Pricing Energy and Ancillary Services in a Day-Ahead Market for a Price-Taker Hydro Generating Company Using a Risk-Constrained Approach

This paper analyzes a price-taker hydro generating company which participates simultaneously in day-ahead energy and ancillary services markets. An approach for deriving marginal cost curves for energy and ancillary services is proposed, taking into consideration price uncertainty and opportunity cost of water, which can later be used to determine hourly bid curves. The proposed approach combines an hourly conditional value-at-risk, probability of occurrence of automatic generation control states and an opportunity cost of water to determine energy and ancillary services marginal cost curves. The proposed approach is in a linear constraint form and is easy to implement in optimization problems. A stochastic model of the hydro-economic river basin is presented, based on the actual Vinodol hydropower system in Croatia, with a complex three-dimensional relationship between the power produced, the discharged water, and the head of associated reservoir.


Introduction
This paper formulates a hydro-economic river basin model (HERBM) of a price-taking hydro generating company (GENCO, Zagreb, Croatia) which participates in a simultaneous energy and ancillary services day-ahead markets (DAMs) [1][2][3][4][5].In order to successfully bid on a market, a hydro

OPEN ACCESS
GENCO first needs to determine energy and ancillary services marginal cost curves.Since the direct cost of energy generation for a hydro GENCO is virtually zero compared to a fossil fuel GENCOs, consequently creating a marginal cost curve is more complicated compared with finding a generation cost curve first derivative for a fossil fuel GENCO.This paper will tackle this problem, and propose easy to implement linear formulation in a form of a constraint.In this paper an opportunity cost of water is viewed as a marginal cost of water.The idea behind the opportunity cost of water is that generating one MWh in a given hour means not being able to generate one MWh in some future hour.Therefore, for a hydro GENCO, the foregone earnings from future generation constitute the opportunity cost of current energy generation.A similar study is done in [6] where a fixed head of a hydro power plant (HPP) is considered.In one of the fundamental works, a marginal cost of water is assumed to be constant over the planning horizon [7,8].For a more realistic approach, a water shadow price should change over the planning horizon every time the reservoir capacity constraints are reached [9,10].This was the case in this research.A lot of work has been done in a field of bidding curve optimization models and a review of these models can be found in [11], with some early works in [12][13][14].
Modeling of ancillary services is specifically addressed in this research using simple linear formulation where automatic generation control (AGC) state (ramp up, ramp down, no ramp) is considered as a random variable.Detailed ancillary services models are important in order to encourage competitiveness and introduction of new ancillary services providers which consequently results in higher penetration of renewable energy sources (RES) particularly wind [15][16][17].
In a short-term planning most of the parameters are usually considered as known and the resulting models are deterministic [18,19].In this paper a stochastic short-term planning model is presented, which enables implementation of a risk measure.It is particularly convenient to use a conditional-value-at-risk (CVaR), a relatively new risk measure introduced in [20].CVaR is a coherent measure of risk [21,22] which is important for the successful representation of real life occurring risks and for implementation in convex optimization problems [23].This paper analyzes the price-taker hydro GENCO which participates simultaneously in energy and ancillary services DAMs.This paper contributes to the previous works [7][8][9][10] with the approach which combines innovative usage of the hourly CVaR, the probability of occurrence of AGC states and the opportunity cost of water to determine the energy and ancillary services marginal cost curves and consequently to determine a bid prices for energy and ancillary services.The proposed approach is in a linear constraint form and is easy to implement in optimization problems.
Relationships between the head of the associated reservoir(s), the discharged water and the generated power is considered.Models describing these relationships are available in [24][25][26][27][28][29].Therefore, in Section 2, the cascade hydropower system (HPS) is carefully modeled as HERBM using the mixed-integer linear programming (MILP) formulation.Afterword's the ancillary services model is formulated and the risk-constrained approach for pricing energy and ancillary services is introduced.Section 3 presents and analyzes the results on the case of HPS Vinodol (Croatia) and in Section 4 conclusions are given.

Problem Description and Formulation
In this chapter, hourly water marginal costs are determined.The profit maximization objective function is formulated.Hourly CVaR as a risk measure is formulated.Ancillary services model is formulated with AGC state as a discrete random variable.Probability of occurrence of AGC states are later used as a rule for redistributing the water marginal cost along the energy, the regulation and the 10 min spinning reserve.Afterwards the energy, regulation and 10 min spinning reserve marginal costs are used for determining an hourly bidding curves.

Hydro Constraints
The continuity equation of the hydro reservoirs is formulated as Equations ( 1) and ( 2).The variable defined after semicolon denotes the associated dual variable: For unit consistency, it should be noted the time periods of 1 h are considered.The τ is the time delay for water transportation between reservoirs i and j.

Head Dependent Power Output
The performance curves (Figure 1) are modeled through a piecewise linear formulation of a Hill chart [30].Non-concavities are modeled through the use of binary variables.Figure 1a shows three-piece linear performance curves.The slope ρ is defined by HPP conversion capabilities (MWs/m 3 ).The five performance curves are associated for the five water contents intervals.Activation of the corresponding performance curve is done by the approach presented in [31].Unlike [32,33] the operation of HPP is not restricted to the local best efficiency points.Activation of the performance curves is modeled with Equations (3)- (6).
The performance curves are activated according to the four binary variables , , and in Equations ( 3)-(6) using the average value of stored water in certain time step t Equation (3) instead of using the final reservoir water content of the time step t [26]: Depending on the binary variables , , and combination the corresponding performance curve will activate.For instance if the binary variables are set to 1, 0, 0, 0 the performance curve Equation ( 8) is activated, if the binary variables are 1, 1, 1, 0 then the curve Equation ( 10) is activated and, etc. Formulation of the performance curves is shown in Equations ( 7)-( 13).

Day-Ahead Market
In this work the hydro GENCO participates in the simultaneous energy and ancillary service DAM.Markets are considered competitive enough to assume the analyzed GENCO is a price-taker.The 10 min spinning reserve and regulation markets are considered as the ancillary service markets.In DAM considered here the independent system operator (ISO) will commit GENCOs in a similar way as ISO did in the vertically integrated structure using security constrained unit commitment (SCUC).Suppliers submit their bids to supply the forecasted daily inelastic demand and ISO uses the bid-in costs submitted by GENCOs for each generating unit (multi-part bids that reflect start-up costs, minimum generation costs, capacity offered, running costs, etc.) to minimize cost of operation and determine which units will be dispatched in how many hours and calculate the corresponding market clearing prices (MCPs) while the system security is retained The prices of energy and ancillary services are set at the shadow price of the market clearing constraints for the corresponding product and as a result, generating units are dispatched to maximize their profits given the prices of energy and ancillary services.Such market design is incentive for a competitive generator to select as a weakly dominating strategy bidding the actual generation cost and operation parameters.DAM market modeled here is assumed to be similar to New York DAM (New York Independent System Operator-NYISO) [34].

Ancillary Services
Besides participating in energy DAM, GENCO participates in the 10 min spinning reserve DAM and the regulation DAM.Ancillary services markets are country-dependent and a review of these markets designs is available in [35].An HPP parameter important for modeling frequency regulation services is the maximum sustain ramp rate parameter (MSR [MW/min]) which defines how much unit can increase (MSR up) or decrease (MSR down) its output in time period of minute.The regulation and spinning reserve services are modeled based on the approach presented in [36], with the improvements Equations ( 16)-( 21) regarding probability of occurrence of particular AGC state.When GENCO participates in the regulation DAM, the following AGC states may occur: • Regulation-up: HPP power output increases.GENCO earns on available regulation capacity ( , ω) priced at π ( , ω) and on electricity produced for regulation service ( , ω) priced at π ( , ω) in a time period t and a price scenario ω.Probability of being in this state is ( , ω).
• Regulation-down: HPP power output decreases.GENCO earns on available regulation capacity ( , ω) priced at π ( , ω), GENCO does not repay nor gets repaid for energy decreased priced at π ( , ω) in a time period t and a price scenario ω.Probability of being in this state is ( , ω).
• No-regulation: There is no change in HPP power output.GENCO earns on available regulation capacity ( , ω) priced at π (t, ω) in a time period t and a price scenario ω.Probability of being in this state is 1 − ( , ω) − ( , ω).
When GENCO commits in the 10 min spinning reserve DAM the following AGC states may occur: • Reserve-up: HPP power output increases.GENCO earns on 10 min spinning reserve capacity ( , ω) available for increase of power output priced at π ( , ω) in time period t and price scenario ω.Also GENCO earns on electricity produced for spinning reserve service ( ) priced at π ( ) in a time period t and a price scenario ω.Probability of being in this state is ( , ω).
• Reserve-neutral: There is no change in HPP power output.GENCO earns on 10 min spinning reserve capacity ( , ω) available for increase of power output priced at π (t, ω) in a time period t and a price scenario ω.Probability of being in this state is 1 − ( , ω).
The operating range of HPP when committed in the regulation and the 10 min spinning reserve DAM is defined with Equation ( 14) and is depicted in Figure 1b.A review of ancillary services models is available in [36][37][38][39]: For the regulation capacity ( , , ω), the value indicates the ability to move both up and down and requires both headroom and footroom, while provision of the spinning reserve ( , , ω) requires headroom but not footroom, meaning HPP cannot provide energy at its maximum output level but can run at minimum output level.Practical implementation of Equation ( 14) is shown in Equation (15) and is depicted in Figure 1: Expected electricity produced when committed in the simultaneous energy, regulation and spinning reserve DAM in a time period t and a price scenario ω Equation ( 16).

Random Variable: Energy and Ancillary Services MCPs
The set Ω will represent future states of knowledge, and the single element ω will define a one future state, meaning ω ∈ Ω.The set Ω will have a mathematical structure of probability space with probability measure p for comparing the likelihood of the future states .
The discrete random variable π as shown in Equation ( 25) will represent random daily price series, and is defined with Π ∶ Ω → , where ⊆ ℝ , since probability won't be defined on each element of ℝ vector space.The probability of particular daily price scenario of market m Π (ω) ∈ is defined as Equation ( 26): Finite number of future states are observed, ω ∈Ω ⊆  .The random variable Equation ( 25) is defined for energy and ancillary service DAM.Moreover, it can be shown that the distribution of hourly prices is approximately lognormal [33].Thus resulting random variable m Π is defined for the each time step t separately Equation ( 27):

Risk Measure
An easy way to incorporate a risk into linear model is to use CVaR [16] as a measure of risk.β-CVaR is defined as an average loss of a 1-β worst scenarios, thus called tail loss (Figure 2a).α-CVaR (Figure 2b) is an average profit in worst 1-α (i.e., 15%) scenarios and α-VaR is minimal profit which GENCO can expect in rest α (i.e., 85%) scenarios.CVaR is a coherent measure of risk, meaning CVaR has properties of convexity, monotonicity, closedness, positive homogeneity, subadditivity and translation invariance [20][21][22][23].In the next section α-VaR and α-CVaR are formulated using the profit PDF as depicted in Figure 2b.
Both α-VaR and α-CVaR can be calculated by solving an elementary optimization problem of convex type in one dimension.For this purpose the special function Equation ( 28) is formulated: where, When maximizing Equation ( 29) over all variables defined in the model, Equations ( 30) and ( 31) are obtained: α-CVaR = (Variables * , ζ * ) (30) where the Variables denotes all the variables defined in nomenclature.The variable ζ is shown separately only to point out its importance.In Equations ( 28)-( 31) a risk measure is expressed as a daily value.Generally, when the daily CVaR is used, the risk will average out throughout the day, and this is not desirable as, especially in a case of a major contingency in particular hour this can result in unpredictable values of the daily CVaR.Consequently, the hourly α-VaR(t) and the hourly α-CVaR(t) are formulated using the hourly special function (Variables , ζ( )) over all (Variables , ζ( )) ∈ ℝ .
× ℝ.Since objective function is already defined, the special function is thus brought in the optimization problem in a form of the constraint Equation ( 32): where the denotes all the variables defined in nomenclature in a time step t.The variable ζ( ) is the hourly value.
When implementing the constraint Equation ( 32), the risk is "shaped" using the profit tolerance minCVaR( ) ∀ ∈ .When the set of profit tolerances minCVaR( , ), ∀ ∈ is introduced in the optimization model and the requirement Equation ( 33) is valid, then Equations ( 34)-( 36) is used to shape the risk with CVaR (Figure 3): The linear formulation of Equation (32), defined using indexed assignment is used in the model Equations ( 34)-( 36): Figure 3. Algorithm for risk-shaping with CVaR.

Objective Function
GENCO's objective function is maximal expected profit (37).GENCO can justify optimization of the expected profit by the law of large numbers because it gives an optimal decision on average.However, because of the variability of the input data, the average of the first few results may be very bad and it does not help that the decisions were optimal on average.For these reasons, quantitative models of risk and risk aversion are needed.
GENCO's profit in a time step t of a price scenario ω (37): Objective function is an expected profit of the planning horizon Equation ( 38): assignments.The presented results were obtained on 3.4 GHz based processor with 8 GB RAM using CPLEX (International Business Machines Corporation, New York, NY, USA) under General Algebraic Modeling System (GAMS) [40].However, in general terms, using the dual information from a mixed integer programming problem may not capture all the terms in the problem related to the binary variables.The Case Study in this research is considered fairly accurate, since optimality tolerance is set to 0%.This means that the solver will stop when it finds a feasible integer solution within 0% of the global optimum.
Figure 5.The approach for pricing energy and ancillary services.

Example and Results
For the needs of this research HPS Vinodol (Croatia) is modeled, which consists of: four reservoirs, one pump station (PS), two pump storage hydropower plants (PSP) and one HPP (Table 1) (Figure 6).Real HPS Vinodol is already used for frequency regulation ancillary services.Price-taker model is considered reasonably accurate since of relatively small size of HPS Vinodol compared to other participants on DAM assumed in this research.The prices are formulated according to Equation (27) where values are normally distributed.Mean and standard deviation of prices in each time step t are statistically obtained according to Equations ( 47) and (48): The energy and ancillary services prices are obtained from New York DAM [34], from 3 January 2011 to 5 December 2011 (337 days are observed) and descriptive statistics are calculated (Table 2).Using pseudo-random number sampling and obtained descriptive statistics ten price scenarios (Ω = 10) are generated for each hour (Figure 7).The probability of occurrence ( , ω) of the each scenario ω in the hour t is 1/Ω .
Observed GENCO owns HPS Vinodol.For clarity and computational effectiveness HPP Vinodol profit is maximized while rest of HPS is used for that objective.There is no significant lose in accuracy when having in mind that in GENCO's profit, HPP Vinodol profit participates with 95% or more.Probabilities are considered stationary, is 100%, is 3%, and are 40% and 35% respectively [35].The confidence level α is set to 0.85.Having in mind ancillary services requirements [41,42] and the HPP Vinodol ramp rate limits (Table 3) the problem is considered sufficiently accurate.The ramp rate limits of HPP Vinodol are as follows; MSR(up) is 255.15MW/min, MSR(down) is 340.2MW/min, ramping capabilities UR and DR between hours is 94.2 MW/h.The profit tolerances are minCVaR(1, ) = 0 , minCVaR(2, ) = 400 , minCVaR(3, ) = 800, and minCVaR(4, ) = 1380 for every t.The perfectly inelastic (vertical) energy marginal cost curves here considered as bid curves are obtained for the ordinate values 38 MWh and 57 MWh (Figure 8).For the price threshold $0 and $1380 the mean shadow prices are $1.9 and $2.4, respectively.This means reducing risk costs, and these costs can be determined.For the regulation DAM (Figure 9), the perfectly inelastic bid curves are obtained for 18 MW and 37 MW and at the same points all other bid curves are ending, reason for this is regulation headroom and footroom.For the price threshold $0 and $1'380 the mean shadow prices are $0.11 and $0.135, respectively.
For spinning reserve DAM (Figure 10), most of the bid curves are perfectly inelastic and located around 38 MW.Optimization showed that GENCO should seek to bid certain amount of capacity.When the price threshold increases from $0 to $1380 the mean spinning reserve prices are $0.065 and $0.075, respectively.
The water shadow prices have relatively low deviation compared to DAM prices and in this research use of the simplified approach with the constant water shadow prices over the planning horizon is justified.The hourly water shadow prices also have relatively low mean value which is a result of having more than enough water with relative relationships maintained so conclusions can be given (e.g., ratio between energy and ancillary services bid prices is similar to the ratio between energy and ancillary services DAM prices).The results were not obtained for lesser amount of water to avoid long computation time.Further bid curve analysis is omitted in this work since its slope, ordinate and aspics position is highly determined by input parameters.When the function is nonincreasing in some intervals two period moving average approach is used for determining the nondecreasing function as depicted on Figure 11.Also, it is assumed the bid curve is perfectly inelastic in last point as depicted on Figure 11 (red line).Thus, GENCO can create perfectly inelastic bid curve consisting of energy bided marginal cost of producing energy or ancillary (Figure 12).
As the profit tolerance rises from $0 to $1380 for the each hour (Table 3), the profit decreases from $98,193 to $96,959, daily α-VaR increases from $0 to $44,037 and daily α-CVaR increases from $0 to $33120.In other words, when GENCO decides the average hourly profit of 15% worst scenarios should be at least $1380 then GENCO can expect daily profit of $33120 in 15% worst scenarios, the minimal daily profit GENCO can expect with the confidence of 85% is $44,037, and the average profit GENCO can expect is $96,959.The obtained results confirm a risk aversion principle, generally GENCO should avoid hours with high price deviation as can be seen for the 13-th and 15-th hour in Figure 12.In this case study GENCO should seek certainty even at cost of reducing its profit.The results obtained using normal PDF for representing DAM prices confirm risk aversion principle, but for higher exactness lognormal probability distribution should be used to the pronounce benefits of this approach.

Conclusions
The approach presented in this work assesses the problem which the price-taker hydro generating company faces when participating in the simultaneous day-ahead energy and ancillary services markets and which results in the risk-constrained approach for pricing energy and ancillary services.
The results of this case study justify the use of a simplified approach where water shadow price was considered constant over planning horizon.It was also shown that perfectly inelastic bid curves are appropriate for a price-taker hydro generating company.It was shown that reducing risk of financial losses costs, and that these cost can be easily measured.
Since the research was done on relatively small number of scenarios, further work should be done to implement computational efficient algorithms which will tackle dimensionality issue and improve computational time.For higher exactness lognormal probability distribution should be used which will pronounce benefits of this risk-constrained approach.

Figure 1 .
Figure 1.(a) The three-dimensional unit performance curves; (b) The operating range of HPP when committed in energy and ancillary service markets.

Figure 2 .
Figure 2. The hypothetical probability density function (PDF) of a (a) loss; (b) profit.

Figure 7 .
Figure 7. Randomly generated daily price scenarios for the energy DAM.
Reservoir X

Table 3 .
Daily profit and daily risk measures α-CVaR and α-VaR for confidence level = 0,85 expressed in dependence of hourly α-CVaR(t).