Maximize ASHP output in dispatch

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 = false # 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: 
@@ -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 = false,
         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,6 +202,20 @@ 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"], maximum(heating_load ./ heating_cf)+maximum(cooling_load ./ cooling_cf))
+        else
+            max_ton = min(defaults["max_ton"], 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)
@@ -244,6 +254,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 +289,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 = false # 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)
@@ -303,6 +315,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 = false,
     heating_cop_reference::Array{<:Real,1} = Real[],
     heating_cf_reference::Array{<:Real,1} = Real[],
     heating_reference_temps_degF::Array{<:Real,1} = Real[],
@@ -323,9 +336,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 +357,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,6 +369,16 @@ 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"], maximum(heating_load)/minimum(heating_cf))
+        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)
@@ -398,6 +414,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,

src/core/reopt_inputs.jl

@@ -953,14 +953,7 @@ function setup_ASHPSpaceHeater_inputs(s, max_sizes, min_sizes, cap_cost_slope, o
     if s.ashp.min_allowable_kw > 0.0
         cap_cost_slope["ASHPSpaceHeater"] = s.ashp.installed_cost_per_kw
         push!(segmented_techs, "ASHPSpaceHeater")
-        #as a reasonable big-M, assume 10 times the size to serve 100% of peak load.
-        ashp_sh_max_size = get_ashp_default_min_allowable_size(s.space_heating_load.loads_kw, 
-            s.ashp.heating_cf, 
-            s.cooling_load.loads_kw_thermal, 
-            s.ashp.cooling_cf,
-            10.0
-        )
-        seg_max_size["ASHPSpaceHeater"] = Dict{Int,Float64}(1 => min(s.ashp.max_kw, ashp_sh_max_size))
+        seg_max_size["ASHPSpaceHeater"] = Dict{Int,Float64}(1 => min(s.ashp.max_kw))
         seg_min_size["ASHPSpaceHeater"] = Dict{Int,Float64}(1 => s.ashp.min_allowable_kw)
         n_segs_by_tech["ASHPSpaceHeater"] = 1
         seg_yint["ASHPSpaceHeater"] = Dict{Int,Float64}(1 => 0.0)
@@ -997,14 +990,7 @@ function setup_ASHPWaterHeater_inputs(s, max_sizes, min_sizes, cap_cost_slope, o
     if s.ashp_wh.min_allowable_kw > 0.0
         cap_cost_slope["ASHPWaterHeater"] = s.ashp_wh.installed_cost_per_kw
         push!(segmented_techs, "ASHPWaterHeater")
-        #as a reasonable big-M, assume 10 times the size to serve 100% of peak load.
-        ashp_wh_max_size = get_ashp_default_min_allowable_size(s.dhw_load.loads_kw, 
-            s.ashp_wh.heating_cf, 
-            Real[], 
-            Real[],
-            10.0
-        )
-        seg_max_size["ASHPWaterHeater"] = Dict{Int,Float64}(1 => min(s.ashp_wh.max_kw, ashp_wh_max_size))
+        seg_max_size["ASHPWaterHeater"] = Dict{Int,Float64}(1 => s.ashp_wh.max_kw)
         seg_min_size["ASHPWaterHeater"] = Dict{Int,Float64}(1 => s.ashp_wh.min_allowable_kw)
         n_segs_by_tech["ASHPWaterHeater"] = 1
         seg_yint["ASHPWaterHeater"] = Dict{Int,Float64}(1 => 0.0)