Adding resistive heating and non-fusion dispatch and UC constraints

Example_Systems/SmallNewEngland/OneZone/Settings/highs_settings.yml

@@ -10,4 +10,4 @@ Method: choose             #HiGHS-specific solver settings # Solver option: "sim
 
 # Run the crossover routine for IPX
 # [type: bool, advanced: true, range: {false, true}, default: true]
-run_crossover: false
+run_crossover: off

Example_Systems/SmallNewEngland/Simple_Test_Case/Generators_data.csv

@@ -1,5 +1,5 @@
-Resource,Zone,THERM,MUST_RUN,STOR,FLEX,HYDRO,VRE,LDS,Num_VRE_Bins,New_Build,Existing_Cap_MW,Existing_Cap_MWh,Existing_Charge_Cap_MW,Max_Cap_MW,Max_Cap_MWh,Max_Charge_Cap_MW,Min_Cap_MW,Min_Cap_MWh,Min_Charge_Cap_MW,Inv_Cost_per_MWyr,Inv_Cost_per_MWhyr,Inv_Cost_Charge_per_MWyr,Fixed_OM_Cost_per_MWyr,Fixed_OM_Cost_per_MWhyr,Fixed_OM_Cost_Charge_per_MWyr,Var_OM_Cost_per_MWh,Var_OM_Cost_per_MWh_In,Heat_Rate_MMBTU_per_MWh,Fuel,Cap_Size,Start_Cost_per_MW,Start_Fuel_MMBTU_per_MW,Up_Time,Down_Time,Ramp_Up_Percentage,Ramp_Dn_Percentage,Hydro_Energy_to_Power_Ratio,Min_Power,Self_Disch,Eff_Up,Eff_Down,Min_Duration,Max_Duration,Max_Flexible_Demand_Advance,Max_Flexible_Demand_Delay,Flexible_Demand_Energy_Eff,Reg_Max,Rsv_Max,Reg_Cost,Rsv_Cost,MinCapTag,MinCapTag_1,MinCapTag_2,MinCapTag_3,MGA,Resource_Type,CapRes_1,ESR_1,ESR_2,region,cluster,LDS
-natural_gas_combined_cycle,1,1,0,0,0,0,0,0,0,1,0,0,0,-1,-1,-1,0,0,0,65400,0,0,10287,0,0,3.55,0,7.43,NG,250,91,2,6,6,0.64,0.64,0,0.468,0,1,1,0,0,0,0,1,0.25,0.5,0,0,2,0,0,0,1,natural_gas_fired_combined_cycle,0.93,0,0,NE,1,0
-solar_pv,1,0,0,0,0,0,1,0,1,1,0,0,0,-1,-1,-1,0,0,0,85300,0,0,18760,0,0,0,0,9.13,None,0,0,0,0,0,1,1,0,0,0,1,1,0,0,0,0,1,0,0,0,0,7,1,0,0,1,solar_photovoltaic,0.8,1,1,NE,1,0
-onshore_wind,1,0,0,0,0,0,1,0,1,1,0,0,0,-1,-1,-1,0,0,0,97200,0,0,43205,0,0,0.1,0,9.12,None,0,0,0,0,0,1,1,0,0,0,1,1,0,0,0,0,1,0,0,0,0,6,0,1,0,1,onshore_wind_turbine,0.8,1,1,NE,1,0
-battery,1,0,0,1,0,0,0,0,0,1,0,0,0,-1,-1,-1,0,0,0,19584,22494,0,4895,5622,0,0.15,0.15,0,None,0,0,0,0,0,1,1,0,0,0,0.92,0.92,1,10,0,0,1,0,0,0,0,12,0,0,1,0,battery_mid,0.95,0,0,NE,0,0
+Resource,Zone,THERM,MUST_RUN,STOR,TS,FLEX,HYDRO,VRE,LDS,Num_VRE_Bins,New_Build,Existing_Cap_MW,Existing_Cap_MWh,Existing_Charge_Cap_MW,Max_Cap_MW,Max_Cap_MWh,Max_Charge_Cap_MW,Min_Cap_MW,Min_Cap_MWh,Min_Charge_Cap_MW,Inv_Cost_per_MWyr,Inv_Cost_per_MWhyr,Inv_Cost_Charge_per_MWyr,Fixed_OM_Cost_per_MWyr,Fixed_OM_Cost_per_MWhyr,Fixed_OM_Cost_Charge_per_MWyr,Var_OM_Cost_per_MWh,Var_OM_Cost_per_MWh_In,Heat_Rate_MMBTU_per_MWh,Fuel,Cap_Size,Start_Cost_per_MW,Start_Fuel_MMBTU_per_MW,Up_Time,Down_Time,Ramp_Up_Percentage,Ramp_Dn_Percentage,Hydro_Energy_to_Power_Ratio,Min_Power,Self_Disch,Eff_Up,Eff_Down,Min_Duration,Max_Duration,Max_Flexible_Demand_Advance,Max_Flexible_Demand_Delay,Flexible_Demand_Energy_Eff,Reg_Max,Rsv_Max,Reg_Cost,Rsv_Cost,MinCapTag,MinCapTag_1,MinCapTag_2,MinCapTag_3,MGA,Resource_Type,CapRes_1,ESR_1,ESR_2,region,cluster,LDS
+natural_gas_combined_cycle,1,1,0,0,1,0,0,0,0,0,1,0,0,0,-1,-1,-1,0,0,0,65400,0,0,10287,0,0,3.55,0,7.43,NG,250,91,2,6,6,0.64,0.64,0,0.468,0,1,1,0,0,0,0,1,0.25,0.5,0,0,2,0,0,0,1,natural_gas_fired_combined_cycle,0.93,0,0,NE,1,0
+solar_pv,1,0,0,0,0,0,0,1,0,1,1,0,0,0,-1,-1,-1,0,0,0,85300,0,0,18760,0,0,0,0,9.13,None,0,0,0,0,0,1,1,0,0,0,1,1,0,0,0,0,1,0,0,0,0,7,1,0,0,1,solar_photovoltaic,0.8,1,1,NE,1,0
+onshore_wind,1,0,0,0,0,0,0,1,0,1,1,0,0,0,-1,-1,-1,0,0,0,97200,0,0,43205,0,0,0.1,0,9.12,None,0,0,0,0,0,1,1,0,0,0,1,1,0,0,0,0,1,0,0,0,0,6,0,1,0,1,onshore_wind_turbine,0.8,1,1,NE,1,0
+battery,1,0,0,1,0,0,0,0,0,0,1,0,0,0,-1,-1,-1,0,0,0,19584,22494,0,4895,5622,0,0.15,0.15,0,None,0,0,0,0,0,1,1,0,0,0,0.92,0.92,1,10,0,0,1,0,0,0,0,12,0,0,1,0,battery_mid,0.95,0,0,NE,0,0

