Medium-Term Power Planning in Electricity Markets with Pool and Bilateral Contracts ✩

Many of the existing electricity markets are of the mixed type, which has pool auction and bilateral contracts between producers and distributors. In this case, the problem faced by a Generation Company (GenCo) is that of maximizing the revenues from participating in the market through the pool auction while honoring the bilateral contracts agreed, for which the revenue is ﬁxed. The extension to mixed markets of a medium-term model, successfully employed for auction-only markets, is presented. It results in a non-convex expected revenue function to be maximized subject to constraints, for which the currently available direct global-optimization solvers prove not to be e ﬃ cient enough. A heuristic procedure based on a sequence of solutions by a nonlinear solver is presented, and numerical results obtained with several realistic cases show satisfactory results. The test cases presented have dispatchable and non-dispatchable renewables and consider medium-term pumping together with conventional units by all GenCos participating in the mixed market. The advantages for GenCos of employing medium-term results as those produced by the model presented, include, among others, the evaluation of the expected proﬁtability of their bilateral contracts.


Introduction
Liberalized mixed electricity markets typically offer two trading systems: the pool and the bilateral contracts (BC). Part of the generation is supplied to customers under bilateral agreements and part is traded in the pool market auctions. The Nordpool in Scandinavia and the MIBEL in Portugal and Spain are examples of mixed trading system.
The responsibility for matching the BC load and for satisfying other medium-term constraints, as the management of hydro resources, is on the GenCos. In mixed systems the GenCos also decide which units will be devoted to supplying their BC load, and during which hours. For every hour the System Operator (SO) in a mixed system market collects the information from the pool auction and from the bilateral exchanges, determines the system load, and checks that the transmission network can safely convey the generations and consumptions agreed. From data released by the SO one can deduce the predicted system load and the total BC load.
The aim of this paper is to present a model for a GenCo participating in a mixed market that wishes to optimize its power planning for a medium term horizon (typically of one to two years, subdivided into several periods). GenCos must also optimally split over the successive periods the available time of their units into time devoted to match BC load, and time to participate in the pool market: while BCs represent a longer-term stability, they also represent loss of opportunities. The procedure for reaching the bilateral agreements (time, amount of supply, and price) is not considered here: in the planning model presented it will be ensured that these agreements are executed, without optimizing the revenues from the BCs, since their supply produces a previously known fixed revenue.
Other authors have approached a similar issue as follows. In [1] and [2] risk-benefit functions are defined and an iterative negotiation is proposed before the agreement on a BC is sealed. A practical process is put forward in which the bargainers take both benefits and risks into account. The authors claim that their process will lead to agreement on a mutually beneficial and risk-tolerable forward bilateral contract. In [3] the authors present the model for optimizing the amount of energy to be sold and bought in a futures market in order to hedge against pool-price risk for a pricetaker company. In [4] a portfolio optimization problem is proposed for allocating the assets of a supplier and in [5] a multi-agent negotiation scheme is proposed to analyze the behavior of a bilateral market.
Another family of models used by GenCos is the resource allocation planning models for medium-term. These models support decisions on fuel procurement, water resource management or strategic planning among other processes. An example of that is in [6], where a mixed model is proposed that optimizes both the water resources allocation and the forward contracts for a small producer that owns only hydro generation. In this formulation, favorable decisions are rewarded and unfavorable ones are punished. Uncertainties, such as water inflows or market price, which is exogenous to the company's generation, are considered. They are represented through scenarios. Forward contracts are used to hedge against price volatility, due to uncertain hydro inflows. A medium-term energy procurement model for a consumer with self-production and access to BCs can be found in [7], where the authors propose a mixed-integer linear programming model. In [8] it is for a large consumer with access to several types of BCs with known tariffs and bounds; the expected value of the procurement cost is minimized while limiting its risk by including risk aversion through the CVaR methodology.
Attention has also been given to pumping together with BCs, as both divert generation resources from satisfying the auctioned load. There are works on the short-term operation of pumped storage units in electricity markets selling its generation and buying energy for pumping based on a previously determined market-price hourly forecast, for a week in [9], and for a day and establishing a BC price in [10].
The work here presented also shows the specific power planning model and presents optimization techniques to solve the problem posed in a mixed system market, where it is necessary to match several load-duration curves in the same period and the objective function becomes non-convex. One of the key issues is how to model the electricity market price. In the context of sufficient information on the GenCos forming the pool and on the market auction, it is possible to estimate a linear relation of the market price with respect to the load duration, leading to quadratic market revenue functions for each period of the expected unit generations [11]. The model and the solution procedure presented in this paper are the extension to a mixed market system of the pure pool model and solution procedure presented in [12], where the price function employed is linear and endogenous, as it depends on the hydro generation level [13]. In this work it is also dependent on the non-dispatchable renewable generation level represented on nodes of a scenario tree as proposed in [14]. The stochastic endogenous model with equilibrium behavior was introduced in [15].
The paper is organized as follows: Section 2 recalls the basic model for a pure pool market with an endogenous price function for cartel and equilibrium behavior; Section 3 analyzes the data referring to BCs available in the MIBEL and presents the mixed-market model; Section 4 is devoted to the optimization techniques employed for solving the quadratic non-convex mixed-market model; Section 5 follows presenting the computational results for several test cases, and Section 6 summarizes the conclusions.

