Application of the Level Method for Computing Locational Convex Hull Prices

Convex hull pricing is a well-documented method for coping with the non-existence of uniform clearing prices in electricity markets with non-convex costs and constraints. We revisit primal and dual methods for computing convex hull prices, and discuss the positioning of existing approximation methods in this taxonomy. We propose a dual decomposition algorithm known as the Level Method and we adapt the basic algorithm to the specificities of convex hull pricing. We benchmark its performance against a column generation algorithm that has recently been proposed in the literature. We provide empirical evidence about the favorable performance of our algorithm on large test instances based on PJM and Central Europe.


A. Pricing Non-Convexities
T HE classical analysis of an economic dispatch problem, together with its dual, provides a fundamental argument for uniform pricing in electricity markets [1] -an optimal dispatch can be supported by a set of competitive equilibrium prices. In other words, even if a central authority cannot effectively control the dispatch of the assets itself, it can provide prices that align the behaviour of selfish profit maximizing agents with social welfare maximization. However, as the argument assumes convexity of the dispatch problem, a fundamental challenge for market efficiency is non-convexity, as the latter implies that it is not guaranteed that a competitive market equilibrium exists.
Non-convexities are at the heart of power system operations [2], in terms of both the network model as well as in the market orders: (i) they are present in the alternating current (AC) power flow equations which characterize the physics of the grid and (ii) in the mixed integer programming (MIP) constraints that describe the market offers. As the day-ahead (DA) markets in Europe and in the US rely on a linear direct current (DC) power Manuscript  flow model of the grid, point (i) is not encountered in these markets. 1 On the other hand, point (ii) is a reality in both US markets that rely on solving a unit commitment (UC) problem, as well as in the EU market which includes integer market ordersthe so-called "block orders". Throughout this paper, we neglect (i) and rather focus on (ii). The inexistence of equilibrium prices in electricity auctions has triggered a long-lasting debate on the choice of an appropriate pricing scheme in the presence of non-convexities. Convex hull pricing (CHP) has arisen as one promising alternative: while being so far mainly debated in the US, it has also recently emerged as a possible option for the EU market [4]. A practical concern of CHP is that its computation can be challenging (e.g. see Issue 7 in [5]). Our paper aims at addressing these computational challenges by putting forward a workable algorithm (the Level Method) for realistic instances subject to network constraints. In the remainder of this section, we sketch the main concepts related to CHP and we discuss the context of non-uniform pricing discussions in the EU. Insofar as the EU market is concerned, we discuss institutional aspects as well as computational issues, which motivate our choice of test instances.

B. Non-Uniform Pricing Schemes
The most widely debated "non-uniform pricing schemes" in the literature include integer programming (IP) pricing proposed by O'Neill [6], convex hull pricing proposed by Gribik and Hogan [7], [8], and "extended LMP" pricing which has been applied early on in the PJM market [9], [10]. They all amount to a convex relaxation of the market clearing problem. These strategies consist of combining a uniform electricity price with discriminatory payments, called uplift payments, which aim at restoring the incentives of market participants for following the market matches. In this framework, the overall market clearing procedure can be described in three steps, which are also followed by our simulations: 1) Solve the primal problem, in order to establish the dispatch and commitment instructions; 2) Solve a pricing problem in order to compute uniform electricity prices; 3) Solve the independent profit maximization problems of all market agents (generators and the network operator) in order to establish uplift payments. Regarding step 3, it is worth noting that uplift payments are often categorized as follows in the literature: 1) Potential Congestion Revenue Shortfall are uplifts associated with the network congestion revenue [3]. 2) Generator side-payments are defined as the difference between the maximum profit achieved by self-scheduling given the market prices and the as-cleared profit. The pricing strategy proposed by O'Neill is a common choice in non-convex settings. We also use it as a benchmark for our simulations. However, it does not attempt to minimize uplift payments, and can therefore possibly lead to high side payments. Uplifts are undesirable, as they can distort the incentives of bidders or create revenue adequacy problems for the market operator that needs to finance them [11].
These concerns motivate Convex Hull Pricing (CHP), the main property of which is to minimize uplifts. Because it is computationally challenging, PJM (and other US ISOs) has recently implemented a new pricing scheme, referred to as "extended LMP" which is more tractable computationally than CHP. For certain forms of simple market orders, it can also be shown to be a reasonable approximation of CHP [9]. We expand on how it relates to the computation of CHP in Section II.