src/load_inputs/load_generators_data.jl

@@ -77,7 +77,6 @@ function load_generators_data!(setup::Dict, path::AbstractString, inputs_gen::Di
 
 	# Add Resource IDs after reading to prevent user errors
 	gen_in[!,:R_ID] = 1:length(collect(skipmissing(gen_in[!,1])))
-
 	# Store DataFrame of generators/resources input data for use in model
 	inputs_gen["dfGen"] = gen_in
 

src/model/resources/thermal_storage/thermal_storage.jl

@@ -33,6 +33,11 @@ function get_nonfus(inputs::Dict)::Vector{Int}
 	dfTS[dfTS.FUS.==0,:R_ID]
 end
 
+function get_resistive_heating(inputs::Dict)::Vector{Int}
+	dfTS = inputs["dfTS"]
+	dfTS[dfTS.RH.==1,:R_ID]
+end
+
 function get_maintenance(inputs::Dict)::Vector{Int}
 	dfTS = inputs["dfTS"]
 	if "MAINT" in names(dfTS)
@@ -102,6 +107,7 @@ function thermal_storage(EP::Model, inputs::Dict, setup::Dict)
 	# Load thermal storage inputs
 	TS = inputs["TS"]
 	dfTS = inputs["dfTS"]
+	RH = get_resistive_heating(inputs)
 
 	by_rid(rid, sym) = by_rid_df(rid, sym, dfTS)
 
@@ -115,6 +121,12 @@ function thermal_storage(EP::Model, inputs::Dict, setup::Dict)
 	vTSCAP[y in TS] >= 0			#thermal storage energy capacity for resource y
 	end)
 
