Core

discharge!(EP::Model, inputs::Dict, setup::Dict)

This module defines the power decision variable $\Theta_{y,t} \forall y \in \mathcal{G}, t \in \mathcal{T}$, representing energy injected into the grid by resource $y$ by at time period $t$. This module additionally defines contributions to the objective function from variable costs of generation (variable O&M) from all resources $y \in \mathcal{G}$ over all time periods $t \in \mathcal{T}$:

\[\begin{aligned} Obj_{Var\_gen} = \sum_{y \in \mathcal{G} } \sum_{t \in \mathcal{T}}\omega_{t}\times(\pi^{VOM}_{y})\times \Theta_{y,t} \end{aligned}\]

investment_discharge!(EP::Model, inputs::Dict, setup::Dict)

This function defines the expressions and constraints keeping track of total available power generation/discharge capacity across all resources as well as constraints on capacity retirements. The total capacity of each resource is defined as the sum of the existing capacity plus the newly invested capacity minus any retired capacity. Note for storage and co-located resources, additional energy and charge power capacity decisions and constraints are defined in the storage and co-located VRE and storage module respectively.

\[\begin{aligned} & \Delta^{total}_{y,z} =(\overline{\Delta_{y,z}}+\Omega_{y,z}-\Delta_{y,z}) \forall y \in \mathcal{G}, z \in \mathcal{Z} \end{aligned}\]

One cannot retire more capacity than existing capacity.

\[\begin{aligned} &\Delta_{y,z} \leq \overline{\Delta_{y,z}} \hspace{4 cm} \forall y \in \mathcal{G}, z \in \mathcal{Z} \end{aligned}\]

For resources where $\overline{\Omega_{y,z}}$ and $\underline{\Omega_{y,z}}$ is defined, then we impose constraints on minimum and maximum power capacity.

\[\begin{aligned} & \Delta^{total}_{y,z} \leq \overline{\Omega}_{y,z} \hspace{4 cm} \forall y \in \mathcal{G}, z \in \mathcal{Z} \\ & \Delta^{total}_{y,z} \geq \underline{\Omega}_{y,z} \hspace{4 cm} \forall y \in \mathcal{G}, z \in \mathcal{Z} \end{aligned}\]

In addition, this function adds investment and fixed O&M related costs related to discharge/generation capacity to the objective function:

\[\begin{aligned} & \sum_{y \in \mathcal{G} } \sum_{z \in \mathcal{Z}} \left( (\pi^{INVEST}_{y,z} \times \overline{\Omega}^{size}_{y,z} \times \Omega_{y,z}) + (\pi^{FOM}_{y,z} \times \overline{\Omega}^{size}_{y,z} \times \Delta^{total}_{y,z})\right) \end{aligned}\]

write_multi_stage_capacities_discharge(outpath::String, settings_d::Dict)

This function writes the file capacities_multi_stage.csv to the Results directory. This file contains starting resource capacities from the first model stage and end resource capacities for the first and all subsequent model stages.

inputs:

• outpath – String which represents the path to the Results directory.
• settings_d - Dictionary containing settings dictionary configured in the multi-stage settings file multi_stage_settings.yml.

write_virtual_discharge(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)

Function for writing the "virtual" discharge of each storage technology. Virtual discharge is used to allow storage resources to contribute to the capacity reserve margin without actually discharging.

non_served_energy!(EP::Model, inputs::Dict, setup::Dict)

This function defines the non-served energy/curtailed demand decision variable $\Lambda_{s,t,z} \forall s \in \mathcal{S}, \forall t \in \mathcal{T}, z \in \mathcal{Z}$, representing the total amount of demand curtailed in demand segment $s$ at time period $t$ in zone $z$. The first segment of non-served energy, $s=1$, is used to denote the cost of involuntary demand curtailment (e.g. emergency load shedding or rolling blackouts), specified as the value of $n_{1}^{slope}$. Additional segments, $s \geq 2$ can be used to specify a segment-wise approximation of a price elastic demand curve, or segments of price-responsive curtailable demands (aka demand response). Each segment denotes a price/cost at which the segment of demand is willing to curtail consumption, $n_{s}^{slope}$, representing the marginal willingness to pay for electricity of this segment of demand (or opportunity cost incurred when demand is not served) and a maximum quantity of demand in this segment, $n_{s}^{size}$, specified as a share of demand in each zone in each time step, $D_{t,z}.$ Note that the current implementation assumes demand segments are an equal share of hourly demand in all zones. This function defines contributions to the objective function from the cost of non-served energy/curtailed demand from all demand curtailment segments $s \in \mathcal{S}$ over all time periods $t \in \mathcal{T}$ and all zones $z \in \mathcal{Z}$:

\[\begin{aligned} Obj_{NSE} = \sum_{s \in \mathcal{S} } \sum_{t \in \mathcal{T}} \sum_{z \in \mathcal{Z}}\omega_{t} \times n_{s}^{slope} \times \Lambda_{s,t,z} \end{aligned}\]

Contributions to the power balance expression from non-served energy/curtailed demand from each demand segment $s \in \mathcal{S}$ are also defined as:

\[\begin{aligned} PowerBal_{NSE} = \sum_{s \in \mathcal{S} } \Lambda_{s,t,z} \hspace{4 cm} \forall s \in \mathcal{S}, t \in \mathcal{T} \end{aligned}\]

Bounds on curtailable demand Demand curtailed in each segment of curtailable demands $s \in \mathcal{S}$ cannot exceed a maximum allowable share of demand:

\[\begin{aligned} \Lambda_{s,t,z} \leq (n_{s}^{size} \times D_{t,z}) \hspace{4 cm} \forall s \in \mathcal{S}, t \in \mathcal{T}, z\in \mathcal{Z} \end{aligned}\]

Additionally, total demand curtailed in each time step cannot exceed total demand:

\[\begin{aligned} \sum_{s \in \mathcal{S} } \Lambda_{s,t,z} \leq D_{t,z} \hspace{4 cm} \forall t \in \mathcal{T}, z\in \mathcal{Z} \end{aligned}\]

load_operational_reserves!(setup::Dict,path::AbstractString, inputs::Dict)

Read input parameters related to frequency regulation and operating reserve requirements

operational_reserves!(EP::Model, inputs::Dict, setup::Dict)

This function sets up reserve decisions and constraints, using the operationalreservescore()and operational_reserves_contingency() functions.

operational_reserves_contingency!(EP::Model, inputs::Dict, setup::Dict)

This function establishes several different versions of contingency reserve requirement expression, $CONTINGENCY$ used in the operationalreservescore() function below.

Contingency operational reserves represent requirements for upward ramping capability within a specified time frame to compensated for forced outages or unplanned failures of generators or transmission lines (e.g. N-1 contingencies).

There are three options for the $Contingency$ expression, depending on user settings: 1. a static contingency, in which the contingency requirement is set based on a fixed value (in MW) specified in the '''Operational_reserves.csv''' input file; 2. a dynamic contingency based on installed capacity decisions, in which the largest 'installed' generator is used to determine the contingency requirement for all time periods; and 3. dynamic unit commitment based contingency, in which the largest 'committed' generator in any time period is used to determine the contingency requirement in that time period.

Note that the two dynamic contigencies are only available if unit commitment is being modeled.

Static contingency Option 1 (static contingency) is expressed by the following constraint:

\[\begin{aligned} Contingency = \epsilon^{contingency} \end{aligned}\]

where $\epsilon^{contingency}$ is static contingency requirement in MWs.

Dynamic capacity-based contingency Option 2 (dynamic capacity-based contingency) is expressed by the following constraints:

\[\begin{aligned} & Contingency \geq \Omega^{size}_{y,z} \times \alpha^{Contingency,Aux}_{y,z} & \forall y \in \mathcal{UC}, z \in \mathcal{Z}\\ & \alpha^{Contingency,Aux}_{y,z} \leq \Delta^{\text{total}}_{y,z} & \forall y \in \mathcal{UC}, z \in \mathcal{Z}\\ & \Delta^{\text{total}}_{y,z} \leq M_y \times \alpha^{Contingency,Aux}_{y,z} & \forall y \in \mathcal{UC}, z \in \mathcal{Z}\\ \end{aligned}\]

where $M_y$ is a `big M' constant equal to the largest possible capacity that can be installed for generation cluster $y$, and $\alpha^{Contingency,Aux}_{y,z} \in [0,1]$ is a binary auxiliary variable that is forced by the second and third equations above to be 1 if the total installed capacity $\Delta^{\text{total}}_{y,z} > 0$ for any generator $y \in \mathcal{UC}$ and zone $z$, and can be 0 otherwise. Note that if the user specifies contingency option 2, and is also using the linear relaxation of unit commitment constraints, the capacity size parameter for units in the set $\mathcal{UC}$ must still be set to a discrete unit size for this contingency to work as intended.

Dynamic commitment-based contingency Option 3 (dynamic commitment-based contingency) is expressed by the following set of constraints:

\[\begin{aligned} & Contingency \geq \Omega^{size}_{y,z} \times Contingency\_Aux_{y,z,t} & \forall y \in \mathcal{UC}, z \in \mathcal{Z}\\ & Contingency\_Aux_{y,z,t} \leq \nu_{y,z,t} & \forall y \in \mathcal{UC}, z \in \mathcal{Z}\\ & \nu_{y,z,t} \leq M_y \times Contingency\_Aux_{y,z,t} & \forall y \in \mathcal{UC}, z \in \mathcal{Z}\\ \end{aligned}\]

where $M_y$ is a `big M' constant equal to the largest possible capacity that can be installed for generation cluster $y$, and $Contingency\_Aux_{y,z,t} \in [0,1]$ is a binary auxiliary variable that is forced by the second and third equations above to be 1 if the commitment state for that generation cluster $\nu_{y,z,t} > 0$ for any generator $y \in \mathcal{UC}$ and zone $z$ and time period $t$, and can be 0 otherwise. Note that this dynamic commitment-based contingency can only be specified if discrete unit commitment decisions are used (e.g. it will not work if relaxed unit commitment is used).

operational_reserves_core!(EP::Model, inputs::Dict, setup::Dict)

This function creates decision variables related to frequency regulation and reserves provision and constraints setting overall system requirements for regulation and operating reserves.