C. Uniform Pricing in the EU
The EU market landscape presents a number of major institutional differences compared to US markets [12]. One such notable difference is that day-ahead energy auctions are operated by for-profit Nominated Electricity Market Operators (NEMOs) while, in the US, it is the (typically non-profit) ISO that operates both the market and the network. One implication of this difference relates to the ability of the market operator to socialize uplift payments. This difference may, in part, justify the currently employed "uniform" pricing scheme that is adopted in Europe, as implemented in Euphemia, the algorithm that clears the pan-European day-ahead auction [13].
In Euphemia parlance, the aforementioned generator side payments can be related to: (i) paradoxically accepted blocks (PAB) -cleared bids actually facing losses, i.e. requiring make-whole payments (as defined in [5]) -and (ii) paradoxically rejected blocks (PRB) -a rejected bid that would have been profitable, i.e. facing a lost opportunity cost. The EU day-ahead market "avoids" uplift payments by (i) constraining the problem by not allowing the acceptance of PABs while (ii) allowing PRBs, but not paying their lost opportunity costs. Ultimately, it does not effectively reduce the uplifts to zero, but it guarantees zero make-whole payments, while increasing the total lost opportunity cost and not paying it. Consequently, this pricing scheme only outputs uniform prices while it does not provide the market participants with any discriminatory payments. This justifies why, in EU NEMO parlance, it is referred to as uniform, in contrast to the three non-uniform pricing schemes that are discussed previously.
This uniform pricing scheme involves "primal-dual" constraints that implicate dispatch and price decisions in a single market clearing model. The solution implemented in Euphemia amounts to an iterative algorithm that matches market orders while aiming to find a feasible price (without PAB). If this is not possible, the algorithm generates a cut in the primal model and repeats the process. In contrast to the non-uniform pricing schemes that work in three steps (dispatch, price, uplifts), the EU uniform pricing scheme works as a single -but iterativestep, and couples dispatch and price problems together.
This makes the problem that Euphemia is called to solve (a mixed integer quadratic program subject to complementarity constraints) computationally challenging. Moreover, the approach deteriorates market welfare, since welfare-enhancing orders can be discarded if no market clearing price can be found to support the aforementioned clearing rule. For these reasons, non-uniform pricing schemes, and in particular convex hull pricing, have recently received consideration by the European NEMOs as a possible option for the European DA energy auction [4]. Considering the aforementioned institutional EU structure, as well as the algorithm implemented in Euphemia, this would constitute a disruptive market design evolution.
2) The market model includes a network of ∼ 40 bidding zones, and its geographic footprint is expected to be further enlarged.
3) The market model is expected to move towards 15-minute granularity in the near future (a horizon of 96 periods). Forty bidding zones for ninety-six periods implies a 3,840dimensional price space. These requirements motivate the considered use cases in Section IV.

D. Contributions and Structure of the Paper
The contribution of the paper is twofold: 1) We propose the Level Method [14] for computing CHP and adapting it to the specificities of our problem. We specifically adapt the algorithm in order to exploit the convexity of the network model. We further introduce a "multi-cut" variant of the Level Method in order to leverage the separability of the sub-problems. Note that two types of approaches have been envisioned in the literature for solving CHP: dual approaches and the primal approaches (we define these in Section II). The Level Method belongs to the former. Primal approaches, and their drawbacks which motivate our choice for a dual approach, are presented in Section II. The review of alternative (tested) dual approaches comes in Section III and motivates our choice of the Level Method. 2) We efficiently solve CHP, using the Level Method, for large instances including a network and a horizon of 96 periods, which anticipates the evolution of the EU market. We conduct a critical comparison of our approach against both primal and dual decomposition approaches.
In particular, we compare it to a notable recent publication by [15], which proposes a Dantzig-Wolfe (D-W) algorithm for computing CHP. The D-W algorithm exhibits favorable performance on a test case without a network and with 24 time periods, as considered in [15]. Given our preoccupation with a market clearing model at the scale of the EU market, the question becomes how the method scales when moving from a 24-dimensional to a 3,840-dimensional price space. When increasing the dimension, the Level Method is empirically shown to attain favorable performance relative to [15]. Our paper is inspired by an older unpublished work [16], and is further motivated by [15]. We describe the mathematical formulation of CHP in Section II. We then introduce the Level Method in Section III. In Section IV, we test the algorithm on multiple large instances and compare the results with D-W. Section V concludes and discusses areas of further research.