The medium-term pure-pool model
The main characteristics of a model used to represent a GenCo operating in a pure-pool market system, where all electricity produced and consumed is traded through the market operator hourly auctions, will be recalled first.

Load and generation in medium term
Although consumers' load in a future medium-term period i (used as supra-index of parameters and variables), of length T i hours, changes randomly over time, it has fairly-well predictable features as its peak p i and base load p i (MW), and its total energy to be supplied e i (MWh). These parameters, together with a well-predictable shape, characterize the probabilistic load distribution function and the associated load-survival function S i ∅ (z), which gives the probability of a load being greater than a non-negative z (MW), i.e., S i ∅ (z) = Prob(demand ≥ z). The load-duration curve (LDC) is an alternative representation of the load-survival function, as the axis are changed and the probabilities are scaled by the period length T i to give the load in decreasing load order for each load duration t. Even though the LDC does not provide a detailed representation of the load dynamics, it represents all the load variability.
Regarding generation units, three parameters are considered in medium term to characterize them: the generation capacity c j (MW), a linear generation price f j ( /MWh) and a failure probability q j for each unit j.

Probabilistic matching of the load
The ensemble of generation units that, over the hours of a certain period, have their generation bids accepted match the LDC of this period. This is a constraint that the expected generations that take place in a market must satisfy, in the same way as it was for regulated electricity systems.
Through the LDC-equivalent load-survival function, the unit outages can be easily taken into account in mediumterm load matching. Considering the load-survival function S i θ (z) of the still unmatched load after using for generation all units in subset θ, the convolution first proposed in [16]: expresses the change to the load-survival function caused by also using unit j after having previously used all units in subset θ ∈ Ω, being Ω the set of all generation units. The expected unsupplied energy s i (θ) after having loaded the units in θ is computed as: where p i is the peak load and T i the duration of period i. Bloom and Gallant [17] formulated the exact probabilistic matching of a LDC as the set of linear inequality constraints: where x i j is the expected energy generated by unit j over period i, e i is the total LDC energy and s i (θ) is calculated as in (2). This set of constraints will be referred to as load-matching constraints (LMCs) and their number is exponential as there are n LMC =2 n u −1 subsets θ of Ω, being n u = |Ω|. It should be noted that not only the number of these constraints matters but also the cumbersome numerical calculation through (2) of the right-hand sides in (3).
Besides the LMCs (3) for each period, there may be other constraints relating the expected generations x i j of a single or of different periods. They will be called non-LMCs. The LMCs and the non-LMCs define a feasible region for the expected generations where a utility function in terms of the expected generations x i j , j ∈ Ω, ∀ i could be maximized. Examples of non-LMCs are single period and multi-period emission caps, take-or-pay contracts of certain fuels, availability limits of hydro generation, etc.
Because in a liberalized market the amount of energy produced by each company is determined daily in the spot market auctions, the only LDC that can be predicted is that of the whole system. Although the model can be used for a specific GenCo, all units participating in the pool must be represented. Therefore Ω must contain all units of all GenCos participating in the market. In practice, in order to have a problem of moderate size, the GenCo that is planning its medium-term generation will consider its units in detail, and those of the rest of GenCos will be merged into several big units of different technologies that will represent the rest of units in the market.
The LMC formulation (3) has been employed with indirect solution procedures not requiring the previous generation of all LMCs [18,19], to the pure pool system. These indirect procedures are not sufficiently efficient for problems of moderate to large size.

The use of the load matching heuristic for the LMCs
Off-the-shelf linearly constrained optimization solvers would require the generation of all LMCs for all periods, (or for all nodes of a scenario tree in case of use of Stochastic Programming), which is not practical. The load-matching heuristic (LMH) for dealing efficiently with the LMCs proposed in [20] is a finite multistage procedure to generate a reduced subset of LMCs that will most likely contain the active LMCs (satisfied as equalities) at the optimizer. This heuristic was employed in [13] and [15]. The heuristic in [20] was modified in [12] to improve its performance for medium-term problems with more complicated utility functions. This modified load-matching heuristic (mLMH), was employed in [14] and is also employed in this work. Alternative procedures of LDC matching are compared in [12] and their associate loss of precision are quantified.
The strategy of this heuristic is to use the same number of sets as LMH uses but with the addition of other sets. All stages in the mLMH are associated to creating a batch of LMCs determined from the current values of unit generations, followed by an optimization solution producing fresh unit generation values.
The all unit but one and the all unit but two sets will be added at first stage of the heuristic (the initialization part) because it will restrict the number of units that are at or close to its maximum value in the cartel solution. The stages are: 1) Initialization: Initialize the list of LMCs with the n u all-but-one-unit LMCs, the n u × (n u − 1) all-but-two-unit LMCs, the all-unit LMC, and the n u single-unit LMCs (upper bounds) c j 0 S ∅ (z)dz and solve the problem.
2) Self-ordering: Based on the previous solution, add LMCs with all combinations of units in the subset χ whose generation in the former solution is at (or close to) its upper bound, expressed as ρ j = x j /x j ≈ 1 . Resolve the problem.
3) Low-load unit ordering: Form the subsets ξ 9 and ξ 7 with the nine and the seven less loaded units, if any, that have ρ j < 1/2 , and the complementary sets σ 9 := Ω \ ξ 9 and σ 7 := Ω \ ξ 7 . Add the LMCs made of the units in σ 9 and some combinations of units in ξ 9 and the LMCs made of the units in σ 7 and some combinations of units in ξ 7 as detailed and justified in [12], and resolve.

