Data-Driven Inverse Optimization for Marginal Offer Price Recovery in Electricity Markets

An inverse optimization (IO) problem seeks to recover unknown parameters in the objective function [27] or constraints [13] of a forward optimization (FO) problem. The FO could be a linear program [3], a conic program [25], a mixed-integer program [2], a linearly constrained separable convex program [51], or more generally, a variational inequality function [8, 48].

The classic IO model assumes that there exists a single set of parameters that make the given solutions optimal. However, in reality, noise in observations can prevent the existence of parameters that make all solutions strictly optimal. This noise can come from various sources, such as (i) measurement errors during the data collection process, (ii) mismatches between the forward model and the actual underlying decision-making process, and (iii) bounded rationality (i.e., when the decision maker settles for sub-optimal results due to cognitive or computational limitations) [5, 32]. Applying the classic IO in noisy settings can lead to infeasibility or uninformative solutions, such as a trivial zero cost vector [14]. Therefore, the goal of IO in noisy settings is to minimize the fitting errors between the model and data to make the given solutions “approximately” optimal [5, 14, 32, 48].

In deterministic settings, IO can precisely recover the parameters of the FO problem based on a single observation. However, in noisy settings, multiple observations are required to minimize the empirical loss of the IO problem. This type of IO, which can handle multiple observations, is also known as data-driven IO [30, 32]. Solving data-driven IO problems can become difficult with an increase in variables as more training data is incorporated. To tackle this challenge, researchers have proposed to use machine learning (ML) methods instead of conventional data-driven IO models [7, 9, 17, 47].

The objective of IO is similar to those of some ML techniques, especially neural networks, as they both aim to estimate the unknown parameters of a model based on available observations. However, the key difference is that the underlying model in IO corresponds to the actual FO problem and its parameters can be interpreted, while in ML, the unknown parameters (such as weights and biases in neural networks) often lack a straightforward explanation. To incorporate physical laws into ML models, researchers have developed physics-informed ML methods [26]. Essentially, if the FO problem with unknown parameters is transformed into a differentiable layer in the neural network, these unknown parameters can be learned through training (similar to other parameters in the network) [1, 4, 12, 19]. For instance, OptNet is a network architecture that integrates a quadratic optimization program as a layer within an end-to-end neural network [4]. In [9], OptNet is used to determine the unknown parameters in a demand response model. The estimation of unknown parameters can also be achieved through other ML methods. For example, [7, 17] reveal the parameters through an online learning process, while [47] presents a deep inverse optimization model that formulates the parameter recovery problem as a deep learning model. The fundamental assumption behind these ML models is that the FO model is differentiable with respect to the unknown parameters, which allows for the computation of gradients and the training based on (stochastic) gradient descent methods.

To compare IO and non-physics-informed ML, [24] evaluated their performance in imputing a convex objective function and found that IO is more suitable for problems with smaller training sets or highly complex objectives/constraints, while ML can handle larger training sets and is less sensitive to noisy training data. However, potential drawbacks of ML, such as its tendency to diverge in extreme situations, restrict its use in critical decision-making. Hence, to gain the trust of decision makers in using ML, formal performance guarantees need to be established for ML models [6].

IO has various applications across fields, including in power systems and power markets. To maximize the profit of market participants in the day-ahead bidding process, IO has been used to recover the piece-wise linear cost function of generators based on historical market-clearing results in [15, 16, 41]. Similarly, [40] also designs a IO problem based on market-clearing results, but the forward problem is an equilibrium model in an oligopolistic electricity market, aimed at determining if the market structure matches the observed outcomes. While [15, 16, 40, 41] assume that all parameters except for offer prices are known, [10] uses IO to recover market structures that are not disclosed to market participants, such as parameters of transmission lines. Alternatively, IO has been used to estimate the parameters of demand response in power systems in studies such as [28, 42, 43]. In [43], IO is recast as a bi-level program and solved heuristically. [42] deals with a non-convex IO model and develops a two-step heuristic approach to solve it. In [28], IO is transformed into a quadratically constrained quadratic program and solved through successive linear programming.

The aim of this paper is to develop a data-driven IO method for recovering the marginal offer prices of generators from historical market-clearing results. The US power systems and markets have the following unique properties that distinguish this IO problem from other application scenarios.