A. Convex Hull Pricing Program
We define the dispatch problem subject to network constraints as follows: Here, G i denotes the set of generators (or market offers) at node i. Each offer is modelled with a total cost c g , a power output p g,t at time t and a set of non-convex constraints X g . The generic variables u g stand for all the binary variables encountered in the generator model. The demand at time t and node i, D i t , appears in the market clearing (MC) constraints (1b). Regarding the network, f l,t stands for the flow on line l, while from(i) is the set of lines originating from i and to(i) the ones directed towards i. No assumption is made on the network constraints F, except that it is a convex set.
Each generator g is assumed to be a selfish agent that maximizes profit, i.e. solves the following program: Here, i(g) stands for the node of generator g, while π i(g) t represents the market price of node i(g) at time t. A fundamental result [7], [8] on CHP establishes that minimizing uplifts amounts to solving the following problem: Here, L(π) denotes the Lagrangian dual function, obtained by relaxing constraints (1b) of problem (1): We recognize in (4b) the profit maximization problems (2) of the generators. As established in [8], using the optimal primal dispatch solution of (1) and injecting it into (4) clarifies why the previous Lagrangian problem does indeed minimize the uplifts. As also pointed out in the literature, the definition (3) of CHP also indicates that the uplifts can be interpreted as the duality gap between (1) and (4).

B. US Versus EU Models
In addition to institutional differences between US and EU markets, another major difference relates to the definition of market products. The US markets follow a unit-bidding model, where each unit is represented explicitly in the market, along with its technical characteristics. On the other hand, the EU day-ahead market follows a portfolio-bidding model (which cannot be subsumed in the unit commitment formulation), where each agent submits multiple generic market orders that represent the portfolio of its assets in an aggregated way. These market orders [13] include convex hourly ordersstepwise and interpolated curves -as well as non-convex orders 2mainly the family of block orders. The latter is a financial order spanning over multiple periods and involving a binary acceptance variable.
Model (1) remains general regarding the bid (generator) constraints (1c), which are simply represented as the non-convex set X g . This implies that the approach outlined in this paper can accommodate all the flavours of unit commitment models as well as the EU-like auctions. This exceeds what a "primal CHP approach" can model.
Finally, model (1) considers a general (but convex) set of network constraints F. Our approach can in fact accommodate any convex representation of the network. In both the US and EU market, F would amount to a set of linear constraints, the main difference being that certain US markets are nodal (larger number of nodes) while the EU market is zonal (roughly one zone per country). We remark in Section III on the specific treatment of the network in our proposed Level Method.

