Dynamic locational marginal emissions via implicit differentiation

TL;DR

This paper introduces a model-agnostic method using implicit differentiation to compute dynamic locational marginal emissions in electricity systems, improving accuracy by accounting for temporal constraints.

Contribution

The authors develop a novel, model-agnostic approach for calculating dynamic LMEs using implicit differentiation, applicable to complex dispatch models with temporal constraints.

Findings

Incorporating dynamic constraints improves model accuracy by 8.2%.

Static and dynamic LMEs differ with a normalized RMS deviation of 28.40%.

The method is validated against state-of-the-art approaches using real U.S. data.

Abstract

Locational marginal emissions rates (LMEs) estimate the rate of change in emissions due to a small change in demand in a transmission network, and are an important metric for assessing the impact of various energy policies or interventions. In this work, we develop a new method for computing the LMEs of an electricity system via implicit differentiation. The method is model agnostic; it can compute LMEs for any convex optimization-based dispatch model, including some of the complex dispatch models employed by system operators in real electricity systems. In particular, this method lets us derive LMEs for dynamic dispatch models, i.e., models with temporal constraints such as ramping and storage. Using real data from the U.S. electricity system, we validate the proposed method against a state-of-the-art merit-order-based method and show that incorporating dynamic constraints improves model accuracy by 8.2%. Finally, we use simulations on a realistic 240-bus model of WECC to demonstrate the flexibility of the tool and the importance of incorporating dynamic constraints. Namely, static LMEs and dynamic LMEs exhibit a normalized average RMS deviation of 28.40%, implying dynamic constraints are essential to accurately modeling emissions rates.

Equations

$$\begin{array}[]{lll}\textrm{minimize}&\sum_{j=1}^{k}f_{j}(g_{j})\\[5.0pt] \textrm{subject to}&g_{j}\in\mathcal{G}_{j},&j\in[1:k],\\ &\mathbf{1}^{T}d_{t}=\mathbf{1}^{T}\tilde{g}_{t},&t\in[1:T],\\ &F(d_{t}-B\tilde{g}_{t})\leq u^{\max},&t\in[1:T],\end{array}$$

$$\mathcal{G}=\{g\mid$$

$$g_t-\rho \le g_{t+1}\le g_t+\rho ,t\in [1:T-1]\}$$

$$\displaystyle\begin{array}[]{rll}0&\leq s\leq C,\quad\quad 0\leq\gamma\leq P,\quad\quad 0\leq\delta\leq P,\\ g_{t}&=\delta_{t}-\gamma_{t},\quad\quad s_{t}=s_{t-1}+\eta\gamma_{t}-(1/\eta)\delta_{t},\end{array}$$

$$J{x}^{\ast }(D_0)=-J_{x}K(D_0,x_0{)}^{-1}{J}_{D}K(D_0,x_0).$$

$$\begin{array}[]{lll}\textrm{minimize}&\sum_{j=1}^{k_{1}}f_{j}(g_{j})\\[5.0pt] \textrm{subject to}&g_{j}\in\mathcal{G}_{j},&j\in[1:k_{1}],\\ &g_{j}=g_{j}^{*},&j\in[k_{1}+1:k],\\ &\mathbf{1}^{T}d_{t}=\mathbf{1}^{T}\tilde{g}_{t},&t\in[1:T],\\ &F(d_{t}-B\tilde{g}_{t})\leq u^{\max},&t\in[1:T],\end{array}$$

Full text

1 Introduction

Policy-makers interested in decarbonizing the electricity grid require reliable emissions data in order to quantify the impact of a particular policy or intervention strategy. Similarly, grid operators conducting generation and transmission expansion studies [1, 2, 3, 4] are increasingly looking to reduce emissions [5, 6, 7, 8, 9] and require detailed information on the relationship between demand and emissions. The need for this data will only grow as various systems, such as transportation and heating, begin to electrify and place additional strain on the grid. Unfortunately, electricity systems depend on complex, interconnected transmission network structures with numerous operational constraints, making it difficult to attribute emissions to energy consumption at a particular time and place [10].

Emissions rates are important metrics that quantify the amount of pollutant emissions, e.g., $\text{CO}_{2}$, $\text{SO}_{2}$, or $\text{NO}_{\text{x}}$, attributable to the consumption of energy. Researchers and decision-makers often examine average emissions rates [11, 12, 13, 14, 15] which measure the average rate of emissions per MWh of energy consumed, and marginal emissions rates [16, 17, 18, 19, 20] (also known as marginal emission factors or marginal emission intensities), which measure the rate at which emissions increase or decrease given a marginal change in energy demand. While average emissions rates quantify emissions over long periods, marginal emissions rates better quantify the emissions impacts of small, local changes in demand, since only a few specific generators are expected to change production in response. This response to marginal changes can be estimated both at the network level—quantifying the aggregate sensitivity across many nodes in the network—or at specific locations. Indeed, these metrics vary on a node-by-node basis, as network constraints and the local energy mix dictate which generators are available to compensate for changes in demand. Hence, locational marginal emissions rates (LMEs) [21, 22] quantify the emission sensitivity at the nodal level, revealing the spatial heterogeneity in marginal emissions rates that emerges from network constraints. LMEs are the emissions-equivalent of locational marginal prices, which have been studied extensively in the power systems community due to their importance to electricity markets [23, 24, 25]. LMEs have been used to quantify the impacts of various policies on carbon emissions, e.g., increasing electric vehicle penetration [26, 27] and changing electric vehicle charging policies [28], and are published live in the PJM Interconnection [29], a major U.S. system operator. LMEs have also been used in transmission expansion studies [30, 31, 9]. In this application, the LMEs define the marginal emissions effect of offsetting demand at a particular node and can be viewed as the gradient of emissions with respect to net load in the planning optimization problem [9].

