Merge pull request #272 from NREL/centralized-ghp

Centralized ghp, hybrid ghp, ghp residual value

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)

src/core/reopt_inputs.jl

@@ -115,6 +115,9 @@ struct REoptInputs{ScenarioType <: AbstractScenario} <: AbstractInputs
     ghp_electric_consumption_kw::Array{Float64,2}  # Array of electric load profiles consumed by GHP
     ghp_installed_cost::Array{Float64,1}  # Array of installed cost for GHP options
     ghp_om_cost_year_one::Array{Float64,1}  # Array of O&M cost for GHP options
+    avoided_capex_by_ghp_present_value::Array{Float64,1} # HVAC upgrade costs avoided
+    ghx_useful_life_years::Array{Float64,1} # GHX useful life years
+    ghx_residual_value::Array{Float64,1} # Residual value of each GHX options
     tech_renewable_energy_fraction::Dict{String, <:Real} # (techs)
     tech_emissions_factors_CO2::Dict{String, <:Real} # (techs)
     tech_emissions_factors_NOx::Dict{String, <:Real} # (techs)
@@ -175,7 +178,8 @@ function REoptInputs(s::AbstractScenario)
     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(s, time_steps, 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(s, time_steps, time_steps_without_grid)
 
     if any(pv.existing_kw > 0 for pv in s.pvs)
         adjust_load_profile(s, production_factor)
@@ -232,6 +236,9 @@ function REoptInputs(s::AbstractScenario)
         ghp_electric_consumption_kw,
         ghp_installed_cost,
         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, 
@@ -974,26 +981,43 @@ function setup_ghp_inputs(s::AbstractScenario, time_steps, time_steps_without_gr
     cooling_thermal_load_reduction_with_ghp_kw = zeros(num, length(time_steps))
     ghp_cooling_thermal_load_served_kw = zeros(num, length(time_steps))        
     ghp_electric_consumption_kw = zeros(num, length(time_steps))
+    avoided_capex_by_ghp_present_value = Vector{Float64}(undef, num)
+    ghx_useful_life_years = Vector{Float64}(undef, num)
+    ghx_residual_value = Vector{Float64}(undef, num)
     if num > 0
         require_ghp_purchase = s.ghp_option_list[1].require_ghp_purchase  # This does not change with the number of options
+
         for (i, option) in enumerate(s.ghp_option_list)
-            ghp_cap_cost_slope, ghp_cap_cost_x, ghp_cap_cost_yint, ghp_n_segments = cost_curve(option, s.financial)
-            ghp_size_ton = option.heatpump_capacity_ton
-            seg = 0
-            if ghp_size_ton <= ghp_cap_cost_x[1]
-                seg = 1
-            elseif ghp_size_ton > ghp_cap_cost_x[end]
-                seg = ghp_n_segments
-            else
-                for n in 2:(ghp_n_segments+1)
-                    if (ghp_size_ton > ghp_cap_cost_x[n-1]) && (ghp_size_ton <= ghp_cap_cost_x[n])
-                        seg = n
-                        break
-                    end
-                end
+            if option.heat_pump_configuration == "WSHP"
+                fixed_cost, variable_cost = get_ghp_installed_cost(option, s.financial, option.heatpump_capacity_ton)
+                ghp_installed_cost[i] = fixed_cost + variable_cost
+
+            elseif option.heat_pump_configuration == "WWHP"
+                temp = option.installed_cost_per_kw
+                option.installed_cost_per_kw = option.wwhp_heating_pump_installed_cost_curve
+                fixed_cost_heating, variable_cost_heating = get_ghp_installed_cost(option, s.financial, option.wwhp_heating_pump_capacity_ton)
+                ghp_installed_cost_heating = 0.5 * fixed_cost_heating + variable_cost_heating
+
+                option.installed_cost_per_kw = option.wwhp_cooling_pump_installed_cost_curve
+                fixed_cost_cooling, variable_cost_cooling = get_ghp_installed_cost(option, s.financial, option.wwhp_cooling_pump_capacity_ton)
+                option.installed_cost_per_kw = temp
+                ghp_installed_cost_cooling = 0.5 * fixed_cost_cooling + variable_cost_cooling
+
+                ghp_installed_cost[i] = ghp_installed_cost_heating + ghp_installed_cost_cooling
             end
-            ghp_installed_cost[i] = ghp_cap_cost_yint[seg-1] + ghp_size_ton * ghp_cap_cost_slope[seg-1]
+
             ghp_om_cost_year_one[i] = option.om_cost_year_one
+            avoided_capex_by_ghp_present_value[i] = option.avoided_capex_by_ghp_present_value
+            ghx_useful_life_years[i] = option.ghx_useful_life_years
+            # ownership guided residual value determination
+            discount_rate = (1 - 1*s.financial.third_party_ownership)*s.financial.offtaker_discount_rate_fraction + s.financial.third_party_ownership*s.financial.owner_discount_rate_fraction
+            ghx_residual_value[i] = option.ghx_only_capital_cost*
+            (
+                (option.ghx_useful_life_years - s.financial.analysis_years)/option.ghx_useful_life_years
+            )/(
+                (1 + discount_rate)^s.financial.analysis_years
+            )
+
             heating_thermal_load = s.space_heating_load.loads_kw + s.dhw_load.loads_kw
             # Using minimum of thermal load and ghp-serving load to avoid small negative net loads
             for j in time_steps
@@ -1021,7 +1045,8 @@ function setup_ghp_inputs(s::AbstractScenario, time_steps, time_steps_without_gr
     return 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
+    ghp_installed_cost, ghp_om_cost_year_one, avoided_capex_by_ghp_present_value,
+    ghx_useful_life_years, ghx_residual_value
 end
 
 function setup_operating_reserve_fraction(s::AbstractScenario, techs_operating_reserve_req_fraction)
@@ -1033,4 +1058,26 @@ function setup_operating_reserve_fraction(s::AbstractScenario, techs_operating_r
     techs_operating_reserve_req_fraction["Wind"] = s.wind.operating_reserve_required_fraction
 
     return nothing
+end
+
+function get_ghp_installed_cost(option::AbstractTech, financial::Financial, ghp_size_ton::Float64)
+
+    ghp_cap_cost_slope, ghp_cap_cost_x, ghp_cap_cost_yint, ghp_n_segments = cost_curve(option, financial)
+    seg = 0
+    if ghp_size_ton <= ghp_cap_cost_x[1]
+        seg = 2
+    elseif ghp_size_ton > ghp_cap_cost_x[end]
+        seg = ghp_n_segments+1
+    else
+        for n in 2:(ghp_n_segments+1)
+            if (ghp_size_ton > ghp_cap_cost_x[n-1]) && (ghp_size_ton <= ghp_cap_cost_x[n])
+                seg = n
+                break
+            end
+        end
+    end
+    fixed_cost = ghp_cap_cost_yint[seg-1] 
+    variable_cost = ghp_size_ton * ghp_cap_cost_slope[seg-1]
+
+    return fixed_cost, variable_cost
 end