+	# resistive heating variables
+	@variables(EP, begin
+		vRH[y in RH, t = 1:T] >= 0 		#electrical energy from grid
+		vRHCAP[y in RH] >= 0			#RH power capacity for resource
+	end)
+
 	### THERMAL CORE CONSTRAINTS ###
 	# Core power output must be <= installed capacity, including hourly capacity factors
 	@constraint(EP, cCPMax[y in TS, t=1:T], vCP[y,t] <= vCCAP[y]*inputs["pP_Max"][y,t])
@@ -145,14 +157,32 @@ function thermal_storage(EP::Model, inputs::Dict, setup::Dict)
 
 	# thermal state of charge balance for interior timesteps:
 	# (previous SOC) - (discharge to turbines) - (turbine startup energy use) + (core power output) - (self discharge)
-	@constraint(EP, cTSocBalInterior[t in INTERIOR_SUBPERIODS, y in TS], (
+
+	# first for resources with no RH
+	@constraint(EP, cTSocBalInterior[t in INTERIOR_SUBPERIODS, y in setdiff(TS, RH)], (
 		vTS[y,t] == vTS[y,t-1]
 		- (1 / dfGen[y, :Eff_Down] * EP[:vP][y,t])
 		- (1 / dfGen[y, :Eff_Down] * dfGen[y, :Start_Fuel_MMBTU_per_MW] * dfGen[y,:Cap_Size] * EP[:vSTART][y,t])
 		+ (dfGen[y,:Eff_Up] * vCP[y,t])
 		- (dfGen[y,:Self_Disch] * vTS[y,t-1]))
 	)
 
+	# then for resources with RH
+	@constraint(EP, cTSocBalInteriorRH[t in INTERIOR_SUBPERIODS, y in intersect(TS, RH)], (
+		vTS[y,t] == vTS[y,t-1]
+		- (1 / dfGen[y, :Eff_Down] * EP[:vP][y,t])
+		- (1 / dfGen[y, :Eff_Down] * dfGen[y, :Start_Fuel_MMBTU_per_MW] * dfGen[y,:Cap_Size] * EP[:vSTART][y,t])
+		+ (dfGen[y,:Eff_Up] * vCP[y,t])
+		- (dfGen[y,:Self_Disch] * vTS[y,t-1])
+		+ (vRH[y, t]))  #100% resistive heating efficiency
+	)
+
+	# add resistive heating to power balance
+	@expression(EP, ePowerBalanceRH[t=1:T, z=1:Z],
+		-1 * sum(vRH[y, t] for y in intersect(RH, dfGen[dfGen[!, :Zone].==z, :R_ID])))
+	EP[:ePowerBalance] += ePowerBalanceRH
+
+
 	# TODO: perhaps avoid recomputing these; instead use sets TS_LONG_DURATION, etc
 	TS_and_LDS, TS_and_nonLDS = split_LDS_and_nonLDS(dfGen, inputs, setup)
 
@@ -161,7 +191,7 @@ function thermal_storage(EP::Model, inputs::Dict, setup::Dict)
 		- (1 / dfGen[y, :Eff_Down] * EP[:vP][y,t])
 		- (1 / dfGen[y, :Eff_Down] * dfGen[y, :Start_Fuel_MMBTU_per_MW] * dfGen[y, :Cap_Size] * EP[:vSTART][y,t])
 		+ (dfGen[y, :Eff_Up] * vCP[y,t]) - (dfGen[y, :Self_Disch] * vTS[y,t + hours_per_subperiod - 1])
-		))
+	))
 
 	if !isempty(TS_and_LDS)
 		REP_PERIOD = inputs["REP_PERIOD"]  # Number of representative periods
@@ -213,6 +243,13 @@ function thermal_storage(EP::Model, inputs::Dict, setup::Dict)
 	@expression(EP, eTotalCFixedTS, sum(eCFixed_TS[y] for y in TS))
 	EP[:eObj] += eTotalCFixedTS
 