Empirical studies on marginal emissions rates in the U.S. and U.K. have used regression-based approaches to estimate emissions rates across large geographical regions [16, 32, 33, 21]. These works leverage publicly available data and fit linear models to changes in emissions and net demand. The main benefit of these methods is that they do not require a model of the underlying electricity system. However, because of their inherent data requirements, these methods are difficult to extend to finer geographic resolutions and hypothetical electricity systems that lack preexisting data.

In contrast, model-based approaches explicitly calculate LMEs using the underlying dispatch model. This calculation has been performed using marginal generator identification in merit-order models [34, 35], LMP-based approximations [36, 37, 38], or, in simple cases, explicit formulae [39, 40, 22]. Model-based methods are promising because dispatch in real-world electricity systems often involves solving an optimization problem. However, the models used in real world systems are highly complex, limiting the applicability of specific derivations in the aforementioned studies and highlighting the need for incorporating and feature specific constraints, e.g., ramping limits or storage, that need to be taken into account when calculating LMEs [35, 41].

For example, merit-order-based models usually neglect the impact of network constraints and only focus on generation cost to identify the marginal generator. LMP-based methods, which rely on matching LMPs with generation cost in order to identify the marginal generator, are not exact [37, 38], and the presence of complex coupling constraints would likely make identification even harder. On the other hand, analytical derivations, while exact, have so far only been conducted on dispatch models with a limited number of constraint types, e.g., transmission constraints. As far as we are aware of, no previous studies have exactly accounted for time dependencies such as energy storage when computing LMEs.

In this work, we address this limitation by developing a method that supports arbitrary convex optimization-based dispatch models.

Contribution and outline

This paper makes three main contributions:

• •

We propose a new method to compute LMEs in arbitrary convex optimization-based dispatch models. This is in contrast to standard regression-based methods that may have significant data requirements and previous model-based methods that have been derived for specific dispatch models. The method we propose generalizes previous analytical derivations [40, 22] and is generic, flexible, and suitable for real systems dispatched by grid operators.

• •

We use the proposed method to derive LMEs in networks with dynamic constraints, such as energy storage. As far as we are aware of, this is the first method to calculate LMEs for dynamic network dispatch models, such as the standard dynamic economic dispatch problem.

• •

We demonstrate the utility of computing LMEs in dynamic models using two different experimental setups. First, using a published dataset [34], we show that dynamic models more accurately represent the real-world relationship between demand and emissions compared to their static counterparts. Second, we use our method to study the impact of dynamic constraints on emissions rates using a realistic model of the Western United States transmission system. The dynamic LMEs are distinct from their static counterparts, demonstrating the importance of accurately including all relevant dynamic constraints when estimating emissions sensitivities.

The paper is structured as follows. In Section 2, we introduce the problem of computing LMEs in dynamic electricity networks. We show that this problem generalizes previous approaches [40, 22] to complex models with temporal constraints. We then show how to solve this problem using implicit differentiation in Section 3. Although we use this technique to compute LMEs for the model specified in Section 2, our technique generalizes to arbitrary convex optimization-based dispatch models, including those used by system operators in real world electricity markets. Lastly, we report simulation results on two datasets in Section 4. In the first experiment, we demonstrate the validity of our approach on real US electricity data and compare our results with an established method [34]. In particular, we show that a dynamic model with unit-commitment constraints more accurately models changes in emissions compared to its static counterpart. Second, we analyze a 240-bus model of the Western United States and show that, in the presence of grid-scale storage, computing LMEs dynamically is essential to accurately quantifying changes in emissions. We conclude and discuss future work in Section 5.

2 Problem formulation

In this section, we formulate the problem of computing the LMEs in a dynamically-constrained electricity system. First, Section 2.1 provides background information on the dynamic dispatch problem, a mathematical model for electricity networks with temporal constraints. We then describe our mathematical model for emissions and marginal emissions in static networks in Section 2.2, where we formally state the dynamic marginal emissions problem. Next, in Section 2.3, we describe two special cases of the dynamic marginal emissions problem that have been solved in previous work. Notably, both special cases are static, i.e., they do not incorporate any dynamic constraints. Finally, Section 2.4 gives three examples of dynamic devices, i.e., devices that cannot be represented in a static model.