4)
Step-by-step order: One by one addition of LMCs by repeatedly including in set χ the most loaded unit (that is, the unit with largest ρ j ) not in χ, adding the new LMC made with the units in χ and resolving, until the size of χ is the total number of units minus two.

Alternatives to using a probabilistic load matching heuristic
The iterative procedure stems from the use of the mLMH, which, as is also the case of the parent LMH [20], converts the single optimization solution of the medium-term generation problem with all LMCs into a succession of optimizations with an increasing number of LMCs considered.
There are three types of alternative methods to using a probabilistic load matching heuristic: • those that employ a multi-block approximation to the LDC (see e.g., [21]), which are presented and computationally compared to the mLMH in [12], • that employed in [22] for matching an LDC with patches of different shapes (rectangular, trapezoidal and triangular) corresponding to different generation technologies placed inside the LDC approximately simulating the successive loading of units in ascending bid price in a day-ahead market auction, and • those that match a given set of hourly loads in a year but do not match LDCs in a probabilistic sense.
Instead of using the available probabilistic procedures of LDC matching by generation units having a failure probability, an approximation to the LDC consisting of a few rectangular blocks with a given power height and a given width of a certain number of hours could be employed. A single optimization is necessary, and the equations to match these blocks with generation units, even when taking into account the unavailability due to random outages, are simple and far less than the exponential number of inequalities necessary to formulate the probabilistic matching. However, the feasible domain represented by the block matching constraints and that represented by the probabilistic matching constraints is in practice sufficiently different to produce quite different results for the less loaded units (see [12]). This is the reason for having preferred the probabilistic matching procedure in this work.
As for the last type of alternative matching, although hourly load profiles could properly represent a given load dynamic, they would be impractical as an inordinate number of scenarios would be necessary only to represent load and generation outage randomness. Moreover, the former scenarios should be combined with those of random generations of renewable generations.
On the other hand, the use of the mLMH yields sufficiently good probabilistic LDC matching and allow for the use of scenario trees of moderate size to represent the randomness of renewable generations only.

The endogenous cartel and the equilibrium models
A realistic model must incorporate the market behavior and how the market players may interact. One approach is to assume that the market price can be predicted with sufficient accuracy (given that the market is an oligopoly with special characteristics) and therefore the decision makers' goal could be to maximize the sum of their profits [11,12,13,14,15].
Another approach is to deem that decision makers react to the decisions taken by the other participants in the pool, as in game theory. In this paper, both behaviors will be considered.
The Nikaido-Isoda relaxation algorithm (NIRA) for calculating the Nash-Cournot equilibrium point of a noncooperative game is described in [23,24]. It is based on a succession of optimizations with changes in parameters of the objective function according to the results of the former optimization. For applying the NIRA procedure to obtain the equilibrium in a liberalized electricity market there should be a formulation of the influence of each GenCo generations on the profits of the rest of GenCos in the market. The first application of the NIRA to a stochastic medium-term generation planning model taking into account the endogenous influence of hydro generation and of other generation types was in [15], where the successive optimizations take place within the stages of the LMH [20]. In [12] there is the application of the NIRA in a medium-term model taking into account the influence of hydro, with optimizations within the stages of the mLMH, and the extension to also using non-dispatchable renewables in a stochastic medium-term model taking into account the endogenous influence of hydro generation and of nondispatchable renewables on market price can be found in [14]. The harm to convergence of the NIRA due to the addition of new LMCs at each mLMH stage is reduced due to the fact that the successive cuts meant by the extra LMCs added are shallower as the final solution is approached.
In [25] the NIRA procedure is employed for finding the short-term (one-day scope) hydro-thermal equilibrium solution through predetermined linear price-demand functions taking into account network constraints, with and without bilateral contracts.
There are also procedures for the direct calculation of an equilibrium solution to the medium term generation planning through a single optimization [21] (extended to a stochastic formulation in [26]). However, this procedure is based on a single load-balance equation per period, whose Lagrange multiplier is the market price. This approach excludes the consideration of the probabilistic load-matching employed in this manuscript, as the probabilistic matching takes into account an exponential number of linear inequality LMCs, of which there is only a small, a priori unknown, active subset. Moreover, [21] assumes the prior knowledge, for each period, of the inelastic demand, of the demand elasticity coefficient with respect to market price, and of the elasticity coefficient of generations of each GenCo in the pool. These assumptions are hardly realistic as compared to the market price assumptions for each period employed here: the market price at peak load, the linear endogenous effect on it due to the renewable generation levels, and the market price decrease slope with load duration, all of them deductible from records of market price, of load levels, and of renewable generations.
In [27] a GenCo short-term equilibrium bidding optimization procedure is described that is also based on a single load-balance equation per period and per transmission network node.