C. Primal and Dual Approaches for Computing CHP
In this section, (i) we locate the Level Method in the perspective of the landscape of all the alternative of approaches for solving CHP and (ii) we motivate the choice of a dual approach in the light of the limitations of the primal approaches.
As noted in Section I, there are two main approaches envisioned for computing convex hull prices -i.e. solving problem (3): (i) the Lagrangian dual approaches, which directly attempt to maximize function L(π) using an iterative algorithm, and (ii) the primal approaches, understood as methods that seek to describe the CH of the non-relaxed constraints (1c)-(1d) by developing tight formulations. Fig. 1 outlines the landscape of approaches for computing convex hull prices. The top problem (A) corresponds to the dispatch problem (1). Below, on the left, we find primal relaxations of (A) while, on the right, we find Lagrangian relaxations of (A) -Lagrangians are indeed a widely employed method for deriving convex relaxations of non-convex programs [17]. The problem (Γ) corresponds to the CHP definition (3), which can be solved by dual decomposition approaches such as the Level Method. Problem (Γ) maps to its primal equivalent in (C). The underlying idea of the primal formulation is that computing the CHP as the Lagrangian multipliers of (3) is equivalent to computing the dual variable π associated to the market clearing constraint (1b) in the primal problem (1), if the latter is expressed on the convex hull of its domain -i.e. conv(X g )∀g ∈ G (see [18], [17] for the general result in Lagrangian relaxation theory or [19] for the specific result related to CHP).
Although (C) is the tightest primal relaxation of (A), there exist looser relaxations, such as (D), which amounts to relaxing the integrality constraints u g,t ∈ {0, 1} to u g,t ∈ [0, 1]. This corresponds to PJM pricing, discussed in the introduction. PJM pricing can be interpreted as a computationally efficient approximation of CHP [9]. In certain cases, relaxing the integrality constraints in X g may provide conv(X g ). In this case, problems (C) and (D) are equivalent and PJM pricing effectively corresponds to convex hull pricing. The fact that relaxation (D) is looser than (C) implies that the duality gap between (A)-(D) will be greater than or equal to the one between (A)-(C).
Interestingly (and to the best of our knowledge, unnoticed in the literature), one can also relate the primal relaxed problem (D) to its Lagrangian dual counterpart (Δ). While CHP is solving the partial Lagrangian dual relaxation (Γ), PJM pricing corresponds to solving the full (looser) Lagrangian dual relaxation (Δ), where all the constraints -and not only the market clearing constraints -are dualized. 3 Regarding the primal CHP problem (C), a way to approach it is to develop a tight -but custom -formulation, specific to the targeted problem (A). Recent researches have embraced this idea: [19] proposes an explicit formulation for the primal model of CHP for classical UC constraints. Madani [24] analyses primal CHP formulations for the constraints of the European day-ahead market clearing model. 4 More recent research further elaborates on the idea, developing tight (custom) formulations for MISO [25] or proposing a network flow model of unit commitment, in order to compute CHP for a broader set of constraints [26]. One value of the primal CHP approaches is to establish the link between CHP theory and the literature dedicated to tight formulations of UC polytopes such as [27]- [35] (see chapter 2 in [16]). Similarly, when including a non-convex network model, the primal CHP approach [3] also establishes the connection between CHP theory and SDP/SOCP relaxations of AC power flow [2].
Nevertheless, as also voiced in [15], there are certain constraints for which the convex hull is not tractable in the sense that it may not be possible to characterize the convex hull with a scalable number of constraints. This already holds for simple ramp constraints [19]. This is also acknowledged by [26], where the authors do not account for them in their network flow model. Instead, [25] needs to combine the proposed tight formulation with an iterative algorithm in order to account for the ramp constraints in a scalable way. It goes without saying that these modelling limitations also hold for more advanced constraints such as multimode CCGT units, detailed battery models, and so on. Thus, since the pricing mechanism becomes dependent on the quality of the primal formulation, the primal approach can be ruined by adding a new constraint -which is particularly concerning, since electricity market models are constantly subject to changes (e.g. triggered by regulatory requirements such as article 40 of EGBL guidelines). These modelling limitations imply that, if the representation of the convex hull is not tractable, the primal approaches are irremediably left with an approximation of convex hull prices, such as the PJM pricing model (D). This is illustrated in our numerical results of Section IV, where a primal method benchmark [19] is included. This motivates our choice for a dual approach. 3 Taylor [2], which inspired Fig. 1, proposes an interesting interpretation of CHP by relating it to the semi-definite programming (SDP) relaxation of problem (1). The proposition is motivated by the well-known SDP relaxation of a nonconvex quadratically constrained program (QCP) [20]- [22] and the fact that a MIP can be expressed as a QCP. However, the above taxonomy reveals an innacuracy in the reasonning: it mixes (Δ) and (Γ), as it omits the fact that CHP relies on a partial (and not complete) Lagrangian relaxation, where only the market clearing constraints are relaxed (i.e. dualizing fewer constraints can only improve the duality gap [23]). 4 Note however that [24] focuses on a subset of the market constraints, ignoring e.g. linked blocks and exclusive groups [13].

A. Review of Existing Algorithms
The appropriate algorithmic scheme for solving (3) is related to the type of function L(π).
Property 1 (Concave): Function L(π) is concave in π. Property 2 (Non-smooth): Function L(π) is a non-smooth (piecewise linear) function, i.e. each facet can be seen as corresponding to a set of binary (commitment) decisions u g .
Property 3 (First-order oracle): A first-order oracle is available, i.e. given a price π, both the function value L(π) as well as its supergradient s ∈∂L(π) can be evaluated.
Property 4 (Supergradient): Let (c * , p * , u * , f * ) be the optimal reactions to π (solving respectively (4b) and (4c)). Then is a supergradient of L in π; i.e. s ∈∂L(π). Each call to the oracle implies solving MIP profit maximization programs (2) for each generator as well as for the network program (4c) -these are thus slave problems. We propose later a special treatment of the network and leverage the separability of the profit maximization problems in order to substantially improve the formulation. Any algorithm tackling this problem would work in three steps: 1) Given a price π k , evaluate L(π k ) and∂L(π k ); 2) Given this information, generate a new price π k+1 ; 3) If the stopping criterion is met, stop. Else, go to step 1.
The main difference between dual decomposition algorithms is in the way that they construct the sequence of iterates {π k } ∞ k=0 : (i) some algorithms simply update the prices based on the latest supergradient information -they are memoryless -; (ii) while other algorithms will keep memory of the sequence of iterates. We briefly summarize three approaches, which were tested (and compared to the Level Method) by the authors in previous work [16].
A well-known scheme belonging to category (i) is the subgradient scheme. Perhaps surprisingly, it is proven to be optimal for general convex non-smooth optimization with arbitrarily high dimension [14]. However, when dealing with problems of "moderate" dimension such as the one presented in our context, there exists more optimistic alternatives.
Indeed, the subgradient scheme for piecewise-linear functions, such as our problem (3), tends to oscillate between the facets of the Lagrangian dual function, around an edge. Therefore, one idea is to "catch the edge" and follow it until the optimum, instead of oscillating from one facet to another, as the subgradient method does. This intuitive reasoning leads to the Extreme-Point Subdifferential (EPSD) algorithm, which is specifically applied in [36], [37] to the CHP problem. However, our experiments in [16] reveal that each iterate of the algorithm is costly, as it requires not only to solve the problems (2) for each generator to optimality, but to enumerate all the optimal solutions.
Unlike these two memoryless schemes, the Analytic Center Cutting Plane Method (ACCPM, see [14], [38] for the theory and [37], [39] for its application to CHPs) is based on the principle of iteratively reducing the search domain: the price domain is initially limited to a box and, at each iterate, the supergradient is used for generating a cut, which shrinks the search domain. The next testing point is chosen as the analytical center of the updated domain.
Our original investigation of these alternative dual approaches (subgradient, EPSD and ACCPM) in [16] concluded that none of them were competitive with the Level Method for computing CHPs.