2.1 Dynamic dispatch model

In electricity systems, a dispatch model is a mathematical model for deciding which electricity generators should be used to meet demand. Dispatch models are often formulated as convex optimization problems, where the variables are the amount of power produced by each generator, and the parameters include both the current electricity demand and the physical constraints on the system. When modeling emissions, past work has often considered static dispatch models, i.e., models that only reflect a single instant in time [40, 22, 15, 34]. However, most real world electricity systems have dynamic constraints—constraints that couple generator outputs across time periods. For example, a generator with ramping constraints can only change its power output by a small amount between successive time periods. In order to effectively model the impact of temporal constraints on emissions, we will study the dynamic optimal power flow problem,111The dynamic (DC) optimal power flow problem has been well studied in the power systems community [42, 43] and is sometimes referred to as the dynamic economic dispatch problem.

[TABLE]

where the variable is $G\in\mathbf{R}^{T\times k}$. In the above, we use $g_{j}$ to refer to $j$-th column of $G$, and $\tilde{g}_{t}$ to refer to the $t$-th row of $G$. The matrix $G$ represents the power output of $k$ devices over $T$ timesteps; each device can supply power, store power, or otherwise interact with the grid. The entry $G_{tj}$ represent the output of device $j$ at time $t$.

Device Costs and Constraints

Each device $j$ has three properties: a convex cost function $f_{j}(g):\mathbf{R}^{T}\rightarrow\mathbf{R}$, a convex constraint set $\mathcal{G}_{j}\subset\mathbf{R}^{T}$, and a location on the network. We model the device locations with a matrix $B\in\{0,1\}^{n\times k}$ that maps each device to a particular node in the network, i.e., $B_{ij}=1$ if device $j$ is located at node $i$, and $B_{ij}=0$ otherwise. The objective of Problem (1) is to minimize the sum of the device costs, and the first constraint states that each device must stay within their constraint set.

Network Constraints

Problem (1) considers an electricity network with $n$ nodes and $m$ edges (transmission lines). Node $i\in[1:n]$ has demand $(d_{t})_{i}$ at time $t\in[1:T]$. The second constraint is that power must be balanced across the network, i.e., $\mathbf{1}^{T}d_{t}=\mathbf{1}^{T}\tilde{g}_{t}$. Finally, the third constraint is that the power flowing across each transmission line is limited by its capacity. We define $F\in\mathbf{R}^{n\times m}$ to be the power flow distribution factor matrix, where $F_{i\ell}$ determines how a power injection at node $i$ (withdrawn at node $n$) affects power flow across line $\ell\in[1:m]$. Because of thermal and voltage phase angle constraints, each line $\ell$ can only transport up to $u^{\max}_{\ell}$ units of power, modeled with the constraints $F(d_{t}-B\tilde{g}_{t})\leq u^{\max}$.

Solution Map

Let $D=(d_{1},\ldots,d_{T})\in\mathbf{R}^{Tn}$ be the concatenated vector of demand schedules and assume the solution to Problem (1) exists and is unique for all $D\in\mathbf{R}^{Tn}$.222We make this assumption with little loss of generality. For example, we can guarantee a solution exists by adding a generator with no capacity limits and extremely high costs to each node (representing curtailed demand). Similarly, we can guarantee uniqueness by adding a small quadratic penalty $\sum_{t,j}G_{tj}^{2}s$ to the objective. Let $G^{*}(D):\mathbf{R}^{Tn}\rightarrow\mathbf{R}^{T\times k}$ denote the optimal choice of $G$ in Problem (1) as a function of $D$. Because we assume the solution to Problem (1) exists uniquely for all $D$, then $G^{*}$ is a well-defined function. We call $G^{*}$ the solution map, and use the vector-valued function $\tilde{g}_{t}^{*}(D):\mathbf{R}^{Tn}\rightarrow\mathbf{R}^{k}$ to denote the $t$-th row of $G^{*}$. As we will see shortly, the solution map will allow us to formalize the relationship between demand and emissions.

2.2 Locational marginal emissions

We model the emissions of generator $i$ as a linear function of power output with rate $c_{i}$, i.e., the total emissions at time $t$ are $c^{T}\tilde{g}_{t}$. Since the generator power outputs are determined by the dispatch model, the emissions at time $t$ generated as a function of demand is $E_{t}(D)=c^{T}\tilde{g}_{t}^{*}(D)$. The total emissions over the entire time horizon are then $E(D)=\sum_{t}E_{t}(D)$. Although we use a linear model throughout the remainder of the paper, it is straightforward to generalize all our results to nonlinear models. For example, if each generator has nonlinear emissions functions $\gamma_{i}(g):\mathbf{R}\rightarrow\mathbf{R}$, then the total emissions at time $t$ is $E_{t}(D)=\sum_{i}\gamma_{i}((\tilde{g}_{t})_{i})$.

Problem statement

