Carbon-Aware Optimal Power Flow

$\mathcal{N}_{i}$

Set of neighbor nodes of node $i$.

$\mathcal{N}_{i}^{+}(\mathcal{N}_{i}^{-})$

Set of neighbor nodes that send power to (receive power from) node $i$.

$\mathcal{G}_{i}$

Set of generators at node $i$.

$\mathcal{L}_{i}$

Set of loads at node $i$.

$\mathcal{T}$

Set of time steps with the time interval $\delta_{t}$.

$w_{i,g}^{\mathrm{G}}$

Generation carbon emission factor of generator $g$ at node $i$.

There have been a number of recent studies [23, 24, 25, 26, 27, 28, 29, 30, 31, 32] that consider carbon emission flow in power system planning and operation. Reference [23] proposes a transmission expansion planning method that defines an index to quantify the equity performance of carbon emission allocation based on the carbon flow model. In [24], a multi-objective power network transition model is built to plan the retirement of aging coal-fired power plants, while one of the objectives is to minimize user-side carbon footprints. Reference [33] studies the low-carbon operation of multiple energy systems and derives the locational energy-carbon integrated price based on the nodal carbon intensities calculated using the carbon flow method. In [25, 26, 27, 28], carbon-aware expansion planning models are established for multi-energy systems under carbon emission constraints on electric devices and energy hubs. Additionally, the carbon flow model and constraints are taken into account in power scheduling [29], energy management [30, 31], and peer-to-peer carbon-electricity trading [32]. Existing studies outlined above focus on specific power system applications. However, it remains unaddressed in establishing a fundamental decision-making methodology necessary for guiding various grid decarbonization decisions and supporting theoretical studies and performance analysis.

Moreover, there are two critical issues in the integration of carbon emission flow into the grid decision models that necessitate to be addressed:

• 1)

  (Power Flow Directions): The carbon flow model [14] needs to pre-determine the power flow directions for all branches to identify the power inflows for each node. However, the directions of branch power flows are typically unknown prior to solving optimal decision models. Most existing works that employ the carbon flow model in grid decision-making either overlook the issue of unknown power flow directions or assume they are predetermined by alternating between grid optimization and carbon flow calculation through iterations [32]. Reference [28] introduces binary indicator variables to handle the unknown power flow directions in the carbon flow model. It results in a mixed-integer nonconvex quadratically constrained optimization problem, which cannot be readily solved using off-the-shelf optimizers due to involving both integer variables and nonconvex constraints. Thus, a tailored heuristic penalty-based iterative algorithm is designed to solve the optimization problem.

• 2)

  (Carbon Footprint Models for Energy Storage): By switching between charging and discharging, energy storage (ES) systems can shift load demand and transfer renewable energy across time, offering substantial potential to reduce power system emissions. Hence, developing carbon footprint models for ES systems is crucial, since it lays the quantitative foundation for making optimal ES planning and operation decisions, such as determining installation location and capacity, and managing charging and discharging. In addition, granular and accurate carbon accounting for ES systems enhances information transparency, supporting the development of low-carbon regulatory policies and market mechanisms for ES. It also incentivizes ES stakeholders to adopt low-carbon practices to reduce their own carbon footprint and overall grid emissions. Carbon accounting for ES systems has attracted considerable recent attention from the industry [34, 35]. References [29, 28, 36] propose different carbon footprint models for ES systems, while they neglect the carbon emission leakage associated with ES energy loss (see Remark 4), and the carbon accounting mechanisms for ES system owners remain unclear.

In this paper, we introduce a generic Carbon-aware Optimal Power Flow (C-OPF) methodology as the fundamental theory for carbon-aware power system decision-making that incorporates the active management of the grid’s carbon footprints. Built upon conventional OPF models, the C-OPF model (see model (7)) further integrates the carbon flow equations and constraints, as well as carbon-related objectives, to co-optimize power flow and carbon flow in the power grid. Essentially, C-OPF is a carbon-aware generalization of OPF, and it produces optimal decisions that satisfy carbon emission constraints and balance the power-related and carbon-related costs. The main contributions of this paper are threefold:

• 1)

  To our knowledge, this is the first work that introduces the generic C-OPF methodology and builds its mathematical model for lossy power networks with formulation examples. In particular, this paper rigorously establishes the conditions that guarantee the feasibility and solution uniqueness of the carbon flow equations (see Theorem 1), and presents the key properties of the C-OPF model.

• 2)

  We propose a reformulation approach to address the issue of unknown power flow directions in the C-OPF model by introducing dual power flow variables and complementarity constraints. This reformulation is exact and eliminates the necessity of knowing power flow directions in advance.

• 3)

  We develop two novel carbon footprint models for ES systems: one precisely models the dynamics of carbon emissions virtually stored in an ES unit and the carbon leakage associated with ES energy loss, while the other treats ES as a load during charging and a carbon-free generator during discharging. The corresponding carbon accounting schemes for ES owners are also introduced.

Furthermore, we develop a carbon-aware economic dispatch model as an example based on the proposed C-OPF method, and demonstrate the effectiveness of C-OPF through numerical experiments in comparison with OPF-based solutions.

The remainder of this paper is organized as follows: Section II introduces the concept and model of carbon flow as well as the use for carbon accounting. Section III presents the C-OPF method with the reformulation approach. Section IV introduces the carbon footprint models for energy storage systems. Numerical experiments are conducted in Section V, and conclusions are drawn in Section VI.

As mentioned above, the location-based method calculates the average emission factor (i.e., the ratio of total carbon emissions to the total generation energy) of an area grid over a defined period for demand-side carbon accounting. As illustrated in Figure 1, this AEF method treats the power system as a large “pool”, assuming that all generations and end-users are connected to one homogeneous and shared infrastructure, without consideration of the power network and power flow. In contrast, the carbon flow method views the carbon emissions from generators as virtual attachments to the power flow. These emissions are considered to be transmitted from generators through power networks and accumulate on the user side, forming carbon flows. The concept of carbon flow is analogous to a “water supply” system, where virtual carbon emissions accompanying power flows can be viewed as invisible particles contained in water flows that are delivered from water sources to users. The carbon emission intensity of electricity is analogous to the particle concentration of water, and a renewable generator is analogous to a pure water source. In this way, the carbon flow method aligns with the physical power grids and underlying power flows, enabling a temporally and spatially granular depiction of the grid’s carbon footprints.

In practice, the grid operators and utilities are expected to implement the carbon flow scheme, as these entities oversee power grid operations and possess detailed power flow profiles. Moreover, the carbon flow information (e.g., real-time nodal carbon intensities) represents useful grid emission signals, which can be sent to end-users (e.g., by displaying on smart meters) to inform user-side decarbonization decisions [37].

Consider a power network described by a connected graph $G(\mathcal{N},\mathcal{E})$, where $\mathcal{N}\!:=\!\{1,\cdots,N\}$ denotes the set of nodes and $\mathcal{E}\subset\mathcal{N}\times\mathcal{N}$ denotes the set of branches. As illustrated in Figure 2, the notations with $P$ denote the active power flow values (in the unit of MW), while the notations with $R$ represent the associated carbon flow rates (in the unit of ton$\mathrm{CO_{2}}$/h). The carbon (emission) intensity is defined as the ratio $w\!=\!\frac{R}{P}$ (in the unit of ton$\mathrm{CO_{2}}$/MWh) that describes the amount of carbon flow associated with one unit of power flow.

The active power flows adhere to the power balance equation (1) at each node $i$:

and $P_{ij}^{i}=P_{ij}^{\mathrm{loss}}+P_{ij}^{j}$ for each branch $ij\in\mathcal{E}$. Then, the carbon flow model is built upon the active power flows and the following two fundamental principles [14]:

• 1)

  Conservation Law of Nodal Carbon Mass: Similar to the power flow balance (1), the total carbon inflows equal the total carbon outflows at each node $i\in\mathcal{N}$, i.e.,

  $\displaystyle\sum_{k\in\mathcal{N}_{i}^{+}}\!\!R_{ki}^{i}+\sum_{g\in\mathcal{G}_{i}}\!R_{i,g}^{\mathrm{G}}=\sum_{j\in\mathcal{N}_{i}^{-}}\!\!R_{ij}^{i}+\sum_{l\in\mathcal{L}_{i}}\!R_{i,l}^{\mathrm{L}}.$

  (2)

  and $R_{ij}^{i}=R_{ij}^{\mathrm{loss}}+R_{ij}^{j}$ for each branch $ij\in\mathcal{E}$. In (2), $R_{i,g}^{\mathrm{G}}\!=\!w_{i,g}^{\mathrm{G}}P_{i,g}^{\mathrm{G}}$ is the generation carbon emission rate of generator-$g$ at node $i$.

• 2)

  Proportional Sharing Principle: At each node $i\in\mathcal{N}$, the allocation of total carbon inflows among all outflows is proportional to their active power flow values, i.e.,

  $\displaystyle R_{ij}^{i}=\frac{R_{i}^{\mathrm{in}}}{P_{i}^{\mathrm{in}}}\cdot P_{ij}^{i},\ \ R_{i,l}^{\mathrm{L}}=\frac{R_{i}^{\mathrm{in}}}{P_{i}^{\mathrm{in}}}\cdot P_{i,l}^{\mathrm{L}},$

  (3)

  where $P_{i}^{\mathrm{in}}\!:=\!\sum_{k\in\mathcal{N}_{i}^{+}}\!P_{ki}^{i}+\sum_{g\in\mathcal{G}_{i}}\!P_{i,g}^{\mathrm{G}}$ and $R_{i}^{\mathrm{in}}\!:=\!\sum_{k\in\mathcal{N}_{i}^{+}}\!R_{ki}^{i}+\sum_{g\in\mathcal{G}_{i}}\!R_{i,g}^{\mathrm{G}}$ denote the total power inflow and total carbon inflow at node $i\in\mathcal{N}$.

The nodal carbon intensity of each node $i$ is calculated as (4):

which is the ratio of total carbon inflow to the total power inflow. Equation (4) indicates that the power flow and power loss of each branch $ij\in\mathcal{E}$ share the same carbon emission intensity that equals to the nodal carbon intensity of the sending node, i.e., $\frac{R_{ij}^{i}}{P_{ij}^{i}}=\frac{R_{ij}^{j}}{P_{ij}^{j}}=\frac{R_{ij}^{\mathrm{loss}}}{P_{ij}^{\mathrm{loss}}}=w_{i}$.

Equation (4) is referred to as the carbon flow model or carbon flow equation, and it can be equivalently reformulated in the matrix form (5). The detailed derivation of (5a) is provided in Appendix A.

Here, $\bm{w}_{\mathrm{N}}\!:=\!(w_{i})_{i\in\mathcal{N}}\!\in\!\mathbb{R}^{N}$ and $\bm{r}_{\mathrm{G}}\!:=\!(\sum_{g\in\mathcal{G}_{i}}\!\!w_{i,g}^{\mathrm{G}}P_{i,g}^{\mathrm{G}})_{i\in\mathcal{N}}\!\in\!\mathbb{R}^{N}$ denote the column vectors that collect the nodal carbon intensities and nodal generation carbon flows, respectively. $\bm{P}_{\mathrm{N}}\!:=\!\text{diag}(P_{i}^{\mathrm{in}})\!\in\!\mathbb{R}^{N\times N}$ is the diagonal matrix whose $i$-th diagonal entry is the nodal active power inflow $P_{i}^{\mathrm{in}}$ to node $i$. $\bm{P}_{\mathrm{B}}\!\in\!\mathbb{R}^{N\times N}$ is the branch power inflow matrix that is built by letting $\bm{P}_{\mathrm{B}}[i,k]\!=\!P_{ki}^{i}$ and $\bm{P}_{\mathrm{B}}[k,i]\!=\!0$ if node $k$ sends power flow $P_{ki}^{i}$ to node $i$. Note that $P_{ki}^{i}$ and $P_{ki}^{k}$ are different due to the line power loss. Given power flow profiles, the linear equations (5a) can be solved to obtain $\bm{w}_{\mathrm{N}}$, e.g., through matrix factorization and backward and forward substitution, rather than directly computing the matrix inverse in (5b).