The extension to a stochastic model
Expected generation of some of the units (and market prices) are influenced by exogenous factors such as water inflows in reservoirs or the wind speed in wind farms. Stochastic programming [28] gives stochastic parameters, such as hydro inflows or wind power (WP) levels, different values in different nodes ν ∈ N of a scenario tree leading to the joint optimization on a number of scenario paths λ ∈ L that go from a root node in the first period to as many leafs in the last period as there are paths.
The treatment of stochasticity of dispatchable renewable resources, as hydro generation, and that of non-dispatchable resources as WP and solar photovoltaic generation (SPV) are described in [14]. The market-price level in each node is influenced by the WP level in the node, as explained in [14].
Scenarios were developed for hydro inflows in reservoirs and final stored hydro energy, and for WP. The procedure followed was the quasi-Monte Carlo technique, first generating a fan tree from a two-dimensional lattice [29] (one dimension for each stochastic variable: hydro and WP) and then reducing it to a tree of the desired size [30].
In short-term planning for the day-ahead auction there are two bidding processes: that of generation bids to match the demand bids, and that for spinning reserve to compensate for mismatches in predictions of consumption and of non-dispatchable generation (wind power and solar photo-voltaic). It is argued that the errors in non-dispatchable generation predictions should be compensated by a non-zero cost for non-dispatchable bids [31], where a procedure is proposed to work out a suitable price for the non-dispatchable generation bids. In the medium-term procedure presented here the stochasticity of non-dispatchable generations is taken into account from the outset by the base unit with low failure probability and the complementary crest unit with high failure probability [14], and by using the probabilistic matching of the LDC, and by the stochastic programming scenario tree employed for hydro, WP and SPV. Should a price be imposed on WP and SPV for spinning reserve compensation, it could be placed on crest-unit generations corresponding to WP and/or SPV scenarios with levels far from the central expected ones.

Pumped-storage hydro stations in medium-term planning
The procedure for modeling a pumped-storage unit in a single medium-term period described in [17] was extended to multi period stochastic problems in [14] and is employed here to account for pumped-storage units.

The medium-term mixed model with bilateral contracts
The analysis of the amount of energy traded through BCs for several months compared to the total energy produced in the MIBEL indicates that it is consistent to consider the modeling of the energy traded through BCs as an LDC (Fig. 1). For future periods, BC duration curves (BCDC) could be forecasted in the same way as the system LDCs are estimated. Two types of BCs are contemplated by the regulations of the MIBEL: base contracts, which span the full length of the period, and peak contracts for several lengths of time. Piling up the peak contracts within a given period, and ordering the hours by decreasing system load, a generally decreasing BCDC will be obtained.

Limits of available information on bilateral contracts
The information on the past bilateral hourly load is available to the participants in the market, and from it, system BCDCs, as that in Fig. 1 right, for future periods can be deducted. Moreover, a specific GenCo (SGC) will know from its own records which are its forecasted BCDCs and, by subtracting them from the forecasted BCDC of the system they will have the BCDCs of the rest of participants (RoP). There are thus two types of BCDCs: that of the SGC and that of the RoP, each of which having its own generation units to satisfy its specific BC load. For notation purposes in the problem formulation, the two-element set B := SGC, RoP will be used, and a tilde on parameter, variable and set symbols will refer to BCs. The set Ω β ⊆ Ω, ∀β ∈ B will contain the units of either the SGC, or RoP, which match its own BCDC.

Time-share hypothesis
In order to ensure that each BCDC is matched by the units of each BC group, the corresponding BCDCs are added to the model. The LMC with all units of a BC group is put as an equality to ensure that BCs are honoured as planned. Then, to compute the revenue of the energy sold at market price (recall that the energy supplied through BCs is paid at a known pre-determined price and it is not subject to optimization), a time-share hypothesis is made to address the problem of a certain unit having the possibility of matching two different LDCs over a given period. This unit may devote all or a part of its generation to honor BCs as long as there is BC load still to be matched. The remaining available generation of that unit may participate in the market and will be rewarded using the market price function.
Suppose that the shaded part of Fig. 2 (left) corresponds to BCs. Let us assume that the solution is given in Fig.  2 (right), where each darker-shaded slice, or part of it, represents the expected generation of a unit satisfying BCs (assuming null outage probability). It is easy to see that the union of the darker areas conforms to the BCDC. As in the model of the pure pool trading system, in order to compute the expected profit, two assumptions are taken: (i) a unit generates at its maximum capacity and (ii) the shape of the contribution is approximated to a rectangle.