The LMEs $\Lambda(D):\mathbf{R}^{Tn}\rightarrow\mathbf{R}^{Tn}$ are the marginal rate of change in total emissions given a marginal change in demand at a specific node and at a given time. In other words, the LMEs are the gradient of emissions with respect to demand, i.e., $\Lambda(D)=\nabla E(D)$. The function $\Lambda(D)$ is vector-valued, since changes in electricity consumption at different nodes and different times may have different impacts on emissions. As an illustration, we report a comparison between total emissions and LMEs for different values of demand at a given node in an arbitrary network (see Fig. 1). Locally, LMEs do indeed provide good approximations to the change in total emissions. It is however clear that those metrics are only locally valid and can sometimes be ill-defined, e.g. at points of non-differentiability of total emissions. The problem we study in this paper is how to compute $\Lambda(D)$ when the solution maps $\tilde{g}_{t}^{*}(D)$ are determined by the dynamic optimal power flow problem. As far as we are aware, no prior published results have shown how to compute LMEs for generic dynamic dispatch models.

2.3 Special Case: Static Generators

When we restrict the devices (the functions $f_{i}$ and sets $\mathcal{G}_{i}$) to be static generators, we recover previous analytical models [39, 40, 22]. The static generator device has constraint set $\mathcal{G}=\{g\mid g^{\min}\leq g\leq g^{\max}\}$, where $g^{\min},g^{\max}\in\mathbf{R}^{T}$ and cost function $f(g)=\sum_{t}ag_{t}^{2}+bg_{t}$. The static generator could represent a traditional dispatchable generator, in which case $g^{\max}=\alpha\mathbf{1}$, or a renewable resource, in which case the entries of $g^{\max}$ may vary. Most importantly, the static generator has no temporal constraints: $g_{t}$ is independent of the choice of $g_{t^{\prime}}$ when $t\neq t^{\prime}$. In a network with only static generator devices, the dynamic problem would simplify to $T$ static optimal power flow problems that can be solved independently. Moreover, if we remove the network constraints by setting $F=0$, we recover the model used in [34] to empirically estimate emissions rates.

2.4 Dynamic Devices

By addressing the dynamic optimal power flow problem in its full generality, our framework allows us to consider dynamic devices as well. These devices have temporal constraints, implying their actions at any given time depend on their actions at other times. We give three examples below.

Ramp-constrained generators

Ramping constraints limit the rate at which a generator can change its output. These generators are modeled with the constraint set,

[TABLE]

where $\rho\in\mathbf{R}^{k}$ is the ramping rate of each generator. Ramp-constrained generator devices have the same cost functions as static generator devices, $f(g)=\sum_{t}ag_{t}^{2}+bg_{t}$. Ramping constraints are particularly useful in dynamic dispatch models with short time intervals, e.g., 15 minutes, and slow dispatching generators, like nuclear and hydro.

Storage devices

Storage devices include pump-hydro resources, grid-scale batteries, and DER aggregators selling storage services to the grid. We define a storage device to have cost function $f(g)=0$ and constraint set $\mathcal{G}$ that is the set of $g\in\mathbf{R}^{T}$ such that there exists $s,\gamma,\delta\in\mathbf{R}^{T}$ satisfying,

[TABLE]

for $t\in[1:T]$, where $s_{0}=0$. The vector $C\in\mathbf{R}$ is the storage capacity, $P\in\mathbf{R}$ is the maximum (dis)charge rate, and $\eta\in(0,1]$ is the storage efficiency.

Unit-commitment-constrained generators

Unit-commitment constraints are integer constraints specifying whether or not a power generator is on. If generator $i$ is on, it must produce a minimum power $g^{\min}_{i}$ and stay on for a specified period of time. We model this by modifying the generator device constraint set to be the set $\mathcal{G}=\mathcal{G}_{\textrm{static}}\cup\{0\}$, where $\mathcal{G}_{\textrm{static}}$ is the equivalent static generator. Although the set $\mathcal{G}$ is not convex, it can be made mixed-integer convex by introducing an integer variable $z\in\{0,1\}$. Although the mixed-integer constraint makes (1) NP-hard, many good heuristics for solving mixed-integer convex programs are readily available through commercial solvers such as Gurobi [44].

3 Implicit differentiation-based LMEs

In previous model-based studies, locational marginal emissions rates are derived by first calculating how generation changes with demand, and then multiplying the change in generation by the emissions rate of each generator. This is a manifestation of the chain rule, which states that the LMEs are $\Lambda(D)=\nabla E(D)=\sum_{t=1}^{T}J\tilde{g}_{t}^{*}(D)^{T}c$, where $Jf(z)\in\mathbf{R}^{k\times n}$ denotes the Jacobian of the function $f:\mathbf{R}^{n}\rightarrow\mathbf{R}^{k}$ evaluated at $z\in\mathbf{R}^{n}$. Therefore, the main technical challenge when computing $\Lambda(D)$ lies in computing an analytical expression for the Jacobians $J\tilde{g}_{t}^{*}(D)$.

