Merge pull request #276 from NREL/develop

v.0.34.0 hybrid GHX and centralized GHP capabilities

CHANGELOG.md

@@ -23,6 +23,18 @@ Classify the change according to the following categories:
     ### Deprecated
     ### Removed
 
+## v0.34.0
+### Added
+- Ability to run hybrid GHX sizing using **GhpGhx.jl** (automatic and fractional sizing)
+- Added financial inputs for **GHP** and updated objective and results to reflect these changes
+- Added central plant **GHP**
+### Fixed
+- Fix output of `get_tier_with_lowest_energy_rate(u::URDBrate)` to return an index and not cartesian coordinates for multi-tier energy rates.
+- Updated **GHP** cost curve calculations so incentives apply to all GHP components
+
+### Changed
+- If a `REoptInputs` object solves with termination status infeasible, altert user and return a dictionary insteadof JuMP model
+
 ## v0.33.0
 ### Added
 - Functionality to evaluate scenarios with Wind can in the ERP (`backup_reliability`)

src/constraints/ghp_constraints.jl

@@ -20,4 +20,13 @@ function add_ghp_constraints(m, p; _n="")
             sum(m[Symbol("binGHP"*_n)][g] for g in p.ghp_options) <= 1
         )
     end
+
+    m[:AvoidedCapexByGHP] = @expression(m,
+        sum(p.avoided_capex_by_ghp_present_value[g] * m[Symbol("binGHP"*_n)][g] for g in p.ghp_options)
+    )
+
+    m[:ResidualGHXCapCost] = @expression(m,
+        sum(p.ghx_residual_value[g] * m[Symbol("binGHP"*_n)][g] for g in p.ghp_options)
+    )
+
 end

src/core/bau_inputs.jl

@@ -102,7 +102,8 @@ function BAUInputs(p::REoptInputs)
     ghp_options, require_ghp_purchase, ghp_heating_thermal_load_served_kw, 
         ghp_cooling_thermal_load_served_kw, space_heating_thermal_load_reduction_with_ghp_kw, 
         cooling_thermal_load_reduction_with_ghp_kw, ghp_electric_consumption_kw, 
-        ghp_installed_cost, ghp_om_cost_year_one = setup_ghp_inputs(bau_scenario, p.time_steps, p.time_steps_without_grid)    
+        ghp_installed_cost, ghp_om_cost_year_one, avoided_capex_by_ghp_present_value,
+        ghx_useful_life_years, ghx_residual_value = setup_ghp_inputs(bau_scenario, p.time_steps, p.time_steps_without_grid)    
 
     # filling export_bins_by_tech MUST be done after techs_by_exportbin has been filled in
     for t in techs.elec
@@ -175,7 +176,10 @@ function BAUInputs(p::REoptInputs)
         cooling_thermal_load_reduction_with_ghp_kw,
         ghp_electric_consumption_kw,
         ghp_installed_cost,
-        ghp_om_cost_year_one,        
+        ghp_om_cost_year_one,
+        avoided_capex_by_ghp_present_value,
+        ghx_useful_life_years,
+        ghx_residual_value,
         tech_renewable_energy_fraction, 
         tech_emissions_factors_CO2, 
         tech_emissions_factors_NOx, 

src/core/electric_tariff.jl

@@ -361,7 +361,7 @@ function get_tier_with_lowest_energy_rate(u::URDBrate)
     """
     #TODO: can eliminate if else if confirm that u.energy_rates is always 2D
     if length(u.energy_tier_limits) > 1
-        return argmin(sum(u.energy_rates, dims=1))
+        return argmin(vec(sum(u.energy_rates, dims=1)))
     else
         return 1
     end

src/core/ghp.jl

@@ -13,6 +13,8 @@ struct with outer constructor:
 ```julia
     require_ghp_purchase::Union{Bool, Int64} = false  # 0 = false, 1 = true
     installed_cost_heatpump_per_ton::Float64 = 1075.0
+    installed_cost_wwhp_heating_pump_per_ton::Float64 = 700.0
+    installed_cost_wwhp_cooling_pump_per_ton::Float64 = 700.0
     heatpump_capacity_sizing_factor_on_peak_load::Float64 = 1.1
     installed_cost_ghx_per_ft::Float64 = 14.0
     installed_cost_building_hydronic_loop_per_sqft = 1.70
@@ -56,11 +58,18 @@ struct with outer constructor:
     om_cost_year_one::Float64 = NaN
 ```
 """
