Add Optional Costs for Existing Boiler and Chiller

CHANGELOG.md

@@ -30,6 +30,14 @@ Classify the change according to the following categories:
 - Memory-clearing commands after each JuMP model instance in `runtests.jl` to avoid memory buildup which were slowing down Actions test job
 - Added back `ubuntu` OS as an additional runner OS for the tests Action job, now that memory buildup is reduced (removed a year ago due to memory crashing the runner)
 
+## hvac-costs
+### Added
+- Add `installed_cost...` for `ExistingBoiler` and `ExistingChiller` which is incurred in the BAU scenario, and may be avoided with other heating and cooling technologies in the Optimal scenario.
+
+# non-hourly-fuel-cost
+### Fixed
+- Fixed handling of non-hourly (e.g. 15-min interval) fuel cost
+
 ## v0.52.0
 ### Added
 - Add **Financial** inputs `min_initial_capital_costs_before_incentives` and `max_initial_capital_costs_before_incentives` which, when provided, provide upper and lower bounds on initial capital costs for all technologies.
@@ -38,10 +46,6 @@ Classify the change according to the following categories:
 ### Fixed
 - Fix implementation of production-based incentives
 
-# non-hourly-fuel-cost
-### Fixed
-- Fixed handling of non-hourly (e.g. 15-min interval) fuel cost
-
 ## v0.51.1
 ### Added
 - Added the following output fields: `year_one_fuel_cost_after_tax` for `ExistingBoiler`, `CHP`, `Generator`, and `Boiler`; `ElectricTariff`: `year_one_bill_after_tax` and `year_one_export_benefit_after_tax`, `Financial`: `capital_costs_after_non_discounted_incentives`, `year_one_total_operating_cost_savings_before_tax`, `year_one_total_operating_cost_savings_after_tax`, `year_one_total_operating_cost_before_tax`, `year_one_total_operating_cost_after_tax`, `year_one_fuel_cost_before_tax`, `year_one_fuel_cost_after_tax`, `year_one_chp_standby_cost_after_tax`, `year_one_chp_standby_cost_after_tax`, `GHP.avoided_capex_by_ghp_present_value`

src/constraints/tech_constraints.jl

@@ -40,4 +40,4 @@ function add_no_curtail_constraints(m, p; _n="")
             fix(m[Symbol("dvCurtail"*_n)][t, ts] , 0.0, force=true)
         end
     end
-end
+end

src/constraints/thermal_tech_constraints.jl

@@ -263,3 +263,35 @@ function no_existing_chiller_production(m, p; _n="")
     end
     fix(m[Symbol("dvSize"*_n)]["ExistingChiller"], 0.0, force=true)
 end
+
+function add_existing_boiler_capex_constraints(m, p; _n="")
+    # @variable(m, binExistingBoiler, Int, lower_bound = 0, upper_bound = 1)  # This is same as below with Bin
+    @variable(m, binExistingBoiler, Bin)
+    # If still using ExistingBoiler in optimal case at all, incur costs (not scaled by size)
+    # Force dvSize["ExistingBoiler] to zero if binExistingBoiler is zero:
+    @constraint(m, ExistingBoilerCostCon, m[Symbol("dvSize"*_n)]["ExistingBoiler"] <= m[Symbol("binExistingBoiler"*_n)] * BIG_NUMBER)
+
+    if p.s.existing_boiler.retire_in_optimal
+        @constraint(m, ExistingBoilerSelect, m[Symbol("binExistingBoiler"*_n)] == 0)
+    else
+        @constraint(m, ExistingBoilerSelect, m[Symbol("binExistingBoiler"*_n)] <= 1)
+    end
+
+    m[:ExistingBoilerCost] = @expression(m, p.third_party_factor *
+        sum(p.s.existing_boiler.installed_cost_dollars * m[Symbol("binExistingBoiler"*_n)])
+    )
+end
+
+function add_existing_chiller_capex_constraints(m, p; _n="")
+    # @variable(m, binExistingChiller, Int, lower_bound = 0, upper_bound = 1)  # This is same as below with Bin
+    @variable(m, binExistingChiller, Bin)
+    # If still using ExistingChiller in optimal case, incur costs (not scaled by size)
+    # Force dvSize["ExistingChiller] to zero if binExistingChiller is zero:
+    @constraint(m, ExistingChillerCostCon, m[Symbol("dvSize"*_n)]["ExistingChiller"] <= m[Symbol("binExistingChiller"*_n)] * BIG_NUMBER)
+
+    @constraint(m, ExistingChillerSelect, m[Symbol("binExistingChiller"*_n)] <= 1)
+
+    m[:ExistingChillerCost] = @expression(m, p.third_party_factor *
+        sum(p.s.existing_chiller.installed_cost_dollars * m[Symbol("binExistingChiller"*_n)])
+    )
+end