Regulation and reserves decisions $f_{y,t,z} \geq 0$ is the contribution of generation or storage resource $y \in Y$ in time $t \in T$ and zone $z \in Z$ to frequency regulation

\[r_{y,t,z} \geq 0\]

is the contribution of generation or storage resource $y \in Y$ in time $t \in T$ and zone $z \in Z$ to operating reserves up

We assume frequency regulation is symmetric (provided in equal quantity towards both upwards and downwards regulation). To reduce computational complexity, operating reserves are only modeled in the upwards direction, as downwards reserves requirements are rarely binding in practice.

Storage resources ($y \in \mathcal{O}$) have two pairs of auxilary variables to reflect contributions to regulation and reserves when charging and discharging, where the primary variables ($f_{y,z,t}$ and $r_{y,z,t}$) becomes equal to sum of these auxilary variables.

Co-located VRE-STOR resources are described further in the reserves function for colocated VRE and storage resources (vre_stor_operational_reserves!()).

Unmet operating reserves

\[unmet\_rsv_{t} \geq 0\]

denotes any shortfall in provision of operating reserves during each time period $t \in T$

There is a penalty $C^{rsv}$ added to the objective function to penalize reserve shortfalls, equal to:

\[\begin{aligned} C^{rvs} = \sum_{t \in T} \omega_t \times unmet\_rsv_{t} \end{aligned}\]

Frequency regulation requirements

Total requirements for frequency regulation (aka primary reserves) in each time step $t$ are specified as fractions of hourly demand (to reflect demand forecast errors) and variable renewable avaialblity in the time step (to reflect wind and solar forecast errors).

\[\begin{aligned} & \sum_{y \in Y, z \in Z} f_{y,t,z} \geq \epsilon^{demand}_{reg} \times \sum_{z \in Z} \mathcal{D}_{z,t} + \epsilon^{vre}_{reg} \times (\sum_{z \in Z} \rho^{max}_{y,z,t} \times \Delta^{\text{total}}_{y,z} \\ & + \sum_{z \in Z} \rho^{max,pv}_{y,z,t} \times \Delta^{\text{total,pv}}_{y,z} + \sum_{z \in Z} \rho^{max,wind}_{y,z,t} \times \Delta^{\text{total,wind}}_{y,z}) \quad \forall t \in T \end{aligned}\]

where $\mathcal{D}_{z,t}$ is the forecasted electricity demand in zone $z$ at time $t$ (before any demand flexibility); $\rho^{max}_{y,z,t}$ is the forecasted capacity factor for standalone variable renewable resources $y \in VRE$ and zone $z$ in time step $t$; $\rho^{max,pv}_{y,z,t}$ is the forecasted capacity factor for co-located solar PV resources $y \in \mathcal{VS}^{pv}$ and zone $z$ in time step $t$; $\rho^{max,wind}_{y,z,t}$ is the forecasted capacity factor for co-located wind resources $y \in \mathcal{VS}^{pv}$ and zone $z$ in time step $t$; $\Delta^{\text{total,pv}}_{y,z}$ is the total installed capacity of co-located solar PV resources $y \in \mathcal{VS}^{pv}$ and zone $z$; $\Delta^{\text{total,wind}}_{y,z}$ is the total installed capacity of co-located wind resources $y \in \mathcal{VS}^{wind}$ and zone $z$; and $\epsilon^{demand}_{reg}$ and $\epsilon^{vre}_{reg}$ are parameters specifying the required frequency regulation as a fraction of forecasted demand and variable renewable generation.

Operating reserve requirements

Total requirements for operating reserves in the upward direction (aka spinning reserves or contingency reserces or secondary reserves) in each time step $t$ are specified as fractions of time step's demand (to reflect demand forecast errors) and variable renewable avaialblity in the time step (to reflect wind and solar forecast errors) plus the largest planning contingency (e.g. potential forced generation outage).

\[\begin{aligned} & \sum_{y \in Y, z \in Z} r_{y,z,t} + r^{unmet}_{t} \geq \epsilon^{demand}_{rsv} \times \sum_{z \in Z} \mathcal{D}_{z,t} + \epsilon^{vre}_{rsv} \times (\sum_{z \in Z} \rho^{max}_{y,z,t} \times \Delta^{\text{total}}_{y,z} \\ & + \sum_{z \in Z} \rho^{max,pv}_{y,z,t} \times \Delta^{\text{total,pv}}_{y,z} + \sum_{z \in Z} \rho^{max,wind}_{y,z,t} \times \Delta^{\text{total,wind}}_{y,z}) + Contingency \quad \forall t \in T \end{aligned}\]

