Long-term Hydrothermal Bid-based Market Simulator

$I$

Set of agents indexed by $i$.

$I_{-i}$

Set of agents excluding agent $i$.

$I^{M}$

Set of price-maker agents, indexed by $i$.

$I^{T}$

Set of price-taker agents, indexed by $i$.

$J^{G}$

Set of thermal plants indexed by $j$.

$J^{H}$

Set of hydro plants indexed by $j$.

$J^{L}$

Set of lags of autoregressive model, indexed by $l$.

$J^{R}$

Set of renewable plants indexed by $j$.

$J^{U}(j)$

Set of hydro plants upstream plant $j$, indexed by $y$.

$J^{V}$

Set of vertices (points) representing a convex hull.

$S$

Set of sampled scenario indices, indexed by $s$.

$T$

Set of time indices, indexed by $t$.

$J_{i}$

Subset of a set $J$ with the indices of elements that belong to agent $i$.

Subscripts related to stage, $t$, and scenario, $s$, will be omitted for simplicity when they are not essential for the reader’s understanding of sampling and chronological relations.

$\mathcal{A}$

Vector of inflows of all hydros for all stages and sampled scenarios.

$C_{j}$

Operating cost of thermal $j$.

${E}^{Q}_{j}$

Energy quantity for an element $j$ of a convex hull.

${E}^{R}_{j}$

Revenue for an element $j$ of a convex hull.

$G_{j}$

Maximum generation of thermal $j$.

$P_{i}$

Price offer of agent $i$.

$\mathcal{P}$

Vector price offer of all agents for all stages and sampled scenarios.

$P^{F}$

Forward contract price.

$Q_{i}$

Quantity offer of agent $i$.

$\mathcal{Q}$

Vector of quantity offer of all agents for all stages and sampled scenarios.

$Q^{F}$

Forward contract quantity.

$\tilde{R}_{j}(\omega)$

Maximum generation of renewable $k$ at scenario $\omega$.

$\mathcal{R}$

Vector of maximum generation of all renewables for all stages and sampled scenarios.

$U_{j}$

Maximum flow through turbine of hydro $j$.

$V_{j}$

Maximum storage of hydro $j$.

$\tilde{\varepsilon}_{i}^{t}(\omega)$

Inflow noise coefficient of hydro $j$, stage $t$ and scenario $\omega$.

$\phi_{j,l}$

Inflow autoregressive coefficient of hydro $j$, lag $l$.

$\pi$

Spot price.

$\mathit{\Pi}$

Vector of spot prices for all stages and sampled scenarios.

$\mathcal{M}$

Vector of transition probabilities matrices, for all stages.

Indexing vectors in calligraphic ($\mathcal{P},\mathcal{Q}$) by $i$ stands for the vector where elements belong to agent $i$ and by $-i$ stands for the vector where elements belong to all agents except $i$.

$a_{j}^{t}$

Inflow at hydro $j$, stage $t$. $a^{[t]}$ stands for the vector of all inflows before stage $t$.

$\textbf{a}^{[t]}$

Vector of inflows at all hydros, for all stages before stage $t$.

$e$

Energy offer.

$g_{j}$

Generation of thermal $j$.

$q$

Bid quantity accepted in a dispatch.

$r_{j}$

Generation of renewable $j$.

$u_{j}$

Turbine flow at hydro $j$.

$v_{j}^{t}$

Storage at hydro $j$, at the beginning of stage $t$, and at the end of stage $t-1$.

$\textbf{v}^{t}$

Vector of storage at all hydros, at the beginning of stage $t$, and the end of stage $t-1$.

$z_{j}$

Spill flow at hydro $j$.

$\lambda_{j}$

Convex hull value for vertex $j$.

Indexing vectors in bold ($\textbf{a},\textbf{v}$) by $i$ stands for the sub-vector where elements belong to agent $i$.

$\mathbb{E}_{\omega}[\cdot]$