+
+
+
 Base.@kwdef mutable struct GHP <: AbstractGHP
     require_ghp_purchase::Union{Bool, Int64} = false  # 0 = false, 1 = true
     installed_cost_heatpump_per_ton::Float64 = 1075.0
+    installed_cost_wwhp_heating_pump_per_ton::Float64 = 700.0
+    installed_cost_wwhp_cooling_pump_per_ton::Float64 = 700.0
     heatpump_capacity_sizing_factor_on_peak_load::Float64 = 1.1
     installed_cost_ghx_per_ft::Float64 = 14.0
+    ghx_useful_life_years::Int = 50
+    ghx_only_capital_cost::Union{Float64, Nothing} = nothing # overwritten afterwards
     installed_cost_building_hydronic_loop_per_sqft = 1.70
     om_cost_per_sqft_year::Float64 = -0.51
     building_sqft::Float64 # Required input
@@ -69,6 +78,12 @@ Base.@kwdef mutable struct GHP <: AbstractGHP
     ghpghx_response::Dict = Dict()
     can_serve_dhw::Bool = false
 
+    aux_heater_type::String = "electric"
+    is_ghx_hybrid::Bool = false
+    aux_heater_installed_cost_per_mmbtu_per_hr::Float64 = 26000.00
+    aux_cooler_installed_cost_per_ton::Float64 = 400.00
+    aux_unit_capacity_sizing_factor_on_peak_load::Float64 = 1.2
+
     macrs_option_years::Int = 5
     macrs_bonus_fraction::Float64 = 0.8
     macrs_itc_reduction::Float64 = 0.5
@@ -91,23 +106,45 @@ Base.@kwdef mutable struct GHP <: AbstractGHP
     cooling_thermal_kw::Vector{Float64} = []
     yearly_electric_consumption_kw::Vector{Float64} = []
     peak_combined_heatpump_thermal_ton::Float64 = NaN
+    heat_pump_configuration::String = ""
 
     # Intermediate parameters for cost processing
     tech_sizes_for_cost_curve::Union{Float64, AbstractVector{Float64}} = NaN
     installed_cost_per_kw::Union{Float64, AbstractVector{Float64}} = NaN
-    heatpump_capacity_ton::Float64 = NaN
+    wwhp_heating_pump_installed_cost_curve::Union{Float64, AbstractVector{Float64}} = NaN
+    wwhp_cooling_pump_installed_cost_curve::Union{Float64, AbstractVector{Float64}} = NaN
+    heatpump_capacity_ton::Float64 = 0
+    wwhp_heating_pump_capacity_ton::Float64 = 0
+    wwhp_cooling_pump_capacity_ton::Float64 = 0
 
     # Process and populate these parameters needed more directly by the model
     om_cost_year_one::Float64 = NaN
+
+    # Account for expenses avoided by addition of GHP.
+    avoided_capex_by_ghp_present_value::Float64 = 0.0
 end
 
 
 function GHP(response::Dict, d::Dict)
     ghp = GHP(; ghpghx_response = response, dictkeys_tosymbols(d)...)
+
+    if !(0 <= ghp.aux_cooler_installed_cost_per_ton <= 1.0e6)
+        @error "out of bounds aux_cooler_installed_cost_per_ton"
+    end
+
+    if !(0 <= ghp.aux_heater_installed_cost_per_mmbtu_per_hr <= 1.0e6)
+        @error "out of bounds aux_heater_installed_cost_per_mmbtu_per_hr"
+    end
+
+    if !(1.0 <= ghp.aux_unit_capacity_sizing_factor_on_peak_load <= 5.0)
+        @error "out of bounds aux_unit_capacity_sizing_factor_on_peak_load"
+    end
+
     # Inputs of GhpGhx.jl, which are still needed in REopt
     ghp.heating_thermal_kw = response["inputs"]["heating_thermal_load_mmbtu_per_hr"] * KWH_PER_MMBTU
     ghp.cooling_thermal_kw = response["inputs"]["cooling_thermal_load_ton"] * KWH_THERMAL_PER_TONHOUR
     # Outputs of GhpGhx.jl