Maximization of stochastic market profits with endogenous price function and bilateral contracts
In order to calculate the profits in the mixed system, the endogenous market price function p ν (t, x h ) is needed, with hydro and WP influence for each node ν of the scenario tree, whose period is denoted by i(ν), where node ν and period i(ν) are indicated as superscripts: where b ν 0 is the node WP-influenced market price at peak period load, decreased by the endogenous hydro generation term (d/T i(ν) ) h∈H x ν h with d being an estimated negative factor, h ∈ H are the hydro pseudo units that represent the hydro basins of set H, and l i(ν) is the estimated negative slope of the market price with respect to the LDC load duration t. A linear price function like (4) with a WP and solar photovoltaic generation-influenced market price at peak load b ν 0 has been employed in [14]. The expected generation x ν j of unit j in node ν may have two parts: x ν j for honoring BCs, and the rest x ν j − x ν j , remunerated at market price. The market revenue is calculated multiplying the capacity of the unit by the integration of the price function during the unit expected generation duration, which starts at x ν j /c j and ends at the total expected duration x ν j /c j , as units matching BC peak load participate in the market in hours that come after the duration x ν j /c j of matched BC load. Subtracting the generation cost from the revenue, the generation unit profit r ν j (x ν j , x ν j ) is obtained (where the revenue from the BCs, which is fixed, is not included).
The revenue function (5) is a non-convex quadratic function of the expected generations x ν j and x ν j .
In each short-term hourly auction, through which supply bids match demand bids, the social welfare of producers and consumers can be calculated from the supply and demand curves. A proposal of its use for clearing the market in an electricity pool is detailed in [27]. In any of the medium-term periods there is neither a single supply nor a single demand bid curve, and the load is represented through a LDC, thus social welfare cannot be calculated; instead, a function of variation of the market price level with the load duration can be predicted, and also its level change with respect to renewable generations during the period. GenCos profits calculated as difference of market revenue and generation costs can be maximized subject to load matching and other operational constraints.
Moreover, the purpose of the procedure presented is to give a planning tool for a GenCo participating in the market. That is why the objective function contemplated is maximizing the expected net profit of the GenCo over a yearly horizon, expressed as the difference between the remuneration from the market (calculated from the linearized endogenous market clearing price curve of each period) minus the variable generation costs.
Market prices over time are influenced by random non-dispatchable (WP and SPV) generation levels and are also influenced by the generation bidding behavior of GenCos, especially regarding hydro generation in the MIBEL market, where the weight of hydro generation over a year is around 15%. Two types of behavior are thus contemplated in the manuscript: endogenous cartel and equilibrium.

External energy
It is assumed that the expected unsupplied energy by the units of the pool x ν 0 will be produced by external units and paid at a high price f 0 . Regarding the BCDCs, it will be considered that their external energy is the lowest probabilistically possible, thus an equality LMC involving all units (6d) will be imposed; its external energy is then a constant x ν 0β = s ν (Ω β ) that needs not be optimized.

Mathematical model of the stochastic mixed market
The complete model for optimizing the medium-term power planning production in a mixed market is presented in (6). Given that the stochastic parameters are modeled with a scenario tree, each variable and parameter refers then to a node ν of the scenario tree, which has an associated probability π ν . (5), the objective function is quadratic and indefinite. Constraints (6f) couple the total generation with the generation devoted to honor BCs. Constraints (6c) and (6e) are the LMCs for each LDC (system and BCDCs LDC). Constraints (6g) represent the single and multi-period non-LMCs, which are usually linear conditions over the total generation of some subsets of units; C λ and D λ are the coefficients, and the set H λ contains all the nodes on path λ. Note that non-LMCs (6g) may require the definition of extra variables, here represented as y ν . Recall also that x ν j refers to the expected generation by unit j over the period of node ν to match the system demand and the symbol is used to refer to the expected generation of units matching BCDCs.

The optimization procedure for the mixed model
Direct methods of global optimization [32] are the most appropriate for solving an indefinite problem such as (6). However, the current implementations of these techniques, such as the available branch-and-reduce based code described in [33] have not provided an acceptable solution in reasonable time for the smallest test cases used in this work. Optimization solvers non specialized for global optimization may have difficulties with problem (6), but there are ways, as those presented in the following subsections, that take advantage of the successive optimizations in the mLMH, leading to satisfactory results. These procedures are based on using different, though equivalent, formulations of (6) and on gradually changing the non-convexity of (6a) from less non-convex to its full non-convexity over the successive optimizations in the mLMH, where batches of LMCs are added.
The successive solutions in the mLMH use as initial point the solution of the former problem (plus default values for new variables associated to LMCs added). Default initial values are used in the first solution.

Solver employed and possible solution outcomes
The solver used is the publicly available Interior Point code Ipopt version 3.10.3 [34]. This code implements a filter method to determine step sizes and allows nonlinear objective functions and nonlinear constraints. Possible outcomes from the optimization of an indefinite problem are • Optimal solution (the normal solution within tolerances) • Acceptable solution (an optimal solution within slightly relaxed tolerances) • Restoration phase failure (a solution with acceptable primal and dual feasibilities that does not satisfy optimality conditions due to the objective function being close to semidefiniteness) • Iteration limit exceeded (no convergence within iteration limit; default limit 3000 iterations) • Infeasibility or other anomalous solution.
The solution outcomes have been checked after the optimizations in the application of the mLMH to the mixed model (6) and it has been observed that the last four outcome types may occur. It must be stressed that it is important that the last optimization in the mLMH (with all LMC batches incorporated) finishes as Optimal, but having had previous non-standard solutions could entail a risk of having wrongly chosen the ensuing LMC batch, which could lead to solutions not matching properly the predicted LDCs and BCDCs. In that regard it has been observed that, with some of the test cases employed, the sequence of plain optimizations in the mLMH led to a series of infeasible outcomes terminating in a spurious solution. That is why the strategies described below are aimed at avoiding the infeasible or anomalous solutions, and at reducing as much as possible the exceeded iterations outcomes.
In a process of successive optimizations, adding a batch of LMCs after each optimization and using the former solution as a starting point, interior point methods appear to lead iterated points more smoothly towards a feasible optimal solution than classical methods. These obtain a feasible point and an active constraint set first, and define directions using the projection of the objective function over the active constraints. Moreover, Ipopt's bi-criteria filter method (using a primal-dual merit function and a quadratic infeasibility function) for determining step sizes in a descent direction appears to be more efficient in obtaining an improved point than alternative solvers with classical single-criteria line search when having to minimize indefinite objective functions. The authors tried to use a classical nonlinear programming solver in the successive optimizations of the application of the mLMH for medium-term generation planning with BCs, and the outcome was that the solver halted due to lack of progress before reaching a solution.
Another example of use of an interior point nonlinear solver for electric generation problems is [22] where longterm generation expansion problems over a horizon of 30 years are addressed. Its procedure involves solving nonconvex optimizations subject to linear equalities and inequalities, and the solver employed is an interior point code based on successive sequential quadratic programming and a combination of classical line search and trust-region steps [35]. The main difference between the problem solved in [22] and that presented in this work is that in [22] the variables are the installed capacities of each generation technology, which change over the years due to the expansion planning, and in our work the variables are the energies generated by each generation unit over each period within a yearly horizon.