Expected value over the random variable $\omega$.

$\tilde{B}(\cdot,\cdot)$

Future cost as a function of states and uncertainty.

$\tilde{\Lambda}(\cdot,\cdot)$

Revenue as a function of energy, $e$ and scenario $\omega$.

$\rho(\cdot)$

Hydro generation as a function of turbine flow. It can also be a function of volume.

$|\cdot|$

Cardinality of a set.

Electricity markets are usually divided into two main groups: centralized markets (also known as cost-based) and liberalized markets (also known as bid-based). In the following, we present two figures to help contrast these two designs. Other market designs are possible [7], and some systems might mix the two (to some extent), like Brazil, which is cost-based but allows for flexible demand bids (small compared to the total load).

The two types of agents introduced in Section II-A2, price takers and price makers (see Figure 2), aim to optimize their bids and maximize profits, but using different approaches and methods.

We present two equilibrium models and associated solution methods that will be instrumental for our developments.

MSO is a framework to handle problems that are coupled in time and depend on random variables. We can formulate the linear programming version through its Bellman recursion, which depicts the problem being solved in each stage and how they are coupled through state variables $\textbf{x}^{t}$:

SDDP is a fundamental algorithm for effectively solving MSO and is applicable as long as they follow two main assumptions: 1) convexity of stage problems (which is trivial in the above linear case); 2) stagewise independence, that is, for all $t$, $\omega^{t}$ is independent of $\omega^{\tau}$ if $\tau<t$.

To apply the SDDP to the above problem, we need to reformulate it. The reformulation consists of replacing the future cost function (FCF) – the second term in the objective function (3) – by its piece-wise linear formulation, which is represented by the epigraph variable, $\alpha$, and cuts (5) as follows:

The above model considers a theoretical representation of the true model, as all the cuts representing the FCF in (5) are not known in advance. Therefore, the SDDP relaxes the set of all cuts, $\mathcal{K}$, to a subset $[K]$ that is iteratively updated. In this case, $\tilde{B}^{t}$ is also replaced by $\tilde{B}^{t}_{K}$ as it will only be a lower approximation of $\tilde{B}^{t}$ (at iteration $K$).

In a very high-level description, the SDDP algorithm starts sampling a scenario, $\omega^{t}$, for each $t\in T$, then solves the stage subproblems in chronological order, generating a feasible solution, this phase is called the forward step. After that, in the so-called backward step, problems are solved in the reverse order of time, generating Benders cuts to improve the representation of the FCF and, thus, propagating information from the future to the present. By averaging multiple feasible solutions for forward passes, we compute a statistical upper bound, while the solution of the first stage problem gives a deterministic lower bound. These steps are repeated until a specified stopping criterion is reached. This process is also known as policy optimization (or training). For algorithm details, the reader is directed to [9, 34]. After convergence is achieved, it is usual to proceed with a final forward pass in which we sample scenarios $s\in S$, and for each of them, we solve all stage subproblems in chronological order. This last procedure, known as simulation, results in a solution for each primal and dual variable. In other words, we finish with values for each stage in $T$ and sample scenario in $S$ for all optimization variables.

In this section, we focused our presentation on the expected value case of MSO for simplicity. However, it is relevant to mention that the methodological developments considered in this paper hold for other risk measures, such as the conditional value at risk, that can be consistently represented in the SDDP scheme (see [35, 36, 37]).

The first version of the single price-maker hydro owner agent was modeled and solved with SDDP in [24], but it did not consider uncertain bids from other agents. Hence, we present an extension of the technique that handles the uncertainty from other agents’ bids with a Markov model. Model (16)–(22), representing the optimization of an individual agent, $i$, is very similar to (9)–(15), the centralized operation model. (16) represents a Bellman recursion analogous to the one of the centralized dispatch. The main difference in the objective function is that there is a new term to describe the agent’s revenue in the given stage $t$ for a random event $\omega^{t}$ as a function of the electricity $e$ produced by the agent. (17) states that $e$ equals the total generation among all resources of a given agent $i$. Finally, (18)–(22) are almost the same as (11)–(15) but only accounting for generators owned by agent $i$.