First, the structures and most parameters of power markets are transparent [35], and the market-clearing results, including schedules and prices, are promptly and accurately published [33, 34]. This allows for the observation of both primal variables (schedules) and some dual variables (prices) in the forward problem, giving our price recovery problem an advantage over other IO problems where only primal solutions can be observed.

Second, the market operation relies on the use of locational marginal prices (LMPs) which represent the incremental costs of providing the next megawatt of energy at different locations within the power system. However, the offer price of a generator is only reflected in the LMPs when it is the marginal generator at that location. Hence, a single observation is insufficient to recover the offer prices of all generators, requiring the use of data-driven IO to explore multiple scenarios where different generators are marginal.

Third, the offer price of a generator tend to be stable over time. In day-ahead markets, generators are paid according to LMPs instead of their offer prices. To secure market participation, they typically set the offer prices close to or slightly above their generation costs. However, the offer price of a generator may vary due to external factors like fuel prices and weather, causing the observations from the market being the optimal solutions to the same forward model but with different cost vectors.

Finally, the noise in market-clearing results is sparse and traceable. According to [46], only a small percentage (less than 0.1%) of real-time intervals have resulted in price correction since 2009, indicating that random errors in published price data are uncommon. Another source of noise is the mismatch between the forward model studied in this paper (i.e., the DC optimal power flow model) and the actual market-clearing model which takes into account more constraints and may involve ad-hoc manipulations by market operators. These noises may affect the accuracy of IO results.

In summary, the transparent market structure and abundant market data make the data-driven IO feasible. However, the noise in historical data and the fluctuating nature of offer prices pose challenges to the IO problem. We plan to tackle these challenges by solving the data-driven IO problem using the gradient descent method to effectively utilize the training data and ensure robustness in noisy settings. Our focus is also on establishing performance guarantees for the proposed approach.

After reviewing the related work on IO in Section 2 and formulating the FO problem in Section 3, our contributions in the rest of the paper are:

• We propose a data-driven IO model for marginal offer price recovery in wholesale power markets. Our method differs from previous studies such as [16, 41] since it incorporates the unknown cost vector directly into the objective function, rather than relying solely on the primal or dual variables of the FO problem. This design results in a computationally manageable single-level optimization problem, compared to the bi-level programs in [5, 32].
  
• Unlike previous work that solves the data-driven IO using off-the-shelf solvers [10, 15, 16, 40, 41], we solve it using gradient descent and derive the corresponding error bound. Our method allows for proving the convergence to the global optimum with a sub-linear rate, whereas the method in [9] only guarantees convergence in the absence of inequality constraints in the FO problem. Furthermore, we rigorously prove the existence and uniqueness of the global optimum to the data-driven IO problem, which are generally missing in other studies.
  
• Different from [10, 15, 16, 41] which accomplish the offer price recovery using IO but neglect the impact of noise in observations, our approach considers potential sources of errors in market-clearing results and evaluates robustness of the data-driven IO in two noisy settings. The validity and robustness of the proposed data-driven IO method are demonstrated through simulations in the IEEE 14-bus system and the 1814-bus NYISO transmission network.
  

While the proposed data-driven IO model focuses on the price recovery problem in day-ahead markets, it can also be adapted to applications in real-time and ancillary service markets by appropriately changing constraints (while preserving linearity) and training data.

The following notations are used in the paper: I_n is an identity matrix of size n × n with ones on the main diagonal and zeros elsewhere, J_{m, n} is an all-ones matrix of size m × n, O_{m, n} is an all-zeros matrix of size m × n, D(·) creates a diagonal matrix from a vector, and ‖ · ‖_p represents the ℓ_p-norm.

Let $x \in \mathbb {R}^n$, $w \in \mathbb {R}^n$, $A \in \mathbb {R}^{m \times n}$, and $b \in \mathbb {R}^m$. We define our linear forward optimization (FO) problem as:

where $\xi \in \mathbb {R}^m$ is the vector of dual variables associated with Ax ≥ b. Let $\mathcal {X}(w)$ be the set of feasible solutions to FO(w) and $\mathcal {X}^{OPT}(w)$ be the set of optimal solutions to FO(w). We assume that $\mathcal {X}(w)$ and $\mathcal {X}^{OPT}(w)$ are non-empty.