+	# Resistive heating investment costs
+	# Fixed costs for resource y
+	@expression(EP, eCFixed_RH[y in RH], by_rid(y, :Fixed_Cost_per_MW_RH) * vRHCAP[y])
+	# Total fixed costs for all resistive heating
+	@expression(EP, eTotalCFixedRH, sum(eCFixed_RH[y] for y in RH))
+	EP[:eObj] += eTotalCFixedRH
+
 	# Parameter Fixing Constraints
 	# Fixed ratio of gross generator capacity to core equivalent gross electric power
 	@constraint(EP, cCPRatMax[y in dfTS[dfTS.Max_Generator_Core_Power_Ratio.>0,:R_ID]],
@@ -228,14 +265,20 @@ function thermal_storage(EP::Model, inputs::Dict, setup::Dict)
 	### FUSION CONSTRAINTS ###
 	FUS =  get_fus(inputs)
 
-	# TODO: define expressions & constraints for NONFUS
-	NONFUS =  get_nonfus(inputs)
-
 	# Use fusion constraints if thermal cores tagged 'FUS' are present
 	if !isempty(FUS)
 		fusion_constraints!(EP, inputs, setup)
 	end
 
+	### NONFUSION CONSTRAINTS ###
+	NONFUS = get_nonfus(inputs)
+
+	# use thermal core constraints for thermal cores not tagged 'FUS' 
+	if !isempty(NONFUS)
+		thermal_core_constraints!(EP, inputs, setup)
+		thermal_core_max_cap_constraint!(EP, inputs, setup)
+	end
+
 	# Capacity Reserves Margin policy
 	if setup["CapacityReserveMargin"] > 0
 		ncap = inputs["NCapacityReserveMargin"]
@@ -252,7 +295,7 @@ function thermal_storage(EP::Model, inputs::Dict, setup::Dict)
 		EP[:eCapResMarBalance] += eCapResMarBalanceFusionAdjustment
 	end
 
-return EP
+	return EP
 end
 
 function fusion_max_cap_constraint!(EP::Model, inputs::Dict, setup::Dict)
@@ -288,6 +331,21 @@ function fusion_max_cap_constraint!(EP::Model, inputs::Dict, setup::Dict)
 	end
 end
 
+# set maximum power capacity constraint for the thermal core (in MWth)
+function thermal_core_max_cap_constraint!(EP::Model, inputs::Dict, setup::Dict)
+
+	dfTS = inputs["dfTS"]
+	by_rid(rid, sym) = by_rid_df(rid, sym, dfTS)
+
+	NONFUS = get_nonfus(inputs)
+
+	#System-wide installed capacity is less than a specified maximum limit
+	HAS_MAX_LIMIT = dfTS[by_rid(NONFUS, :Max_Core_Power_Capacity) .> 0, :R_ID]
+	@constraint(EP, cCoreMaxCapacity[y in HAS_MAX_LIMIT], EP[:vCCAP][y] <= by_rid(y, :Max_Core_Power_Capacity) / by_rid(y, :Eff_Down))
+
+end
+
+
 @doc raw"""
     fusion_constraints!(EP::Model, inputs::Dict)
 
@@ -410,6 +468,131 @@ function fusion_constraints!(EP::Model, inputs::Dict, setup::Dict)
 	EP[:ePowerBalance] += ePowerBalanceRecircFus
 end
 
+# TODO make compatible with reserves
+function thermal_core_constraints!(EP::Model, inputs::Dict, setup::Dict)
+
+	dfGen = inputs["dfGen"]
+	dfTS = inputs["dfTS"]
+	by_rid(rid, sym) = by_rid_df(rid, sym, dfTS)
+
+	T = inputs["T"]     # Number of time steps (hours)
+	Z = inputs["Z"]     # Number of zones
+	G = inputs["G"]     # Number of resources
+	NONFUS = get_nonfus(inputs) # non fusion thermal cores
+
+	hours_per_subperiod = inputs["hours_per_subperiod"] #total number of hours per subperiod
+	START_SUBPERIODS = inputs["START_SUBPERIODS"]
+	INTERIOR_SUBPERIODS = inputs["INTERIOR_SUBPERIODS"]
+
+	COMMIT = intersect(inputs["THERM_COMMIT"], NONFUS)
+	NON_COMMIT = intersect(inputs["THERM_NO_COMMIT"])
+
+	# constraints for generators not subject to UC
+	if !isempty(NON_COMMIT)
+
+		# ramp up and ramp down rates
+		@constraints(EP, begin
+
+			# ramp up, start
+			[y in NON_COMMIT, t in START_SUBPERIODS], EP[:vCP][y, t] - EP[:vCP][y, t+hours_per_subperiod-1] <= by_rid(y, :Ramp_Up_Percentage) * EP[:vCCAP][y]
+
+			# ramp up, interior
+			[y in NON_COMMIT, t in INTERIOR_SUBPERIODS], EP[:vCP][y, t] - EP[:vCP][y, t-1] <= by_rid(y, :Ramp_Up_Percentage) * EP[:vCCAP][y]
+
+			# ramp dn, start
+			[y in NON_COMMIT, t in START_SUBPERIODS], EP[:vCP][y,t+hours_per_subperiod-1] - EP[:vCP][y,t] <= by_rid(y, :Ramp_Dn_Percentage) * EP[:vCCAP][y]
+
+			# ramp dn, interior
+			[y in NON_COMMIT, t in INTERIOR_SUBPERIODS], EP[:vCP][y,t-1] - EP[:vCP][y,t] <= by_rid(y, :Ramp_Dn_Percentage) * EP[:vCCAP][y]
+		end)
+
+		# minimum stable power (assumes capacity factor of 1, so max power already implemented)
+		@constraint(EP, [y in NON_COMMIT, t=1:T], EP[:vCP][y,t] >= by_rid(y, :Min_Power)* EP[:vCCAP][y])
+	end
+
+	# constraints for generatiors subject to UC
+	if !isempty(COMMIT)
+
+		### Decision variables for unit commitment  ###
+		# commitment state variable
+		@variable(EP, vCCOMMIT[y in COMMIT, t=1:T] >= 0)
+		# startup event variable
+		@variable(EP, vCSTART[y in COMMIT, t=1:T] >= 0)
+		# shutdown event variable
+		@variable(EP, vCSHUT[y in COMMIT, t=1:T] >= 0)
+
+		### TODO: STARTUP COSTS ???????? ###
+
+		## Declaration of integer/binary variables
+		if setup["UCommit"] == 1 # Integer UC constraints
+			for y in COMMIT
+				set_integer.(vCCOMMIT[y,:])
+				set_integer.(vCSTART[y,:])
+				set_integer.(vCSHUT[y,:])
+				set_integer(EP[:vCCAP][y])
+			end
+		end
+
+		### Capacitated limits on unit commitment decision variables 
+		@constraints(EP, begin
+			[y in COMMIT, t=1:T], vCCOMMIT[y,t] <= EP[:vCCAP][y] / by_rid(y,:Cap_Size)
+			[y in COMMIT, t=1:T], vCSTART[y,t] <= EP[:vCCAP][y] / by_rid(y,:Cap_Size)
+			[y in COMMIT, t=1:T], vCSHUT[y,t] <= EP[:vCCAP][y] / by_rid(y,:Cap_Size)
+		end)
+
+		# Commitment state constraint linking startup and shutdown decisions (Constraint #4)
+		@constraints(EP, begin
+			# For Start Hours, links first time step with last time step in subperiod
+			[y in COMMIT, t in START_SUBPERIODS], vCCOMMIT[y,t] == vCCOMMIT[y,(t+hours_per_subperiod-1)] + vCSTART[y,t] - vCSHUT[y,t]
+			# For all other hours, links commitment state in hour t with commitment state in prior hour + sum of start up and shut down in current hour
+			[y in COMMIT, t in INTERIOR_SUBPERIODS], vCCOMMIT[y,t] == vCCOMMIT[y,t-1] + vCSTART[y,t] - vCSHUT[y,t]
+		end)
+
+		#ramp up, start
+		@constraint(EP,[y in COMMIT, t in START_SUBPERIODS],
+			EP[:vCP][y,t]-EP[:vCP][y,(t+hours_per_subperiod-1)] <= by_rid(y,:Ramp_Up_Percentage)*by_rid(y,:Cap_Size)*(vCCOMMIT[y,t]-vCSTART[y,t])
+			+ min(1, max(by_rid(y,:Min_Power), by_rid(y,:Ramp_Up_Percentage)))*by_rid(y,:Cap_Size)*vCSTART[y,t]
+			- by_rid(y,:Min_Power)*by_rid(y,:Cap_Size)*vCSHUT[y,t])
+
+		#ramp up, interior
+		@constraint(EP,[y in COMMIT, t in INTERIOR_SUBPERIODS],
+			EP[:vCP][y,t]-EP[:vCP][y,(t-1)] <= by_rid(y,:Ramp_Up_Percentage)*by_rid(y,:Cap_Size)*(vCCOMMIT[y,t]-vCSTART[y,t])
+				+ min(1,max(by_rid(y,:Min_Power), by_rid(y,:Ramp_Up_Percentage)))*by_rid(y,:Cap_Size)*vCSTART[y,t]
+				- by_rid(y,:Min_Power)*by_rid(y,:Cap_Size)*vCSHUT[y,t])
+
+		#ramp down, start
+		@constraint(EP,[y in COMMIT, t in START_SUBPERIODS],
+			EP[:vCP][y,(t+hours_per_subperiod-1)]-EP[:vCP][y,t] <= by_rid(y,:Ramp_Dn_Percentage)*by_rid(y,:Cap_Size)*(vCCOMMIT[y,t]-vCSTART[y,t])
+			- by_rid(y,:Min_Power)*by_rid(y,:Cap_Size)*vCSTART[y,t]
+			+ min(1,max(by_rid(y,:Min_Power), by_rid(y,:Ramp_Dn_Percentage)))*by_rid(y,:Cap_Size)*vCSHUT[y,t])
+
+		#ramp down, interior
+		@constraint(EP,[y in COMMIT, t in INTERIOR_SUBPERIODS],
+			EP[:vCP][y,(t-1)]-EP[:vCP][y,t] <= by_rid(y,:Ramp_Dn_Percentage)*by_rid(y,:Cap_Size)*(vCCOMMIT[y,t]-vCSTART[y,t])
+			- by_rid(y,:Min_Power)*by_rid(y,:Cap_Size)*vCSTART[y,t]
+			+ min(1,max(by_rid(y,:Min_Power), by_rid(y,:Ramp_Dn_Percentage)))*by_rid(y,:Cap_Size)*vCSHUT[y,t])
+
+		# minimum and maximum stable power (assumes capacity factor of 1, so max power already implemented)
+		@constraints(EP, begin
+			[y in COMMIT, t=1:T], EP[:vCP][y,t] >= by_rid(y,:Min_Power)*by_rid(y,:Cap_Size)*vCCOMMIT[y,t]
+		end)
+
+		### Minimum up and down times (Constraints #9-10)
+		p = hours_per_subperiod
+		Up_Time = zeros(Int, nrow(dfGen))
+		Up_Time[COMMIT] .= Int.(floor.(dfTS[COMMIT,:Up_Time]))
+		@constraint(EP, [y in COMMIT, t in 1:T],
+			vCCOMMIT[y,t] >= sum(vCSTART[y, hoursbefore(p, t, 0:(Up_Time[y] - 1))])
+		)
+		Down_Time = zeros(Int, nrow(dfGen))
+		Down_Time[COMMIT] .= Int.(floor.(dfTS[COMMIT,:Down_Time]))
+		@constraint(EP, [y in COMMIT, t in 1:T],
+			EP[:vCCAP][y]/by_rid(y,:Cap_Size)-vCCOMMIT[y,t] >= sum(vCSHUT[y, hoursbefore(p, t, 0:(Down_Time[y] - 1))])
+		)
+
+	end
+end
+
 function maintenance_constraints!(EP::Model, inputs::Dict, setup::Dict)
 
 	println("Fusion Core Maintenance Module")