+    ghp.heat_pump_configuration = response["outputs"]["heat_pump_configuration"]
     ghp.yearly_electric_consumption_kw = response["outputs"]["yearly_total_electric_consumption_series_kw"]
     ghp.peak_combined_heatpump_thermal_ton = response["outputs"]["peak_combined_heatpump_thermal_ton"]
 
@@ -139,34 +176,71 @@ function setup_installed_cost_curve!(ghp::GHP, response::Dict)
     big_number = 1.0e10
     # GHX and GHP sizing metrics for cost calculations
     total_ghx_ft = response["outputs"]["number_of_boreholes"] * response["outputs"]["length_boreholes_ft"]
-    heatpump_peak_ton = response["outputs"]["peak_combined_heatpump_thermal_ton"]
+
+    if ghp.heat_pump_configuration == "WSHP"
+        heatpump_peak_ton = response["outputs"]["peak_combined_heatpump_thermal_ton"]
+    elseif ghp.heat_pump_configuration == "WWHP"
+        wwhp_heating_pump_peak_ton = response["outputs"]["peak_heating_heatpump_thermal_ton"]
+        wwhp_cooling_pump_peak_ton = response["outputs"]["peak_cooling_heatpump_thermal_ton"]
+    end
 
     # Use initial cost curve to leverage existing incentives-based cost curve method in data_manager
     # The GHX and hydronic loop cost are the y-intercepts ([$]) of the cost for each design
-    ghx_cost = total_ghx_ft * ghp.installed_cost_ghx_per_ft
     hydronic_loop_cost = ghp.building_sqft * ghp.installed_cost_building_hydronic_loop_per_sqft
 
+    if isnothing(ghp.ghx_only_capital_cost)
+        ghp.ghx_only_capital_cost = total_ghx_ft * ghp.installed_cost_ghx_per_ft
+    else
+        @info "Using user provided GHX costs, please validate that this is intentional"
+    end
+
+    aux_heater_cost = 0.0
+    aux_cooler_cost = 0.0
+    if ghp.is_ghx_hybrid
+        aux_heater_cost = ghp.aux_heater_installed_cost_per_mmbtu_per_hr*
+        response["outputs"]["peak_aux_heater_thermal_production_mmbtu_per_hour"]*
+        ghp.aux_unit_capacity_sizing_factor_on_peak_load
+
+        aux_cooler_cost = ghp.aux_cooler_installed_cost_per_ton*
+        response["outputs"]["peak_aux_cooler_thermal_production_ton"]*
+        ghp.aux_unit_capacity_sizing_factor_on_peak_load
+    end
+
     # The DataManager._get_REopt_cost_curve method expects at least a two-point tech_sizes_for_cost_curve to
     #   to use the first value of installed_cost_per_kw as an absolute $ value and
     #   the initial slope is based on the heat pump size (e.g. $/ton) of the cost curve for
     #   building a rebate-based cost curve if there are less-than big_number maximum incentives
     ghp.tech_sizes_for_cost_curve = [0.0, big_number]
-    ghp.installed_cost_per_kw = [ghx_cost + hydronic_loop_cost, 
-                                        ghp.installed_cost_heatpump_per_ton]
+
+    if ghp.heat_pump_configuration == "WSHP"
+        # Use this with the cost curve to determine absolute cost
+        ghp.heatpump_capacity_ton = heatpump_peak_ton * ghp.heatpump_capacity_sizing_factor_on_peak_load
+    elseif ghp.heat_pump_configuration == "WWHP"
+        ghp.wwhp_heating_pump_capacity_ton = wwhp_heating_pump_peak_ton * ghp.heatpump_capacity_sizing_factor_on_peak_load
+        ghp.wwhp_cooling_pump_capacity_ton = wwhp_cooling_pump_peak_ton * ghp.heatpump_capacity_sizing_factor_on_peak_load
+    end
 
     # Using a separate call to _get_REopt_cost_curve in data_manager for "ghp" (not included in "available_techs")