Let $x^0 \in \mathbb {R}^n$ denote an observed solution to FO(w) in (1). In a deterministic setting, there exists a cost vector w such that x⁰ is an optimal solution to FO(w). Therefore, given A, b and x⁰, the goal is to find w such that $x^0 \in \mathcal {X}^{opt}(w)$. Note that there may exist multiple w that can achieve this goal [27]. For example, if $x^0 \in \mathcal {X}^{opt}(w)$, then $x^0 \in \mathcal {X}^{opt}(G(w))$ also holds, where G(·) is any convex increasing function. The uniqueness of w can be guaranteed by normalization and regulation, e.g., by requiring ‖w‖ = 1. Further, if FO(w) is a linear program, its inverse problem IO(x⁰) is also a linear program [3].

One way to evaluate the optimality of the estimated w is via the optimality conditions of FO(w), which can be satisfied exactly by x⁰ when the optimal w is found. Therefore the loss function of FO(w) is trivial and can be set to 0. The optimality primal-dual (PD) conditions of FO(w) can be formulated as [3]:

where (2b) ensures a unique optimal solution, (2c) embodies primal feasibility, (2d) and (2e) denote dual feasibility, and (2f) represents strong duality. Note that (2c) can be omitted since it does not depend on any primal or dual variables.

Alternatively, the optimality conditions of FO(w) can also be formulated as the Karush-Kuhn-Tucker (KKT) conditions [27]:

where (3b) is complementary slackness. Note that (2d) is referred to as the stationarity constraint in the KKT constraints. The IO-PD and IO-KKT formulations in (2) and (3) are equivalent, i.e.,

In noisy settings, it is assume that the observed solution x⁰ is feasible but not necessarily optimal in FO(w), i.e., $x^0 \in \mathcal {X}(w)$. Thus, the optimality conditions in (2) or (3) may not be exactly met, and the objective of IO becomes finding the vector w such that x⁰ is an “ϵ -optimal” solution, i.e., the choice of w enables x⁰ to approximately satisfy the optimality conditions of FO(w). This IO problem is referred to as generalized IO (GIO) in [14] and there are two GIO models reported in the literature.

Based on the IO-PD formulation, GIO can be formulated by adding a slack variable ϵ_pd to relax the constraints of primal feasibility (2c) and strong duality (2f) [8, 14, 40], as given by:

It is proved in [14] that (4b) can be omitted since it is automatically satisfied by the optimal solution to GIO-PD(x⁰) without (4b).

Alternatively, GIO can also be formulated based on the IO-KKT formulation by relaxing the stationarity conditions in (2d) and complementary slackness in (3b) with slack variables ϵ_stat and ϵ_comp [27, 43], respectively:

If the observed solution is not feasible, i.e., $x^0 \notin \mathcal {X}(w)$, we can add slack variables ϵ_prim and ϵ_dual to the primal feasibility (2c) and the dual feasibility (2e), respectively, instead of forcing them to hold.

Both forms of GIO are typically straightforward to solve. For example, [14] derives the closed-form solutions to GIO-PD under different forms of loss functions, while [27] proves that GIO-KKT in (5) is a finite-dimensional convex optimization problem which can be solved efficiently with off-the-shelf solvers.

The study in [5] shows that GIO-PD and GIO-KKT are statistically inconsistent, meaning that their results may converge to incorrect values even with a large amount of training data. In response, [22] introduces the nominal IO (NIO), in which the objective is to minimize the difference between the observed solution x⁰ and the optimal solution of FO(w).

Similarly, the NIO based on KKT constraints can be formulated. The distinction between the GIO and NIO formulations lies in the choice of loss functions that measure the discrepancy between predictions and observations. In the GIO formulation, the loss is measured by the slackness required to make the observed solution satisfy the approximate optimality conditions under the estimated parameters. In the NIO formulation, the loss is the distance between the observed solution and the optimal solution to FO under the estimated parameters. Moreover, the NIO formulation is statistical consistent [5] and more robust to noise [48], but it is a bi-level program and thus NP-hard [5], while the GIO formulations are relatively easier to solve.

