Merge pull request #387 from NREL/develop

v. 0.46.0

Author: zolanaj

.github/workflows/CI.yml

@@ -13,7 +13,8 @@ jobs:
       matrix:
         julia-version: ['1.8']
         julia-arch: [x64]
-        os: [ubuntu-latest, windows-latest, macOS-11]
+        # os: [ubuntu-latest, windows-latest, macOS-11]
+        os: [windows-latest, macOS-11]
 
     steps:
       - uses: actions/checkout@v2

CHANGELOG.md

@@ -23,8 +23,43 @@ Classify the change according to the following categories:
     ### Deprecated
     ### Removed
 
+## v. 0.46.0
+### Added 
+- In `src/core/absorption_chiller.jl` struct, added field **heating_load_input** to the AbsorptionChiller struct
+- Added new variables **dvHeatToStorage** and **dvHeatFromStorage** which are indexed on `p.heating_loads` and added reconciliation constraints so that **dvProductionToStorage** and **dvDischargeFromStorage** maintain their relationship to state of charge for Hot thermal energy storage.
+- In `src/constraints/thermal_tech_constraints.jl`, added function **no_existing_boiler_production** which prevents ExistingBoiler from producing heat in optimized (non-BAU) scenarios 
+- for all heating techs and CHP, added fields **can_serve_space_heating**, **can_serve_dhw**, and **can_serve_process_heat** in core structs and added new results fields **thermal_to_dhw_load_series_mmbtu_per_hour**, **thermal_to_space_heating_load_series_mmbtu_per_hour**, and **thermal_to_process_heat_load_series_mmbtu_per_hour**
+- In `src/core/techs.jl`, added new sets **ghp_techs**, **cooling_techs**, **techs_can_serve_space_heating**, **techs_can_serve_dhw**, and **techs_can_serve_process_heat**
+- In `src/core/reopt_inputs.jl`, added new fields **heating_loads**, **heating_loads_kw**, **heating_loads_served_by_tes**, and **absorption_chillers_using_heating_load** to the REoptInputs and BAUInputs structs. in the math, new set `p.heating_loads` has index q (to represent "qualities" of heat).
+- In `src/core/heating_cooling_loads.jl`, added new struct **ProcessHeatLoad**
+- In `src/core/scenario.jl`, added new field **process_heat_load**
+- In `src/mpc/inputs.jl`, added new field **heating_loads**
+- In `src/core/existing_boiler.jl`, added field **retire_in_optimal** to the ExistingBoiler struct
+- Info to user including name of PV and/or temperature datasource used and distance from site location to datasource location
+- Warning to user if data is not from NSRDB or if data is more than 200 miles away
+- In `results/heating_cooling_load.jl`, added new fields **process_heat_thermal_load_series_mmbtu_per_hour**, **process_heat_boiler_fuel_load_series_mmbtu_per_hour**, **annual_calculated_process_heat_thermal_load_mmbtu**, and **annual_calculated_process_heat_boiler_fuel_load_mmbtu** to HeatingLoad results, with sum heating loads now including process heat 
+### Changed
+- Change the way we determine which dataset to utilize in the PVWatts API call. Previously, we utilized defined lat-long bounds to determine if "nsrdb" or "intl" data should be used in PVWatts call. Now, we call the Solar Dataset Query API (v2) (https://developer.nrel.gov/docs/solar/data-query/v2/) to determine the dataset to use, and include "tmy3" as an option, as this is currently the best-available data for many locations in Alaska. 
+- Refactored **dvThermalProduction** to be separated in **dvCoolingProduction** and **dvHeatingProduction** with **dvHeatingProduction** now indexed on `p.heating_loads`
+- Refactored heating load balance constraints so that a separate flow balance is reconciled for each heating load in `p.heating_loads`
+- Renamed **dvThermalProductionYIntercept** to **dvHeatingProductionYIntercept**
+- Divided **ThermalStorage** into **HotThermalStorage** and **ColdThermalStorage** as the former now has attributes related to the compatible heat loads as input or output.
+- Changed technologies included **dvProductionToWaste** to all heating techs.  NOTE: this variable is forced to zero to allow steam turbine tests to pass, but I believe that waste heat should be allowed for the turbine.  A TODO is in place to review this commit (a406cc5df6e4a27b56c92815c35d04815904e495).
+- Changed test values and tolerances for CHP Sizing test.
+- Updated test sets "Emissions and Renewable Energy Percent" and "Minimize Unserved Load" to decrease computing time.
+- Test for tiered TOU demand rates in `test/runtests.jl`
+- Updated `pop_year` and `income_year` used in call to EASIUR data (`get_EASIUR2005`) each to 2024, from 2020. 
+- Updated usd conversion used for EASIUR health cost calcs from USD_2010_to_2020 = 1.246 to USD_2010_to_2024 = 1.432
+### Fixed  
+- Added a constraint in `src/constraints/steam_turbine_constraints.jl` that allows for heat loads to reconcile when thermal storage is paired with a SteamTurbine. 
+- Fixed a bug in which net-metering system size limits could be exceeded while still obtaining the net-metering benefit due to a large "big-M".
+- Fixed a reshape call in function `parse_urdb_tou_demand` that incorrectly assumed row major instead of column major ordering
+- Fixed a loop range in function `parse_urdb_tou_demand` that incorrectly started at 0 instead of 1
+- Added the missing tier index when accessing `p.s.electric_tariff.tou_demand_rates` in function `add_elec_utility_expressions`
+
 ## v0.45.0
 ### Fixed 
+- Fixed bug in call to `GhpGhx.jl` when sizing hybrid GHP using the fractional sizing method
 - Added `export_rate_beyond_net_metering_limit` to list of inputs to be converted to type Real, to avoid MethodError if type is vector of Any. 
 - Fix blended CRB processing when one or more load types have zero annual energy
 - When calculating CHP fuel intercept and slope, use 1 for the HHV because CHP fuel measured in units of kWh, instead of using non-existent **CHP.fuel_higher_heating_value_kwh_per_gal**

Project.toml

@@ -1,7 +1,7 @@
 name = "REopt"
 uuid = "d36ad4e8-d74a-4f7a-ace1-eaea049febf6"
 authors = ["Nick Laws", "Hallie Dunham <hallie.dunham@nrel.gov>", "Bill Becker <william.becker@nrel.gov>", "Bhavesh Rathod <bhavesh.rathod@nrel.gov>", "Alex Zolan <alexander.aolan@nrel.gov>", "Amanda Farthing <amanda.farthing@nrel.gov>"]
-version = "0.45.0"
+version = "0.46.0"
 
 [deps]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"

data/absorption_chiller/absorption_chiller_defaults.json

@@ -39,7 +39,8 @@
             20.0,
             18.0
         ],
