Merge pull request #473 from NREL/develop

Jan 2025 Updates: Load Year Fixes, ASHP Max Dispatch

CHANGELOG.md

@@ -26,6 +26,17 @@ Classify the change according to the following categories:
     ### Removed
 
 
+## develop
+### Added
+- New parameter `force_dispatch` in the **ASHPSpaceHeater** and **ASHPWaterHeater** technologies (default = `true`).  When kept at `true`, the ASHP's thermal output will be the minimum of the site load(s) served and the system size (adjusted for timestep-specific capacity factor) in each period. If set to `false`, ASHP will do economic dispatch considering COP and CF along with electricity prices.
+### Fixed
+- Align heating and cooling load profiles based on electric load year input, if using custom electric load profile with simulated (CRB or schedule-based flatloads) heating/cooling loads
+- Handling of leap years for `ElectricLoad.loads_kw` inputs to align with URDB rate structures
+### Changed
+- Make `year` input required with any custom load profile input (e.g. `ElectricLoad.loads_kw`, `SpaceHeatingLoad.fuel_loads_mmbtu_per_hour`)
+- Shift and adjust CRB load profiles (i.e. with `doe_reference_name` input) based on the `year` input
+
+
 ## v0.49.1
 ### Changed
 - Swap an error for a warning with inconsistent load-year between electric and heating; soon to 

src/constraints/thermal_tech_constraints.jl

@@ -107,11 +107,75 @@ end
 
 
 function add_ashp_force_in_constraints(m, p; _n="")
