CO2 module and fuel module

src/configure_settings/configure_settings.jl

@@ -24,6 +24,7 @@ function default_settings()
         "IncludeLossesInESR" => 0,
         "EnableJuMPStringNames" => false,
         "HydrogenHourlyMatching" => 0,
+        "PiecewiseHeatRate" => 0
     )
 end
 

src/load_inputs/load_fuels_data.jl

@@ -1,6 +1,5 @@
 @doc raw"""
-    load_fuels_data!(setup::Dict, path::AbstractString, inputs::Dict)
-
+load_fuels_data!(setup::Dict, path::AbstractString, inputs::Dict)
 Read input parameters related to fuel costs and CO$_2$ content of fuels
 """
 function load_fuels_data!(setup::Dict, path::AbstractString, inputs::Dict)
@@ -32,9 +31,11 @@ function load_fuels_data!(setup::Dict, path::AbstractString, inputs::Dict)
     scale_factor = setup["ParameterScale"] == 1 ? ModelScalingFactor : 1
 
     for i = 1:length(fuels)
+            # fuel cost is in $/MMBTU w/o scaling, $/Billon BTU w/ scaling
             fuel_costs[fuels[i]] = costs[:,i] / scale_factor
-            # fuel_CO2 is kton/MMBTU with scaling, or ton/MMBTU without scaling.
-            fuel_CO2[fuels[i]] = CO2_content[i] / scale_factor
+            # fuel_CO2 is ton/MMBTU w/o scaling, or kton/Billion BTU w/ scaling.
+            # No need to scale fuel_CO2
+            fuel_CO2[fuels[i]] = CO2_content[i] 
     end
 
     inputs["fuels"] = fuels

src/load_inputs/load_generators_data.jl