In previous studies that only consider static dispatch models, i.e., $T=1$, one only needs to derive a single expression for $J\tilde{g}_{1}^{*}(D)\in\mathbf{R}^{k\times n}$. In the general setting, the situation is much more complex—one must derive $T$ Jacobians $J\tilde{g}_{t}^{*}(D)$ of size $k\times Tn$. Although deriving an analytical expression might be possible, we take a simpler and more powerful approach in this paper: we use the implicit function theorem to compute the Jacobians $J_{D}\tilde{g}_{t}^{*}(D)$. Our approach essentially generalizes the analytical derivations of [40, 22] to arbitrary convex optimization-based dispatch models, producing identical in the simpler static setting.

Theorem 1** (Implicit Function Theorem, [45]).**

Suppose $K:\mathbf{R}^{n}\times\mathbf{R}^{r}\rightarrow\mathbf{R}^{k}$ is strictly differentiable at $(D_{0},x_{0})\in\mathbf{R}^{n}\times\mathbf{R}^{r}$ and $K(D_{0},x_{0})=0$. Moreover, suppose $J_{x}K(D_{0},x_{0})$ is nonsingular, where $J_{z}f(z,y)\in\mathbf{R}^{k\times n}$ denotes the partial Jacobian of the function $f:\mathbf{R}^{n}\times\mathbf{R}^{r}\rightarrow\mathbf{R}^{k}$ with respect to $z$ evaluated at $(z,y)$. Then the solution mapping $x^{*}(D)=\{x\in\mathbf{R}^{r}\mid K(D,x)=0\}$ is single-valued in a neighborhood around $(D_{0},x_{0})$ and strictly differentiable at $(D_{0},x_{0})$ with Jacobian

[TABLE]

The implicit function theorem states that if a differentiable system of equations $K(D,x)=0$ has a solution at $(D_{0},x_{0})$, and the corresponding partial Jacobian $J_{x}K(D_{0},x_{0})$ is non-singular, then the solution map $x^{*}(D)$ is a locally well-defined function with Jacobian given by Theorem 1.

In our setting, the solution map $G^{*}(D)$ is not the solution to a system of equations, but is rather the solution to a convex optimization problem. From convex analysis, we know that $G$ solves Problem 1 if and only if it solves the Karush-Kuhn-Tucker (KKT) equations $K(D,G)=0$ [46, 45]. Therefore, we can apply the implicit function theorem to the KKT equations to compute $JG^{*}(D)$. If we assume that the device objectives $f_{j}$ are twice differentiable and that the device constraints $\mathcal{G}_{j}$ can be parametrized via a set of twice differentiable inequalities, then the KKT equations are strictly differentiable, and $(D,G^{*}(D))$ satisfy the conditions of Theorem 1 (for more details on the implicit function theorem and its applications to optimization, we refer the reader to [45] and the references therein). By combining this with the chain rule, we can then derive marginal emissions rates for the dynamic optimal power flow problem specified by Problem (1). To summarize, we compute LMEs in two steps. First, we compute the gradient of emissions with respect to the device outputs, i.e. $dE/d\tilde{g}_{t}=c$ in the case of linear emissions functions. Then, we multiply this by the Jacobian $JG^{*}(D)$, which is computed using the implicit function theorem. The resulting metrics indicate the changes in emissions resulting from a marginal change in demand, as illustrated in Fig. 1.

Critically, this derivation works for any dynamic dispatch model that fits the form in Problem 1, i.e., regardless of the choice of device cost functions $f_{j}$ and constraint sets $\mathcal{G}_{j}$ (as long as they are convex and twice differentiable). When calculating the LMEs, different choices of device characteristics in Problem (1) only change the KKT operator $K(D,G)$ and its corresponding partial Jacobians $J_{G}K(D,G)$ and $J_{D}K(D,G)$; however, the general derivation using the implicit function theorem remains unchanged. In practice, these Jacobians can either be derived analytically or computed using an automatic differentiation library, such as [47], given a programmatic specification $f_{j}$ and $\mathcal{G}_{j}$.

Remark: Since the dynamic OPF problem includes the static OPF problem and the economic dispatch problem as special cases, this work generalizes the derivations in [40, 22] and [34]. However, our method is by no means constrained to the dynamic dispatch model used in this paper; implicit differentiation can be used to derive the marginal emissions rates for any convex-optimization based dispatch model. Importantly, this includes many of the dispatch models used by system operators in practice, a point we revisit in Section 5.

3.1 Complexity and Software Implementation

We implement our method in Julia [48] for all the aforementioned dispatch models and constraints. Our implementation is publicly available on GitHub in the package DynamicMarginalEmissions.jl. We use Convex.jl [49] to model the problem and solve it with external solvers, such as ECOS [50] or Gurobi [44]. For the large-scale network discussed in Sections 4.2 and 4.3 (i.e., $n=240$ nodes, $k=136$ generators, $m=448$ lines) and $T=1$ time periods, our software package solves the dispatch problem and computes the resulting LMEs in just under a second on a modern laptop with a 2.3 GHz quad-core processor. For the same problem with $T=24$ time periods, the same machine takes about two minutes to solve the problem and compute the LMEs.