B. Kelley's Approach
The Kelley algorithm [14] forms the basis for the proposed Level Method. It is based on the idea of iteratively constructing a model (upper approximation) of the Lagrangian function L(π), using its supergradients. Definition 1 (model function): Let Q be the initial domain of our problem (i.e. a box limiting the prices, which can be economically interpreted as price caps) and let {π k } ∞ k=0 be a sequence in Q. Let s k be the supergradient at iterate π k . Then is a model for the Lagrangian function L(π), such thatL(π, k) ≥ L(π).
In order words, the piecewise linear function L(π) is upperapproximated at each iterate by a model functionL(π, k) consisting of supporting hyperplanes. At iteration 0, this is a single hyperplane. Then, as the iterate count k is increasing, the model functionL(π, k) is becoming increasingly accurate.

Definition 2 (master program):
The maximization of the model function yields the master program at iterate k: Here, s j are the "cut coefficients" (as defined in Property 4) and b j = L(π j ) − s j , π j are the "cut constants". This is a computationally tractable linear program.
Having the upper-approximation functionL(π, k) at hand, one needs to decide the rule for building the sequence of iterates {π k } ∞ k=0 . The more intuitive way to pick the next iterate is: i.e. the solution of the master program (6). This defines Kelley's cutting plane method. One of its benefits is that it explicitly provides an upper bound as well as a lower bound at each iterate k: a lower bound is defined as LB k = max j=0..k L(π j ), while an upper bound is UB k = max πL (π, k). Note that the sequence of upper bounds {UB j } k j=0 is decreasing, as the definition of the model function implies thatL(π, k + 1) ≤L(π, k). The upper and lower bounds can be combined to define the relative gap, which is used as a stopping criterion for the Kelley (and Level) Method:

C. Level Stabilization
Kelley's algorithm is finite, because each iterate adds a new hyperplane and the number of hyperplanes supporting the function is finite. Nevertheless, despite its simplicity and its good behaviour in low dimension, it tends to be unstable in higher dimension. This is due to the unstable nature of piecewise linear functions: adding a new supporting hyperplane can move the optimum far from the previous point (i.e. to a corner of the box Q). This well-known drawback [40] justifies why multiple stabilization approaches have been proposed in the literature, including the Level Method [14], [40].
The underlying idea of the Level Method is to update prices more smoothly: instead of using the optimum of the model function as the next iterate, the algorithm chooses π k+1 such that it is "better" than π k (as evaluated by the model function L(π k+1 , k)) without being optimal at all costs. We observe in Section IV that this stabilization has a major influence on the practical performance of the algorithm.
A graphical illustration in 1-D is presented in Fig. 2. The cuts, the LB and the UB are obtained as in Kelley's method, by solving the master program (6). However, unlike in Kelley's method, the next price candidate is selected by solving a projection program.
Regarding the calibration of α, α = 1 corresponds to the classic Kelley method, while α = 0 implies that the iterate simply does not move. We note that a theoretically optimal α exists [14] for general convex non-smooth functions, but that a calibration to the specific problem can still be meaningful. Our empirical tests on the CHP problem reveal that, for the high-dimensional instances that we are interested in, the approach is largely insensitive to the choice of α. This is shown later in Table III, where any value of α between 0.2 and 0.7 exhibits similar performances. Following [16], the value α = 0.2 is chosen for all of our experiments in the present work.
Regarding the choice of the box Q, experimental evidences show that the Level Method is not too sensitive to its exact value, despite it impacts the quality of the UB estimate. In all of our experiments, Q is initially set to ±300$/M W h and is then progressively shrunk after 10, 20 and 30 iterates to ±25$/M W h around the latest price candidate. This is justified by an analysis of the volatility of the price iterates, which rapidly reach a price close to the CHP.