@@ -185,44 +185,49 @@ function load_generators_data!(setup::Dict, path::AbstractString, inputs_gen::Di
         end
     end
 
-	if setup["UCommit"]>=1
-		# Fuel consumed on start-up (million BTUs per MW per start) if unit commitment is modelled
-		start_fuel = convert(Array{Float64}, gen_in[!,:Start_Fuel_MMBTU_per_MW])
-		# Fixed cost per start-up ($ per MW per start) if unit commitment is modelled
+	if setup["UCommit"] >= 1
 		start_cost = convert(Array{Float64}, gen_in[!,:Start_Cost_per_MW])
 		inputs_gen["C_Start"] = zeros(Float64, G, inputs_gen["T"])
-		gen_in[!,:CO2_per_Start] = zeros(Float64, G)
 	end
 
-	# Heat rate of all resources (million BTUs/MWh)
-	heat_rate = convert(Array{Float64}, gen_in[!,:Heat_Rate_MMBTU_per_MWh])
-	# Fuel used by each resource
-	fuel_type = gen_in[!,:Fuel]
-	# Maximum fuel cost in $ per MWh and CO2 emissions in tons per MWh
-	inputs_gen["C_Fuel_per_MWh"] = zeros(Float64, G, inputs_gen["T"])
-	gen_in[!,:CO2_per_MWh] = zeros(Float64, G)
+	# The fuel (including start up fuel) costs and CO2 emissions will be separately tracked in Fuel.jl and CO2.jl
+	# No need to determine C_Fuel_per_MWh and CO2_per_MWh, and C_Start
+	# Only scale the start up cost here
 	for g in 1:G
-		# NOTE: When Setup[ParameterScale] =1, fuel costs are scaled in fuels_data.csv, so no if condition needed to scale C_Fuel_per_MWh
-		inputs_gen["C_Fuel_per_MWh"][g,:] = fuel_costs[fuel_type[g]].*heat_rate[g]
-		gen_in[g,:CO2_per_MWh] = fuel_CO2[fuel_type[g]]*heat_rate[g]
-		gen_in[g,:CO2_per_MWh] *= scale_factor
-		# kton/MMBTU * MMBTU/MWh = kton/MWh, to get kton/GWh, we need to mutiply 1000
 		if g in inputs_gen["COMMIT"]
 			# Start-up cost is sum of fixed cost per start plus cost of fuel consumed on startup.
-			# CO2 from fuel consumption during startup also calculated
-
-			inputs_gen["C_Start"][g,:] = gen_in[g,:Cap_Size] * (fuel_costs[fuel_type[g]] .* start_fuel[g] .+ start_cost[g])
-			# No need to re-scale C_Start since Cap_size, fuel_costs and start_cost are scaled When Setup[ParameterScale] =1 - Dharik
-			gen_in[g,:CO2_per_Start]  = gen_in[g,:Cap_Size]*(fuel_CO2[fuel_type[g]]*start_fuel[g])
-			gen_in[g,:CO2_per_Start] *= scale_factor
-			# Setup[ParameterScale] =1, gen_in[g,:Cap_Size] is GW, fuel_CO2[fuel_type[g]] is ktons/MMBTU, start_fuel is MMBTU/MW,
-			#   thus the overall is MTons/GW, and thus gen_in[g,:CO2_per_Start] is Mton, to get kton, change we need to multiply 1000
-			# Setup[ParameterScale] =0, gen_in[g,:Cap_Size] is MW, fuel_CO2[fuel_type[g]] is tons/MMBTU, start_fuel is MMBTU/MW,
-			#   thus the overall is MTons/GW, and thus gen_in[g,:CO2_per_Start] is ton
+			inputs_gen["C_Start"][g,:] .= gen_in[g,:Cap_Size] * ( start_cost[g])
 		end
 	end
 
 	load_vre_stor_data!(inputs_gen, setup, path)
+
+
+    # write zeros if col names are not in the gen_in dataframe
+	missing_cols = ["Incremental_Heat_Rate_Segment2", "Incremental_Heat_Rate_Segment3", "Intercept_Fuel_Consumption_Segment1","Intercept_Fuel_Consumption_Segment2", "Intercept_Fuel_Consumption_Segment3",
+	"Biomass", "CO2_Capture_Rate", "CO2_Capture_Cost_per_Metric_Ton"]
+	write_zeros_if_not_exist!(gen_in, missing_cols)
+
+
+    # Scale CO2_Capture_Cost_per_Metric_Ton for CCS units 
+	gen_in.CO2_Capture_Cost_per_Metric_Ton /= scale_factor
+
+    # Scale Intercept of fuel consumption of segments
+	# Users should at least provide Incremental_Heat_Rate_Segment1 and Intercept_Fuel_Consumption_Segment1 
+	# if Users didn't provide data for but turn on piecewiseheatrate, we will set the Incremental_Heat_Rate_Segment to be the same as conventional heat rate
+	if setup["UCommit"] > 0
+        if setup["PiecewiseHeatRate"] == 1
+			gen_in.Incremental_Heat_Rate_Segment1 = ("Incremental_Heat_Rate_Segment1" in names(gen_in)) ? gen_in.Incremental_Heat_Rate_Segment1 : gen_in.Heat_Rate_MMBTU_per_MWh
+
+            # no need to scale incremental heat rate, but the intercept of fuel consumption in each segment needs to be scaled.
+			gen_in.Intercept_Fuel_Consumption_Segment1 /= scale_factor
+			gen_in.Intercept_Fuel_Consumption_Segment2 /= scale_factor
+			gen_in.Intercept_Fuel_Consumption_Segment3 /= scale_factor
+		end
+	end
+
+
+
 	println(filename * " Successfully Read!")
 end
 
@@ -501,4 +506,15 @@ function load_vre_stor_data!(inputs_gen::Dict, setup::Dict, path::AbstractString
 		inputs_gen["dfVRE_STOR"] = DataFrame()
 	end
 	summarize_errors(error_strings)
-end
+end
+
+
+function write_zeros_if_not_exist!(dfGen::DataFrame, col_names:: Vector{String})
+	for col_name in col_names
+		if !(col_name in names(dfGen))
+			dfGen[!, col_name] = zeros(Int, nrow(dfGen))
+		end
+	end
+	return dfGen
+end
+

src/model/core/co2.jl

@@ -0,0 +1,107 @@
+@doc raw""" 
+
+co2!(EP::Model, inputs::Dict)
+
+This function creates expression to account for CO2 emissions and 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. This module will displace the emissions module.
+
+***** Expressions *****
+
+For thermal generators that use fuels that contain CO2 content (e.g., coal, natural gas, and biomass), the CO2 emissions are a function of fuel consumption, CO2 capture rate, and whether the feedstock is biomass. Biomass (e.g., wastes or agriculture resides) derived energy is typically considered to be carbon-neutral because the carbon in the biomass is originated from the atmosphere. When bioenergy is coupled with carbon capture and storage (CCS), it creates negative emissions.
+
+Here we create a column called Biomass in the Generator data file (1 or 0), which determines if a generator $y$ uses biomass or not. The CO2 emissions from a generator should be zero without CCS and negative with CCS.
+
+The CO2 emissions from the generator $y$ at time $t$, denoted by $eEmissionsByPlant_{y,t}$, is determined by total fuel consumption (MMBTU, including startup fuel) multiplied by the CO2 content of the fuel (t CO2/MMBTU), then times (1 - Biomass - CO2 capture rate). In short, the CO2 emissions depend on total CO2 content from fuel consumption, the CO2 capture rate, and whether the generators use biomass.
+
+```math
+\begin{aligned}
+eEmissionsByPlant_{g,t} = (1-Biomass_y-  CO2\_Capture\_Rate_y) * (vFuel_{y,t} + 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 that determines if the generator $y$ uses biomass (Biomass = 1) or not (Biomass = 0), $CO2\_Capture\_Rate_y$ represents a fraction (between 0 - 1) for CO2 capture rate.
+
+In addition to CO2 emissions, for generators with non-zero CO2 capture rate, we also determine the amount of CO2 being captured and sequestrated. The CO2 emissions from the generator $y$ at time $t$, denoted by $eEmissionsCaptureByPlant_{g,t}$, is determined by total fuel consumption (MMBTU, including startup fuel) multiplied by the $CO_2$ content of the fuel (t CO2/MMBTU), then times CO2 capture rate.
+
+```math
+\begin{aligned}
+eEmissionsCaptureByPlant_{g,t} = CO2\_Capture\_Rate_y * (vFuel_{y,t} + eStartFuel_{y,t}) * CO2_{content}
+\hspace{1cm} \forall y \in G, \forall t \in T
+\end{aligned}
+```
+
+
+
+"""
+function co2!(EP::Model, inputs::Dict)
+
+    println("C02 Module")
+
+    dfGen = inputs["dfGen"]
+    G = inputs["G"]     # Number of resources (generators, storage, DR, and DERs)
+    T = inputs["T"]     # Number of time steps (hours)
+    Z = inputs["Z"]     # Number of zones
+
+    dfGen.Biomass = "Biomass" in names(dfGen) ? dfGen.Biomass : zeros(Int, nrow(dfGen))
+    dfGen.CO2_Capture_Rate = "CO2_Capture_Rate" in names(dfGen) ? dfGen.CO2_Capture_Rate : zeros(Int, nrow(dfGen))
+
+    ### Expressions ###
+    # CO2 emissions from power plants in "Generator_data.csv"
+    # if all the CO2 capture rates from generator data are zeros, the CO2 emissions from thermal generators are determined by fuel consumptiono times CO2 content per MMBTU 
+    if all(x -> x == 0, dfGen.CO2_Capture_Rate)
+        @expression(EP, eEmissionsByPlant[y=1:G, t=1:T], 
+            ((1-dfGen[y, :Biomass]) *(EP[:vFuel][y, t] + EP[:eStartFuel][y, t]) * inputs["fuel_CO2"][dfGen[y,:Fuel]]))
+    else 
+        @expression(EP, eEmissionsByPlant[y=1:G, t=1:T],
+            (1-dfGen[y, :Biomass] - dfGen[y, :CO2_Capture_Rate]) * 
+            ((EP[:vFuel][y, t] + EP[:eStartFuel][y, t]) * 
+                inputs["fuel_CO2"][dfGen[y,:Fuel]]))
+
+        # CO2  captured from power plants in "Generator_data.csv"
+        @expression(EP, eEmissionsCaptureByPlant[y=1:G, t=1:T],
+            (dfGen[y, :CO2_Capture_Rate]) * 
+            ((EP[:vFuel][y, t] + EP[:eStartFuel][y, t]) * 
+                inputs["fuel_CO2"][dfGen[y,:Fuel]]))
+
+
+        #************************************* 
+        @expression(EP, eEmissionsCaptureByPlantYear[y=1:G], 
+            sum(inputs["omega"][t] * eEmissionsCaptureByPlant[y, t] 
+                for t in 1:T))
+        @expression(EP, eEmissionsCaptureByZone[z=1:Z, t=1:T], 
+            sum(eEmissionsCaptureByPlant[y, t] 
+                for y in dfGen[(dfGen[!, :Zone].==z), :R_ID]))
+        @expression(EP, eEmissionsCaptureByZoneYear[z=1:Z], 
+            sum(eEmissionsCaptureByPlantYear[y] 
+                for y in dfGen[(dfGen[!, :Zone].==z), :R_ID]))
+        #*************************************
+
+        # add CO2 sequestration cost to objective function
+        # when scale factor is on tCO2/MWh = > kt CO2/GWh
+        @expression(EP, ePlantCCO2Sequestration[y=1:G], 
+            sum(inputs["omega"][t] * eEmissionsCaptureByPlant[y, t] * 
+                dfGen[y, :CO2_Capture_Cost_per_Metric_Ton] for t in 1:T))
+
+        @expression(EP, eZonalCCO2Sequestration[z=1:Z], 
+            sum(ePlantCCO2Sequestration[y] 
+                for y in dfGen[(dfGen[!, :Zone].==z), :R_ID]))
+
+        @expression(EP, eTotaleCCO2Sequestration, 
+            sum(eZonalCCO2Sequestration[z] for z in 1:Z))
+
+        add_to_expression!(EP[:eObj], EP[:eTotaleCCO2Sequestration])
+    end
+
+    @expression(EP, eEmissionsByPlantYear[y = 1:G], 
+        sum(inputs["omega"][t] * eEmissionsByPlant[y, t] for t in 1:T))
+
+    @expression(EP, eEmissionsByZone[z = 1:Z, t = 1:T], 
+        sum(eEmissionsByPlant[y, t] for y in dfGen[(dfGen[!, :Zone].==z), :R_ID]))
+
+    @expression(EP, eEmissionsByZoneYear[z = 1:Z], 
+        sum(inputs["omega"][t] * eEmissionsByZone[z, t] for t in 1:T))
+
+
+    return EP
+
+end

src/model/core/discharge/discharge.jl

@@ -27,9 +27,8 @@ function discharge!(EP::Model, inputs::Dict, setup::Dict)
 
 	## Objective Function Expressions ##
 
-	# Variable costs of "generation" for resource "y" during hour "t" = variable O&M plus fuel cost
-	@expression(EP, eCVar_out[y=1:G,t=1:T], (inputs["omega"][t]*(dfGen[y,:Var_OM_Cost_per_MWh]+inputs["C_Fuel_per_MWh"][y,t])*vP[y,t]))
-	#@expression(EP, eCVar_out[y=1:G,t=1:T], (round(inputs["omega"][t]*(dfGen[y,:Var_OM_Cost_per_MWh]+inputs["C_Fuel_per_MWh"][y,t]), digits=RD)*vP[y,t]))
+	# Variable costs of "generation" for resource "y" during hour "t" = variable O&M, the fuel costs will be determined in fuel.jl
+	@expression(EP, eCVar_out[y=1:G,t=1:T], (inputs["omega"][t]*(dfGen[y,:Var_OM_Cost_per_MWh]*vP[y,t])))
 	# Sum individual resource contributions to variable discharging costs to get total variable discharging costs
 	@expression(EP, eTotalCVarOutT[t=1:T], sum(eCVar_out[y,t] for y in 1:G))
 	@expression(EP, eTotalCVarOut, sum(eTotalCVarOutT[t] for t in 1:T))

src/model/core/fuel.jl

@@ -0,0 +1,117 @@
+
+@doc raw"""
+fuel!(EP::Model, inputs::Dict, setup::Dict)
+
+This function creates expression to account for total fuel consumption (e.g., coal, natural gas, hydrogen, etc). It also has the capability to model the piece-wise fuel consumption in part load (if data is available)
+
+***** Expressions ******
+
+The fuel consumption for power generation $vFuel_{y,t}$ is determined by power generation ($vP_{y,t}$) mutiplied by the corresponding heat rate ($Hear\_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$, is determined by fuel consumption ($vFuel_{y,t}$ and $eStartFuel$) multiplied by the fuel costs (\$/MMBTU)
+
+From above formulations, thermal generators are expected to have the same fuel consumption per generating 1 MWh electricity, regardless of minimum load or full load. However, thermal generators tend to have decreased efficiency when operating at part load, leading to higher fuel consumption per generating the same amount of electricity. To have more precise representation of fuel consumption at part load, the piecewise-linear fitting of heat input can be introduced. 
+
+```math
+\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 incremental heat rate of a thermal generator $y$ in segment $x$ [MMBTU/MWh] and $f_{y,x}$ represents intercept of fuel consumption of a thermal generator $y$ in segment $x$ [MMBUT], and $U_{y,t}$ represents the commit status of a thermal generator $y$ at time $t$. We include at most three segments to represent the piecewise heat consumption. 
+
+Since fuel consumption has a positive value, the optimization will optimize the fuel consumption by enforcing the inequity to equal to the highest piecewise segment. 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, and users should provide Incremental_Heat_Rate_Segmenti and Intercept_Fuel_Consumption_Segmenti for at least one segment. 
+
+If users only want to model piecewise heat rate for some of thermal generators, then set the Incremental_Heat_Rate_Segment1 for thoese plants to be the same as conventional heat rate, and then set Incremental_Heat_Rate_Segment2/3 and Intercept_Fuel_Consumption_Segment1/2/3 to be zero.
+"""
+
+function fuel!(EP::Model, inputs::Dict, setup::Dict)
+    println("Fuel Module")
+    dfGen = inputs["dfGen"]
+    T = inputs["T"]     # Number of time steps (hours)
+    Z = inputs["Z"]     # Number of zones
+    G = inputs["G"]
+    THERM_COMMIT = inputs["THERM_COMMIT"]
+    FUEL = length(inputs["fuels"])
+    ALLGEN = collect(1:G)
+    # create variable for fuel consumption for output
+    @variable(EP, vFuel[y in 1:G, t = 1:T] >= 0)
+
+    ### Expressions ####
+    # Fuel consumed on start-up (MMBTU or kMMBTU (scaled)) 
+    # if unit commitment is modelled
+    @expression(EP, eStartFuel[y in 1:G, t = 1:T],
+        if y in THERM_COMMIT
+            (dfGen[y,:Cap_Size] * EP[:vSTART][y, t] * 
+                dfGen[y,:Start_Fuel_MMBTU_per_MW])
+        else
+            EP[:vZERO]
+        end)
+
+    # fuel_cost is in $/MMBTU (M$/billion BTU if scaled)
+    # vFuel and eStartFuel is MMBTU (or billion BTU if scaled)
+    # Therefore eCFuel_start or eCFuel_out is $ or Million$)
+    # Separately track the start up fuel and fuel consumption for power generation
+
+    # Start up fuel cost
+    @expression(EP, eCFuelStart[y = 1:G, t = 1:T], 
+        (inputs["fuel_costs"][dfGen[y,:Fuel]][t] * EP[:eStartFuel][y, t]))
+    # plant level start-up fuel cost for output
+    @expression(EP, ePlantCFuelStart[y = 1:G], 
+        sum(inputs["omega"][t] * EP[:eCFuelStart][y, t] for t in 1:T))
+    # zonal level total fuel cost for output
+    @expression(EP, eZonalCFuelStart[z = 1:Z], 
+        sum(EP[:ePlantCFuelStart][y] for y in dfGen[dfGen[!, :Zone].==z, :R_ID]))
+
+    # Fuel cost for power generation
+    @expression(EP, eCFuelOut[y = 1:G, t = 1:T], 
+        (inputs["fuel_costs"][dfGen[y,:Fuel]][t] * EP[:vFuel][y, t]))
+    # plant level start-up fuel cost for output
+    @expression(EP, ePlantCFuelOut[y = 1:G], 
+        sum(inputs["omega"][t] * EP[:eCFuelOut][y, t] for t in 1:T))
+    # zonal level total fuel cost for output
+    @expression(EP, eZonalCFuelOut[z = 1:Z], 
+        sum(EP[:ePlantCFuelOut][y] for y in dfGen[dfGen[!, :Zone].==z, :R_ID]))
+
+
+    # system level total fuel cost for output
+    @expression(EP, eTotalCFuelOut, sum(eZonalCFuelOut[z] for z in 1:Z))
+    @expression(EP, eTotalCFuelStart, sum(eZonalCFuelStart[z] for z in 1:Z))
+
+
+    add_to_expression!(EP[:eObj], EP[:eTotalCFuelOut] + EP[:eTotalCFuelStart])
+
+    #fuel consumption (MMBTU or Billion BTU)
+    @expression(EP, eFuelConsumption[f in 1:FUEL, t in 1:T],
+        sum(EP[:vFuel][y, t] + EP[:eStartFuel][y,t]
+            for y in dfGen[dfGen[!,:Fuel] .== string(inputs["fuels"][f]) ,:R_ID]))
+
+    @expression(EP, eFuelConsumptionYear[f in 1:FUEL],
+        sum(inputs["omega"][t] * EP[:eFuelConsumption][f, t] for t in 1:T))
+
+
+    ### Constraint ###
+    @constraint(EP, FuelCalculation[y in setdiff(ALLGEN, THERM_COMMIT), t = 1:T],
+        EP[:vFuel][y, t] - EP[:vP][y, t] * dfGen[y, :Heat_Rate_MMBTU_per_MWh] == 0)
+    if !isempty(THERM_COMMIT)
+        if setup["PiecewiseHeatRate"] == 1
+            # Piecewise heat rate UC
+            @constraint(EP, First_segment[y in THERM_COMMIT, t = 1:T],
+                EP[:vFuel][y, t] >= (EP[:vP][y, t] * dfGen[y, :Incremental_Heat_Rate_Segment1] + 
+                    EP[:vCOMMIT][y, t] * dfGen[y, :Intercept_Fuel_Consumption_Segment1]))
+            @constraint(EP, Second_segment[y in THERM_COMMIT, t = 1:T],
+                EP[:vFuel][y, t] >= (EP[:vP][y, t] * dfGen[y, :Incremental_Heat_Rate_Segment2] + 
+                    EP[:vCOMMIT][y, t] * dfGen[y, :Intercept_Fuel_Consumption_Segment2]))
+            @constraint(EP, Third_segment[y in THERM_COMMIT, t = 1:T],
+                EP[:vFuel][y, t] >= (EP[:vP][y, t] * dfGen[y, :Incremental_Heat_Rate_Segment3] + 
+                    EP[:vCOMMIT][y, t] * dfGen[!, :Intercept_Fuel_Consumption_Segment3][y]))
+        else
+            @constraint(EP, FuelCalculationCommit[y in THERM_COMMIT, t = 1:T],
+                EP[:vFuel][y, t] - EP[:vP][y, t] * dfGen[y, :Heat_Rate_MMBTU_per_MWh] == 0)
+        end
+    end
+
+    return EP
+end