Although renewable energy generation is still considered stagewise independent and the inflow dependency is still handled by the auto-regressive process of Section II-D1, the revenue function uncertainty is not restricted to the RHS. Hence, as highlighted in Section II-D2, the solution of the above model will require the Markovian SDDP, in fact, a mixed auto-regressive (for inflows) and a Markovian (for the revenue) SDDP. The uncertainty representation will be detailed in the next paragraphs.

A second fundamental challenge in this procedure is requiring all hydro plants in the same cascade to belong to a single agent. This hypothesis is also assumed in previous works like [25]. Handling different owners in the same river system can be done by considering a water wholesale market besides the electricity market as done in [45]. Alternatively, other market designs can be used. Two possible alternative designs are: 1) hydro slicing, in which agents own fractions of the cascade and, consequently, can optimize their strategies as if they did not share the cascade, and 2) virtual hydro reservoirs, in which all hydro plants of a system are aggregated in a single reservoir which is then split proportionally among all agents. The interested reader is referred to [28] for more details.

In the case of price-taker agents, their operation can be optimized by considering the revenue function, as in [29]:

where $\pi(\omega)$ is a scenario of spot prices that might be obtained from any model. In particular, we could consider spot price scenarios $\mathit{\Pi}$ produced by the simulation from Section III.

On the other hand, price-maker agents have the ability to affect spot prices with their energy offers, as in the Stackelberg Equilibrium of Section II-C1. In other words, $\pi$ is no longer an exogenous input data, it is a function of the energy offer of agent $i$, i.e., $e$. In this case, price and quantity bids from other agents are considered random variables depending on the uncertainty $\omega$, and are given $\{(P_{j}(\omega),Q_{j}(\omega)),j\in I_{-i}\}$. With those bids, we follow the same rationale of [24, 25]. Therefore, we write the following model that expresses a simple bid-based dispatch problem:

This is analogous to a market clearing problem, in which a system operator selects the optimal amounts of energy to attain the demand in (25), under the limits (26)–(27) given by the offers of the current agent and the other players. The main result is the spot prices (the dual variable of (25)). Now we define:

However, this strategic agent revenue function is a non-convex piece-wise linear discontinuous function as detailed in [24, 25]. Thus, to satisfy the SDDP convexity requirements, we follow the method proposed in [24] and represent the convex hull of $\tilde{\Lambda}(e,\omega)$ with respect to $e$ for a fixed $\omega$:

where each pair $(E^{Q}_{j}(\omega),E^{R}_{j}(\omega))$ represents a vertex (in the set of vertices $J^{V}$) of the convex hull of the hypo-graph of ${\tilde{\Lambda}}(e,\omega)$. We keep the vertices indexed by $\omega$ to make it clear that they depend on the random variables of the problem. Theoretically, it would be possible to exactly represent the non-convex revenue function with methods like SDDiP [46]. However, this would come with the cost of a considerably larger computational burden. Because we are accounting for many computationally intensive features in this work, we adopted the convex hull representation.

As one of the key mechanisms to mitigate market power, we must also be able to represent forward contracts [20]. This was not described in the previous MSO models solved by SDDP. However, it is simple to modify the revenue function to consider two additional terms as follows:

the first term is the fixed revenue of the forward contract, the second represents the energy that must be delivered due to the contract, and the third is the previously represented earnings from the spot market. The constants $P^{F}$ and $Q^{F}$ are input data. Consequently, the first term is constant and does not affect the optimization of a risk-neutral (strategic or not) agent. On the other hand, agents are discouraged from generating below contracted quantities (where $e-Q^{F}<0$) as they will have to buy energy to match their obligations, thereby losing money. This has an effect that opposes the willingness of price-maker agents to withhold energy (reducing $e$) to increase spot prices.