where $\mathcal{D}_{z,t}$ is the forecasted electricity demand in zone $z$ at time $t$ (before any demand flexibility); $\rho^{max}_{y,z,t}$ is the forecasted capacity factor for standalone variable renewable resources $y \in VRE$ and zone $z$ in time step $t$; $\rho^{max,pv}_{y,z,t}$ is the forecasted capacity factor for co-located solar PV resources $y \in \mathcal{VS}^{pv}$ and zone $z$ in time step $t$; $\rho^{max,wind}_{y,z,t}$ is the forecasted capacity factor for co-located wind resources $y \in \mathcal{VS}^{wind}$ and zone $z$ in time step $t$; $\Delta^{\text{total}}_{y,z}$ is the total installed capacity of standalone variable renewable resources $y \in VRE$ and zone $z$; $\Delta^{\text{total,pv}}_{y,z}$ is the total installed capacity of co-located solar PV resources $y \in \mathcal{VS}^{pv}$ and zone $z$; $\Delta^{\text{total,wind}}_{y,z}$ is the total installed capacity of co-located wind resources $y \in \mathcal{VS}^{wind}$ and zone $z$; and $\epsilon^{demand}_{rsv}$ and $\epsilon^{vre}_{rsv}$ are parameters specifying the required contingency reserves as a fraction of forecasted demand and variable renewable generation. $Contingency$ is an expression defined in the operational_reserves_contingency!() function meant to represent the largest N-1 contingency (unplanned generator outage) that the system operator must carry operating reserves to cover and depends on how the user wishes to specify contingency requirements.

function dcopf_transmission!(EP::Model, inputs::Dict, setup::Dict)

The addtional constraints imposed upon the line flows in the case of DC-OPF are as follows: For the definition of the line flows, in terms of the voltage phase angles:

\[\begin{aligned} & \Phi_{l,t}=\mathcal{B}_{l} \times (\sum_{z\in \mathcal{Z}}{(\varphi^{map}_{l,z} \times \theta_{z,t})}) \quad \forall l \in \mathcal{L}, \; \forall t \in \mathcal{T}\\ \end{aligned}\]

For imposing the constraint of maximum allowed voltage phase angle difference across lines:

\[\begin{aligned} & \sum_{z\in \mathcal{Z}}{(\varphi^{map}_{l,z} \times \theta_{z,t})} \leq \Delta \theta^{\max}_{l} \quad \forall l \in \mathcal{L}, \forall t \in \mathcal{T}\\ & \sum_{z\in \mathcal{Z}}{(\varphi^{map}_{l,z} \times \theta_{z,t})} \geq -\Delta \theta^{\max}_{l} \quad \forall l \in \mathcal{L}, \forall t \in \mathcal{T}\\ \end{aligned}\]

Finally, we enforce the reference voltage phase angle constraint:

\[\begin{aligned} \theta_{1,t} = 0 \quad \forall t \in \mathcal{T} \end{aligned}\]

function investment_transmission!(EP::Model, inputs::Dict, setup::Dict)

This function model transmission expansion and adds transmission reinforcement or construction costs to the objective function. Transmission reinforcement costs are equal to the sum across all lines of the product between the transmission reinforcement/construction cost, $pi^{TCAP}_{l}$, times the additional transmission capacity variable, $\bigtriangleup\varphi^{cap}_{l}$.

\[\begin{aligned} & \sum_{l \in \mathcal{L}}\left(\pi^{TCAP}_{l} \times \bigtriangleup\varphi^{cap}_{l}\right) \end{aligned}\]

Note that fixed O\&M and replacement capital costs (depreciation) for existing transmission capacity is treated as a sunk cost and not included explicitly in the GenX objective function. Accounting for Transmission Between Zones Available transmission capacity between zones is set equal to the existing line's maximum power transfer capacity, $\overline{\varphi^{cap}_{l}}$, plus any transmission capacity added on that line (for lines eligible for expansion in the set $\mathcal{E}$).

\[\begin{aligned} &\varphi^{cap}_{l} = \overline{\varphi^{cap}_{l}} , &\quad \forall l \in (\mathcal{L} \setminus \mathcal{E} ),\forall t \in \mathcal{T}\\ % trasmission expansion &\varphi^{cap}_{l} = \overline{\varphi^{cap}_{l}} + \bigtriangleup\varphi^{cap}_{l} , &\quad \forall l \in \mathcal{E},\forall t \in \mathcal{T} \end{aligned}\]

The additional transmission capacity, $\bigtriangleup\varphi^{cap}_{l} $, is constrained by a maximum allowed reinforcement, $\overline{\bigtriangleup\varphi^{cap}_{l}}$, for each line $l \in \mathcal{E}$.

\[\begin{aligned} & \bigtriangleup\varphi^{cap}_{l} \leq \overline{\bigtriangleup\varphi^{cap}_{l}}, &\quad \forall l \in \mathcal{E} \end{aligned}\]

transmission!(EP::Model, inputs::Dict, setup::Dict)

This function establishes decisions, expressions, and constraints related to transmission power flows between model zones and associated transmission losses (if modeled).

Power flow and transmission loss terms are also added to the power balance constraint for each zone:

\[\begin{aligned} & - \sum_{l\in \mathcal{L}}{(\varphi^{map}_{l,z} \times \Phi_{l,t})} - \frac{1}{2} \sum_{l\in \mathcal{L}}{(\varphi^{map}_{l,z} \times \beta_{l,t}(\cdot))} \end{aligned}\]