D. Refinements of the Level Method in the Context of CHP
We now propose adjustments to the basic algorithm which exploit the structure of our problem. We specifically leverage the fact that: (i) the network model is convex and (ii) the profit maximization programs of the generators are separable.
In our development so far, we have been treating the convex network term (4c) identically to the non-convex generators, i.e. by solving the network profit maximization given a price π, and generating a supergradient. We illustrate below the treatment of the convex parts of the primal program by focusing our discussion on the network. The idea applies identically to convex generators (e.g. the convex orders in Euphemia, which are numerous), a convex pumped-storage model, etc. (see section 3.6 and appendix A in [16] for a treatment of these cases).
For the sake of illustration, let us assume that the network constraints F correspond to the DC (voltage angle) power flow. Term (4c) then reads as follows: Here, B l stands for the susceptance of line l, and F l and F l are its max and min capacity, while or(l) and dest(l) denote the origin and destination nodes of line l. The dual of (10) can be expressed as: Problem (11) can now be injected into (4) as a substitute for (4c), meaning that the network dual variables (μ, ν, λ) would explicitly be variables of the master (and projection) program. Secondly, the classical Kelley and Level Methods add a single cut at each iterate, namely one single cut for all the generators. Nevertheless, the dual function is separable with respect to the generators. We therefore propose a multi-cut Level Method, whereby we compute one cut (one lower approximation) for each generator profit maximization subproblem. Our experiments reveal that this adaptation can deliver substantial computational benefits. Generating more cuts makes the model function more accurate, which enables the algorithm to converge faster. Note that multi-cut versions of other approaches have been applied successfully in different contexts, such as for two-stages stochastic programs [41], [42].
To summarise, after the inclusion of both the network dual and the multi-cut approach, the master program (6) at iterate k becomes: Here, {p j g } k j=0 , corresponds to the sequence of generator g power output for iterates j = 0..k. These parameters are also cut coefficients for generator g. On the other hand, {c j g } k j=0 , which corresponds to the sequence of generator g cost for iterates j = 0..k, are the cut constants. The translation of the projection program (9) is applied as discussed previously.
In the classical Kelley/Level Methods, estimating the lower bound (evaluating (4) at a given π) follows directly from the resolution of the slave subproblems. The inclusion of the network into the master program, as described above, complicates the process. Indeed, the network contribution in the dual function (4c) is not solved explicitly anymore, but now comes in the master objective (12a), together with constraints (12c) and (12d) that should not be violated. Therefore, estimating the value of L(π) after having retrieved the cuts from the slaves (for the same π) amounts to solving the master (linear) program (12) with the variables π fixed. The overall procedure is described schematically in Fig. 3. Note that the resolution of the two master programs (with π fixed and variable) can be parallelized.

IV. SIMULATION RESULTS
This section presents the numerical results of the (multi-cut) Level Method on instances of realistic scale. The Level Method has been benchmarked against other dual approaches in earlier work by the authors [16]. It is chosen as the most promising method for computing CHPs among all tested alternatives. In the present section, we therefore focus on its comparison with a recent work [15] which employs a D-W column generation algorithm [43] (i.e. the dual of Kelley) for iteratively building the convex hull of the dispatch problem, i.e. D-W gradually discovers the corners of the primal formulation. As in the case of the Level Method, it can be applied to any UC formulation. We use it as a performance benchmark in our analysis, due to its favourable empirical performance. We also include O'Neill pricing, discussed in the introduction, as another benchmark in our analysis, as well as PJM pricing (discussed in Section II-C) as a primal method benchmark.
Unlike other computational researches on CHP [15], [37] which are mainly concerned about the number of generators in the problem, we rather focus our investigations on the sensitivity of the algorithms with respect to the dimension of the price space. Indeed, although the number of generators is surely relevant, since the ultimate goal is to compute prices by optimizing L(π), the price-space dimension is expected to have a significant impact on the performance of any tested method. Therefore, we first present results without a network, with a horizon of 24 periods, and then introduce network constraints and extend the time horizon to 96 periods.
For all our test cases, the comprehensive market procedure for computing the prices and measuring uplifts follows the steps that are described in Section I. Concretely, there are three steps: dispatch, price, and uplift computation. The Level Method and D-W differ with respect to the second step. Both approaches have been implemented in Julia (JuMP) and all the tests are run on a personal computer (Intel Core i5, 2.6 GHz with 8 GB of RAM) using Gurobi 9.1.1.