Our software package offers a flexible interface and can be used to analyze networks with different physical characteristics (e.g., locations of generators, transmission line ratings) and constraints (e.g., ramping rates). After specifying the network parameters, one can compute the LMEs with a single function call. Because our implementation is open-source, reasonably fast, and easy to use, we believe it is of value to the broader community.

In general, we expect our method to scale well to realistic, large problems encountered by grid operators. Specifically, let $z=T\cdot\max(m,n,k)$. Solving the dispatch problem itself requires $O(z^{4})$ operations. Once the dispatch problem is solved, constructing and inverting the Jacobian to compute LMEs only requires $O(z^{3})$ operations, which is negligible compared to the complexity of solving the dispatch problem itself. Since most grid operators must solve dispatch problems at regular intervals, e.g., every 15 minutes to clear the real-time market, computing LMEs can be viewed as a post-processing step requiring little additional compute time.

4 Simulation Results

In this section, we illustrate the applicability and utility of the suggested approach using two different simulation setups. First, in Section 4.1, we compute the LMEs for a static model and a dynamic model with unit-commitment constraints using real demand and generator data from the U.S. Western Interconnection [34]. We compare each model’s LMEs to real changes in emissions and to estimates from a merit-order-based model [34]. Second, in Sections 4.2,  4.3 and 4.4, we illustrate the methodology on a recent reduced network model of the same region [51]. Using the original dataset, we highlight the geographic variance of marginal emissions across the network. Then, we investigate the potential impacts of hypothetically large renewable penetration of storage and renewable generation. We conclude by comparing the LMEs obtained from a static approximation to those of the dynamic model.

4.1 Economic Dispatch in the Western United States

In our first experiment, we analyze electricity data from the U.S. Western Interconnection system in 2016. The Western Interconnection dataset is compiled in [34] and contains weekly generator capacities, generator costs, and generator emission rates for large (above 25 MW capacity) fossil fuel generators in the Western Interconnection, as well as hourly total demand and total carbon emissions. Because no transmission data is available, we consider models without transmission constraints. The LMEs for a static model, a dynamic model with unit commitment constraints, and a state-of-the-art static merit-order method are compared to the real hourly rate of change in emissions.

Models

We analyze two models, which we compare to a baseline. First, we analyze the results of the simple economic dispatch model (1), with linear costs $f_{i}(g_{i})=b_{i}g_{i}$, where $b_{i}$ is the marginal operation cost of generator $i$. Second, we analyze a dynamic economic dispatch model with unit commitment constraints, over a time horizon of $T=24$. The unit-commitment constraints are only applied to coal generators, all of which are given a minimum output of $g^{\min}_{i}=0.4g^{\max}_{i}$. We benchmark our results against the reduced-order dispatch model (RODM) described in [34]. The core of the RODM is a merit order-based dispatch process: generators are dispatched in ascending order of cost. After dispatching generators via the merit-order, the marginal generator—the generator next in line to modify its output to meet an increase in demand—is identified to find the marginal emissions rate of the system. In [34], post-processing steps are applied to generate the marginal emission rates. Notably, when no post-processing is applied, the RODM is identical to the economic dispatch model in (1) with linear costs $f_{i}(g_{i})=b_{i}g_{i}$.

Results

After generating LMEs $\lambda_{t}$ for every hour $t=1,\ldots,T$ of the year, where $T=8760$, we compare the resulting LMEs to the actual hourly changes in emissions. Specifically, we compute the change in demand $\Delta d_{t}$ and change in emissions $\Delta E_{t}$ for every hour of the year. Each model’s estimated change in emissions is given by $\Delta\hat{E}_{t}=\lambda_{t}\Delta d_{t}$. In order to compare the models, we compute the absolute error $|\Delta\hat{E}_{t}-\Delta E_{t}|/Z$ of each model’s estimate against the true hourly change in emissions at each timepoint, where errors are normalized by the mean absolute change in emissions $Z=(1/T)\sum_{t=1}^{T}|\Delta E_{t}|$. We use absolute error, instead of square error, to minimize the effect of outliers.

A violin plot of absolute errors is displayed in Figure 2, Panel A. As expected, the economic dispatch model and the merit-order model from [34] perform similarly—the merit-order model only differs from the economic dispatch model in its post-processing. Notably, the unit-commitment model better models hourly changes in emissions than both the economic dispatch model and the merit-order model, reducing the mean absolute error by 8.2%. We attribute this to the fact that the unit-commitment model accurately represents dynamic constraints that appear in real-world dispatch processes, namely that coal generators cannot rapidly turn on and off again.