Figure 3 shows revenue curves for a simple case where $P^{F}=0$$/MWh, $D=40$MWh, and we consider $3$ price-quantity offers from the other agents: (3$/MWh,10MWh), (2$/MWh,15MWh), (1$/MWh,20MWh). Note that the total quantity from the sum of the other agents’ offers is $45$MWh, which is higher than the demand, and, hence, no deficit occurs even if $e=0$MWh. This function is also non-convex, but we can represent its convex hull in the optimization problem in the exact same way as the previously described case with no contracts. For instance, we consider $Q^{F}=0$MWh for which we only have the last term of (31). In this case, the revenue is a linear function of $e$, with a slope equal to the spot price, $\pi$. The most expensive accepted bid sets the spot price, then a bid $e<5$MWh for the strategic hydro leads to $\pi=3$$/MWh. When the strategic hydro offers more than $5$MWh, it displaces the bid previously defining the spot price, and the price drops to the variable cost immediately below, i.e., $2$$/MWh. So, from this point on, the curve has a smaller slope, and the first break appears. The other curves and breaks follow the same rationale but for different values of $Q^{F}$ and $e$. For more details, we refer the reader to [24], which only considers the last term in (31). We emphasize that for different values of $Q^{F}$, the first term in (31) does not affect this example; as for simplicity, we considered $P^{F}=0$$/MWh. However, the second term subtracts a piece-wise constant function from the revenue function as $Q^{F}$ is constant for each curve and $\pi(e,\omega)$ is a piece-wise constant function. Because $P^{F}=0$$/MWh, the blue curve, which considers no contract, is the highest. In more realistic cases, these curves can cross each other.

We remark that in both [24] and [25], the revenue function is not stochastic since it only considers thermal bids in the form of thermal installed capacities and operating costs and does not require the Markovian SDDP from Section II-D2.

In contrast, we allow stochastic and time-dependent uncertainty in the two cases: 1) spot prices in the price-taker case and 2) energy bids from other agents in the price-maker case. In both cases, we assume we have scenarios of either spot prices, $\mathit{\Pi}$, or price and quantity bids, $(\mathcal{P},\mathcal{Q})$. Such sets of scenarios are jointly distributed with the inflow, $\mathcal{A}$, and renewable energy, $\mathcal{R}$ scenarios if they were generated by SDDP simulations.

Given all these jointly distributed scenarios, we apply the same Expectation Maximization (EM) algorithm as [47] to obtain the transition probabilities of the Markov model for all stages, $\mathcal{M}$, and the classification of each sample, $s$, of each stage, $t$ in a Markov state.

To start the diagonalization-based algorithm, we need an initial set of actions for all agents. Such actions will be obtained from a simulation of the centralized audited-cost-based market design of the power system through the typical cost-minimization model, (9)–(15), from Section III solved with SDDP from Sections II-D and II-D1. Besides the physical description of the system, the main input data for this model are the inflows and renewable generation scenarios, $(\mathcal{A},\mathcal{R})$. The main results from the simulation procedure are, at a given stage $t$ and scenario $s$: 1) the generation decisions of all units of a generation agent $i$, which lead quantity part of the bid, $Q_{i,t,s}$, and 2) the spot prices of the system ($\pi_{s,t}$) that will be to the price parcel of the bids. We denote the set of price and quantity bids of an agent $i$ as $(\mathcal{P}_{i},\mathcal{Q}_{i})=\{P_{i,t,s},Q_{i,t,s}\}_{t\in T,s\in S}$ and the spot prices as $\mathit{\Pi}$. We will refer to this procedure as $CentralizedOperation(\mathcal{A},\mathcal{R})$.