A critical fundamental question is whether the carbon flow equations (4) or (5) for a power network are feasible and admit a unique solution. Based on the matrix form (5), it translates to the invertibility of the carbon flow matrix $\bm{P}_{\mathrm{C}}:=\bm{P}_{\mathrm{N}}-\bm{P}_{\mathrm{B}}$, which implies the feasibility and solution uniqueness. Hence, this subsection studies the properties of the matrix $\bm{P}_{\mathrm{C}}$.

By definition, the matrix $\bm{P}_{\mathrm{C}}$ is diagonally dominant [38], and we define the set $\mathcal{J}$:

In practice, $\mathcal{J}\neq\emptyset$, i.e., there exist some rows $i\in\mathcal{J}$ of $\bm{P}_{\mathrm{C}}$ that are strictly diagonally dominant; such rows correspond to the nodes with generation power injections $P_{i,g}^{\mathrm{G}}$. Then, we establish the invertibility property of $\bm{P}_{\mathrm{C}}$ in Theorem 1.

Theorem 1 indicates that under the condition in Theorem 1, the matrix $\bm{P}_{\mathrm{C}}$ is invertible, and thus the carbon flow equations are feasible and have a unique solution regarding the nodal carbon intensities and other carbon flow values.

Given power flow results, one can solve the set of carbon flow equations (4) or (5) to obtain the nodal carbon intensities $\bm{w}_{\mathrm{N}}$, and then calculate all carbon flow values. According to the GHG Protocol [5, 3], the carbon accounting rules are:

• •

  Generator-$g$ at node $i$ shall account for the (Scope 1) direct carbon emission rate $R_{i,g}^{\mathrm{G}}=w_{i,g}^{\mathrm{G}}P_{i,g}^{\mathrm{G}}$;

• •

  Load-$l$ at node $i$ shall account for the (Scope 2) attributed carbon emission rate $R_{i,l}^{\mathrm{L}}=w_{i}P_{i,l}^{\mathrm{L}}$;

• •

  The power network owner shall account for the (Scope 2) attributed carbon emission rate associated with the network power loss, i.e., $\sum_{ij\in\mathcal{E}}\!R_{ij}^{\mathrm{loss}}=\sum_{ij\in\mathcal{E}}\!w_{i}P_{ij}^{\mathrm{loss}}$.

As the generation fuel mix and power flows change over time, the carbon flow profiles are time-varying, and the carbon footprints over a time period $\mathcal{T}$ (e.g., one day or one year) are the time accumulation of $R$. For example, the carbon footprint of load-$l$ at node $i$ is calculated as $\hat{E}_{i,l,\mathcal{T}}:=\sum_{t\in\mathcal{T}}\delta_{t}w_{i,t}P_{i,l,t}^{\mathrm{L}}$, where $\delta_{t}$ is the time interval (e.g., 1 hour or 15 minutes).

To enable the optimal management of carbon footprints in power system decision-making, we propose the generic C-OPF model (7) to co-optimize power flow and carbon flow.

Here, $\bm{x}$ denotes the decision variables subject to the feasible set $\mathcal{X}$. The specification of $\bm{x}$ and $\mathcal{X}$ depends on the practical applications. For instance, $\bm{x}$ denotes the generation decisions of various generators and $\mathcal{X}$ represents the generation capacity limits and ramping constraints in economic dispatch [41]. $\bm{x}$ can also be the load adjustment decisions in demand response, or the site and size decisions of new renewable generators in grid planning. $\bm{y}$ denotes the power flow-related variables, such as network voltage profiles and branch power flows. $\bm{z}$ represents the carbon flow-related variables, such as nodal carbon intensities $\bm{w}_{\mathrm{N}}$ and carbon flow values.