A. FERC (US) Test Cases
The first test cases in our analysis are based on FERC datasets [34], [44]. The test sets are publicly available, together with the associated UC model, and are also used by [15]. These test cases consist of a detailed UC model. The only adaptations in our work are the removal of reserve and the netting out of renewable supply from the load. The UC model includes, among others, min up and down time constraints, ramp constraints (including start-up and shut-down ramp rates), variable start-up costs which depend on how long a unit has been off, no-load costs, and piecewise linear production costs. The model has no network, but gathers > 930 generators. This corresponds to an instance of realistic size, barring for the absence of the network. As in [15], we conduct our analysis on a 24-period horizon with hourly time step.
Table I presents the average results over 11 FERC instances, while Fig. 4 illustrates the convergence behaviour of both approaches on one of the instances. The 11 instances essentially correspond to 11 different load profiles, with slight changes in the production fleet, which varies from 934 to 978 generators. The stopping criterion of the Level Method (8) is set to 0.01%. The number of iterates reported in Table I for D-W corresponds  TABLE I  RESULTS OF THE LEVEL METHOD AND THE DANTZIG-WOLFE ALGORITHM ON  FERC DATASETS (AVERAGE OVER 11 INSTANCES) a (·) denotes the average time per iterate for solving the "master programs" (i.e. master plus projection in the case of the Level Method).  (Table I).
It should be noted that this number of iterates is a reasonable measure for comparing the performance of both approaches. Concretely, both methods have to solve the same subproblems and mainly differ in the other computations that they are required to perform. Whereas the Level Method has to solve both a linear master and a quadratic projection, D-W is only required to solve the linear (master) extended formulation. On the other hand, the extended formulation solved by D-W is larger than the Level master program, as illustrated in Fig. 5. Overall, this results in a similar run time per iterate, as reported in Table I which shows both the average run time per iterate as well as, between parentheses, the average run time spent in the master programs (master plus projection programs for the Level Method). This implies that the number of iterations (the analytical complexity: the number of calls of the oracle to reach a reliability target) is a reasonable measure for comparing performance. It also has the benefit of being less dependent on the specific machine or on the implementation details. Note that, for both approaches, the slave subproblems can be parallelized.
Furthermore, there is a gain in robustness: the Level Method exhibits a more stable performance, as observed in Fig. 4. Indeed, Fig. 4 suggests that it does not seem possible to stop the D-W algorithm long before its termination, since uplifts remain high for a large number of iterations (we also refer the reader to Fig. 7 of the next use case, which shows how the convergence of uplift over iterates translates to the distance of prices from CHPs). Instead, the Level Method reaches near-optimal prices in fewer iterations. This is an inherent advantage of the Level Method, which is by design a stabilization approach.
Finally, we comment on the primal method benchmark. The FERC model exceeds what a primal CHP approach such as [19] can model, since it includes ramp constraints and time-dependent startup costs. The integer relaxation is therefore expected to lead to an approximation of CHP. The quality of the primal method largely depends on the tightness of the formulation. In this respect, the FERC model is derived from a careful review of the literature dedicated to tight formulations of the unit commitment model [27], [35]. The quality of the model is discussed in [34], where it is accompanied by computational experiments of its tightness. As observed in Table I, the primal method turns out to provide a close approximation of CHP on these FERC instances. Nonetheless, this is not always guaranteed, as we observe in the next test case (Table IV), where the primal method leads to an average uplift which is ∼60, 000€ higher than CHP, for a market of comparable dispatch cost.
The test cases analysed so far suggest a promising performance for the Level Method. Nevertheless, even if these FERC instances are of realistic scale insofar as the number of power plants are concerned, we are interested in computing prices. This suggests that it is the dimension of the price space that matters the most. There are essentially two ways 5 to increase the price dimension: (i) augmenting the time horizon -the horizon of future EU markets will be 96 periods of 15 minutes -and (ii) adding a network -which is unavoidable in both the EU and the US markets. This motivates the next test cases.