-    if "ASHPSpaceHeater" in p.techs.ashp && p.s.ashp.force_into_system
-        for t in setdiff(p.techs.can_serve_space_heating, ["ASHPSpaceHeater"])
-            for ts in p.time_steps
-                fix(m[Symbol("dvHeatingProduction"*_n)][t,"SpaceHeating",ts], 0.0, force=true)
-                fix(m[Symbol("dvProductionToWaste"*_n)][t,"SpaceHeating",ts], 0.0, force=true)
+    if "ASHPSpaceHeater" in p.techs.ashp 
+        if p.s.ashp.force_into_system
+            for t in setdiff(p.techs.can_serve_space_heating, ["ASHPSpaceHeater"])
+                for ts in p.time_steps
+                    fix(m[Symbol("dvHeatingProduction"*_n)][t,"SpaceHeating",ts], 0.0, force=true)
+                    fix(m[Symbol("dvProductionToWaste"*_n)][t,"SpaceHeating",ts], 0.0, force=true)
+                end
+            end
+        elseif p.s.ashp.force_dispatch
+            dv = "binASHPSHSizeExceedsThermalLoad"*_n
+            m[Symbol(dv)] = @variable(m, [p.time_steps], binary=true, base_name=dv)
+            dv = "dvASHPSHSizeTimesExcess"*_n
+            m[Symbol(dv)] = @variable(m, [p.time_steps], lower_bound=0, base_name=dv)
+            if p.s.ashp.can_serve_cooling
+                max_sh_size_bigM = 2*max(p.max_sizes["ASHPSpaceHeater"], maximum(p.heating_loads_kw["SpaceHeating"] ./ p.heating_cf["ASHPSpaceHeater"])+maximum(p.s.cooling_load.loads_kw_thermal ./ p.cooling_cf["ASHPSpaceHeater"]))
+                @constraint(m, [ts in p.time_steps],
+                    m[Symbol("binASHPSHSizeExceedsThermalLoad"*_n)][ts] >= (
+                        m[Symbol("dvSize"*_n)]["ASHPSpaceHeater"]
+                        - (p.heating_loads_kw["SpaceHeating"][ts] / p.heating_cf["ASHPSpaceHeater"][ts]) 
+                        - (p.s.cooling_load.loads_kw_thermal[ts] / p.cooling_cf["ASHPSpaceHeater"][ts])  
+                        ) / max_sh_size_bigM
+                )
+                @constraint(m, [ts in p.time_steps],
+                    m[Symbol("binASHPSHSizeExceedsThermalLoad"*_n)][ts] <= 1 - (
+                        (p.heating_loads_kw["SpaceHeating"][ts] / p.heating_cf["ASHPSpaceHeater"][ts]) 
+                        + (p.s.cooling_load.loads_kw_thermal[ts] / p.cooling_cf["ASHPSpaceHeater"][ts]) 
+                        - m[Symbol("dvSize"*_n)]["ASHPSpaceHeater"]
+                        ) / max_sh_size_bigM
+                )
+                # set dvASHPSHSizeTimesExcess = binASHPSHSizeExceedsThermalLoad * dvSize
+                # big-M is min CF times heat load
+                @constraint(m, [ts in p.time_steps],
+                    m[Symbol("dvASHPSHSizeTimesExcess"*_n)][ts] >= m[Symbol("dvSize"*_n)]["ASHPSpaceHeater"] - max_sh_size_bigM * (1-m[Symbol("binASHPSHSizeExceedsThermalLoad"*_n)][ts])  
+                )
+                @constraint(m, [ts in p.time_steps],
+                    m[Symbol("dvASHPSHSizeTimesExcess"*_n)][ts] <= m[Symbol("dvSize"*_n)]["ASHPSpaceHeater"]
+                )
+                @constraint(m, [ts in p.time_steps],
+                    m[Symbol("dvASHPSHSizeTimesExcess"*_n)][ts] <= max_sh_size_bigM * m[Symbol("binASHPSHSizeExceedsThermalLoad"*_n)][ts]
+                )
+                #Enforce dispatch: output = system size - (overage)
+                @constraint(m, [ts in p.time_steps],
+                    m[Symbol("dvHeatingProduction"*_n)]["ASHPSpaceHeater","SpaceHeating",ts] / p.heating_cf["ASHPSpaceHeater"][ts] + m[Symbol("dvCoolingProduction"*_n)]["ASHPSpaceHeater",ts] / p.cooling_cf["ASHPSpaceHeater"][ts] >= m[Symbol("dvSize"*_n)]["ASHPSpaceHeater"] - m[Symbol("dvASHPSHSizeTimesExcess"*_n)][ts] + (p.heating_loads_kw["SpaceHeating"][ts] / p.heating_cf["ASHPSpaceHeater"][ts] + p.s.cooling_load.loads_kw_thermal[ts] / p.cooling_cf["ASHPSpaceHeater"][ts] ) * m[Symbol("binASHPSHSizeExceedsThermalLoad"*_n)][ts]
+                )
+            else
+                # binary variable enforcement for size >= load
+                max_sh_size_bigM = 2*max(p.max_sizes["ASHPSpaceHeater"], maximum(p.heating_loads_kw["SpaceHeating"] ./ p.heating_cf["ASHPSpaceHeater"]))
+                @constraint(m, [ts in p.time_steps],
+                    m[Symbol("binASHPSHSizeExceedsThermalLoad"*_n)][ts] >= (m[Symbol("dvSize"*_n)]["ASHPSpaceHeater"] - p.heating_loads_kw["SpaceHeating"][ts] / p.heating_cf["ASHPSpaceHeater"][ts]) / max_sh_size_bigM
+                )
+                @constraint(m, [ts in p.time_steps],
+                    m[Symbol("binASHPSHSizeExceedsThermalLoad"*_n)][ts] <= 1 - (p.heating_loads_kw["SpaceHeating"][ts] / p.heating_cf["ASHPSpaceHeater"][ts] - m[Symbol("dvSize"*_n)]["ASHPSpaceHeater"]) / max_sh_size_bigM
+                )
+                # set dvASHPSHSizeTimesExcess = binASHPSHSizeExceedsThermalLoad * dvSize
+                # big-M is min CF times heat load
+
+                @constraint(m, [ts in p.time_steps],
+                    m[Symbol("dvASHPSHSizeTimesExcess"*_n)][ts] >= p.heating_cf["ASHPSpaceHeater"][ts]*m[Symbol("dvSize"*_n)]["ASHPSpaceHeater"] - max_sh_size_bigM * (1-m[Symbol("binASHPSHSizeExceedsThermalLoad"*_n)][ts])  
+                )
+                @constraint(m, [ts in p.time_steps],
+                    m[Symbol("dvASHPSHSizeTimesExcess"*_n)][ts] <= p.heating_cf["ASHPSpaceHeater"][ts]*m[Symbol("dvSize"*_n)]["ASHPSpaceHeater"]
+                )
+                @constraint(m, [ts in p.time_steps],
+                    m[Symbol("dvASHPSHSizeTimesExcess"*_n)][ts] <= max_sh_size_bigM * m[Symbol("binASHPSHSizeExceedsThermalLoad"*_n)][ts]
+                )
+                #Enforce dispatch: output = system size - (overage)
+                @constraint(m, [ts in p.time_steps],
+                    m[Symbol("dvHeatingProduction"*_n)]["ASHPSpaceHeater","SpaceHeating",ts] >= p.heating_cf["ASHPSpaceHeater"][ts]*m[Symbol("dvSize"*_n)]["ASHPSpaceHeater"] - m[Symbol("dvASHPSHSizeTimesExcess"*_n)][ts] + p.heating_loads_kw["SpaceHeating"][ts] * m[Symbol("binASHPSHSizeExceedsThermalLoad"*_n)][ts]
+                )
             end
         end
     end