In this step, the actions of each agent will be updated, while the actions of the other ones will be fixed until the convergence criterion is satisfied. The main goal is to simulate the power market in an agent-based fashion. If convergence is strictly attained, we might have reached a Nash equilibrium. However, such an equilibrium might not even exist in this setting. The convergence criterion considered in this work is assuring that the maximum absolute variation of the price and quantity bids of all agents between two subsequent iterations vary less than 1% of their values in the centralized operation.

This system was constructed from real data from the Brazilian system expansion scenario for $2025$ (based on the Brazilian market operator data [48]) and contains $46$ thermal plants and $21$ hydro plants with a simplified system topology. Note that this is already more than the current $13$ aggregated reservoirs considered in the official model. The overall hydro installed capacity represents $70\%$ of the system’s installed capacity. To illustrate the algorithm’s functioning, first, we consider three main hydro generation companies, one considered as a price-taker agent, and the other two as price makers. The $21$ hydro plants are split into the $3$ aforementioned hydro generation companies, $7$ plants for each one. The three different agents do not share cascades as mentioned in Section IV. All thermal plants are considered individual price-taker agents. We considered a single load block for the sake of simplicity, that is, a constant demand for the entire month. The following simulations were carried out in a $5$-year monthly horizon, that is, $60$ stages. We considered $1000$ sampled scenarios. The large-scale simulations were performed in a PSRCloud cluster of $8$ servers, each one with $64$ cores and $128$GB of RAM. Each simulation took approximately $3$ hours. All algorithms were coded in the Julia language [49], and the optimization models were coded with JuMP [50] and solved with the Xpress solver[51].

In this section, we present five controlled simulation experiments to understand the effect of different market concentrations on market power abuse. We adjust the hydro plants so that we have different databases defined by tuples $({share}_{1}\%,{share}_{2}\%)$. In each case, ${share}_{1}$ and ${share}_{2}$ represent the shares (in percentage) of the hydro system belonging to each of the two price-maker companies (agents). The number of plants of each hydro agent is fixed as mentioned in Section VI-A. We multiplied and installed capacities and reservoir sizes by scalars to adjust the hydro system size of each agent in each database. The remainder of the hydro system is assumed to be of price-taker agents.

Figure 6 presents a typical convergence profile of the proposed iterative method. These profiles were taken from the case study with price-maker shares of $(25\%,50\%)$ and no contracts. Figure 6 shows both average absolute and relative differences between two consecutive iterations. First, absolute (or relative) differences are computed for each stage and scenario, and then, they are averaged.

Figure 5 contains results from this first study and from the next study for better understanding. For now, we shall focus on the centralized cases (in green) and the multiple cases with no contracts (all in red). In Figure 5a, we show the average spot prices as a function of the market distributions in a bar plot. We add the first bar with the spot price of the centralized dispatch. As expected, the spot prices rapidly grow as the concentration increases. Figure 5b depicts the captured price, i.e., the average value of energy sold, which is the total revenue of each agent divided by its total generated energy. We present results in such normalized forms to compare the results of different-sized agents. Figures 5c and 5d show additional results with average reservoir levels of the price-maker agents throughout the study period and the average spillage. Interestingly, market power abuse is mostly characterized by excess spillage, whereas, on average, the total reservoir levels do not vary much compared to the benchmark (centralized dispatch with audited costs).

In this section, we repeat the above analyses but with contracts. To come up with monthly contract quantities, we take average generation values from the centralized dispatch and prices from the average spot prices in the centralized dispatch. This will stimulate the agents to produce energy and reduce the spot prices since being short in the contract together with high spot prices will lead the agent to have expenses due to the second term in (31).

Now we look again to Figures 5a, 5b, 5c and 5d described in the previous section but we also will focus on the cases of agents that are $75\%$ (purple) and $100\%$ (blue) contracted, respectively. In the case studies, we can clearly see that the contracts almost completely eradicated the market power, and the resulting spot prices are very close to the ones obtained in the centralized dispatch. At the same time, the captured revenues and spillage levels were shrunk to values close to the ones observed on the centralized dispatch with audited costs. By analyzing Figure 5d, we can see that spillage is a relevant marker for market power abuse.

