Hydrogen

src/REopt.jl

@@ -17,6 +17,10 @@ export
     run_mpc,
     build_mpc!, 
     backup_reliability,
+    shift_gen_storage_prob_matrix!,
+    backup_reliability_reopt_inputs,
+    get_maximum_generation,
+    storage_leakage!,
     get_chp_defaults_prime_mover_size_class,
     get_steam_turbine_defaults_size_class,
     simulated_load,
@@ -118,9 +122,14 @@ include("core/wind.jl")
 include("core/energy_storage/storage.jl")
 include("core/energy_storage/electric_storage.jl")
 include("core/energy_storage/thermal_storage.jl")
+include("core/energy_storage/hydrogen_storage.jl")
+include("core/electrolyzer.jl")
+include("core/compressor.jl")
+include("core/fuel_cell.jl")
 include("core/generator.jl")
 include("core/doe_commercial_reference_building_loads.jl")
 include("core/electric_load.jl")
+include("core/hydrogen_load.jl")
 include("core/existing_boiler.jl")
 include("core/boiler.jl")
 include("core/existing_chiller.jl")
@@ -159,6 +168,7 @@ include("constraints/ghp_constraints.jl")
 include("constraints/steam_turbine_constraints.jl")
 include("constraints/renewable_energy_constraints.jl")
 include("constraints/emissions_constraints.jl")
+include("constraints/hydrogen_constraints.jl")
 
 include("mpc/structs.jl")
 include("mpc/scenario.jl")
@@ -189,6 +199,11 @@ include("results/ghp.jl")
 include("results/steam_turbine.jl")
 include("results/electric_heater.jl")
 include("results/heating_cooling_load.jl")
+include("results/hydrogen_load.jl")
+include("results/hydrogen_storage.jl")
+include("results/electrolyzer.jl")
+include("results/compressor.jl")
+include("results/fuel_cell.jl")
 
 include("core/reopt.jl")
 include("core/reopt_multinode.jl")

src/constraints/electric_utility_constraints.jl