@@ -124,13 +188,44 @@ function add_ashp_force_in_constraints(m, p; _n="")
         end
     end
 
-    if "ASHPWaterHeater" in p.techs.ashp && p.s.ashp_wh.force_into_system
-        for t in setdiff(p.techs.can_serve_dhw, ["ASHPWaterHeater"])
-            for ts in p.time_steps
-                fix(m[Symbol("dvHeatingProduction"*_n)][t,"DomesticHotWater",ts], 0.0, force=true)
-                fix(m[Symbol("dvProductionToWaste"*_n)][t,"DomesticHotWater",ts], 0.0, force=true)
+    if "ASHPWaterHeater" in p.techs.ashp 
+        if p.s.ashp_wh.force_into_system
+            for t in setdiff(p.techs.can_serve_dhw, ["ASHPWaterHeater"])
+                for ts in p.time_steps
+                    fix(m[Symbol("dvHeatingProduction"*_n)][t,"DomesticHotWater",ts], 0.0, force=true)
+                    fix(m[Symbol("dvProductionToWaste"*_n)][t,"DomesticHotWater",ts], 0.0, force=true)
+                end
             end
-        end
+        elseif p.s.ashp_wh.force_dispatch
+            dv = "binASHPWHSizeExceedsThermalLoad"*_n
+            m[Symbol(dv)] = @variable(m, [p.time_steps], binary=true, base_name=dv)
+            dv = "dvASHPWHSizeTimesExcess"*_n
+            m[Symbol(dv)] = @variable(m, [p.time_steps], lower_bound=0, base_name=dv)
+            # binary variable enforcement for size >= load
+            max_wh_size_bigM = 2*max(p.max_sizes["ASHPWaterHeater"], maximum(p.heating_loads_kw["DomesticHotWater"] ./ p.heating_cf["ASHPWaterHeater"]))
+            @constraint(m, [ts in p.time_steps],
+                m[Symbol("binASHPWHSizeExceedsThermalLoad"*_n)][ts] >= (m[Symbol("dvSize"*_n)]["ASHPWaterHeater"] - p.heating_loads_kw["DomesticHotWater"][ts] / p.heating_cf["ASHPWaterHeater"][ts]) / max_wh_size_bigM
+            )
+            @constraint(m, [ts in p.time_steps],
+                m[Symbol("binASHPWHSizeExceedsThermalLoad"*_n)][ts] <= 1 - (p.heating_loads_kw["DomesticHotWater"][ts] / p.heating_cf["ASHPWaterHeater"][ts] - m[Symbol("dvSize"*_n)]["ASHPWaterHeater"]) / max_wh_size_bigM
+            )
+            # set dvASHPWHSizeTimesExcess = binASHPWHSizeExceedsThermalLoad * dvSize
+            # big-M is min CF times heat load
+
+            @constraint(m, [ts in p.time_steps],
+                m[Symbol("dvASHPWHSizeTimesExcess"*_n)][ts] >= p.heating_cf["ASHPWaterHeater"][ts]*m[Symbol("dvSize"*_n)]["ASHPWaterHeater"] - max_wh_size_bigM * (1-m[Symbol("binASHPWHSizeExceedsThermalLoad"*_n)][ts])  
+            )
+            @constraint(m, [ts in p.time_steps],
+                m[Symbol("dvASHPWHSizeTimesExcess"*_n)][ts] <= p.heating_cf["ASHPWaterHeater"][ts]*m[Symbol("dvSize"*_n)]["ASHPWaterHeater"]
+            )
+            @constraint(m, [ts in p.time_steps],
+                m[Symbol("dvASHPWHSizeTimesExcess"*_n)][ts] <= max_wh_size_bigM * m[Symbol("binASHPWHSizeExceedsThermalLoad"*_n)][ts]
+            )
+            #Enforce dispatch: output = system size - (overage)
+            @constraint(m, [ts in p.time_steps],
+                m[Symbol("dvHeatingProduction"*_n)]["ASHPWaterHeater","DomesticHotWater",ts] >= p.heating_cf["ASHPWaterHeater"][ts]*m[Symbol("dvSize"*_n)]["ASHPWaterHeater"] - m[Symbol("dvASHPWHSizeTimesExcess"*_n)][ts] + p.heating_loads_kw["DomesticHotWater"][ts] * m[Symbol("binASHPWHSizeExceedsThermalLoad"*_n)][ts]
+            )
+        end 
     end
 end
 