-        "cop_thermal": 0.74
+        "cop_thermal": 0.74,
+        "heating_load_input": "DomesticHotWater"
     },
     "steam":{
         "installed_cost_per_ton": [
@@ -78,6 +79,7 @@
             23.0,
             20.0
         ],
-        "cop_thermal":1.42
+        "cop_thermal":1.42,
+        "heating_load_input": "DomesticHotWater"
     }
 }

docs/src/reopt/inputs.md

@@ -147,6 +147,11 @@ REopt.DomesticHotWaterLoad
 REopt.SpaceHeatingLoad
 ```
 
+## ProcessHeatLoad
+```@docs
+REopt.ProcessHeatLoad
+```
+
 ## FlexibleHVAC
 ```@docs
 REopt.FlexibleHVAC
@@ -166,3 +171,8 @@ REopt.GHP
 ```@docs
 REopt.SteamTurbine
 ```
+
+## ElectricHeater
+```@docs
+REopt.ElectricHeater
+```

src/constraints/chp_constraints.jl

@@ -47,36 +47,37 @@ function add_chp_thermal_production_constraints(m, p; _n="")
     thermal_prod_slope = (thermal_prod_full_load - thermal_prod_half_load) / (1.0 - 0.5)  # [kWt/kWe]
     thermal_prod_intercept = thermal_prod_full_load - thermal_prod_slope * 1.0  # [kWt/kWe_rated
 
-    # Conditionally add dvThermalProductionYIntercept if coefficient p.s.chpThermalProdIntercept is greater than ~zero
+
+    # Conditionally add dvHeatingProductionYIntercept if coefficient p.s.chpThermalProdIntercept is greater than ~zero
     if abs(thermal_prod_intercept) > 1.0E-7