The difference of convex function decomposition of the objective function in the endogenous bilateral model
Using that, given two variables v and w, their product vw is (u + w) 2 − (u − w) 2 /4, which is a difference of convex quadratic functions, the indefinite objective function (6a) can be rewritten as a difference of convex quadratic (DcxQ) functions: s.t.: LMCs (6b-6e) and non- LMCs (6f-6i) where, as d and l i(ν) are negative, (7a) is a positive definite quadratic function and (7b) is a concave quadratic form.
There is an alternative equivalent formulation [32] to the former problem that uses an extra variable v RC : LMCs (6b-6e) and non-LMCs (6f-6i) where (8b) is a reverse convex constraint (RCC) because it makes the convex domain that it contains infeasible. Both the DcxQ problem (7) and the RCC problem (8) are non-convex and hard to solve, especially from initial points far away from the solution.
Let us also consider an easier problem: s.t.: LMCs (6b-6e) and non-LMCs (6f-6i) where the objective function (9a) is the positive definite quadratic part (7a) of the objective function of the DcxQ problem (7). Problem (9) is convex, easy to solve and has a unique solution.

Equilibrium solution of the medium-term bilateral planning
The NIRA procedure within the application of the mLMH is initialized during its stage 3), described in §2.3, and followed during the iterations of stage 4). The objective function of the NIRA optimization, based on the NIRA function [12] will be here extended to the stochastic formulation for the two observable participants: SGC and RoP. Using the constant terms calculated from the results x ν h ∀h ∈ H of a former solution, the equilibrium objective function, expressed as a DcxQ function (eDcxQ), and the equilibrium problem are: s.t.: LMCs (6b-6e) and non-LMCs (6f-6i) Comparing the DcxQ problem (7) and the eDcxQ problem (11) it can be noted that the later has more linear terms in the quadratic function and less quadratic terms both in the quadratic function and in the quadratic form, therefore the equilibrium problem is less non-convex and should be easier to solve than the endogenous problem.
An alternative formulation to problem (11) using a RCC and an extra variable v RC would be: LMCs (6b-6e) and non-LMCs (6f-6i)

Changing the non-convexity of the objective function
The degree of non-convexity of the DcxQ and of the eDcxQ objective functions can be easily changed by substituting the 4 in the denominator of the fraction d/(4T i(ν) ) in (7a) and in (11a) by a parameter ρ and giving to ρ values different from 4: the lower the value, the less non-convex will be the DcxQ and eDcxQ functions, (7a+7b) and (11a+11b). Starting with ρ = 2 in the early stages of the mLMH and gradually changing it to 4 in the later stages of it, the non-convexity of the problem is reduced during the first stages of the mLMH and recovers its proper degree towards the final stages.
There exists a value of ρ for which the objective function is just semidefinite, and falling sufficiently close to this value in some stage of the gradual change of ρ causes the solver Ipopt to abort the solution indicating the restoration phase failure outcome, mentioned in §4.1, with no important effect on the convergence of the solution strategy.