Power flows, $\Phi_{l,t}$, on each line $l$ into or out of a zone (defined by the network map $\varphi^{map}_{l,z}$), are considered in the demand balance equation for each zone. By definition, power flows leaving their reference zone are positive, thus the minus sign is used for this term. Losses due to power flows increase demand, and one-half of losses across a line linking two zones are attributed to each connected zone. The losses function $\beta_{l,t}(\cdot)$ will depend on the configuration used to model losses (see below). Accounting for Transmission Between Zones Power flow, $\Phi_{l,t}$, on each line (or more likely a `path' aggregating flows across multiple parallel lines) is constrained to be less than or equal to the line's power transfer capacity, $\varphi^{cap}_{l}$, plus any transmission capacity added on that line (for lines eligible for expansion in the set $\mathcal{E}$). The additional transmission capacity, $\bigtriangleup\varphi^{cap}_{l} $, is constrained by a maximum allowed reinforcement, $\overline{\bigtriangleup\varphi^{cap}_{l}}$, for each line $l \in \mathcal{E}$.

\[\begin{aligned} % trasmission constraints &-\varphi^{cap}_{l} \leq \Phi_{l,t} \leq \varphi^{cap}_{l} , &\quad \forall l \in \mathcal{L},\forall t \in \mathcal{T}\\ \end{aligned}\]

Accounting for Transmission Losses Transmission losses due to power flows can be accounted for in three different ways. The first option is to neglect losses entirely, setting the value of the losses function to zero for all lines at all hours. The second option is to assume that losses are a fixed percentage, $\varphi^{loss}_{l}$, of the magnitude of power flow on each line, $\mid \Phi_{l,t} \mid$ (e.g., losses are a linear function of power flows). Finally, the third option is to calculate losses, $\ell_{l,t}$, by approximating a quadratic-loss function of power flow across the line using a piecewise-linear function with total number of segments equal to the size of the set $\mathcal{M}$.

\[\begin{aligned} %configurable losses formulation & \beta_{l,t}(\cdot) = \begin{cases} 0 & \text{if~} \text{losses.~0} \\ \\ \varphi^{loss}_{l}\times \mid \Phi_{l,t} \mid & \text{if~} \text{losses.~1} \\ \\ \ell_{l,t} &\text{if~} \text{losses.~2} \end{cases}, &\quad \forall l \in \mathcal{L},\forall t \in \mathcal{T} \end{aligned}\]

For the second option, an absolute value approximation is utilized to calculate the magnitude of the power flow on each line (reflecting the fact that negative power flows for a line linking nodes $i$ and $j$ represents flows from node $j$ to $i$ and causes the same magnitude of losses as an equal power flow from $i$ to $j$). This absolute value function is linearized such that the flow in the line must be equal to the subtraction of the auxiliary variable for flow in the positive direction, $\Phi^{+}_{l,t}$, and the auxiliary variable for flow in the negative direction, $\Phi^{+}_{l,t}$, of the line. Then, the magnitude of the flow is calculated as the sum of the two auxiliary variables. The sum of positive and negative directional flows are also constrained by the line flow capacity.

\[\begin{aligned} % trasmission losses simple &\Phi_{l,t} = \Phi^{+}_{l,t} - \Phi^{-}_{l,t}, &\quad \forall l \in \mathcal{L}, \forall t \in \mathcal{T}\\ &\mid \Phi_{l,t} \mid = \Phi^{+}_{l,t} + \Phi^{-}_{l,t}, &\quad \forall l \in \mathcal{L}, \forall t \in \mathcal{T}\\ &\Phi^{+}_{l,t} + \Phi^{-}_{l,t} \leq \varphi^{cap}_{l}, &\quad \forall l \in \mathcal{L}, \forall t \in \mathcal{T} \end{aligned}\]

If discrete unit commitment decisions are modeled, ``phantom losses'' can be observed wherein the auxiliary variables for flows in both directions ($\Phi^{+}_{l,t}$ and $\Phi^{-}_{l,t}$) are both increased to produce increased losses so as to avoid cycling a thermal generator and incurring start-up costs or opportunity costs related to minimum down times. This unrealistic behavior can be eliminated via inclusion of additional constraints and a set of auxiliary binary variables, $ON^{+}_{l,t} \in {0,1} \forall l \in \mathcal{L}$. Then the following additional constraints are created:

\[\begin{aligned} \Phi^{+}_{l,t} \leq TransON^{+}_{l,t}, &\quad \forall l \in \mathcal{L}, \forall t \in \mathcal{T}\\ \Phi^{-}_{l,t} \leq \varphi^{cap}_{l} -TransON^{+}_{l,t}, &\quad \forall l \in \mathcal{L}, \forall t \in \mathcal{T} \end{aligned}\]

where $TransON^{+}_{l,t}$ is a continuous variable, representing the product of the binary variable $ON^{+}_{l,t}$ and the expression, $\varphi^{cap}_{l}$. This product cannot be defined explicitly, since it will lead to a bilinear expression involving two variables. Instead, we enforce this definition via the Glover's Linearization as shown below (also referred McCormick Envelopes constraints for bilinear expressions, which is exact when one of the variables is binary).

\[\begin{aligned} TransON^{+}_{l,t} \leq (\overline{varphi^{cap}_{l}} + \overline{\bigtriangleup\varphi^{cap}_{l}}) \times TransON^{+}_{l,t}, &\quad \forall l \in \mathcal{L}, \forall t \in \mathcal{T} \\ TransON^{+}_{l,t} \leq \varphi^{cap}_{l}, &\quad \forall l \in \mathcal{L}, \forall t \in \mathcal{T} \\ TransON^{+}_{l,t} \geq \varphi^{cap}_{l} - (\overline{\varphi^{cap}_{l}} + \overline{\bigtriangleup\varphi^{cap}_{l}}) \times(1- TransON^{+}_{l,t}), &\quad \forall l \in \mathcal{L}, \forall t \in \mathcal{T} \\ \end{aligned}\]

These constraints permit only the positive or negative auxiliary flow variables to be non-zero at a given time period, not both. For the third option, losses are calculated as a piecewise-linear approximation of a quadratic function of power flows. In order to do this, we represent the absolute value of the line flow variable by the sum of positive stepwise flow variables $(\mathcal{S}^{+}_{m,l,t}, \mathcal{S}^{-}_{m,l,t})$, associated with each partition of line losses computed using the corresponding linear expressions. This can be understood as a segmentwise linear fitting (or first order approximation) of the quadratic losses function. The first constraint below computes the losses a the accumulated sum of losses for each linear stepwise segment of the approximated quadratic function, including both positive domain and negative domain segments. A second constraint ensures that the stepwise variables do not exceed the maximum size per segment. The slope and maximum size for each segment are calculated as per the method in \cite{Zhang2013}.

\[\begin{aligned} & \ell_{l,t} = \frac{\varphi^{ohm}_{l}}{(\varphi^{volt}_{l})^2}\bigg( \sum_{m \in \mathcal{M}}( S^{+}_{m,l}\times \mathcal{S}^{+}_{m,l,t} + S^{-}_{m,l}\times \mathcal{S}^{-}_{m,l,t}) \bigg), &\quad \forall l \in \mathcal{L}, \forall t \in \mathcal{T} \\ & \text{\quad Where:} \\ & \quad S^{+}_{m,l} = \frac{2+4 \times \sqrt{2}\times (m-1)}{1+\sqrt{2} \times (2 \times M-1)} (\overline{\varphi^{cap}_{l}} + \overline{\bigtriangleup\varphi^{cap}_{l}}) &\quad \forall m \in [1 \colon M], l \in \mathcal{L} \\ & \quad S^{-}_{m,l} = \frac{2+4 \times \sqrt{2}\times (m-1)}{1+\sqrt{2} \times (2 \times M-1)} (\overline{\varphi^{cap}_{l}} + \overline{\bigtriangleup\varphi^{cap}_{l}}) &\quad \forall m \in [1 \colon M], l \in \mathcal{L}\\ & \\ & \mathcal{S}^{+}_{m,l,t}, \mathcal{S}^{-}_{m,l,t} <= \overline{\mathcal{S}_{m,l}} &\quad \forall m \in [1:M], l \in \mathcal{L}, t \in \mathcal{T} \\ & \text{\quad Where:} \\ & \quad \overline{S_{l,z}} = \begin{cases} \frac{(1+\sqrt{2})}{1+\sqrt{2} \times (2 \times M-1)} (\overline{\varphi^{cap}_{l}} + \overline{\bigtriangleup\varphi^{cap}_{l}}) & \text{if~} m = 1 \\ \frac{2 \times \sqrt{2} }{1+\sqrt{2} \times (2 \times M-1)} (\overline{\varphi^{cap}_{l}} + \overline{\bigtriangleup\varphi^{cap}_{l}}) & \text{if~} m > 1 \end{cases} \end{aligned}\]

Next, a constraint ensures that the sum of auxiliary segment variables ($m \geq 1$) minus the "zero" segment (which allows values to go into the negative domain) from both positive and negative domains must total the actual power flow across the line, and a constraint ensures that the sum of negative and positive flows do not exceed the flow capacity of the line.

\[\begin{aligned} &\sum_{m \in [1:M]} (\mathcal{S}^{+}_{m,l,t}) - \mathcal{S}^{+}_{0,l,t} = \Phi_{l,t}, &\quad \forall l \in \mathcal{L}, \forall t \in \mathcal{T}\\ &\sum_{m \in [1:M]} (\mathcal{S}^{-}_{m,l,t}) - \mathcal{S}^{-}_{0,l,t} = - \Phi_{l,t} \end{aligned}\]

As with losses option 2, this segment-wise approximation of a quadratic loss function also permits 'phantom losses' to avoid cycling thermal units when discrete unit commitment decisions are modeled. In this case, the additional constraints below are also added to ensure that auxiliary segments variables do not exceed maximum value per segment and that they are filled in order; i.e., one segment cannot be non-zero unless prior segment is at its maximum value. Binary constraints deal with absolute value of power flow on each line. If the flow is positive, $\mathcal{S}^{+}_{0,l,t}$ must be zero; if flow is negative, $\mathcal{S}^{+}_{0,l,t}$ must be positive and takes on value of the full negative flow, forcing all $\mathcal{S}^{+}_{m,l,t}$ other segments ($m \geq 1$) to be zero. Conversely, if the flow is negative, $\mathcal{S}^{-}_{0,l,t}$ must be zero; if flow is positive, $\mathcal{S}^{-}_{0,l,t}$ must be positive and takes on value of the full positive flow, forcing all $\mathcal{S}^{-}_{m,l,t}$ other segments ($m \geq 1$) to be zero. Requiring segments to fill in sequential order and binary variables to ensure variables reflect the actual direction of power flows are both necessary to eliminate ``phantom losses'' from the solution. These constraints and binary decisions are ommited if the model is fully linear.

\[\begin{aligned} &\mathcal{S}^{+}_{m,l,t} <= \overline{\mathcal{S}_{m,l}} \times ON^{+}_{m,l,t}, &\quad \forall m \in [1:M], \forall l \in \mathcal{L}, \forall t \in \mathcal{T}\\ &\mathcal{S}^{-}_{m,l,t} <= \overline{\mathcal{S}_{m,l}} \times ON^{-}_{m,l,t}, &\quad \forall m \in[1:M], \forall l \in \mathcal{L}, \forall t \in \mathcal{T}\\ &\mathcal{S}^{+}_{m,l,t} \geq ON^{+}_{m+1,l,t} \times \overline{\mathcal{S}_{m,l}}, &\quad \forall m \in [1:M], \forall l \in \mathcal{L}, \forall t \in \mathcal{T}\\ &\mathcal{S}^{-}_{m,l,t} \geq ON^{-}_{m+1,l,t} \times \overline{\mathcal{S}_{m,l}} , &\quad \forall m \in [1:M], \forall l \in \mathcal{L}, \forall t \in \mathcal{T}\\ &\mathcal{S}^{+}_{0,l,t} \leq \varphi^{max}_{l} \times (1- ON^{+}_{1,l,t}), &\quad \forall l \in \mathcal{L}, \forall t \in \mathcal{T}\\ &\mathcal{S}^{-}_{0,l,t} \leq \varphi^{max}_{l} \times (1- ON^{-}_{1,l,t}), &\quad \forall l \in \mathcal{L}, \forall t \in \mathcal{T} \end{aligned}\]

ucommit!(EP::Model, inputs::Dict, setup::Dict)

This function creates decision variables and cost expressions associated with thermal plant unit commitment or start-up and shut-down decisions (cycling on/off)

Unit commitment decision variables:

This function defines the following decision variables:

\[\nu_{y,t,z}\]

designates the commitment state of generator cluster $y$ in zone $z$ at time $t$; $\chi_{y,t,z}$ represents number of startup decisions in cluster $y$ in zone $z$ at time $t$; $\zeta_{y,t,z}$ represents number of shutdown decisions in cluster $y$ in zone $z$ at time $t$.

Cost expressions:

The total cost of start-ups across all generators subject to unit commitment ($y \in UC$) and all time periods, t is expressed as:

\[\begin{aligned} C^{start} = \sum_{y \in UC, t \in T} \omega_t \times start\_cost_{y,t} \times \chi_{y,t} \end{aligned}\]

The sum of start-up costs is added to the objective function.

co2!(EP::Model, inputs::Dict)

This function creates expressions to account for CO2 emissions as well as captured and sequestrated CO2 from thermal generators. It also has the capability to model the negative CO2 emissions from bioenergy with carbon capture and storage.

***** Expressions *****

For thermal generators which combust fuels (e.g., coal, natural gas, and biomass), the net CO2 emission to the environment is a function of fuel consumption, CO2 emission factor, CO2 capture fraction, and whether the feedstock is biomass. Biomass is a factor in this equation because biomass generators are assumed to generate zero net CO2 emissions, or negative net CO2 emissions in the case that the CO2 they emit is captured and sequestered underground.

If a user wishes to represent a generator that combusts biomass, then in the resource .csv files, the "Biomass" column (boolean, 1 or 0), which represents if a generator $y$ uses biomass or not, should be set to 1. The CO2 emissions from such a generator will be assumed to be zero without CCS and negative with CCS.

The CO2 emissions from generator $y$ at time $t$ are determined by total fuel consumption (MMBTU) multiplied by the CO2 content of the fuel (tCO2/MMBTU), and by (1 - Biomass [0 or 1] - CO2 capture fraction [a fraction, between 0 - 1]). The CO2 capture fraction could be differernt during the steady-state and startup events (generally startup events have a lower CO2 capture fraction), so we use distinct CO2 capture fractions to determine the emissions. In short, the CO2 emissions for a generator depend on the CO2 emission factor from fuel combustion, the CO2 capture fraction, and whether the generator uses biomass.

\[\begin{aligned} eEmissionsByPlant_{g,t} = (1-Biomass_y- CO2\_Capture\_Fraction_y) * vFuel_{y,t} * CO2_{content} + (1-Biomass_y- CO2\_Capture\_Fraction\_Startup_y) * eStartFuel_{y,t} * CO2_{content} \hspace{1cm} \forall y \in G, \forall t \in T, Biomass_y \in {{0,1}} \end{aligned}\]

Where $Biomass_y$ represents a binary variable (1 or 0) that determines if the generator $y$ uses biomass, and $CO2\_Capture\_Fraction_y$ represents a fraction for CO2 capture rate.

In addition to CO2 emissions, for generators with a non-zero CO2 capture rate, we also determine the amount of CO2 being captured and sequestered. The CO2 emissions from generator $y$ at time $t$, denoted by $eEmissionsCaptureByPlant_{g,t}$, are determined by total fuel consumption (MMBTU) multiplied by the $CO_2$ content of the fuel (tCO2/MMBTU), times CO2 capture rate.

\[\begin{aligned} eEmissionsCaptureByPlant_{g,t} = CO2\_Capture\_Fraction_y * vFuel_{y,t} * CO2_{content} + CO2\_Capture\_Fraction\_Startup_y * eStartFuel_{y,t} * CO2_{content} \hspace{1cm} \forall y \in G, \forall t \in T \end{aligned}\]

write_co2(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)

Function for reporting time-dependent CO2 emissions by zone.

fuel!(EP::Model, inputs::Dict, setup::Dict)

This function creates expressions to account for total fuel consumption (e.g., coal, natural gas, hydrogen, etc). It also has the capability to model heat rates that are a function of load via a piecewise-linear approximation.

***** Expressions ****** Users have two options to model the fuel consumption as a function of power generation: (1). Use a constant heat rate, regardless of the minimum load or maximum load; and (2). Use the PiecewiseFuelUsage-related parameters to model the fuel consumption via a piecewise-linear approximation of the heat rate curves. By using this option, users can represent the fact that most generators have a decreasing heat rate as a function of load.

(1). Constant heat rate. The fuel consumption for power generation $vFuel_{y,t}$ is determined by power generation ($vP_{y,t}$) mutiplied by the corresponding heat rate ($Heat\_Rate_y$). The fuel costs for power generation and start fuel for a plant $y$ at time $t$, denoted by $eCFuelOut_{y,t}$ and $eFuelStart$, are determined by fuel consumption ($vFuel_{y,t}$ and $eStartFuel$) multiplied by the fuel costs ($/MMBTU) (2). Piecewise-linear approximation With this formulation, the heat rate of generators becomes a function of load. In reality this relationship takes a nonlinear form, but we model it through a piecewise-linear approximation:

\[\begin{aligned} vFuel_{y,t} >= vP_{y,t} * h_{y,x} + U_{g,t}* f_{y,x} \hspace{1cm} \forall y \in G, \forall t \in T, \forall x \in X \end{aligned}\]

Where $h_{y,x}$ represents the heat rate slope for generator $y$ in segment $x$ [MMBTU/MWh], $f_{y,x}$ represents the heat rate intercept (MMBTU) for a generator $y$ in segment $x$ [MMBTU], and $U_{y,t}$ represents the commitment status of a generator $y$ at time $t$. These parameters are optional inputs to the resource .csv files. When Unit commitment is on, if a user provides slope and intercept, the standard heat rate (i.e., HeatRateMMBTUperMWh) will not be used. When unit commitment is off, the model will always use the standard heat rate. The user should determine the slope and intercept parameters based on the CapSize of the plant. For example, when a plant is operating at the full load (i.e., power output equal to the CapSize), the fuel usage determined by the effective segment divided by Cap_Size should be equal to the heat rate at full-load.

Since fuel consumption and fuel costs are postive, the optimization will force the fuel usage to be equal to the highest fuel usage segment for any given value of vP. When the power output is zero, the commitment variable $U_{g,t}$ will bring the intercept to be zero such that the fuel consumption is zero when thermal units are offline.

In order to run piecewise fuel consumption module, the unit commitment must be turned on (UC = 1 or 2), and users should provide PWFUSlope* and PWFUIntercept* for at least one segment.

To enable resources to use multiple fuels during both startup and normal operational processes, three additional variables were added: fuel $i$ consumption by plant $y$ at time $t$ ($vMulFuel_{y,i,t}$); startup fuel consumption for single-fuel plants ($vStartFuel_{y,t}$); and startup fuel consumption for multi-fuel plants ($vMulStartFuel_{y,i,t}$). By making startup fuel consumption variables, the model can choose the startup fuel to meet the constraints.

For plants using multiple fuels:

During startup, heat input from multiple startup fuels are equal to startup fuel requirements in plant $y$ at time $t$: $StartFuelMMBTUperMW$ $\times$ $Capsize$.

\[\begin{aligned} \sum_{i \in \mathcal{I} } vMulStartFuels_{y, i, t}= CapSize_{y} \times StartFuelMMBTUperMW_{y} \times vSTART_{y,t} \end{aligned}\]

During normal operation, the sum of fuel consumptions from multiple fuels dividing by the correspnding heat rates, respectively, is equal to $vPower$ in plant $y$ at time $t$.

\[\begin{aligned} \sum_{i \in \mathcal{I} } \frac{vMulFuels_{y, i, t}} {HeatRate_{i,y} } = vPower_{y,t} \end{aligned}\]

There are also constraints on how much heat input each fuel can provide, which are specified by $MinCofire$ and $MaxCofire$.

\[\begin{aligned} vMulFuels_{y, i, t} >= vPower_{y,t} \times MinCofire_{i} \end{aligned} \begin{aligned} vMulFuels_{y, i, t} <= vPower_{y,t} \times MaxCofire_{i} \end{aligned}\]