-    #    and then use the value below for heat pump capacity to calculate the final absolute cost for GHP
+    #    and then use the value above for heat pump capacity to calculate the final absolute cost for GHP
+
+    if ghp.heat_pump_configuration == "WSHP"
+        ghp.installed_cost_per_kw = [0, (ghp.ghx_only_capital_cost + hydronic_loop_cost + aux_cooler_cost + aux_heater_cost) / 
+                                        ghp.heatpump_capacity_ton + ghp.installed_cost_heatpump_per_ton]
+    elseif ghp.heat_pump_configuration == "WWHP"
+        # Divide by two to avoid double counting non-heatpump costs
+        ghp.wwhp_heating_pump_installed_cost_curve = [0, (ghp.ghx_only_capital_cost + aux_cooler_cost + aux_heater_cost) / 2 /
+                                                          ghp.wwhp_heating_pump_capacity_ton + ghp.installed_cost_wwhp_heating_pump_per_ton]
+        ghp.wwhp_cooling_pump_installed_cost_curve = [0, (ghp.ghx_only_capital_cost + aux_cooler_cost + aux_heater_cost) / 2 /
+                                                          ghp.wwhp_cooling_pump_capacity_ton + ghp.installed_cost_wwhp_cooling_pump_per_ton]
+    end
 
-    # Use this with the cost curve to determine absolute cost
-    ghp.heatpump_capacity_ton = heatpump_peak_ton * ghp.heatpump_capacity_sizing_factor_on_peak_load
 end
 
 function setup_om_cost!(ghp::GHP)
     # O&M Cost
     ghp.om_cost_year_one = ghp.building_sqft * ghp.om_cost_per_sqft_year
 end
 
-
 function assign_thermal_factor!(d::Dict, heating_or_cooling::String)
     if heating_or_cooling == "space_heating"
         name = "space_heating_efficiency_thermal_factor"

src/core/reopt.jl

@@ -144,7 +144,7 @@ function run_reopt(ms::AbstractArray{T, 1}, p::REoptInputs) where T <: JuMP.Abst
 		Threads.@threads for i = 1:2
 			rs[i] = run_reopt(inputs[i])
 		end
-		if typeof(rs[1]) <: Dict && typeof(rs[2]) <: Dict
+		if typeof(rs[1]) <: Dict && typeof(rs[2]) <: Dict && rs[1]["status"] != "error" && rs[2]["status"] != "error"
 			# TODO when a model is infeasible the JuMP.Model is returned from run_reopt (and not the results Dict)
 			results_dict = combine_results(p, rs[1], rs[2], bau_inputs.s)
 			results_dict["Financial"] = merge(results_dict["Financial"], proforma_results(p, results_dict))
@@ -153,7 +153,7 @@ function run_reopt(ms::AbstractArray{T, 1}, p::REoptInputs) where T <: JuMP.Abst
 			end
 			return results_dict
 		else
-			return rs
+			throw(@error("REopt scenarios solved either with errors of non-optimal solutions."))
 		end
 	catch e
 		if isnothing(e) # Error thrown by REopt
@@ -244,7 +244,9 @@ function build_reopt!(m::JuMP.AbstractModel, p::REoptInputs)
 	m[:TotalCHPStandbyCharges] = 0
 	m[:OffgridOtherCapexAfterDepr] = 0.0
     m[:GHPCapCosts] = 0.0
-    m[:GHPOMCosts] = 0.0   
+    m[:GHPOMCosts] = 0.0
+	m[:AvoidedCapexByGHP] = 0.0
+	m[:ResidualGHXCapCost] = 0.0
 
 	if !isempty(p.techs.all)
 		add_tech_size_constraints(m, p)
@@ -446,7 +448,10 @@ function build_reopt!(m::JuMP.AbstractModel, p::REoptInputs)
 		p.s.financial.offgrid_other_annual_costs * p.pwf_om * (1 - p.s.financial.owner_tax_rate_fraction) +
 
 		# Additional capital costs, depreciable (only applies when `off_grid_flag` is true)
-		m[:OffgridOtherCapexAfterDepr]
+		m[:OffgridOtherCapexAfterDepr] -
+
+		# Subtract capital expenditures avoided by inclusion of GHP and residual present value of GHX.
+		m[:AvoidedCapexByGHP] - m[:ResidualGHXCapCost]
 
 	);
 	if !isempty(p.s.electric_utility.outage_durations)