Solution strategies
A solution strategy means here a combination of successive problem modes, as described in subsections 4.2 and 4.3 in the application of the mLMH, along with changes in the non-convexity parameter ρ introduced in subsection 4.4. Solvers non -specialized in global optimization, as Ipopt is, may experience less trouble with convex objective functions subject to non-convex constraints, than with non-convex objective functions subject to linear constraints, This is the reason for changing problem modes depending on the former solution outcome obtained.
The purpose of a strategy is to avoid the risk of not generating LMCs that should be active in the optimal solution due to having landed, through an abnormal solution, on infeasible or spurious points.
Many strategies have been tried. The one presented next uses the same type of problem modes for the endogenous and for the equilibrium solutions and the same variation of parameter ρ in all test cases, and leads to satisfactory results.
The outline of the strategy applied in each stage (see 2. 3) of the mLMH in the endogenous problem is: 1) Let ρ := 2, substitute the 4 in (7a) by ρ, and solve (9).
In each loop of addition of LMCs change gradually ρ from 2 to 4, and solve current problem mode; if in mode (7) the solution outcome is Iteration limit exceeded or Infeasibility, change mode to (8); if in mode (8) the solution outcome is Optimal or Acceptable, change mode to (7).
In the equilibrium problem, stages 3) and 4) are of the same type as in the endogenous problem but different, due to the different objective function: 3) Solve (8) and compute terms (10).
In each loop of addition of LMCs, change gradually ρ from 2 to 4, solve current problem mode and update terms (10); if in mode (11) the solution outcome is Iteration limit exceeded or Infeasibility, change mode to (12); if in mode (12) the solution outcome is Optimal or Acceptable, change mode to (11).
The convex problem (9) has been used in the first two stages of the application of the mLMH for endogenous and for the equilibrium mixed-market problem. This is to ensure that stages 1) and 2) will terminate with feasible points satisfying the basic LMCs incorporated.
According to the authors' experience with using the Ipopt solver, the solution after having added a large batch of LMCs, as in stage 3) of the mLMH, is numerically more stable for modes with RCC (8,12) than with a DcxQ (7,11), this being the reason for employing the RCC modes in stage 3) and to start stage 4) and for maintaining a low value for parameter ρ. At each loop of stage 4), only a small batch of LMCs is added after each solution. Thus the DcxQ modes are employed and the non convexity is restored by increasing ρ to 4. The Restoration phase failure outcome only occurs occasionally when changes in the convexity of the objective function are forced.

Risk of profit loss aversion
The inclusion of risk-aversion terms in the objective function in order to minimize the value at risk of the market participants is implemented in practice by considering new variables and by adding value at risk constraints where a profit loss expression is employed [36]. Given that the losses are minus the profits and the profits are here an indefinite function, the value at risk constraints would also be indefinite, which would add further difficulty in solving the mixed-market problem.
The inclusion of the indefinite value at risk constraints is an area for further research now considered to be beyond the scope of this work.

Modeling language and server employed
Models and heuristics have been implemented as an AMPL [37] script. Programs were run on a Fujitsu Primergy RX200 S6 with two processors Xeon X5680 Six Core/12Threads (3,33 GHz), 12 MB cache and 96 GB RAM.

Characteristics of case studies
Realistic data from the MIBEL have been employed. The modeling fundamentals and details for stochasticity, features of code developed, periods employed, and hydro, WP, and pumping models are described in [14]. The planning horizon is one year decomposed into 6 periods of durations in days: 7, 23, 61, 92, 92 and 90, the starting date being 4 October. Loads (of 2010) to be satisfied amount to 157425 GWh, and generation units are considered in two degrees of disaggregation for different generation technologies. Hydro generation is represented by a simplified model of several basins and has 13,153 GWh of expected inflows. A 1350 MW pumped storage scheme, whose upper reservoir can store up to 340200 MWh, which is assumed to be half full both at the beginning and at the end of the yearly planning horizon, is considered in several test cases. As for BCs, the units of the SGC considered have about 47% of system capacity.
The numbers of generation units of cases are 18, and 24, where, in some cases, two of the units are related to pumping (a compensating unit and the generation one [14]). Three scenario-tree sizes with 21, 59 and 75 scenarios have been considered. Twenty-four possible problems (endogenous or equilibrium, with or without pumping, for 18 or for 24 units, using 21, 59, or 75 scenarios) have been successfully solved. A sample of these is presented in Tables 1 and 2. For comparison of computational requirements, two extended cases in generation units and in scenario paths have been also prepared, solved, and presented. Table 1 has the fixed dimensions of each case. The number of units, followed by a letter "n" for absence, or "P" for presence of pumping and terminated by the number of scenarios identify the case. Problem size is specified with numbers n u of units, n ν of scenario tree nodes, n var of variables, and number of equality and of inequality non-LMCs. Case 18n115 is the extension of case 18n75 with 115 scenario paths instead of 75, and case 36P75 is the extension of case 24P75 with 36 units. Given that some units considered in 24P75 were mergers of several units of the same generation technology, splitting some mergers into separate units yields case 36P75, which should produce similar but not equivalent results.