LMEs as a function of demand are also reported in Panel B of Figure 2. Historical LMEs are computed as $\lambda_{t}=\Delta E_{t}/(\Delta d_{t}+\epsilon)$, where $\epsilon=0.5\ \text{MWh}$ is a small value used to avoid unreasonably large LMEs when $\Delta d_{t}$ is small. Following a similar procedure to [34], the LMEs for the data and for each model are smoothed using the mean of a rolling window of 20% of the data. Shaded regions representing the interquartile range (IQR) of the data are also plotted to better understand the variance of each model. After averaging, the LMEs produced by the economic dispatch model most closely resemble the data. However, the variation is significantly reduced in the unit-commitment model, and the IQR most closely resembles that of the data compared to both other models.

4.2 240-bus Network Model

In this experiment, we study LMEs using a recent 240-bus network that models the Western United States in 2018 [51]. The dataset includes generator capacities, costs, and ramping rates; hourly demand and renewable generation for $T=8760$ hours; and the efficiencies and capacities of four pumped hydro storage units. We solve the dispatch problem and compute hourly LMEs using a dynamic model with storage devices as described in Section 2. We use this experiment to demonstrate the impact of network constraints on the LMEs. Since the network has relatively little storage and only a few generators with ramping constraints, we do not analyze the impact of storage and dynamic constraints in this experiment.

We report the distribution of daily total emissions in Figure 3, Panel A (left frame) and the distribution of LMEs at 6pm in August (left frame) in Panel B. We observe that, on average, the distribution of nodal LMEs is narrowly concentrated around its mode. However, we also note that the distribution of LMEs has relatively long tails; although most of the LMEs are close to the mode, a few have drastically different values. In Panel C, we display a map of LMEs at 6pm in August (averaged over days of the month), illustrating the large geographic variation in emissions rates. Panel C demonstrates how transmission constraints can create large discrepancies in LMEs even within small geographic regions. For example, despite their geographic proximity, the California Bay Area differs significantly from the Sacramento area. The local diversity of the LMEs emphasizes the importance of modeling the network when computing emissions rates: a small geographic change can cause large differences in emissions rates.

4.3 High Renewable Scenario in the 240-bus Network

To illustrate the impact of grid-level changes on emissions, a high-renewable version of the 240-bus network is presented in this section. Specifically, we uniformly scale renewable generation in the original 2018 model so that renewable generators meet 27% of total demand (compared to 13.5% originally). We also add 15 GW of 4-hour battery storage distributed among the ten largest renewable nodes proportional to generation capacity. These batteries have a round trip efficiency of $89.8\%$ with symmetric charging and discharging efficiencies and are constrained to start and end each day with $50\%$ of their total capacity. As in Section 2, we assume the grid operator manages these batteries to minimize the overall cost of electricity. The right frame of Panels A and B in Figure 3 show the distribution of total emissions and LMEs in an identical manner to the 2018 case. Similarly, Panel D displays a map of average nodal LMEs akin to Panel C. The 2018 and the high renewable scenarios differ in several ways. First, as expected, total emissions decrease significantly in the high renewable scenario. The changes to the locational marginal emissions rates, on the other hand, vary significantly from node to node. For example, adding renewable generation and storage causes LMEs to decrease at nodes in southern California, but to increase at nodes in Oregon and Washington. In general, the distribution of nodal LMEs exhibits high variance and displays two modes, in contrast with the 2018 case.333We observe these differences for most hours of the day, but only display results for 6pm for concision. Overall, the changes in LMEs are complex and unintuitive: because of the presence of storage, nodal LMEs depend not only on non-local generation and transmission constraints, but also on that of every other time period. We believe this complexity is one reason grid operators should use an exact method for calculating LMEs (instead of relying, for example, on network-wide heuristics).

4.4 Comparison Between Static and Dynamic LMEs

In order to demonstrate the value of explicitly integrating dynamic constraints, we compare the true dynamic LMEs to the analogous “static LMEs” that arise from eliminating dynamic constraints. Specifically, we consider how the LMEs would differ between static and dynamic models with the exact same loads and device outputs. Since static models cannot incorporate dynamic devices, we first solve the dynamic dispatch problem, then fix the battery charging and discharging schedules and consider them as parameters to a series of independent static problems, eliminating any dynamic constraints between subsequent timesteps. We then compute the LMEs of the resulting model, which now only has static devices and constraints.

More formally, consider the dynamic optimal power flow problem in (1), with the devices ordered so that the first $k_{1}$ devices are static and the remaining $k_{2}$ devices are dynamic. After solving the dynamic problem with $k=k_{1}+k_{2}$ devices to obtain device outputs $G^{*}$ and LMEs $\Lambda^{*}$, we solve the static problem,

[TABLE]

where the variable is again $G\in\mathbf{R}^{T\times k}$. Since the schedules of the $k_{2}$ dynamic devices are fixed, $\tilde{g}_{t}$ is independent of $\tilde{g}_{t^{\prime}}$ for $t\neq t^{\prime}$, and Problem (2) can be decomposed into $T$ independent optimization problems, if desired. We compute the resulting ‘static’ LMEs $\Lambda^{\textrm{static}}$ from solving Problem (2), and compare them to $\Lambda^{*}$. In theory, the difference between the LMEs of a dynamic model and its static approximation can be arbitrarily large, as seen in the following example.