B. EU Test Cases
We now extend our analysis to use cases with a network. The EU dataset that we utilize is the one used in [45]. The network data is based on [46], and is constructed among others from an ENTSO-E database. The market suppliers are modelled as a slightly simpler version of the UC model than the FERC test case, essentially simplifying the cost structure: there is a single start-up cost, instead of the variable start-up costs of FERC, and the marginal production cost is constant. All the cases are simulated over 6 different load profiles. As we are interested in studying the scalability of the Level Method and D-W algorithm with respect to the network and the time horizon, the data has been aggregated into two test cases: BE and BE-NL, which are described in Table II. As detailed in Section I, Euphemia, the EU market clearing algorithm, currently computes prices for ∼40 bidding zones, and is expected to move to 15-minute granularity (96 time periods) in the near future. This makes our two tests cases with 96 periods very relevant proxies of the evolving EU context with respect to price dimensionality.
The final results are obtained with the stopping criterion set to 0.01%, as for the FERC cases. Table III shows the sensitivity of the Level Method towards parameter α, previously discussed in Section III-C.   Fig. 6. We observe that, within 6 iterates, it already reaches a price that achieves lower uplifts than those of O'Neill pricing. Fig. 7 also presents the convergence of both algorithms on the same instance in terms of price distance to the optimum. Being capable to reach quickly decent price candidates is an attractive feature for EU implementation of CHP, recalling from Section I that Euphemia is currently granted 12 minutes for computing the EU day-ahead market matchings and prices.
As far as the network size is concerned, Table V presents the sensitivity with respect to the two use cases. Perhaps surprisingly, neither of the methods is strongly affected by the size of the network, rather the contrary. On the instance with 24 periods, the benefits of the Level Method are similar as in the FERC cases. On the 96-period instances, the Level Method moves from five to four times faster than D-W on the BE and BE-NL cases, in terms of iteration count. Overall, D-W seems much more affected by the increase in the time horizon rather than the presence of a network, to which the test cases suggest D-W is rather robust. This is possibly due to the fact that D-W is required to explore in the space of promising power plant schedules -and these schedules become more and more numerous when increasing the horizon -while the network size does not affect immediately the number of schedules.
It should be stressed that the aforementioned computational gains can make a difference for the practical implementation of CHP, keeping in mind the 12-minute run time limit of Euphemia. From Table IV, we observe that the Level Method requires less than 5 minutes on average for solving a 96-period instance. The D-W algorithm requires 27 minutes.
The computational times reported in our results may of course not be representative of the implementation of the EU NEMOs, as solving the slaves in parallel and increasing the computational power would reduce the run time. Assuming an idealized parallelization of the slaves -which is very optimistic considering the NEMOs currently run Euphemia on 8 threads [47] -, the run time per iterate would be lower-bounded by the time for solving the master programs (master plus projection programs for the Level Method, as reported between brackets in the tables). As an example, the "most difficult" BE-instance was solved in 266 iterates by D-W, with 2.3 sec/iter for solving the master program. This implies a lower bound of more than 10 minutes for obtaining the CHP. On the same instance, the Level Method required 37 iterates, with 1.4 sec/iter for solving the masters, which amounts to a total of less than 1 minute. Furthermore, whereas the price dimension of our test cases has been selected so as to be comparable to the EU market, the number of generators (or market bids) is well below the value that occurs in practice. As an order of magnitude, Euphemia currently solves instances with around 160,000 hourly orders (convex) and 4,000 block orders (non-convex) [4]. This suggests that the time for solving the master programs would likely be higher on the real instances of Euphemia.

V. CONCLUSION
Our paper proposes a (known) bundle stabilization approach for efficiently solving convex hull pricing. We demonstrate that the Level Method is able to converge within few iterations to a certain target gap, while exhibiting a stable behaviour, on large instances which, in terms of price space dimension, are comparable to the size of the EU day-ahead auction.
It is likely that the choice of the best algorithm for solving CHP will depend on the specific use-case: the dimension of the network, the time horizon, the complexity of the unit commitment/market orders, the run time that is afforded to the algorithm, etc. Although no method can conceivably provide an ultimate solution for computing CHP in an arbitrarily complex setting, the Level Method indicates the promising behaviour of a family of "bundle approaches". This suggests areas of future research on alternative bundle approaches, such as the Proximal Stabilization method, the Doubly-Stabilized Bundle Method [40] or the Boxstep method [48], which appear to be well suited for solving the CHP Lagrangian relaxation.
Another question for future research relates to how the proposed approach can be adapted in case one of the following assumptions is relaxed: the convexity of the grid model and the separability of the generators profit maximization problems.
Having scalable algorithms capable to compute CHP on large instances also enables more extensive quantitative analysis of its economical behaviour. As far as the EU market is concerned, we are interested in (i) expanding tests on realistic instances of Euphemia -our preliminary tests show that the Level Method can solve the 4MMC run of Euphemia [13] in ∼1 minute -, (ii) examining the effects of non-uniform pricing on enhancing welfare in the EU day-ahead market, and (iii) understanding distributional effects of non-uniform pricing as well as gaming effects.