The IO formulations presented in (2)–(7) are based on a single observation x⁰, but they can be extended to the data-driven form to accommodate multiple observations. For instance, given N observed solutions to FO $\lbrace x^0_i\rbrace _{i \in \mathcal {I}}$ and the corresponding parameters $\lbrace A_i,b_i\rbrace _{i \in \mathcal {I}}$, the NIO in (7) can be adapted as follows:

In the US, the operation of power systems and power markets are managed by the designated independent system operators (ISOs). ISOs ensure the safe and reliable operation of power systems and regulate markets for various products such as capacity, energy, and ancillary services [44]. Fig. 1 illustrates the typical two-settlement process in energy markets, including the day-ahead market (DAM) and real-time market (RTM). The clearing of DAM is based on security-constrained unit commitment (SCUC) and security-constrained economic dispatch (SCED), which determines the least-cost generator schedules and hourly LMPs for the next day based on day-ahead generation offers, load bids, as well as the forecast of load and renewable energy source (RES) generation. The clearing of RTM is based on real-time commitment (RTC) and real-time dispatch (RTD). The commitment of units is predominantly determined in the DAM, making it a crucial aspect for power suppliers. However, the actual scheduling may diverge from the DAM schedule due to changes in operating conditions, influences of additional RT generation offers, and variations in actual load. This paper focuses on recovering the DA offer prices of generators from the DA market-clearing results.

The actual SCUC model employed by ISOs is complex and involves multiple solution steps, as described in [37], and it is simplified to a representative SCUC model in some papers [21]. Following similar simplifications, we formulate the following UC model for a system with n generators, m nodes and l lines:

where (9a) minimizes system operation cost using a linear cost function, (9b) ensures the power balance in the system, (9d) and (9c) limit the output of generators to their technical bounds x^min and x^max, (9f) and (9e) limit the power flow on each line to the maximum transmission capacity f^max. Other notations used in (9) are explained in Table 1.

Model (9) is a mixed-integer linear program (MILP) due to the presence of binary variables u. Although non-convex, it can be solved by modern solvers such as CPLEX or Gurobi. To obtain the dual variables for marginal price calculation, we need to convert (9) into an equivalent convex linear program (LP) by replacing binary variables to real variable $u\in \mathbb {R}^n$ and setting u = u^*, where u^* is the optimal commitment results obtained by solving (9) with a MILP solver. This approach follows the results in [36], where the authors show that the optimal solution of the MILP is equal to the optimal solution of the LP when u = u^* is enforced in LP. The resulting LP is the following DC optimal power flow (DCOPF) model:

Greek letters in parentheses in (10b)–(10f) denote Lagrangian multipliers (dual variables) of the respective constraints. Specifically, λ denotes the energy price at the reference node in the power system. The KKT optimality conditions of (10) are:

The Lagrangian function of (10) is:

The use of LMPs is the cornerstone of market operation [46]. The marginal energy price, represented by the dual variable λ in (10b), is the LMP at all nodes when the transmission capacity is unlimited. It reflects the cost of serving the next increment of load in the most efficient way. Based on (11a), we can express λ as:

We start from the single-observation case in a deterministic setting, where the observed solutions λ⁰, ω⁰ and x⁰ are assumed to be the optimal solutions of (10) and therefore satisfy the KKT conditions in (11). From (11a), we know:

Then, we extend (21) to the multi-observation case in a deterministic setting. Recall that the offer price of a generator may vary at different times due to external factors. In this case, while the observed solutions still satisfy the KKT conditions in (11), calculating c⁰ based on different sets of λ⁰ and ω⁰ would yield different offer prices for the same generator. To address this, we propose the following data-driven IO formulation. Given multiple observed solutions $\lbrace \lambda ^0_i,\omega ^0_i\rbrace _{i \in \mathcal {I}}$ where $\mathcal {I}$ is the set of training data indexed by i, the objective of data-driven IO is:

Note that the data-driven IO in (22) is a single-level problem, so it has lower computational complexity compared to the bi-level data-driven IO in (16) and (17). This advantage comes from the availability of a closed-form relationship between the unknown parameter c⁰ and the observed solutions λ⁰ and ω⁰ as specified in (21). Although this relationship is insightful in the deterministic setting, it does not hold in some noisy settings, which will be discussed in section 4.4.

Since only the offer prices of free generators can be recovered from (21), $c_i^0$ is a sparse vector. The zero elements in vector $c_i^0$ do not provide any additional information, so they should be eliminated from the loss function in (22) as:

We can use the gradient descent (GD) method to iteratively update the values of $\hat{c}$. The update rule for $\hat{c}_k$ is:

Algorithm 1 outlines the procure of solving the proposed data-driven IO model based on GD. To enhance computational performance, stochastic gradient descent (SGD) can be used when dealing with large training data sets or when the price recovery task becomes an online learning problem.

In this subsection, we prove that the loss function in (23) converges to its optimal value with a sub-linear rate using the gradient descent (GD) method. The established convergence theorem for GD on convex, Lipschitz continuous functions can be found in Appendix A.1. We first discuss the convexity and Lipschitz continuity of $\left\Vert x \right\Vert _p^p$, which are essential to prove the convergence of the proposed data-driven IO formulation.

The proof of Lemma 4.1 can be found in Appendix B.1. Based on Theorem A.1 and Lemma 4.1, we can propose the convergence theorem for the data-driven IO in Section 4.1:

Specifically, if ℓ₁-norm is used in l(c), i.e., when p = 1, the error bound in (25) becomes $\epsilon \le n \bar{c} /\sqrt {T}$. If the loss function is strongly convex, for example when ℓ₂-norm is used, the convergence rate could be greater [39]. To satisfy additional limits on the value of $\hat{c}_k$, inequality constraints can be introduced into the unconstrained optimization model in (23). The resulting problem can then be solved using projected gradient descent. However, for the convergence theorem to hold, it is crucial that the feasible set of $\hat{c}_k$, as defined by the inequality constraints, be both closed and convex [49].

This subsection focuses on demonstrating the existence and uniqueness of the global optimal solution to the data-driven IO formulation in (23). The established theorems for the existence and uniqueness of optimal solutions to convex problems can be found in Appendix A.2. First we prove the strict convexity of $\left\Vert x \right\Vert _p^p$, which is essential to prove the uniqueness of global optimum to the proposed data-driven IO formulation.

The proof of Lemma 4.3 can be found in Appendix B.2. Based on Theorems A.2– A.3 and Lemma 4.3, we can propose the following theorem for the data-driven IO formulation in (23):

Inspecting Theorems 4.2 and 4.4, we notice two differences on the assumptions of l(c): Theorem 4.2 assumes p ≥ 1 and $0\le |c_k|\le \bar{c}$ for ∀k, while Theorem 4.4 assumes p > 1 and $\underline{c} \le c_k\le \overline{c}$ for ∀k. Variable c_k represents the offer price of generator k, which is greater than its generation cost and capped due to market power regulations. Hence, the assumption of $\underline{c} \le c_k\le \overline{c}$ is realistic. We remark that in some exceptional cases, under heavy production credits, renewable generators can submit negative offer prices; however, this practice is expected to gradually phase out.

Theorems 4.2 and 4.4 only apply to generators that have been marginal at least once in the training data set. The offer prices of non-marginal generators are eliminated from the loss function l(c) in (23), so the estimations of these prices remain unchanged from their initial values, rather than converging to an optimal value.

We first evaluate the performance of the data-driven IO in the IEEE 14-bus system (details in Appendix C.1). We assume that each generator sets a piece-wise linear price baseline with five equally divided blocks based on the quadratic power generation cost, as shown in Fig. 2. The actual offer prices fluctuate around the baseline with deviations following a normal distribution N(μ, σ²), where μ = 0 and σ = 2 (unit: $/MWh) are assumed in this case.

5.1.1 Performance in Deterministic Settings. To recover as many offer prices as possible, we create 200 training data points and validate that each of the five generators is marginal at least once at each block. The initial values of the five blocks are set to 10, 20, 30, 40, and 50, respectively. Fig. 3 shows the variations of four estimated offer prices at each iteration and compares the convergence rates based on the ℓ₁ and ℓ₂-norms. In all four cases, the estimated prices converge to the true values, with faster convergence observed when using the ℓ₂-norm due to its strong convexity.