-        dv = "dvThermalProductionYIntercept"*_n
+        dv = "dvHeatingProductionYIntercept"*_n
         m[Symbol(dv)] = @variable(m, [p.techs.chp, p.time_steps], base_name=dv)
 
         #Constraint (2a-1): Upper Bounds on Thermal Production Y-Intercept
         @constraint(m, CHPYInt2a1Con[t in p.techs.chp, ts in p.time_steps],
-            m[Symbol("dvThermalProductionYIntercept"*_n)][t,ts] <= thermal_prod_intercept * m[Symbol("dvSize"*_n)][t]
+            m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts] <= thermal_prod_intercept * m[Symbol("dvSize"*_n)][t]
         )
         # Constraint (2a-2): Upper Bounds on Thermal Production Y-Intercept
         @constraint(m, CHPYInt2a2Con[t in p.techs.chp, ts in p.time_steps],
-            m[Symbol("dvThermalProductionYIntercept"*_n)][t,ts] <= thermal_prod_intercept * p.s.chp.max_kw 
+            m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts] <= thermal_prod_intercept * p.s.chp.max_kw 
             * m[Symbol("binCHPIsOnInTS"*_n)][t,ts]
         )
         #Constraint (2b): Lower Bounds on Thermal Production Y-Intercept
         @constraint(m, CHPYInt2bCon[t in p.techs.chp, ts in p.time_steps],
-            m[Symbol("dvThermalProductionYIntercept"*_n)][t,ts] >= thermal_prod_intercept * m[Symbol("dvSize"*_n)][t] 
+            m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts] >= thermal_prod_intercept * m[Symbol("dvSize"*_n)][t] 
             - thermal_prod_intercept * p.s.chp.max_kw * (1 - m[Symbol("binCHPIsOnInTS"*_n)][t,ts])
         )
         # Constraint (2c): Thermal Production of CHP
         # Note: p.HotWaterAmbientFactor[t,ts] * p.HotWaterThermalFactor[t,ts] removed from this but present in math
         @constraint(m, CHPThermalProductionCon[t in p.techs.chp, ts in p.time_steps],
-            m[Symbol("dvThermalProduction"*_n)][t,ts] ==
+        sum(m[Symbol("dvHeatingProduction"*_n)][t,q,ts] for q in p.heating_loads) ==
             thermal_prod_slope * p.production_factor[t,ts] * m[Symbol("dvRatedProduction"*_n)][t,ts] 
-            + m[Symbol("dvThermalProductionYIntercept"*_n)][t,ts] +
+            + m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts] +
             m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts]
         )
     else
         @constraint(m, CHPThermalProductionConLinear[t in p.techs.chp, ts in p.time_steps],
-            m[Symbol("dvThermalProduction"*_n)][t,ts] ==
+            sum(m[Symbol("dvHeatingProduction"*_n)][t,q,ts] for q in p.heating_loads) ==
             thermal_prod_slope * p.production_factor[t,ts] * m[Symbol("dvRatedProduction"*_n)][t,ts] +
             m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts]
         )        
@@ -99,7 +100,7 @@ function add_chp_supplementary_firing_constraints(m, p; _n="")
     # Constrain upper limit of dvSupplementaryThermalProduction, using auxiliary variable for (size * useSupplementaryFiring)
     @constraint(m, CHPSupplementaryFireCon[t in p.techs.chp, ts in p.time_steps],
                 m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] <=
-                (p.s.chp.supplementary_firing_max_steam_ratio - 1.0) * p.production_factor[t,ts] * (thermal_prod_slope * m[Symbol("dvSupplementaryFiringSize"*_n)][t] + m[Symbol("dvThermalProductionYIntercept"*_n)][t,ts])
+                (p.s.chp.supplementary_firing_max_steam_ratio - 1.0) * p.production_factor[t,ts] * (thermal_prod_slope * m[Symbol("dvSupplementaryFiringSize"*_n)][t] + m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts])
                 )
     if solver_is_compatible_with_indicator_constraints(p.s.settings.solver_name)
         # Constrain lower limit of 0 if CHP tech is off
@@ -110,7 +111,7 @@ function add_chp_supplementary_firing_constraints(m, p; _n="")
         #There's no upper bound specified for the CHP supplementary firing, so assume the entire heat load as a reasonable maximum that wouldn't be exceeded (but might not be the best possible value). 
         max_supplementary_firing_size = maximum(p.s.dhw_load.loads_kw .+ p.s.space_heating_load.loads_kw)
         @constraint(m, NoCHPSupplementaryFireOffCon[t in p.techs.chp, ts in p.time_steps],
-                m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] <= (p.s.chp.supplementary_firing_max_steam_ratio - 1.0) * p.production_factor[t,ts] * (thermal_prod_slope * max_supplementary_firing_size + m[Symbol("dvThermalProductionYIntercept"*_n)][t,ts])
+                m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] <= (p.s.chp.supplementary_firing_max_steam_ratio - 1.0) * p.production_factor[t,ts] * (thermal_prod_slope * max_supplementary_firing_size + m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts])
                 )
     end
 end

src/constraints/electric_utility_constraints.jl

@@ -65,8 +65,21 @@ function add_export_constraints(m, p; _n="")
                     !binNEM => {sum(m[Symbol("dvSize"*_n)][t] for t in NEM_techs) <= p.s.electric_utility.interconnection_limit_kw}
                 )
             else