src/core/ashp.jl

@@ -23,6 +23,7 @@ ASHPSpaceHeater has the following attributes:
     cooling_cf::Array{<:Real,1} # ASHP's cooling capacity factor curves
     can_serve_cooling::Bool # If ASHP can supply heat to the cooling load
     force_into_system::Bool # force into system to serve all space heating loads if true
+    force_dispatch::Bool # force ASHP to meet load or maximize output if true
     back_up_temp_threshold_degF::Real # Degree in F that system switches from ASHP to resistive heater 
     avoided_capex_by_ashp_present_value::Real # avoided capital expenditure due to presence of ASHP system vs. defaults heating and cooling techs
     max_ton::Real # maximum allowable thermal power (tons) 
@@ -48,6 +49,7 @@ struct ASHP <: AbstractThermalTech
     can_serve_process_heat::Bool
     can_serve_cooling::Bool
     force_into_system::Bool
+    force_dispatch::Bool
     back_up_temp_threshold_degF::Real
     avoided_capex_by_ashp_present_value::Real
     max_ton::Real
@@ -81,6 +83,7 @@ function ASHPSpaceHeater(;
     macrs_bonus_fraction::Real = 0.0, # Fraction of upfront project costs to depreciate under MACRS
     can_serve_cooling::Union{Bool, Nothing} = nothing # If ASHP can supply heat to the cooling load
     force_into_system::Bool = false # force into system to serve all space heating loads if true
+    force_dispatch::Bool = true # force ASHP to meet load or maximize output if true
     avoided_capex_by_ashp_present_value::Real = 0.0 # avoided capital expenditure due to presence of ASHP system vs. defaults heating and cooling techs
 
     #The following inputs are used to create the attributes heating_cop and heating cf: 
@@ -103,7 +106,7 @@ function ASHPSpaceHeater(;
 """
 function ASHPSpaceHeater(;
         min_ton::Real = 0.0,
-        max_ton::Real = BIG_NUMBER,
+        max_ton::Union{Real, Nothing} = nothing,
         min_allowable_ton::Union{Real, Nothing} = nothing,
         min_allowable_peak_capacity_fraction::Union{Real, Nothing} = nothing, 
         sizing_factor::Union{Real, Nothing} = nothing, 
@@ -114,6 +117,7 @@ function ASHPSpaceHeater(;
         avoided_capex_by_ashp_present_value::Real = 0.0,
         can_serve_cooling::Union{Bool, Nothing} = nothing,
         force_into_system::Bool = false,
+        force_dispatch::Bool = true,
         heating_cop_reference::Array{<:Real,1} = Real[],
         heating_cf_reference::Array{<:Real,1} = Real[],
         heating_reference_temps_degF::Array{<:Real,1} = Real[],
@@ -144,9 +148,6 @@ function ASHPSpaceHeater(;
     if isnothing(back_up_temp_threshold_degF)
         back_up_temp_threshold_degF = defaults["back_up_temp_threshold_degF"]
     end
-    if isnothing(max_ton)
-        max_ton = defaults["max_ton"]
-    end
     if isnothing(sizing_factor)
         sizing_factor = defaults["sizing_factor"]
     end
@@ -177,11 +178,6 @@ function ASHPSpaceHeater(;
     can_serve_process_heat = defaults["can_serve_process_heat"]
 
 
-    # Convert max sizes, cost factors from mmbtu_per_hour to kw
-    min_kw = min_ton * KWH_THERMAL_PER_TONHOUR
-    max_kw = max_ton * KWH_THERMAL_PER_TONHOUR
-
-
     installed_cost_per_kw = installed_cost_per_ton / KWH_THERMAL_PER_TONHOUR
     om_cost_per_kw = om_cost_per_ton / KWH_THERMAL_PER_TONHOUR
 
@@ -206,14 +202,34 @@ function ASHPSpaceHeater(;
         cooling_cf = Float64[]
     end
 
+    if isnothing(max_ton)
+        if can_serve_cooling 
+            #these are in kW rathern than tons so will be larger as a measure of conservatism
+            max_ton = min(defaults["max_ton"], 10*(maximum(heating_load ./ heating_cf)+maximum(cooling_load ./ cooling_cf)))
+        else
+            max_ton = min(defaults["max_ton"], 10*maximum(heating_load ./ heating_cf))
+        end
+    else
+        if can_serve_cooling 
+            max_ton = min(max_ton, 10*(maximum(heating_load ./ heating_cf)+maximum(cooling_load ./ cooling_cf)))
+        else
+            max_ton = min(max_ton, 10*maximum(heating_load ./ heating_cf)) 
+        end
+        max_ton = max(max_ton, min_ton)
+    end
+
+    # Convert max sizes, cost factors from mmbtu_per_hour to kw
+    min_kw = min_ton * KWH_THERMAL_PER_TONHOUR
+    max_kw = max_ton * KWH_THERMAL_PER_TONHOUR    
+
     if !isnothing(min_allowable_ton) && !isnothing(min_allowable_peak_capacity_fraction)
         throw(@error("at most one of min_allowable_ton and min_allowable_peak_capacity_fraction may be input."))
     elseif !isnothing(min_allowable_ton)
         min_allowable_kw = min_allowable_ton * KWH_THERMAL_PER_TONHOUR
         @warn("user-provided minimum allowable ton is used in the place of the default; this may provided very small sizes if set to zero.")
     else
         if isnothing(min_allowable_peak_capacity_fraction)
-            min_allowable_peak_capacity_fraction = 0.5
+            min_allowable_peak_capacity_fraction = 0.25
         end
         min_allowable_kw = get_ashp_default_min_allowable_size(heating_load, heating_cf, cooling_load, cooling_cf, min_allowable_peak_capacity_fraction)
     end
@@ -244,6 +260,7 @@ function ASHPSpaceHeater(;
         can_serve_process_heat,
         can_serve_cooling,
         force_into_system,
+        force_dispatch,
         back_up_temp_threshold_degF,
         avoided_capex_by_ashp_present_value,
         max_ton,
@@ -278,6 +295,7 @@ function ASHPWaterHeater(;
     can_supply_steam_turbine::Union{Bool, nothing} = nothing # If the boiler can supply steam to the steam turbine for electric production
     avoided_capex_by_ashp_present_value::Real = 0.0 # avoided capital expenditure due to presence of ASHP system vs. defaults heating and cooling techs
     force_into_system::Bool = false # force into system to serve all hot water loads if true
+    force_dispatch::Bool = true # force ASHP to meet load or maximize output if true
     
     #The following inputs are used to create the attributes heating_cop and heating cf: 
     heating_cop_reference::Array{<:Real,1}, # COP of the heating (i.e., thermal produced / electricity consumed)
@@ -293,7 +311,7 @@ function ASHPWaterHeater(;
 """
 function ASHPWaterHeater(;
     min_ton::Real = 0.0,
-    max_ton::Real = BIG_NUMBER,
+    max_ton::Union{Real, Nothing} = nothing,
     min_allowable_ton::Union{Real, Nothing} = nothing,
     min_allowable_peak_capacity_fraction::Union{Real, Nothing} = nothing, 
     sizing_factor::Union{Real, Nothing} = nothing, 
@@ -303,6 +321,7 @@ function ASHPWaterHeater(;
     macrs_bonus_fraction::Real = 0.0,
     avoided_capex_by_ashp_present_value::Real = 0.0,
     force_into_system::Bool = false,
+    force_dispatch::Bool = true,
     heating_cop_reference::Array{<:Real,1} = Real[],
     heating_cf_reference::Array{<:Real,1} = Real[],
     heating_reference_temps_degF::Array{<:Real,1} = Real[],
@@ -323,9 +342,6 @@ function ASHPWaterHeater(;
     if isnothing(back_up_temp_threshold_degF)
         back_up_temp_threshold_degF = defaults["back_up_temp_threshold_degF"]
     end
-    if isnothing(max_ton)
-        max_ton = defaults["max_ton"]
-    end
     if isnothing(sizing_factor)
         sizing_factor = defaults["sizing_factor"]
     end
@@ -347,10 +363,6 @@ function ASHPWaterHeater(;
      can_serve_process_heat = defaults["can_serve_process_heat"]
      can_serve_cooling = defaults["can_serve_cooling"]
 
-    # Convert max sizes, cost factors from mmbtu_per_hour to kw
-    min_kw = min_ton * KWH_THERMAL_PER_TONHOUR
-    max_kw = max_ton * KWH_THERMAL_PER_TONHOUR
-
     installed_cost_per_kw = installed_cost_per_ton / KWH_THERMAL_PER_TONHOUR
     om_cost_per_kw = om_cost_per_ton / KWH_THERMAL_PER_TONHOUR
 
@@ -363,14 +375,26 @@ function ASHPWaterHeater(;
 
     heating_cf[heating_cop .== 1] .= 1
 
+    if isnothing(max_ton)
+        #these are in kW rathern than tons so will be larger as a measure of conservatism
+        max_ton = min(defaults["max_ton"], 10*maximum(heating_load ./ heating_cf))
+    else
+        max_ton = min(max_ton, 10*maximum(heating_load ./ heating_cf))
+    end
+    max_ton = max(max_ton, min_ton)
+
+    # Convert max sizes, cost factors from mmbtu_per_hour to kw
+    min_kw = min_ton * KWH_THERMAL_PER_TONHOUR
+    max_kw = max_ton * KWH_THERMAL_PER_TONHOUR  
+
     if !isnothing(min_allowable_ton) && !isnothing(min_allowable_peak_capacity_fraction)
         throw(@error("at most one of min_allowable_ton and min_allowable_peak_capacity_fraction may be input."))
     elseif !isnothing(min_allowable_ton)
         min_allowable_kw = min_allowable_ton * KWH_THERMAL_PER_TONHOUR
         @warn("user-provided minimum allowable ton is used in the place of the default; this may provided very small sizes if set to zero.")
     else
         if isnothing(min_allowable_peak_capacity_fraction)
-            min_allowable_peak_capacity_fraction = 0.5
+            min_allowable_peak_capacity_fraction = 0.25
         end
         min_allowable_kw = get_ashp_default_min_allowable_size(heating_load, heating_cf, Real[], Real[], min_allowable_peak_capacity_fraction)
     end
@@ -398,6 +422,7 @@ function ASHPWaterHeater(;
         can_serve_process_heat,
         can_serve_cooling,
         force_into_system,
+        force_dispatch,
         back_up_temp_threshold_degF,
         avoided_capex_by_ashp_present_value,
         max_ton,