This dataset was created from the original data of the Brazilian system obtained from the Brazilian market operator webpage [48]. The dataset contains $137$ thermal plants representing $18\%$ of the installed capacity, $364$ renewable energy plants (wind and solar) representing $19\%$ of the installed capacity and $32$ hydropower plants representing the remaining $63\%$ of the installed capacity. Here, we improve the load duration curve accuracy and consider $5$ load blocks [52, 53] instead of just one (from the previous section) to better represent peak demand hours, $5$ is more than the $3$ currently used in the Brazilian official planning tools [48]. Like in previous studies, we consider $5$ years ($60$ monthly stages), $1000$ scenarios, and the same PSRCloud cluster configuration.

In this study, we consider an approximation of the real market concentration in the system. We represent three companies as price-maker agents, with respectively $32\%$, $9\%$, and $7\%$ of the resources in terms of firm energy representing the three largest generation owners. The firm energy metric is used as hydro and renewable capacity is not fully available all the time, so this is a proxy of their average available generation, see [54]. Meanwhile, the remaining $55\%$ of the generators are distributed among one price-taker hydro company, and each thermoelectric generator is considered another small price-taker agent. For completeness, we computed the Herfindahl-Hirschman Index (HHI) of the system, which led to the value of $0.125$, which is below $0.15$, a number that typically represents a small concentration. However, as will be empirically demonstrated in the next study, there is clear room for market power in the system in question. We believe that it is yet another indication of the limitations of HHI that were already highlighted by [55, 20].

Each simulation from this study took approximately $11$ hours using the same code and computational infrastructure of the previous study.

First, we simulated the centralized operation and the bid-based market without contracts. Then, we considered agents with multiple contracting levels in the bid-based market. In particular, we analysed the cases of: $25\%$, $50\%$, $75\%$ and $100\%$ of contracting level. The key results from each of the $6$ simulations are depicted in Figures 7 and 8.

Figure 7 shows the spot price under the above $6$ conditions. We can see the progression from the case with the highest spot prices, where no contracts are considered, until the fully contracted case, where the average spot price is much closer to the one obtained from the centralized operation with audited costs. Given the number of resources allocated to the price-maker agents, $45\%$, we note that a $75\%$ contracting level effectively reduces the gap to the centralized operation. This closely follows what we have seen in the previous section with the southeastern system. In such case, $75\%$ was also effective in containing market power in the case where the first price-maker agent had $15\%$ and the second price maker had $30\%$ of the resources, resulting in a $45\%$ share of the system in the hands of price-maker agents. Finally, we note that even low contracting levels, such as $25\%$, significantly reduce the final average spot prices. So, we highlight the empirical evidence in our case study that the typical contracting levels observed in practice, ranging from $75\%$ to $100\%$, can be seen as a relevant instrument to prevent market power abuse in the Brazilian power system.

Figure 8 presents the overall revenue, including operation and contract costs and spot and contract revenue throughout the $5$-year horizon for each of the three price-maker agents. Overall, the figure reproduces information already verified in Figure 7 with excessive revenues for all agents for the bid-based cases with no contracts. Notably, the $100\%$ contracting case leads to revenues that are very similar to the ones from the centralized dispatch. Interestingly, $75\%$ seems to be a relevant threshold value as it induces behaviors resulting in revenues and spot prices relatively close to the ones obtained in the centralized case. This observation is also relevant because when generators are risk-averse, it is expected that their optimal forward involvement should be slightly below the 100% case (see [56] and references therein).

Figure 9 presents the seasonal average spot price profile. Spot price changes are driven by both withheld and spilled hydro water. Noticeable, the spot-price peaks caused by the market power abuse are concentrated at the end of the dry season, when reservoirs are at their lowest levels (September-December), and in the month with the highest inflow (March), in which strategic agents can easily spill.