@@ -264,7 +264,9 @@ function add_simultaneous_export_import_constraint(m, p; _n="")
         @constraint(m, NoGridPurchasesBinary[ts in p.time_steps],
             m[Symbol("binNoGridPurchases"*_n)][ts] => {
                 sum(m[Symbol("dvGridPurchase"*_n)][ts, tier] for tier in 1:p.s.electric_tariff.n_energy_tiers) +
-                sum(m[Symbol("dvGridToStorage"*_n)][b, ts] for b in p.s.storage.types.elec) <= 0
+                sum(m[Symbol("dvGridToStorage"*_n)][b, ts] for b in p.s.storage.types.elec +
+                m[Symbol("dvGridToElectrolyzer"*_n)][ts] +
+                m[Symbol("dvGridToCompressor"*_n)][ts]) <= 0
             }
         )
         @constraint(m, ExportOnlyAfterSiteLoadMetCon[ts in p.time_steps],
@@ -276,7 +278,9 @@ function add_simultaneous_export_import_constraint(m, p; _n="")
         bigM_hourly_load = maximum(p.s.electric_load.loads_kw)+maximum(p.s.space_heating_load.loads_kw)+maximum(p.s.process_heat_load.loads_kw)+maximum(p.s.dhw_load.loads_kw)+maximum(p.s.cooling_load.loads_kw_thermal)
         @constraint(m, NoGridPurchasesBinary[ts in p.time_steps],
             sum(m[Symbol("dvGridPurchase"*_n)][ts, tier] for tier in 1:p.s.electric_tariff.n_energy_tiers) +
-            sum(m[Symbol("dvGridToStorage"*_n)][b, ts] for b in p.s.storage.types.elec) <= bigM_hourly_load*(1-m[Symbol("binNoGridPurchases"*_n)][ts])
+            sum(m[Symbol("dvGridToStorage"*_n)][b, ts] for b in p.s.storage.types.elec + 
+            m[Symbol("dvGridToElectrolyzer"*_n)][ts] +
+            m[Symbol("dvGridToCompressor"*_n)][ts]) <= bigM_hourly_load*(1-m[Symbol("binNoGridPurchases"*_n)][ts])
         )
         @constraint(m, ExportOnlyAfterSiteLoadMetCon[ts in p.time_steps],
             sum(m[Symbol("dvProductionToGrid"*_n)][t,u,ts] for t in p.techs.elec, u in p.export_bins_by_tech[t]) <= bigM_hourly_load * m[Symbol("binNoGridPurchases"*_n)][ts]

src/constraints/emissions_constraints.jl

@@ -22,28 +22,39 @@ end
 
 function add_yr1_emissions_calcs(m,p)
 	# Components:
-	m[:yr1_emissions_onsite_fuel_lbs_CO2], m[:yr1_emissions_onsite_fuel_lbs_NOx], 
-	m[:yr1_emissions_onsite_fuel_lbs_SO2], m[:yr1_emissions_onsite_fuel_lbs_PM25] = 
+	yr1_emissions_onsite_fuel_series_lbs_CO2, yr1_emissions_onsite_fuel_series_lbs_NOx, 
+	yr1_emissions_onsite_fuel_series_lbs_SO2, yr1_emissions_onsite_fuel_series_lbs_PM25 = 
 		calc_yr1_emissions_from_onsite_fuel(m,p; tech_array=p.techs.fuel_burning)
 
-	m[:yr1_emissions_from_elec_grid_lbs_CO2], m[:yr1_emissions_from_elec_grid_lbs_NOx], 
-	m[:yr1_emissions_from_elec_grid_lbs_SO2], m[:yr1_emissions_from_elec_grid_lbs_PM25] = 
+	yr1_emissions_from_elec_grid_series_lbs_CO2, yr1_emissions_from_elec_grid_series_lbs_NOx, 
+	yr1_emissions_from_elec_grid_series_lbs_SO2, yr1_emissions_from_elec_grid_series_lbs_PM25 = 
 		calc_yr1_emissions_from_elec_grid_purchase(m, p)
 
-	yr1_emissions_offset_from_elec_exports_lbs_CO2, 
-	yr1_emissions_offset_from_elec_exports_lbs_NOx, 
-	yr1_emissions_offset_from_elec_exports_lbs_SO2, 
-	yr1_emissions_offset_from_elec_exports_lbs_PM25 = 
+	yr1_emissions_offset_from_elec_exports_series_lbs_CO2, 
+	yr1_emissions_offset_from_elec_exports_series_lbs_NOx, 
+	yr1_emissions_offset_from_elec_exports_series_lbs_SO2, 
+	yr1_emissions_offset_from_elec_exports_series_lbs_PM25 = 
 		calc_yr1_emissions_offset_from_elec_exports(m, p)
 
-	m[:yr1_emissions_from_elec_grid_net_if_selected_lbs_CO2] = (m[:yr1_emissions_from_elec_grid_lbs_CO2] - 
-		yr1_emissions_offset_from_elec_exports_lbs_CO2)
-	m[:yr1_emissions_from_elec_grid_net_if_selected_lbs_NOx] = (m[:yr1_emissions_from_elec_grid_lbs_NOx] - 
-		yr1_emissions_offset_from_elec_exports_lbs_NOx)
-	m[:yr1_emissions_from_elec_grid_net_if_selected_lbs_SO2] = (m[:yr1_emissions_from_elec_grid_lbs_SO2] - 
-		yr1_emissions_offset_from_elec_exports_lbs_SO2)
-	m[:yr1_emissions_from_elec_grid_net_if_selected_lbs_PM25] = (m[:yr1_emissions_from_elec_grid_lbs_PM25] - 
-		yr1_emissions_offset_from_elec_exports_lbs_PM25)
+	m[:yr1_emissions_from_elec_grid_net_if_selected_series_lbs_CO2] = @expression(m, [ts in p.time_steps], 
+		yr1_emissions_from_elec_grid_series_lbs_CO2[ts] - yr1_emissions_offset_from_elec_exports_series_lbs_CO2[ts])
+	m[:yr1_emissions_from_elec_grid_net_if_selected_series_lbs_NOx] = @expression(m, [ts in p.time_steps], 
+		yr1_emissions_from_elec_grid_series_lbs_NOx[ts] - yr1_emissions_offset_from_elec_exports_series_lbs_NOx[ts])
+	m[:yr1_emissions_from_elec_grid_net_if_selected_series_lbs_SO2] = @expression(m, [ts in p.time_steps], 
+		yr1_emissions_from_elec_grid_series_lbs_SO2[ts] - yr1_emissions_offset_from_elec_exports_series_lbs_SO2[ts])
+	m[:yr1_emissions_from_elec_grid_net_if_selected_series_lbs_PM25] = @expression(m, [ts in p.time_steps], 
+		yr1_emissions_from_elec_grid_series_lbs_PM25[ts] - yr1_emissions_offset_from_elec_exports_series_lbs_PM25[ts])
+
+	# Summed variables
+	m[:yr1_emissions_onsite_fuel_lbs_CO2] = @expression(m, sum(yr1_emissions_onsite_fuel_series_lbs_CO2[ts] for ts in p.time_steps))
+	m[:yr1_emissions_onsite_fuel_lbs_NOx] = @expression(m, sum(yr1_emissions_onsite_fuel_series_lbs_NOx[ts] for ts in p.time_steps))
+	m[:yr1_emissions_onsite_fuel_lbs_SO2] = @expression(m, sum(yr1_emissions_onsite_fuel_series_lbs_SO2[ts] for ts in p.time_steps))
+	m[:yr1_emissions_onsite_fuel_lbs_PM25] = @expression(m, sum(yr1_emissions_onsite_fuel_series_lbs_PM25[ts] for ts in p.time_steps))
+
+	m[:yr1_emissions_from_elec_grid_net_if_selected_lbs_CO2] = @expression(m, sum(m[:yr1_emissions_from_elec_grid_net_if_selected_series_lbs_CO2][ts] for ts in p.time_steps))
+	m[:yr1_emissions_from_elec_grid_net_if_selected_lbs_NOx] = @expression(m, sum(m[:yr1_emissions_from_elec_grid_net_if_selected_series_lbs_NOx][ts] for ts in p.time_steps))
+	m[:yr1_emissions_from_elec_grid_net_if_selected_lbs_SO2] = @expression(m, sum(m[:yr1_emissions_from_elec_grid_net_if_selected_series_lbs_SO2][ts] for ts in p.time_steps))
+	m[:yr1_emissions_from_elec_grid_net_if_selected_lbs_PM25] = @expression(m, sum(m[:yr1_emissions_from_elec_grid_net_if_selected_series_lbs_PM25][ts] for ts in p.time_steps))
 
 	m[:EmissionsYr1_Total_LbsCO2] = m[:yr1_emissions_onsite_fuel_lbs_CO2] + m[:yr1_emissions_from_elec_grid_net_if_selected_lbs_CO2]
 	m[:EmissionsYr1_Total_LbsNOx] = m[:yr1_emissions_onsite_fuel_lbs_NOx] + m[:yr1_emissions_from_elec_grid_net_if_selected_lbs_NOx]
@@ -65,22 +76,22 @@ Function to calculate annual emissions from onsite fuel consumption.
 	account for expected operations during modeled outages (time_steps_without_grid is empty)
 """
 function calc_yr1_emissions_from_onsite_fuel(m,p; tech_array=p.techs.fuel_burning) # also run this with p.techs.boiler
-	yr1_emissions_onsite_fuel_lbs_CO2 = @expression(m,p.hours_per_time_step*
-		sum(m[:dvFuelUsage][t,ts]*p.tech_emissions_factors_CO2[t] for t in tech_array, ts in p.time_steps))
+	yr1_emissions_onsite_fuel_series_lbs_CO2 = @expression(m, [ts in p.time_steps], p.hours_per_time_step*
+		sum(m[:dvFuelUsage][t,ts]*p.tech_emissions_factors_CO2[t] for t in tech_array))
 
-	yr1_emissions_onsite_fuel_lbs_NOx = @expression(m,p.hours_per_time_step*
-		sum(m[:dvFuelUsage][t,ts]*p.tech_emissions_factors_NOx[t] for t in tech_array, ts in p.time_steps))
+	yr1_emissions_onsite_fuel_series_lbs_NOx = @expression(m, [ts in p.time_steps], p.hours_per_time_step*
+		sum(m[:dvFuelUsage][t,ts]*p.tech_emissions_factors_NOx[t] for t in tech_array))
 
-	yr1_emissions_onsite_fuel_lbs_SO2 = @expression(m,p.hours_per_time_step*
-		sum(m[:dvFuelUsage][t,ts]*p.tech_emissions_factors_SO2[t] for t in tech_array, ts in p.time_steps))
+	yr1_emissions_onsite_fuel_series_lbs_SO2 = @expression(m, [ts in p.time_steps], p.hours_per_time_step*
+		sum(m[:dvFuelUsage][t,ts]*p.tech_emissions_factors_SO2[t] for t in tech_array))
 
-	yr1_emissions_onsite_fuel_lbs_PM25 = @expression(m,p.hours_per_time_step*
-		sum(m[:dvFuelUsage][t,ts]*p.tech_emissions_factors_PM25[t] for t in tech_array, ts in p.time_steps))
+	yr1_emissions_onsite_fuel_series_lbs_PM25 = @expression(m, [ts in p.time_steps], p.hours_per_time_step*
+		sum(m[:dvFuelUsage][t,ts]*p.tech_emissions_factors_PM25[t] for t in tech_array))
 
-	return yr1_emissions_onsite_fuel_lbs_CO2, 
-		   yr1_emissions_onsite_fuel_lbs_NOx, 
-		   yr1_emissions_onsite_fuel_lbs_SO2, 
-		   yr1_emissions_onsite_fuel_lbs_PM25
+	return yr1_emissions_onsite_fuel_series_lbs_CO2, 
+		   yr1_emissions_onsite_fuel_series_lbs_NOx, 
+		   yr1_emissions_onsite_fuel_series_lbs_SO2, 
+		   yr1_emissions_onsite_fuel_series_lbs_PM25
 end
 
 """
@@ -96,54 +107,54 @@ Function to calculate annual emissions from grid electricity consumption.
 	account for expected operations during modeled outages (time_steps_without_grid is empty)
 """
 function calc_yr1_emissions_from_elec_grid_purchase(m,p)
-	yr1_emissions_from_elec_grid_lbs_CO2 = @expression(m,p.hours_per_time_step*
-		sum(m[:dvGridPurchase][ts, tier]*p.s.electric_utility.emissions_factor_series_lb_CO2_per_kwh[ts] for ts in p.time_steps, tier in 1:p.s.electric_tariff.n_energy_tiers))
+	yr1_emissions_from_elec_grid_series_lbs_CO2 = @expression(m, [ts in p.time_steps], p.hours_per_time_step*
+		sum(m[:dvGridPurchase][ts, tier]*p.s.electric_utility.emissions_factor_series_lb_CO2_per_kwh[ts] for tier in 1:p.s.electric_tariff.n_energy_tiers))
 
-	yr1_emissions_from_elec_grid_lbs_NOx = @expression(m,p.hours_per_time_step*
-		sum(m[:dvGridPurchase][ts, tier]*p.s.electric_utility.emissions_factor_series_lb_NOx_per_kwh[ts] for ts in p.time_steps, tier in 1:p.s.electric_tariff.n_energy_tiers))
+	yr1_emissions_from_elec_grid_series_lbs_NOx = @expression(m, [ts in p.time_steps], p.hours_per_time_step*
+		sum(m[:dvGridPurchase][ts, tier]*p.s.electric_utility.emissions_factor_series_lb_NOx_per_kwh[ts] for tier in 1:p.s.electric_tariff.n_energy_tiers))
 
-	yr1_emissions_from_elec_grid_lbs_SO2 = @expression(m,p.hours_per_time_step*
-		sum(m[:dvGridPurchase][ts, tier]*p.s.electric_utility.emissions_factor_series_lb_SO2_per_kwh[ts] for ts in p.time_steps, tier in 1:p.s.electric_tariff.n_energy_tiers))
+	yr1_emissions_from_elec_grid_series_lbs_SO2 = @expression(m, [ts in p.time_steps], p.hours_per_time_step*
+		sum(m[:dvGridPurchase][ts, tier]*p.s.electric_utility.emissions_factor_series_lb_SO2_per_kwh[ts] for tier in 1:p.s.electric_tariff.n_energy_tiers))
 
-	yr1_emissions_from_elec_grid_lbs_PM25 = @expression(m,p.hours_per_time_step*
-		sum(m[:dvGridPurchase][ts, tier]*p.s.electric_utility.emissions_factor_series_lb_PM25_per_kwh[ts] for ts in p.time_steps, tier in 1:p.s.electric_tariff.n_energy_tiers))
+	yr1_emissions_from_elec_grid_series_lbs_PM25 = @expression(m, [ts in p.time_steps], p.hours_per_time_step*
+		sum(m[:dvGridPurchase][ts, tier]*p.s.electric_utility.emissions_factor_series_lb_PM25_per_kwh[ts] for tier in 1:p.s.electric_tariff.n_energy_tiers))
 
-	return yr1_emissions_from_elec_grid_lbs_CO2, 
-		   yr1_emissions_from_elec_grid_lbs_NOx, 
-		   yr1_emissions_from_elec_grid_lbs_SO2, 
-		   yr1_emissions_from_elec_grid_lbs_PM25
+	return yr1_emissions_from_elec_grid_series_lbs_CO2, 
+		   yr1_emissions_from_elec_grid_series_lbs_NOx, 
+		   yr1_emissions_from_elec_grid_series_lbs_SO2, 
+		   yr1_emissions_from_elec_grid_series_lbs_PM25
 end
 
 
 function calc_yr1_emissions_offset_from_elec_exports(m, p)
 	if !(p.s.site.include_exported_elec_emissions_in_total)
 		return 0.0, 0.0, 0.0, 0.0
 	end
-	yr1_emissions_offset_from_elec_exports_lbs_CO2 = @expression(m, p.hours_per_time_step *
+	yr1_emissions_offset_from_elec_exports_series_lbs_CO2 = @expression(m, [ts in p.time_steps], p.hours_per_time_step *
 		sum(m[:dvProductionToGrid][t,u,ts] * (p.s.electric_utility.emissions_factor_series_lb_CO2_per_kwh[ts])
-		for t in p.techs.elec, ts in p.time_steps, u in p.export_bins_by_tech[t])
+		for t in p.techs.elec, u in p.export_bins_by_tech[t])
 	)
 		# if battery ends up being able to discharge to grid, need to incorporate here- might require complex tracking of what's charging battery
 
-	yr1_emissions_offset_from_elec_exports_lbs_NOx = @expression(m, p.hours_per_time_step *
+	yr1_emissions_offset_from_elec_exports_series_lbs_NOx = @expression(m, [ts in p.time_steps], p.hours_per_time_step *
 		sum(m[:dvProductionToGrid][t,u,ts] * (p.s.electric_utility.emissions_factor_series_lb_NOx_per_kwh[ts])
-		for t in p.techs.elec, ts in p.time_steps, u in p.export_bins_by_tech[t])
+		for t in p.techs.elec, u in p.export_bins_by_tech[t])
 	)
 
-	yr1_emissions_offset_from_elec_exports_lbs_SO2 = @expression(m, p.hours_per_time_step *
+	yr1_emissions_offset_from_elec_exports_series_lbs_SO2 = @expression(m, [ts in p.time_steps], p.hours_per_time_step *
 		sum(m[:dvProductionToGrid][t,u,ts] * (p.s.electric_utility.emissions_factor_series_lb_SO2_per_kwh[ts])
-		for t in p.techs.elec, ts in p.time_steps, u in p.export_bins_by_tech[t])
+		for t in p.techs.elec, u in p.export_bins_by_tech[t])
 	)
 
-	yr1_emissions_offset_from_elec_exports_lbs_PM25 = @expression(m, p.hours_per_time_step *
+	yr1_emissions_offset_from_elec_exports_series_lbs_PM25 = @expression(m, [ts in p.time_steps], p.hours_per_time_step *
 		sum(m[:dvProductionToGrid][t,u,ts] * (p.s.electric_utility.emissions_factor_series_lb_PM25_per_kwh[ts])
-		for t in p.techs.elec, ts in p.time_steps, u in p.export_bins_by_tech[t])
+		for t in p.techs.elec, u in p.export_bins_by_tech[t])
 	)
 
-	return yr1_emissions_offset_from_elec_exports_lbs_CO2, 
-		   yr1_emissions_offset_from_elec_exports_lbs_NOx, 
-		   yr1_emissions_offset_from_elec_exports_lbs_SO2, 
-		   yr1_emissions_offset_from_elec_exports_lbs_PM25
+	return yr1_emissions_offset_from_elec_exports_series_lbs_CO2, 
+		   yr1_emissions_offset_from_elec_exports_series_lbs_NOx, 
+		   yr1_emissions_offset_from_elec_exports_series_lbs_SO2, 
+		   yr1_emissions_offset_from_elec_exports_series_lbs_PM25
 end