+                #leverage max system sizes for interconnect limit size, alternate is max monthly fully-electrified load in kWh
+                #assume electric heater with COP of 1 for conversion of heat to electricity
+                max_interconnection_size = minimum([
+                    p.s.electric_utility.interconnection_limit_kw, 
+                    sum(p.max_sizes[t] for t in NEM_techs),
+                    p.hours_per_time_step * maximum([sum((
+                        p.s.electric_load.loads_kw[ts] + 
+                        p.s.cooling_load.loads_kw_thermal[ts]/p.cop["ExistingChiller"] + 
+                        (p.s.space_heating_load.loads_kw[ts] + p.s.dhw_load.loads_kw[ts] + p.s.process_heat_load.loads_kw[ts]) 
+                    ) for ts in p.s.electric_tariff.time_steps_monthly[m]) for m in p.months
+                    ])
+                ])
+                
                 @constraint(m,
-                    sum(m[Symbol("dvSize"*_n)][t] for t in NEM_techs) <= p.s.electric_utility.interconnection_limit_kw - (p.s.electric_utility.interconnection_limit_kw - p.s.electric_utility.net_metering_limit_kw)*binNEM 
+                    sum(m[Symbol("dvSize"*_n)][t] for t in NEM_techs) <= max_interconnection_size - (max_interconnection_size - p.s.electric_utility.net_metering_limit_kw)*binNEM 
                 )
             end
 
@@ -371,7 +384,7 @@ function add_elec_utility_expressions(m, p; _n="")
 
     if !isempty(p.s.electric_tariff.tou_demand_rates)
         m[Symbol("DemandTOUCharges"*_n)] = @expression(m, 
-            p.pwf_e * sum( p.s.electric_tariff.tou_demand_rates[r] * m[Symbol("dvPeakDemandTOU"*_n)][r, tier] 
+            p.pwf_e * sum( p.s.electric_tariff.tou_demand_rates[r, tier] * m[Symbol("dvPeakDemandTOU"*_n)][r, tier] 
             for r in p.ratchets, tier in 1:p.s.electric_tariff.n_tou_demand_tiers)
         )
     else

src/constraints/emissions_constraints.jl

@@ -1,17 +1,21 @@
 # REopt®, Copyright (c) Alliance for Sustainable Energy, LLC. See also https://github.com/NREL/REopt.jl/blob/master/LICENSE.
 
 function add_emissions_constraints(m,p)
-	if !isnothing(p.s.site.CO2_emissions_reduction_min_fraction)
-		@constraint(m, MinEmissionsReductionCon, 
-			m[:Lifecycle_Emissions_Lbs_CO2] <= 
-			(1-p.s.site.CO2_emissions_reduction_min_fraction) * m[:Lifecycle_Emissions_Lbs_CO2_BAU]
-		)
-	end
-	if !isnothing(p.s.site.CO2_emissions_reduction_max_fraction)
-		@constraint(m, MaxEmissionsReductionCon, 
-			m[:Lifecycle_Emissions_Lbs_CO2] >= 
-			(1-p.s.site.CO2_emissions_reduction_max_fraction) * m[:Lifecycle_Emissions_Lbs_CO2_BAU]
-		)
+	if !isnothing(p.s.site.bau_emissions_lb_CO2_per_year)
+		if !isnothing(p.s.site.CO2_emissions_reduction_min_fraction)
+			@constraint(m, MinEmissionsReductionCon, 
+				m[:Lifecycle_Emissions_Lbs_CO2] <= 
+				(1-p.s.site.CO2_emissions_reduction_min_fraction) * m[:Lifecycle_Emissions_Lbs_CO2_BAU]
+			)
+		end
+		if !isnothing(p.s.site.CO2_emissions_reduction_max_fraction)
+			@constraint(m, MaxEmissionsReductionCon, 
+				m[:Lifecycle_Emissions_Lbs_CO2] >= 
+				(1-p.s.site.CO2_emissions_reduction_max_fraction) * m[:Lifecycle_Emissions_Lbs_CO2_BAU]
+			)
+		end
+	else
+		@warn "No emissions reduction constraints added, as BAU emissions have not been calculated."
 	end
 end
 
@@ -32,14 +36,19 @@ function add_yr1_emissions_calcs(m,p)
 	yr1_emissions_offset_from_elec_exports_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_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[: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]
+	m[:EmissionsYr1_Total_LbsSO2] = m[:yr1_emissions_onsite_fuel_lbs_SO2] + m[:yr1_emissions_from_elec_grid_net_if_selected_lbs_SO2]
+	m[:EmissionsYr1_Total_LbsPM25] = m[:yr1_emissions_onsite_fuel_lbs_PM25] + m[:yr1_emissions_from_elec_grid_net_if_selected_lbs_PM25]
 	nothing
 end