In this numerical experiment, we study the NYISO system with 1814 buses, 2207 lines, 362 generators and 33 wind farms (details in Appendix C.2). The structure and parameters of this NYISO system are mined and estimated from publicly available data sources [35]. The NYISO publishes daily LMPs at [34] and other market-clearing results every three months at [33]. We use this information to generate hourly market-clearing results over a three-month period (February, April, and August 2018), yielding a total of 2136 training data points. Note that in this case we assume that the baseline price of each generator is equally divided into ten price blocks, and the upper and lower bounds of each block are known. However, if the block setting rules are unpublished, we need to estimate the bounds of each block first based on historical data. Subsequently, we can recover the offer price of each block based on the estimated bounds of the blocks.

5.2.1 Performance in Deterministic Settings. As previously mentioned in Section 4.1, a training data point is valid for the price recovery of a generator if the generator is committed and marginal in that scenario. Fig. 4 summarizes the number of valid training data points for all the generators in the system. As shown in Fig. 4, most of the generators only have valid training data at a few steps of their offer prices, thus only a portion of the prices can be recovered. The results indicate that 44.03% of all offer prices can be recovered from the valid training data, while 81.21% of these recovered prices are based on less than 5 training data points. Meanwhile, for 85.35% of the generators, at least one block of price can be recovered. Note that there is an theoretical upper limit on the proportion of recoverable prices because an offer price can only be recovered from the LMPs if it affected the values of LMPs. Therefore, other approaches, such as NIO, also face the same limitation with regards to the recovery rate.

In total, 53 generators (out of 362) are never found to be marginal in the 2163 data points we analyzed, so their offer prices cannot be estimated from the training data. This result does not indicate a failure of the data-driven IO for large systems, as these non-marginal generators do not typically compete with others. For example, due to a relatively high generation cost, generator No.14 is never committed in this three months. Furthermore, incorporating more training data can decrease the number of unrecoverable generators. As more market-clearing data posted by market operators, the portion of recovered offer prices will increase, which is important for dealing with noisy data.

5.2.2 Performance with Model Mismatch. To assess the performance of the data-driven IO approach in the presence of model mismatches, we replicate the SCUC model used by NYISO in its day-ahead scheduling process, as outlined in [37]. This SCUC model, as detailed in [31], includes additional constraints such as ramping limits for generators, which may challenge the accuracy of recovered prices when those constraints are binding. Fig. 5 shows the binding frequency of two sets of constraints, including the power limits in (26d) and (26e) and the ramping limits in (27) and (28). It also compares the results across three months, namely February (with low demand), April (with medium demand), and August (with high demand). As depicted in Fig. 5 (b), the occurrence of only the ramping limits being binding is rare. As a result, the frequency of both sets of constraints being binding is approximately the same as the frequency of only the power limits being binding.

Fig. 5 conveys a crucial message that the ramping limits of free generators (whose power limits are not binding) are typically not binding either. As explained in Section 5.2.2, the ramping limits will only affect the recovered prices when they are binding. Hence, the impact of ramping limits is mostly eliminated from the valid training data set, which only includes the data of free generators.

With this model mismatch, we are able to recover 36.28% of all offer prices with an average relative error of 3.47%. Fig. 6 shows the errors in recovered prices at the seventh block where 171 generators are marginal. The errors are compared in the deterministic setting and two noisy settings, where noisy setting I includes additional Gaussian noise to 1% LMPs data with a mean value of 50 $/MWh and a standard deviation of 5 $/MWh, and noise setting II involves the model mismatch discussed above. The deterministic setting features small and random errors with an average close to zero, while in noisy setting I, sparse noise in the training data results in significant deviations of the recovered prices of some generators from their true values. Noisy setting II features higher errors compared to the other two cases, however, they are still within an acceptable level compared to the offer prices. Note that although the errors in Fig. 6 (c) tend to be mostly positive, it is not a universal conclusion since the values of errors in this noisy setting depend on the values of dual multipliers of binding constraints in the SCUC model, which can change over time.