As mentioned in the introduction section and indicated by the set $\mathcal{N}_{i}^{+}$, the carbon flow equation $\eqref{eq:mcieq}$ requires the pre-determination of branch power flow directions to identify the power inflows for each node. However, the directions of branch power flows are generally unknown prior to solving the C-OPF problem. To address this issue, we introduce two non-negative power flow variables $\hat{P}_{ji}^{i}\!\geq\!0$ and $\hat{P}_{ij}^{i}\!\geq\!0$ for each branch $ij\in\mathcal{E}$ with $P_{ij}^{i}=\hat{P}_{ij}^{i}-\hat{P}_{ji}^{i}$. Specifically, $\hat{P}_{ji}^{i}$ and $\hat{P}_{ij}^{i}$ denote the power flow components from node $j$ to node $i$ and from node $i$ to node $j$, respectively, both of which are measured on the side of node $i$. Then, we can equivalently reformulate the carbon flow equation (III-A3) as (17), and also need to replace $P_{ij}^{i}$ with $\hat{P}_{ij}^{i}-\hat{P}_{ji}^{i}$ in the power flow equations (9) (or the DC power flow model (11)) and constraints (10a).

In (17a), we replace $\mathcal{N}_{i}^{+}$ (the set of neighbor nodes that send power to node $i$) by $\mathcal{N}_{i}$ (the set of all neighbor nodes of node $i$). In addition, the complementarity constraint (17c) is added to ensure that either $\hat{P}_{ji}^{i}$ or $\hat{P}_{ij}^{i}$ must be zero for each branch $ij$. Here, a useful trick to facilitate the nonlinear optimization is to replace the complementarity constraint (17c) by the relaxed constraint (18) [43], and this relaxation is exact due to (17b).

Alternatively, we can introduce a binary variable $\gamma_{ij}$ for $ij\in\mathcal{E}$ and linearize the complementarity constraint (17c) with

Note that both reformulations via (17) or (19) are equivalent to the original carbon flow equation (III-A3), but the branch flow directions do not need to be known in advance.

In [28], binary indicator variables are introduced to handle the unknown power flow directions in the carbon flow model. It results in a mixed-integer nonconvex quadratically constrained optimization problem, and a tailored penalty-based iterative algorithm is designed in [28] to solve the optimization problem through a number of iterations. In contrast, our proposed dual power flow variables reformulation method with the complementarity constraints (18) renders the C-OPF model a standard nonconvex optimization problem without integer variables. Therefore, the C-OPF model can be directly solved using off-the-shelf nonlinear optimizers such as IPOPT [44], without the need for designing ad hoc solution algorithms.

The “water tank” model of ES is conceptually aligned with the analogy of a “water supply” system used to explain the carbon flow model in Section II-A. Analogous to a water tank that stores both water and invisible particles, an ES system is viewed to store both electric energy $e_{i,t}^{\mathrm{es}}$ (in the unit of MWh) and virtual carbon emissions $E_{i,t}^{\mathrm{es}}$ (in the unit of ton$\mathrm{CO_{2}}$). We then define the internal ES carbon intensity $w_{i,t}^{\mathrm{es}}$ as (21):

Note that $e_{i,t}^{\mathrm{es}}$ should not be zero to make $w_{i,t}^{\mathrm{es}}$ well-defined. This can be achieved by setting the lower energy bound $\underline{e}_{i}$ in (20d) to be a small positive value rather than zero. It can also prevent numerical issues during the optimization process.

Based on the ES power model (20), we develop the dynamical carbon footprint model of the ES unit as (22) for $t\in\mathcal{T}$:

Model (22) implies that (virtual) carbon emissions are injected into the ES unit when it charges with electricity in the nodal carbon intensity $w_{i,t}$; and it releases the stored emissions back to the grid when it discharges with electricity in the ES carbon intensity $w_{i,t}^{\mathrm{es}}$. In particular, (22) models the carbon emission leakage associated with the energy loss during the storage, charging, and discharging processes of the ES unit. Figure 4 illustrates the carbon flow and power flow of the ES unit.

