Create write_operating_reserve_price_revenue.jl

CHANGELOG.md

@@ -24,6 +24,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Maintenance formulation for thermal-commit plants (#556).
 - Add new tests for GenX: three-zone, multi-stage, electrolyzer, VRE+storage, 
   piecewise_fuel+CO2, and TDR (#563 and #578).
+- Added write_operating_reserve_price_revenue.jl to compute annual operating reserve and regulation revenue (PR # 611)
+- Added the operating reserve and regulation revenue to net revenue
 
 
 ### Fixed
@@ -38,8 +40,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Fix access of eELOSSByZone expr before initialization (#541)
 - Correctly write unmet reserves (in reg_dn.csv) (#575)
 - Correctly scale total reserves column (in reg_dn.csv) (#594)
+- Fixes issue #46
 - Add validation for `Reg_Max` and `Rsv_Max` columns in `Generators_data.csv` when `MUST_RUN` is set to 1 (#576)
 - Fix scaling of transmission losses in write_transmission_losses.jl (#621)
+- Fixes issue #46 (#611)
 
 ### Changed
 - Use add_to_expression! instead of the += and -= operators for memory performance improvements (#498).

docs/src/write_outputs.md

@@ -92,6 +92,11 @@ Pages = ["write_energy_revenue.jl"]
 Modules = [GenX]
 Pages = ["write_subsidy_revenue.jl"]
 ```
+## Write Operating Reserve and Regulation Revenue
+```@autodocs
+Modules = [GenX]
+Pages = ["write_operating_reserve_price_revenue.jl"]
+```
 
 ## Write Capacity Revenue
 ```@autodocs

src/load_inputs/load_cap_reserve_margin.jl

@@ -19,6 +19,7 @@ function load_cap_reserve_margin!(setup::Dict, path::AbstractString, inputs::Dic
     mat = extract_matrix_from_dataframe(df, "CapRes")
     inputs["dfCapRes"] = mat
     inputs["NCapacityReserveMargin"] = size(mat, 2)
+    print("Number of CRM types is ",inputs["NCapacityReserveMargin"])
 
     println(filename * " Successfully Read!")
 end

src/write_outputs/reserves/write_operating_reserve_price_revenue.jl

@@ -0,0 +1,68 @@
+@doc raw"""
+	write_operating_reserve_price_revenue(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
+
+Function for reporting the operating reserve and regulation revenue earned by generators listed in the input file.
+    GenX will print this file only when operating reserve and regulation are modeled and the shadow price can be obtained from the solver.
+    The revenues are calculated as the operating reserve and regulation contributions in each time step multiplied by the corresponding shadow price, and then the sum is taken over all modeled time steps.
+    The last column is the total revenue received from all operating reserve and regulation constraints.
+    As a reminder, GenX models the operating reserve and regulation at the time-dependent level, and each constraint either stands for an overall market or a locality constraint.
+"""
+function write_operating_reserve_regulation_revenue(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
+  	scale_factor = setup["ParameterScale"] == 1 ? ModelScalingFactor : 1
+	dfGen = inputs["dfGen"]
+	RSV = inputs["RSV"]
+	REG = inputs["REG"]
+
+	dfOpRsvRevenue = DataFrame(Region = dfGen[RSV, :region], Resource = dfGen[RSV, :Resource], Zone = dfGen[RSV, :Zone], Cluster = dfGen[RSV, :cluster], AnnualSum = Array{Float64}(undef, length(RSV)),)
+	dfOpRegRevenue = DataFrame(Region = dfGen[REG, :region], Resource = dfGen[REG, :Resource], Zone = dfGen[REG, :Zone], Cluster = dfGen[REG, :cluster], AnnualSum = Array{Float64}(undef, length(REG)),)
+
+	weighted_reg_price = operating_regulation_price(EP, inputs, setup)
+	weighted_rsv_price = operating_reserve_price(EP, inputs, setup)
+
+	rsvrevenue = value.(EP[:vRSV][RSV, :].data).* transpose(weighted_rsv_price)
+	regrevenue = value.(EP[:vREG][REG, :].data) .* transpose(weighted_reg_price)
+
+
+	rsvrevenue *= scale_factor
+	regrevenue *= scale_factor
+
+	dfOpRsvRevenue.AnnualSum .= rsvrevenue * inputs["omega"]
+	dfOpRegRevenue.AnnualSum .= regrevenue * inputs["omega"]
+	write_simple_csv(joinpath(path, "OperatingReserveRevenue.csv"), dfOpRsvRevenue)
+	write_simple_csv(joinpath(path, "OperatingRegulationRevenue.csv"), dfOpRegRevenue)
+	return dfOpRegRevenue, dfOpRsvRevenue
+end
+
+@doc raw"""
+    operating_regulation_price(EP::Model,
+                                  inputs::Dict,
+                                  setup::Dict)::Vector{Float64}
+
+Operating regulation price for each time step.
+This is equal to the dual variable of the regulation requirement constraint.
+
+    Returns a vector, with units of $/MW
+"""
+
+function operating_regulation_price(EP::Model, inputs::Dict, setup::Dict)::Vector{Float64}
+    ω = inputs["omega"]
+    scale_factor = setup["ParameterScale"] == 1 ? ModelScalingFactor : 1
+    return dual.(EP[:cReg]) ./ ω * scale_factor
+end
+
+@doc raw"""
+    operating_reserve_price(EP::Model,
+                                  inputs::Dict,
+                                  setup::Dict)::Vector{Float64}
+
+Operating reserve price for each time step.
+This is equal to the dual variable of the reserve requirement constraint.
+
+    Returns a vector, with units of $/MW
+"""
+
+function operating_reserve_price(EP::Model, inputs::Dict, setup::Dict)::Vector{Float64}
+    ω = inputs["omega"]
+    scale_factor = setup["ParameterScale"] == 1 ? ModelScalingFactor : 1
+    return dual.(EP[:cRsvReq]) ./ ω * scale_factor
+end

src/write_outputs/write_energy_revenue.jl

@@ -11,8 +11,8 @@ function write_energy_revenue(path::AbstractString, inputs::Dict, setup::Dict, E
 	NONFLEX = setdiff(collect(1:G), FLEX)
 	dfEnergyRevenue = DataFrame(Region = dfGen.region, Resource = inputs["RESOURCES"], Zone = dfGen.Zone, Cluster = dfGen.cluster, AnnualSum = Array{Float64}(undef, G),)
 	energyrevenue = zeros(G, T)
-    price = locational_marginal_price(EP, inputs, setup)
-    energyrevenue[NONFLEX, :] = value.(EP[:vP][NONFLEX, :]) .* transpose(price)[dfGen[NONFLEX, :Zone], :]
+	price = locational_marginal_price(EP, inputs, setup)
+	energyrevenue[NONFLEX, :] = value.(EP[:vP][NONFLEX, :]) .* transpose(price)[dfGen[NONFLEX, :Zone], :]
 	if !isempty(FLEX)
 		energyrevenue[FLEX, :] = value.(EP[:vCHARGE_FLEX][FLEX, :]).data .* transpose(price)[dfGen[FLEX, :Zone], :]
 	end

src/write_outputs/write_net_revenue.jl

@@ -1,13 +1,15 @@
 @doc raw"""
-	write_net_revenue(path::AbstractString, inputs::Dict, setup::Dict, EP::Model, dfCap::DataFrame, dfESRRev::DataFrame, dfResRevenue::DataFrame, dfChargingcost::DataFrame, dfPower::DataFrame, dfEnergyRevenue::DataFrame, dfSubRevenue::DataFrame, dfRegSubRevenue::DataFrame, dfVreStor::DataFrame)
+	write_net_revenue(path::AbstractString, inputs::Dict, setup::Dict, EP::Model, dfCap::DataFrame, dfESRRev::DataFrame, dfResRevenue::DataFrame, dfChargingcost::DataFrame, dfPower::DataFrame, dfEnergyRevenue::DataFrame, dfSubRevenue::DataFrame, dfRegSubRevenue::DataFrame, dfVreStor::DataFrame, dfOpRegRevenue::DataFrame, dfOpRsvRevenue::DataFrame)
 
 Function for writing net revenue of different generation technologies.
 """
-function write_net_revenue(path::AbstractString, inputs::Dict, setup::Dict, EP::Model, dfCap::DataFrame, dfESRRev::DataFrame, dfResRevenue::DataFrame, dfChargingcost::DataFrame, dfPower::DataFrame, dfEnergyRevenue::DataFrame, dfSubRevenue::DataFrame, dfRegSubRevenue::DataFrame, dfVreStor::DataFrame)
+function write_net_revenue(path::AbstractString, inputs::Dict, setup::Dict, EP::Model, dfCap::DataFrame, dfESRRev::DataFrame, dfResRevenue::DataFrame, dfChargingcost::DataFrame, dfPower::DataFrame, dfEnergyRevenue::DataFrame, dfSubRevenue::DataFrame, dfRegSubRevenue::DataFrame, dfVreStor::DataFrame, dfOpRegRevenue::DataFrame, dfOpRsvRevenue::DataFrame)
 	dfGen = inputs["dfGen"]
 	T = inputs["T"]     			# Number of time steps (hours)
 	Z = inputs["Z"]     			# Number of zones
-	G = inputs["G"]     			# Number of generators
+	G = inputs["G"]
+	RSV = inputs["RSV"]
+	REG = inputs["REG"]     			# Number of generators
 	COMMIT = inputs["COMMIT"]		# Thermal units for unit commitment
 	STOR_ALL = inputs["STOR_ALL"]
 	VRE_STOR = inputs["VRE_STOR"]
@@ -53,7 +55,7 @@ function write_net_revenue(path::AbstractString, inputs::Dict, setup::Dict, EP::
 	# Add operations and maintenance cost to the dataframe
 	dfNetRevenue.Fixed_OM_cost_MW = dfGen[!,:Fixed_OM_Cost_per_MWyr] .* dfCap[1:G,:EndCap]
  	dfNetRevenue.Fixed_OM_cost_MWh = dfGen[!,:Fixed_OM_Cost_per_MWhyr] .* dfCap[1:G,:EndEnergyCap]
-    dfNetRevenue.Fixed_OM_cost_charge_MW = dfGen[!, :Fixed_OM_Cost_Charge_per_MWyr] .* dfCap[1:G, :EndChargeCap]
+    	dfNetRevenue.Fixed_OM_cost_charge_MW = dfGen[!, :Fixed_OM_Cost_Charge_per_MWyr] .* dfCap[1:G, :EndChargeCap]
 
  	dfNetRevenue.Var_OM_cost_out = (dfGen[!,:Var_OM_Cost_per_MWh]) .* dfPower[1:G,:AnnualSum]
 	if !isempty(VRE_STOR)
@@ -117,33 +119,41 @@ function write_net_revenue(path::AbstractString, inputs::Dict, setup::Dict, EP::
 	end
 	# Add charge cost to the dataframe
 	dfNetRevenue.Charge_cost = zeros(nrow(dfNetRevenue))
-	if has_duals(EP) == 1
+	if has_duals(EP)
 		dfNetRevenue.Charge_cost = dfChargingcost[1:G,:AnnualSum] # Unit is confirmed to be US$
 	end
 
 	# Add energy and subsidy revenue to the dataframe
 	dfNetRevenue.EnergyRevenue = zeros(nrow(dfNetRevenue))
 	dfNetRevenue.SubsidyRevenue = zeros(nrow(dfNetRevenue))
-	if has_duals(EP) == 1
+	if has_duals(EP)
 		dfNetRevenue.EnergyRevenue = dfEnergyRevenue[1:G,:AnnualSum] # Unit is confirmed to be US$
 	 	dfNetRevenue.SubsidyRevenue = dfSubRevenue[1:G,:SubsidyRevenue] # Unit is confirmed to be US$
 	end
 
+	# Add energy and subsidy revenue to the dataframe
+	dfNetRevenue.OperatingReserveRevenue = zeros(nrow(dfNetRevenue))
+	dfNetRevenue.OperatingRegulationRevenue = zeros(nrow(dfNetRevenue))
+	if setup["Reserves"] > 0 && has_duals(EP)
+		dfNetRevenue.OperatingReserveRevenue[RSV] = dfOpRsvRevenue.AnnualSum # Unit is confirmed to be US$
+	 	dfNetRevenue.OperatingRegulationRevenue[REG] = dfOpRegRevenue.AnnualSum # Unit is confirmed to be US$
+	end
+
 	# Add capacity revenue to the dataframe
 	dfNetRevenue.ReserveMarginRevenue = zeros(nrow(dfNetRevenue))
- 	if setup["CapacityReserveMargin"] > 0 && has_duals(EP) == 1 # The unit is confirmed to be $
+ 	if setup["CapacityReserveMargin"] > 0 && has_duals(EP) # The unit is confirmed to be $
  		dfNetRevenue.ReserveMarginRevenue = dfResRevenue[1:G,:AnnualSum]
  	end
 
 	# Add RPS/CES revenue to the dataframe
 	dfNetRevenue.ESRRevenue = zeros(nrow(dfNetRevenue))
- 	if setup["EnergyShareRequirement"] > 0 && has_duals(EP) == 1 # The unit is confirmed to be $
+ 	if setup["EnergyShareRequirement"] > 0 && has_duals(EP) # The unit is confirmed to be $
  		dfNetRevenue.ESRRevenue = dfESRRev[1:G,:Total]
  	end
 
 	# Calculate emissions cost
 	dfNetRevenue.EmissionsCost = zeros(nrow(dfNetRevenue))
-	if setup["CO2Cap"] >=1 && has_duals(EP) == 1
+	if setup["CO2Cap"] >=1 && has_duals(EP)
 		for cap in 1:inputs["NCO2Cap"]
 			co2_cap_dual = dual(EP[:cCO2Emissions_systemwide][cap])
 			CO2ZONES = findall(x->x==1, inputs["dfCO2CapZones"][:,cap])
@@ -166,11 +176,11 @@ function write_net_revenue(path::AbstractString, inputs::Dict, setup::Dict, EP::
 
 	# Add regional technology subsidy revenue to the dataframe
 	dfNetRevenue.RegSubsidyRevenue = zeros(nrow(dfNetRevenue))
-	if setup["MinCapReq"] >= 1 && has_duals(EP) == 1 # The unit is confirmed to be US$
+	if setup["MinCapReq"] >= 1 && has_duals(EP)# The unit is confirmed to be US$
 		dfNetRevenue.RegSubsidyRevenue = dfRegSubRevenue[1:G,:SubsidyRevenue]
 	end
 
-	dfNetRevenue.Revenue = dfNetRevenue.EnergyRevenue .+ dfNetRevenue.SubsidyRevenue .+ dfNetRevenue.ReserveMarginRevenue .+ dfNetRevenue.ESRRevenue .+ dfNetRevenue.RegSubsidyRevenue
+	dfNetRevenue.Revenue = dfNetRevenue.EnergyRevenue .+ dfNetRevenue.SubsidyRevenue .+ dfNetRevenue.ReserveMarginRevenue .+ dfNetRevenue.ESRRevenue .+ dfNetRevenue.RegSubsidyRevenue .+ dfNetRevenue.OperatingReserveRevenue .+ dfNetRevenue.OperatingRegulationRevenue
 	dfNetRevenue.Cost = (dfNetRevenue.Inv_cost_MW
                       .+ dfNetRevenue.Inv_cost_MWh
                       .+ dfNetRevenue.Inv_cost_charge_MW

src/write_outputs/write_outputs.jl

@@ -155,17 +155,17 @@ function write_outputs(EP::Model, path::AbstractString, setup::Dict, inputs::Dic
 		end
 
 		elapsed_time_time_weights = @elapsed write_time_weights(path, inputs)
-	  println("Time elapsed for writing time weights is")
-	  println(elapsed_time_time_weights)
+	  	println("Time elapsed for writing time weights is")
+	  	println(elapsed_time_time_weights)
 		dfESR = DataFrame()
 		dfESRRev = DataFrame()
-		if setup["EnergyShareRequirement"]==1 && has_duals(EP) == 1
+		if setup["EnergyShareRequirement"]==1 && has_duals(EP)
 			dfESR = write_esr_prices(path, inputs, setup, EP)
 			dfESRRev = write_esr_revenue(path, inputs, setup, dfPower, dfESR, EP)
 		end
 		dfResMar = DataFrame()
 		dfResRevenue = DataFrame()
-		if setup["CapacityReserveMargin"]==1 && has_duals(EP) == 1
+		if setup["CapacityReserveMargin"]==1 && has_duals(EP)
 			dfResMar = write_reserve_margin(path, setup, EP)
 			elapsed_time_rsv_margin = @elapsed write_reserve_margin_w(path, inputs, setup, EP)
 			dfVirtualDischarge = write_virtual_discharge(path, inputs, setup, EP)
@@ -179,26 +179,34 @@ function write_outputs(EP::Model, path::AbstractString, setup::Dict, inputs::Dic
 				dfResMar_slack = write_reserve_margin_slack(path, inputs, setup, EP)
 			end		  
 		end
-		if setup["CO2Cap"]>0 && has_duals(EP) == 1
+
+		if setup["Reserves"]==1 && has_duals(EP)
+			dfOpRegRevenue, dfOpRsvRevenue = write_operating_reserve_regulation_revenue(path, inputs, setup, EP)
+			elapsed_time_op_res_rev = @elapsed write_operating_reserve_regulation_revenue(path, inputs, setup, EP)
+		  	println("Time elapsed for writing oerating reserve and regulation revenue is")
+		  	println(elapsed_time_op_res_rev)
+		end
+
+		if setup["CO2Cap"]>0 && has_duals(EP)
 			dfCO2Cap = write_co2_cap(path, inputs, setup, EP)
 		end
-		if setup["MinCapReq"] == 1 && has_duals(EP) == 1
+		if setup["MinCapReq"] == 1 && has_duals(EP)
 			dfMinCapReq = write_minimum_capacity_requirement(path, inputs, setup, EP)
 		end
 
-		if setup["MaxCapReq"] == 1 && has_duals(EP) == 1
+		if setup["MaxCapReq"] == 1 && has_duals(EP)
 			dfMaxCapReq = write_maximum_capacity_requirement(path, inputs, setup, EP)
 		end
 
-		if !isempty(inputs["ELECTROLYZER"]) && has_duals(EP) == 1
+		if !isempty(inputs["ELECTROLYZER"]) && has_duals(EP)
 			dfHydrogenPrice = write_hydrogen_prices(path, inputs, setup, EP)
 			if setup["HydrogenHourlyMatching"] == 1
 				dfHourlyMatchingPrices = write_hourly_matching_prices(path, inputs, setup, EP)
 			end
 		end
 
 
-		elapsed_time_net_rev = @elapsed write_net_revenue(path, inputs, setup, EP, dfCap, dfESRRev, dfResRevenue, dfChargingcost, dfPower, dfEnergyRevenue, dfSubRevenue, dfRegSubRevenue, dfVreStor)
+		elapsed_time_net_rev = @elapsed write_net_revenue(path, inputs, setup, EP, dfCap, dfESRRev, dfResRevenue, dfChargingcost, dfPower, dfEnergyRevenue, dfSubRevenue, dfRegSubRevenue, dfVreStor, dfOpRegRevenue, dfOpRsvRevenue)
 	  	println("Time elapsed for writing net revenue is")
 	  	println(elapsed_time_net_rev)
 	end