src/core/existing_boiler.jl

@@ -7,6 +7,9 @@ const existing_boiler_efficiency_defaults = Dict(
 struct ExistingBoiler <: AbstractThermalTech  # useful to create AbstractHeatingTech or AbstractThermalTech?
     max_kw::Real
     production_type::String
+    max_thermal_factor_on_peak_load::Real    
+    installed_cost_per_kw::Real
+    installed_cost_dollars::Real 
     efficiency::Real
     fuel_cost_per_mmbtu::Union{<:Real, AbstractVector{<:Real}}
     fuel_type::String
@@ -29,6 +32,8 @@ end
     max_heat_demand_kw::Real=0, # Auto-populated based on SpaceHeatingLoad and DomesticHotWaterLoad inputs
     production_type::String = "hot_water", # Can be "steam" or "hot_water"
     max_thermal_factor_on_peak_load::Real = 1.25,
+    installed_cost_per_mmbtu_per_hour::Real = 0.0  # Represents needed CapEx in BAU, assuming net present value basis; cost is scaled to the size of boiler needed
+    installed_cost_dollars::Real = 0.0  # Represents needed CapEx in BAU, assuming net present cost basis; also incurred in Optimal case if still using at all
     efficiency::Real = NaN, # Existing boiler system efficiency - conversion of fuel to usable heating thermal energy. See note below.
     fuel_cost_per_mmbtu::Union{<:Real, AbstractVector{<:Real}} = [], # REQUIRED. Can be a scalar, a list of 12 monthly values, or a time series of values for every time step
     fuel_type::String = "natural_gas", # "restrict_to": ["natural_gas", "landfill_bio_gas", "propane", "diesel_oil"]
@@ -71,6 +76,8 @@ function ExistingBoiler(;
     max_heat_demand_kw::Real=0,
     production_type::String = "hot_water",
     max_thermal_factor_on_peak_load::Real = 1.25,
+    installed_cost_per_mmbtu_per_hour::Real = 0.0,
+    installed_cost_dollars::Real = 0.0,
     efficiency::Real = NaN,
     fuel_cost_per_mmbtu::Union{<:Real, AbstractVector{<:Real}} = [], # REQUIRED. Can be a scalar, a list of 12 monthly values, or a time series of values for every time step
     fuel_type::String = "natural_gas", # "restrict_to": ["natural_gas", "landfill_bio_gas", "propane", "diesel_oil"]
@@ -97,11 +104,21 @@ function ExistingBoiler(;
         efficiency = existing_boiler_efficiency_defaults[production_type]
     end
 
-    max_kw = max_heat_demand_kw * max_thermal_factor_on_peak_load
+    max_kw = max_heat_demand_kw * max_thermal_factor_on_peak_load  # This is really the **actual** size in BAU
+
+    installed_cost_per_kw = 0.0
+    if !(installed_cost_per_mmbtu_per_hour == 0.0) && (installed_cost_dollars == 0.0)
+        installed_cost_per_kw = installed_cost_per_mmbtu_per_hour / KWH_PER_MMBTU * max_thermal_factor_on_peak_load
+    elseif !(installed_cost_per_mmbtu_per_hour == 0.0) && !(installed_cost_dollars == 0.0)
+        throw(@error("A non-zero value for both installed_cost_per_mmbtu_per_hour and installed_cost_dollars was input for ExistingBoiler; only provide one or the other"))
+    end
 
     ExistingBoiler(
         max_kw,
         production_type,
+        max_thermal_factor_on_peak_load,
+        installed_cost_per_kw,
+        installed_cost_dollars,
         efficiency,
         fuel_cost_per_mmbtu,
         fuel_type,

src/core/existing_chiller.jl

@@ -5,11 +5,14 @@
     loads_kw_thermal::Vector{<:Real},
     cop::Union{Real, Nothing} = nothing,
     max_thermal_factor_on_peak_load::Real=1.25
+    installed_cost_per_kw::Real = 0.0  # Represents needed CapEx in BAU, assuming net present value basis based on current size; also incurred in Optimal case if still using at all
+    installed_cost_dollars::Real = 0.0  # Represents needed CapEx in BAU, assuming net present cost basis; also incurred in Optimal case if still using at all
     retire_in_optimal::Bool = false  # Do NOT use in the optimal case (still used in BAU)
 ```
 
 !!! note "Max ExistingChiller size" 
-    The maximum size [kW] of the `ExistingChiller` will be set based on the peak thermal load as follows:
+    The maximum size [kW] of the `ExistingChiller` will be set based on the peak thermal load as follows, and
+       this is really the **actual** estimated size of the existing chiller at the site:
     ```julia 
     max_kw = maximum(loads_kw_thermal) * max_thermal_factor_on_peak_load
     ```
@@ -18,6 +21,8 @@ struct ExistingChiller <: AbstractThermalTech
     max_kw::Real
     cop::Union{Real, Nothing}
     max_thermal_factor_on_peak_load::Real
+    installed_cost_per_kw::Real
+    installed_cost_dollars::Real
     retire_in_optimal::Bool
 end
 
@@ -26,13 +31,25 @@ function ExistingChiller(;
         loads_kw_thermal::Vector{<:Real},
         cop::Union{Real, Nothing} = nothing,
         max_thermal_factor_on_peak_load::Real=1.25,
+        installed_cost_per_ton::Real = 0.0,
+        installed_cost_dollars::Real = 0.0,
         retire_in_optimal::Bool = false
     )
-    max_kw = maximum(loads_kw_thermal) * max_thermal_factor_on_peak_load
+    max_kw = maximum(loads_kw_thermal) * max_thermal_factor_on_peak_load  # This is really the **actual** size in BAU
+
+    installed_cost_per_kw = 0.0
+    if !(installed_cost_per_ton == 0.0) && (installed_cost_dollars == 0.0)
+        installed_cost_per_kw = installed_cost_per_ton / KWH_THERMAL_PER_TONHOUR * max_thermal_factor_on_peak_load
+    elseif !(installed_cost_per_ton == 0.0) && !(installed_cost_dollars == 0.0)
+        throw(@error("A non-zero value for both installed_cost_per_ton and installed_cost_dollars was input for ExistingChiller; only provide one or the other"))
+    end
+
     ExistingChiller(
         max_kw,
         cop,
         max_thermal_factor_on_peak_load,
+        installed_cost_per_kw,
+        installed_cost_dollars,
         retire_in_optimal
     )
 end

src/core/reopt.jl

@@ -266,6 +266,8 @@ function build_reopt!(m::JuMP.AbstractModel, p::REoptInputs)
 	m[:AvoidedCapexByASHP] = 0.0
 	m[:ResidualGHXCapCost] = 0.0
 	m[:ObjectivePenalties] = 0.0
+	m[:ExistingBoilerCost] = 0.0
+	m[:ExistingChillerCost] = 0.0
 
 	if !isempty(p.techs.all) || !isempty(p.techs.ghp)
 		if !isempty(p.techs.all)
@@ -309,10 +311,16 @@ function build_reopt!(m::JuMP.AbstractModel, p::REoptInputs)
             add_boiler_tech_constraints(m, p)
 			m[:TotalPerUnitProdOMCosts] += m[:TotalBoilerPerUnitProdOMCosts]
 			m[:TotalFuelCosts] += m[:TotalBoilerFuelCosts]
+			if ("ExistingBoiler" in p.techs.boiler) && (p.s.existing_boiler.installed_cost_dollars > 0.0)
+				add_existing_boiler_capex_constraints(m, p)
+			end			
         end
 
 		if !isempty(p.techs.cooling)
             add_cooling_tech_constraints(m, p)
+			if ("ExistingChiller" in p.techs.cooling) && (p.s.existing_chiller.installed_cost_dollars > 0.0)
+				add_existing_chiller_capex_constraints(m, p)
+			end
         end
 
 		# Zero out ExistingChiller production if retire_in_optimal; setdiff avoids zeroing for BAU 
@@ -554,6 +562,14 @@ function build_reopt!(m::JuMP.AbstractModel, p::REoptInputs)
 			for s in p.s.electric_utility.scenarios, tz in p.s.electric_utility.outage_start_time_steps)
 	end
 
+	if "ExistingBoiler" in p.techs.all && (p.s.existing_boiler.installed_cost_dollars > 0.0)
+		add_to_expression!(Costs, m[:ExistingBoilerCost])
+	end
+
+	if "ExistingChiller" in p.techs.all && (p.s.existing_chiller.installed_cost_dollars > 0.0)
+		add_to_expression!(Costs, m[:ExistingChillerCost])
+	end
+
 	# Set model objective 
 	@objective(m, Min, m[:Costs] + m[:ObjectivePenalties] )
 

src/core/reopt_inputs.jl

@@ -723,7 +723,7 @@ function setup_existing_boiler_inputs(s::AbstractScenario, max_sizes, min_sizes,
     max_sizes["ExistingBoiler"] = s.existing_boiler.max_kw
     min_sizes["ExistingBoiler"] = 0.0
     existing_sizes["ExistingBoiler"] = 0.0
-    cap_cost_slope["ExistingBoiler"] = 0.0
+    cap_cost_slope["ExistingBoiler"] = s.existing_boiler.installed_cost_per_kw
     boiler_efficiency["ExistingBoiler"] = s.existing_boiler.efficiency
     # om_cost_per_kw["ExistingBoiler"] = 0.0
     tech_renewable_energy_fraction["ExistingBoiler"] = s.existing_boiler.fuel_renewable_energy_fraction
@@ -790,7 +790,7 @@ function setup_existing_chiller_inputs(s::AbstractScenario, max_sizes, min_sizes
     max_sizes["ExistingChiller"] = s.existing_chiller.max_kw
     min_sizes["ExistingChiller"] = 0.0
     existing_sizes["ExistingChiller"] = 0.0
-    cap_cost_slope["ExistingChiller"] = 0.0
+    cap_cost_slope["ExistingChiller"] = s.existing_chiller.installed_cost_per_kw
     cooling_cop["ExistingChiller"] .= s.existing_chiller.cop
     cooling_cf["ExistingChiller"]  = ones(8760*s.settings.time_steps_per_hour)
     # om_cost_per_kw["ExistingChiller"] = 0.0

src/core/reopt_multinode.jl

@@ -86,6 +86,12 @@ function add_variables!(m::JuMP.AbstractModel, ps::AbstractVector{REoptInputs{T}
 		m[Symbol(ex_name)] = 0
 
 		add_elec_utility_expressions(m, p; _n=_n)
+
+		# Existing HVAC Cost NOT hooked up for multi-node; this is a patch to avoid errors
+		ex_name = "ExistingBoilerCost"*_n
+		m[Symbol(ex_name)] = 0
+		ex_name = "ExistingChillerCost"*_n
+		m[Symbol(ex_name)] = 0
 
 		#################################  Objective Function   ########################################
 		m[Symbol("Costs"*_n)] = @expression(m,

src/core/scenario.jl

@@ -214,7 +214,7 @@ function Scenario(d::Dict; flex_hvac_from_json=false)
         generator = Generator(; max_kw=0)
     end
 
-    max_heat_demand_kw = 0.0
+    total_heating_load_series_kw = zeros(8760 * settings.time_steps_per_hour)
 
     if haskey(d, "DomesticHotWaterLoad") && !haskey(d, "FlexibleHVAC")
         add_doe_reference_names_from_elec_to_thermal_loads(d["ElectricLoad"], d["DomesticHotWaterLoad"])
@@ -227,7 +227,7 @@ function Scenario(d::Dict; flex_hvac_from_json=false)
                                         time_steps_per_hour=settings.time_steps_per_hour,
                                         existing_boiler_efficiency = existing_boiler_efficiency
                                         )
-        max_heat_demand_kw = maximum(dhw_load.loads_kw)
+        total_heating_load_series_kw .+= dhw_load.loads_kw
     else
         dhw_load = HeatingLoad(;
             load_type = "domestic_hot_water", 
@@ -249,7 +249,7 @@ function Scenario(d::Dict; flex_hvac_from_json=false)
                                             time_steps_per_hour=settings.time_steps_per_hour,
                                             existing_boiler_efficiency = existing_boiler_efficiency
                                             )
-        max_heat_demand_kw = maximum(space_heating_load.loads_kw .+ max_heat_demand_kw)
+        total_heating_load_series_kw .+= space_heating_load.loads_kw
     else
         space_heating_load = HeatingLoad(; 
             load_type = "space_heating",        
@@ -271,7 +271,7 @@ function Scenario(d::Dict; flex_hvac_from_json=false)
                                             existing_boiler_efficiency = existing_boiler_efficiency                                           
                                             )
 
-        max_heat_demand_kw = maximum(process_heat_load.loads_kw .+ max_heat_demand_kw)
+        total_heating_load_series_kw .+= process_heat_load.loads_kw
     else
         process_heat_load = HeatingLoad(;
                 load_type = "process_heat",                
@@ -354,6 +354,7 @@ function Scenario(d::Dict; flex_hvac_from_json=false)
         end
     end
 
+    max_heat_demand_kw = maximum(total_heating_load_series_kw)
     if max_heat_demand_kw > 0 && !haskey(d, "FlexibleHVAC")  # create ExistingBoiler
         boiler_inputs = Dict{Symbol, Any}()
         boiler_inputs[:max_heat_demand_kw] = max_heat_demand_kw

src/results/existing_boiler.jl

@@ -20,7 +20,9 @@
 """
 function add_existing_boiler_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
     r = Dict{String, Any}()
-    r["size_mmbtu_per_hour"] = round(value(m[Symbol("dvSize"*_n)]["ExistingBoiler"]) / KWH_PER_MMBTU, digits=3)
+    max_prod_kw = maximum(value.(sum(m[Symbol("dvHeatingProduction"*_n)]["ExistingBoiler",q,:] for q in p.heating_loads)))
+    size_actual_mmbtu_per_hour = max_prod_kw / KWH_PER_MMBTU * p.s.existing_boiler.max_thermal_factor_on_peak_load
+    r["size_mmbtu_per_hour"] = round(size_actual_mmbtu_per_hour, digits=3)
 	r["fuel_consumption_series_mmbtu_per_hour"] = 
         round.(value.(m[:dvFuelUsage]["ExistingBoiler", ts] for ts in p.time_steps) ./ KWH_PER_MMBTU, digits=5)
     r["annual_fuel_consumption_mmbtu"] = round(sum(r["fuel_consumption_series_mmbtu_per_hour"]), digits=5)

src/results/existing_chiller.jl

@@ -11,6 +11,8 @@
 function add_existing_chiller_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
     r = Dict{String, Any}()
 
+	r["size_ton"] = round(value(m[Symbol("dvSize"*_n)]["ExistingChiller"]) * p.s.existing_chiller.max_thermal_factor_on_peak_load / KWH_THERMAL_PER_TONHOUR, digits=3)
+
 	@expression(m, ELECCHLtoTES[ts in p.time_steps],
 		sum(m[:dvProductionToStorage][b,"ExistingChiller",ts] for b in p.s.storage.types.cold)
     )