Alternatively, we can plug in the ES carbon intensity (21) into (22) to eliminate the variable $E_{i,t}^{\mathrm{es}}$ and reformulate the ES carbon footprint model equivalently as (23):²²2In (23b), the term $\frac{\delta_{t}}{\eta_{i}^{\mathrm{dc}}}P_{i,t}^{\mathrm{dc}}$ is dropped from the definition of $\lambda_{i,t}$, because $P_{i,t}^{\mathrm{dc}}=0$ when the ES unit charges and $\lambda_{i,t}=1$ when it discharges, which are not affected by this term.

which describes the dynamics of the ES carbon intensity $w_{i,t}^{{\mathrm{es}}}$. Interestingly, equation (23a) indicates an intuitive feature that the ES carbon intensity $w_{i,t+1}^{\mathrm{es}}$ at the next time $t\!+\!1$ is a convex combination of the ES carbon intensity $w_{i,t}$ at time $t$ and the nodal carbon intensity $w_{i,t}$ with the weight coefficient $\lambda_{i,t}\in[0,1]$. When charging, the ES unit behaves as if mixing the internally stored electricity with the newly charged electricity; when it discharges, $\lambda_{i,t}=1$ and the ES carbon intensity remains unchanged, i.e., $w_{i,t+1}^{\mathrm{es}}=w_{i,t}^{\mathrm{es}}$.

Accordingly, the carbon flow equation (17a) in the C-OPF model is modified as (IV-A) to incorporate the ES system.

From (IV-A), it is seen that an ES unit affects the nodal carbon intensities and carbon flow only when it discharges.

The proposed “water tank” model above requires continuous monitoring of the virtually stored carbon emissions within an ES system. To facilitate practical implementation, we propose an alternative model termed “load/carbon-free generator” to characterize the carbon footprints of an ES unit. Specifically, this model directly treats an ES unit as a load during charging and as a carbon-free clean generator during discharging. Accordingly, the carbon flow equation (17a) in the C-OPF model is modified as (IV-A) with $w_{i,t}^{\mathrm{es}}\equiv 0$ for all $t$, since the ES unit is regarded as a carbon-free generator during discharging.

In terms of carbon accounting, under the “load/carbon-free generator” model, the ES owner shall account for the (Scope 2) attributed carbon emissions $\hat{E}_{i,\mathcal{T}}^{\mathrm{es}}$ calculated by (26):

It implies that the virtual carbon emissions absorbed from the grid during the charging period accumulate locally at the ES unit and are not released back to the grid. Thus, the ES owner incurs higher Scope-2 carbon emissions compared with the carbon accounting scheme (25) under the “water tank” model. Nevertheless, under the “load/carbon-free generator” model, ES owners can make profits by acting as clean energy suppliers in the carbon-electricity market, e.g., selling renewable energy certificates (RECs) [9], since the discharging power is considered carbon-free. As a result, ES owners are incentivized to charge their ES units when the grid is clean with low carbon intensity to reduce their own carbon footprints (26); and they are incentivized to discharge when the grid is in high emission to make more profit, as the price of clean electricity or RECs is expected to be higher at that time. In this way, the “load/carbon-free generator” model is easy to use and naturally incentivizes the carbon-aware operation of ES systems. The simulation comparison between the “water tank” and “load/carbon-free generator” ES carbon footprint models are provided in Section V-E.

We note that alongside the two proposed carbon footprint models for ES systems, there could be other valid models. The use of different carbon models can lead to discrepancies in carbon accounting results, operational decision-making, and market design for ES systems. Hence, it is crucial to assess the impact and select an appropriate model that aligns with practical objectives and specific requirements.

Based on the C-OPF method (7), we develop a carbon-aware economic dispatch (C-ED) model (27) as a specific application example. It involves the multi-period optimal power scheduling of various generators and ES systems, while considering a complete AC power flow model, network operational constraints, and carbon flow equations and constraints. The C-ED model (27) is formulated as a nonconvex optimization problem and we solve it using the solver IPOPT [44].