Results with the endogenous cartel and the equilibrium models
The 24 problem set (excluding the extended cases) solved through the solution strategy presented in Section 4.5 meant 454 optimizations, of which there were 12 acceptable outcomes, 4 exceeded iterations and 3 restoration phase  Table 2 also includes in the column LMCs the number of these generated by the mLMH. The Ipopt exe column contains the number of different executions of the solver (in the application of the mLMH), columns "ac", "ei" and "rf" stand for the number of executions that terminate abnormally as an acceptable, excess iteration, or restoration failure, and the "ite" column indicates the number of iterations. The CPU sec. column contains the AMPL plus solver execution times. The seven first lines of Table 2 contain the endogenous cartel results of problem (7) and the last ten lines contain the results of the equilibrium planning solutions of (11). The profits of the equilibrium solution are lower than those of the endogenous solution of the same case, and having pumping means extra profits [14].
The profits at an equilibrium solution are necessarily lower than those of an endogenous cartel solution because the equilibrium point is an endogenous cartel solution satisfying an extra condition: that the Nikaido-Isoda function of the problem be zero at the equilibrium point. The equilibrium point obtained through the application of the NIRA satisfies this condition and has thus lower profits than the endogenous cartel solution. The practical meaning of having zero valued Nikaido-Isoda function is that no market participant can increase its own profits by unilaterally changing the generation of its units. This can be appreciated in Table 2 by comparing the optimal profit results of each endogenous case and the equilibrium solution of the same case.
The results of the extended cases indicate that the number of units considered matter a lot in the CPU time required, as shown for case 36P75 that is seven times longer than case 24P75. Enlarging the number of scenarios, only dedicated to the stochasticity of renewables in the procedure presented, is not worth, as the results obtained for case 18n115 are practically equivalent to those with less scenario paths.
The results for each generation setting and number of units are always given for the three scenario tree sizes employed. For perfectly size-reduced scenario trees from the same generated tree, the expected profit should be almost the same for the same case with different reduced tree sizes.
The abnormal solutions that occur during the application of the proposed solution strategy appear not to trouble significantly the final results. The differences in expected profit for each case with different tree size present little variation, which can be caused by a non absolutely perfect tree reduction. Having obtained widely different expected profits would have meant convergence to different local minimizers. However, this not being so, is no guarantee that better local optimizers may not exist. Figure 3 shows three aspects of the solution to equilibrium case 24P21: the evolution of the energy kept in the upper reservoir of the pumped storage scheme along the nodes of each scenario at the top, the expected power productions of different technologies in the middle, and the evolution of the expected market price at the bottom.

Analysis of results of equilibrium case 24P21
As for the expected power productions, the WP generations of each period are totally determined by the WP scenarios generated [14]. The rest of generations are the result of the optimization. Given that the length of the periods is their duration in days, the areas of the rectangles in each period correspond to the energies generated by each technology. These results could be individualized for each generation unit in the data set used.
Regarding market prices, the results at the bottom of Fig. 3 show the evolution over the periods of the expected market price (black dashed line) and the range of the market price at peak load (in orange), and of market price at base load (in red). It is important to stress that these prices result from the influence of all stochastic parameters (hydro and WP), from the endogenous influence of all optimal unit generation policies (i.e., from the decisions of the participants in the market), and from the influence of all technical constraints, whether LMCs or non-LMCs, considered in the model. This is in contrast with existing procedures for determining short-term and medium-term hourly spot prices based solely on records of past spot prices and on future forward prices, and on the different treatment of the observed hourly price regimes of base, of upward spikes and of downward spikes [38]. Figure 4 shows the mean yearly energies of load, pumping, and of all technologies of the generation mix for the market and for the BCs in each scenario of the solution to equilibrium case 24P21. One important feature is that WP energy is split in two parts between market and BCs, the same as with thermal generation. Technologies with lower capacities, as hydro and pumped hydro generation, go almost entirely either to the market or to BCs in the solution obtained.

Conclusions
A stochastic model has been presented for the medium-term planning of GenCos participating in a mixed pool and BC market including pumped storage under the behavioral principles of endogenous cartel and equilibrium. The mLMH was employed for LDC matching and the NIRA for obtaining an equilibrium solution.
The objective function of the medium-term mixed-market problem is non-convex, but its solution using a current implementation of global-optimization direct methods proved not to be practical. A procedure based on the use of an interior-point nonlinear programming solver using a filter step-size reduction technique has been presented, where a solution strategy was employed to circumvent the likely trouble caused by non-convexity in the sequence of optimizations of the mLMH. Computational results show that the results are satisfactory. Figure 3: Solution to equilibrium case 24P21 showing the evolution over the periods of the energy stored in the pumped storage reservoir for each scenario (above), evolution of expected power productions of different technologies: "thr" thermal, "hyd" hydro, "phy" generation from pumped hydro, "wp" WP, and "pmp" pumping, in negative (middle) and evolution of expected market price (dashed line) with maximum price range in orange and minimum price range in red (below). Abscissas: time in days. Figure 4: Yearly expected power of pumping "pmp", in negative, for each scenario in solution to equilibrium case 24P21, with system load (1st column), and generation mix by technologies "thr" thermal, "hyd" hydro, "phy" pumped hydro generation and WP, in market (2nd column) and in BCs (3rd column).
The model put forward is a useful tool for medium-term mixed-market generation planning. Expected profits for the SGC and for the RoP are obtained from the results, together with the optimal hydro and pumping use over the horizon considered. From the optimal expected energies for each unit in each period the medium-term fuel procurement decisions can be readily calculated. The results also provide the split between market and BCs in expected generation and duration of each unit. The stochastic medium-term planning described can be also employed as a reliable way to measure the effect of a given penetration of non-dispatchable renewables (wind power) on the market share by other generation technologies.
The models presented and demonstrated are useful both for a price-maker participant in the mixed market and for a price-taker; they are also useful for the Market Operator. Price-makers can determine their generation planning and the splitting of generations for market auction and for BC honoring in each period. Price-takers will learn which units and for how long to use for satisfying BCs, and how much energy can they expect to produce for the market in each period. The Market Operator will be able to find, from the results obtained, whether the behavior of the market participants is closer to an endogenous cartel or to an equilibrium.