Example: Consider a single node network with $T=2$ timesteps and $k=3$ devices. The first device is a gas generator with constraint set $\mathcal{G}_{1}=[0,10]\times[0,10]$ (i.e., the generator has capacity 10 in both time periods), cost $f_{1}(g)=1$, and emissions rate $c_{1}=500$. The second device is a solar panel with constraint set $\mathcal{G}_{2}=[0,10]\times\{0\}$ (i.e., the generator has capacity 10 in the first period, and no capacity in the second period), cost $f_{3}(g)=0.1$, and emissions rate $c_{3}=0$. Finally, the third device is a battery with constraint set $\mathcal{G}_{3}$ specified by Section 2.4, with capacity $C=10$, charging rate $P=10$, and efficiency $\eta=1$. Assume a constant demand schedule $d=(1,1)$. The economic dispatch will result in the following device outputs: $g_{1}=(0,0)$, $g_{2}=(2,0)$, and $g_{3}=(-1,1)$, i.e., the solar panel will produce two units of power, storing some of it in the battery to serve the second time period, and curtail the remaining 8 units. The dynamic LMEs are naturally $\Lambda(D)=(0,0)$: if we had slightly higher demand in either period, the solar panel would curtail its output less to meet demand. Now, if we fix the battery charging schedule to $g_{3}=(-1,1)$, solving the static problem gives the same device outputs $g_{1}=(0,0)$ and $g_{2}=(2,0)$. However, the resulting static LMEs are $\Lambda(D)=(0,500)$, a drastically different result. This is because the static approximation fixes the battery schedule, so changes in demand during period two are met by a change in the gas generator output.

The toy example above demonstrates that dynamic and static LMEs can differ significantly in theory. We verify that this occurs in practice using the 240-bus network from Section 4.4, where we use the same procedure to compute dynamic LMEs and their static approximations.

We report these differences in Figure 4, Panel A, where we display the distribution across all nodes and all days of the year of the root mean squared (RMS) deviation between the vector of daily emissions rates for the static model, $\lambda_{\textrm{static}}\in\mathbf{R}^{24}$, and the dynamic model $\lambda_{\textrm{dynamic}}\in\mathbf{R}^{24}$. The average RMS deviation (normalized by the median LME) is 28.40%, indicating that the static and dynamic models yield significantly different results. In Panel B, we illustrate the static and dynamic LMEs for three randomly sampled days. While static LMEs are very good approximations in some instances (e.g., morning hours in top-left), they deviate significantly in others. These results suggest that ignoring dynamic constraints and simply computing static LMEs is not sufficient to model emissions in dynamic networks: explicitly computing dynamic LMEs is essential to understanding emissions rates in systems with significant dynamic constraints, such as large grid storage capacity.

5 Conclusion

In this paper, we introduce a novel method for computing locational marginal emissions rates using implicit differentiation. We use this method to compute the LMEs of dynamic dispatch models, i.e., dispatch problems containing temporal constraints.

Using real WECC electricity and emissions data, we find that incorporating these dynamic constraints improves model accuracy by 8.2%. Finally, we observe that dynamic LMEs are difficult to approximate with their static counterparts: in a synthetic approximation of the WECC network, static and dynamic LMEs have a normalized average RMS deviation of 28.40%. Since flexible loads and energy storage are expected to play a large role in future grids, we believe incorporating dynamic constraints will be essential to accurately modeling LMEs.

The method presented in this paper generalizes previous methods to arbitrary convex optimization-based dispatch models. Since many system operators use convex optimization-based dispatch models in day-ahead and real-time electricity markets [52], they could use this method to publish marginal emissions factors in real time. Although these models can be notably more complex than those analyzed in academic research, the proposed method can compute marginal emissions factors for any such model, as long as they can be represented as convex optimization programs. Moreover, by leveraging automatic differentiation software and optimization modeling languages [53], the system operator would only need to specify the objective and constraints of their dispatch problem. LMEs could then be published alongside LMPs to provide real time emissions information to electricity market participants. This could be helpful, for example, to a large internet company choosing to reduce emissions by directing internet traffic to servers in low emitting regions, a problem considered in [54], or more generally to operators wanting to define optimal load management strategies [38].

Finally, we comment on three directions for future work. First, our experimental results indicate that LMEs in dynamic models often display complex behaviors and are difficult to interpret due to temporal and network constraints. Deciphering the mechanisms underlying the structure of the LMEs in different settings would be useful in communicating these results and translating them into grid planning or policy decisions. Second, we note that computing LMEs in large networks could be computationally intensive. Exploiting network structure and using distributed computation could yield significant performance gains. Third, our paper shows how to compute LMEs when the full network model is available. In some cases, however, the network model may be unavailable to the interested party. Understanding how to estimate the parameters of the electricity network from publicly available data (using the methods developed in [55], for example) and then deriving marginal emissions factors from the learned model is an interesting area of research.

Acknowledgements

The authors thank Liang Min and Inês Azevedo for their valuable comments and suggestions.

Disclaimer: This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.