Allow presized ghp/ghx

.github/workflows/CI.yml

@@ -22,5 +22,5 @@ jobs:
           version: ${{ matrix.julia-version }}
       - uses: julia-actions/julia-buildpkg@latest
       # - uses: mxschmitt/action-tmate@v3  # for interactive debugging
-      - run: julia --project=. -e 'using Pkg; Pkg.activate("test"); Pkg.rm("Xpress"); Pkg.activate("."); using TestEnv; TestEnv.activate(); cd("test"); include("runtests.jl")'
+      - run: julia --project=. -e 'using Pkg; Pkg.activate("test"); Pkg.rm("Xpress"); Pkg.activate("."); using TestEnv; TestEnv.activate(); Pkg.add(PackageSpec(name="GhpGhx", url="https://github.com/NREL/GhpGhx.jl.git")); using GhpGhx; cd("test"); include("runtests.jl")'
         shell: bash

CHANGELOG.md

@@ -100,6 +100,7 @@ Classify the change according to the following categories:
 
 ## v0.48.2
 ### Added
+- Add new optional parameter **max_ton** to GHP module to allow user to size GHP smaller than peak load
 - Battery residual value if choosing replacement strategy for degradation
 - Add new **ElectricStorage** parameters **max_duration_hours** and **min_duration_hours** to bound the energy duration of battery storage
 ### Changed

src/REopt.jl

@@ -58,11 +58,12 @@ function __init__()
 end
 
 const EXISTING_BOILER_EFFICIENCY = 0.8
-const GAL_PER_M3 = 264.172  # [gal/m^3]
-const KWH_PER_GAL_DIESEL = 40.7  # [kWh/gal_diesel] higher heating value of diesel
-const KWH_PER_MMBTU = 293.07107  # [kWh/mmbtu]
+const GAL_PER_M3 = 264.172              # [gal/m^3]
+const KWH_PER_GAL_DIESEL = 40.7         # [kWh/gal_diesel] higher heating value of diesel
+const KWH_PER_MMBTU = 293.07107         # [kWh/mmbtu]
 const KWH_THERMAL_PER_TONHOUR = 3.51685
-const TONNE_PER_LB = 1/2204.62  # [tonne/lb]
+const TONNE_PER_LB = 1/2204.62          # [tonne/lb]
+const TONNE_PER_MMBTU_HOUR = 0.012      # [tonne/mmbtu-hour]
 const FUEL_TYPES = ["natural_gas", "landfill_bio_gas", "propane", "diesel_oil"]
 const BIG_NUMBER = 1.0e10  #used for max size.  TODO use this number elsewhere.
 const PRIME_MOVERS = ["recip_engine", "micro_turbine", "combustion_turbine", "fuel_cell"]  #TODO replace `prime_movers` references in CHP code

src/core/ghp.jl

@@ -19,11 +19,14 @@ struct with outer constructor:
     installed_cost_ghx_per_ft::Float64 = 14.0
     installed_cost_building_hydronic_loop_per_sqft = 1.70
     om_cost_per_sqft_year::Float64 = -0.51
-    building_sqft::Float64 # Required input
+    building_sqft::Float64                                  # Required input
     space_heating_efficiency_thermal_factor::Float64 = NaN  # Default depends on building and location
-    cooling_efficiency_thermal_factor::Float64 = NaN # Default depends on building and location
+    cooling_efficiency_thermal_factor::Float64 = NaN        # Default depends on building and location
     ghpghx_response::Dict = Dict()
     can_serve_dhw::Bool = false
+    max_ton::Real                                           # Maximum heat pump capacity size. Default at a big number
+    max_number_of_boreholes::Real                           # Maximum GHX size
+    load_served_by_ghp::String                              # "scaled" or "nonpeak"
 
     macrs_option_years::Int = 5
     macrs_bonus_fraction::Float64 = 0.6
@@ -77,6 +80,9 @@ Base.@kwdef mutable struct GHP <: AbstractGHP
     can_serve_space_heating::Bool = true
     can_serve_process_heat::Bool = false
     can_supply_steam_turbine::Bool = false
+    max_ton::Real = BIG_NUMBER
+    max_number_of_boreholes::Real = BIG_NUMBER
+    load_served_by_ghp::String = "nonpeak"
 
     aux_heater_type::String = "electric"
     is_ghx_hybrid::Bool = false
@@ -154,7 +160,7 @@ function GHP(response::Dict, d::Dict)
     end
     # incentives = IncentivesNoProdBased(**d_mod)
 
-    setup_installed_cost_curve!(ghp, response)
+    setup_installed_cost_curve!(d, ghp, response)
 
     setup_om_cost!(ghp)
 
@@ -168,10 +174,10 @@ function GHP(response::Dict, d::Dict)
 end
 
 """
-    setup_installed_cost_curve!(response::Dict, ghp::GHP)
+    setup_installed_cost_curve!(d::Dict, response::Dict, ghp::GHP)
 
 """
-function setup_installed_cost_curve!(ghp::GHP, response::Dict)
+function setup_installed_cost_curve!(d::Dict, 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"]

src/core/scenario.jl

@@ -480,7 +480,7 @@ function Scenario(d::Dict; flex_hvac_from_json=false)
     space_heating_thermal_load_reduction_with_ghp_kw = zeros(8760 * settings.time_steps_per_hour)
     cooling_thermal_load_reduction_with_ghp_kw = zeros(8760 * settings.time_steps_per_hour)
     eval_ghp = false
-    get_ghpghx_from_input = false    
+    get_ghpghx_from_input = false
     if haskey(d, "GHP") && haskey(d["GHP"],"building_sqft")
         eval_ghp = true
         if haskey(d["GHP"], "ghpghx_responses") && !isempty(d["GHP"]["ghpghx_responses"])
@@ -620,13 +620,103 @@ function Scenario(d::Dict; flex_hvac_from_json=false)
             end
 
             ghpghx_results = Dict()
+
             try
                 # Call GhpGhx.jl to size GHP and GHX
                 @info "Starting GhpGhx.jl"
                 # Call GhpGhx.jl to size GHP and GHX
+                # If user provides undersized GHP, calculate load to send to GhpGhx.jl, and load to send to REopt for backup
+                thermal_load_ton = ghpghx_inputs["heating_thermal_load_mmbtu_per_hr"] * 1/TONNE_PER_MMBTU_HOUR
+                if haskey(ghpghx_inputs,"cooling_thermal_load_ton")
+                    cooling_load_ton = ghpghx_inputs["cooling_thermal_load_ton"]
+                    thermal_load_ton .+= cooling_load_ton
+                end
+                peak_thermal_load_ton = maximum(thermal_load_ton)
+                if haskey(d["GHP"],"max_ton") && peak_thermal_load_ton > d["GHP"]["max_ton"]
+                    @info "User entered undersized GHP. Calculating load that can be served by user specified undersized GHP"
+                    # When user specifies undersized GHP, calculate the load to be served by GHP and send the rest to REopt
+                    if !haskey(d["GHP"], "load_served_by_ghp")
+                        d["GHP"]["load_served_by_ghp"] = "nonpeak"
+                    end
+                    # If user choose to scale down total load (load_served_by_ghp="scaled"), calculate the ratio of the udersized GHP size and peak load
+                    if d["GHP"]["load_served_by_ghp"] == "scaled"
+                        @info "GHP served scaled down of total thermal load"
+                        peak_ratio = d["GHP"]["max_ton"] / peak_thermal_load_ton
+                        # Scale the total load profile down by the peak_ratio and use this scaled down load to rerun GhpGhx.jl
+                        heating_load_mmbtu = ghpghx_inputs["heating_thermal_load_mmbtu_per_hr"]
+                        heating_load_mmbtu = heating_load_mmbtu .* peak_ratio
+                        ghpghx_inputs["heating_thermal_load_mmbtu_per_hr"] = heating_load_mmbtu
+                        if haskey(ghpghx_inputs,"cooling_thermal_load_ton")
+                            ghpghx_inputs["cooling_thermal_load_ton"] = cooling_load_ton .* peak_ratio
+                        end
+                    elseif d["GHP"]["load_served_by_ghp"] == "nonpeak"
+                        @info "GHP serves all thermal load below peak"
+                        heating_load_mmbtu = ghpghx_inputs["heating_thermal_load_mmbtu_per_hr"]
+                        # if cooling load is included, cut down total thermal load and send as much heating load to GhpGhx.jl as possible
+                        if haskey(ghpghx_inputs,"cooling_thermal_load_ton")
+                            thermal_load_ton = heating_load_mmbtu .* 1/TONNE_PER_MMBTU_HOUR .+ cooling_load_ton
+                            peak_thermal_load_ton = maximum(thermal_load_ton)
+                            # If total thermal load (heating + cooling) is more than user-defined GHP size, 
+                            # first reduce heating load as much as possible while keeping cooling load the same
+                            if peak_thermal_load_ton > d["GHP"]["max_ton"]
+                                thermal_load_ton[thermal_load_ton .>= d["GHP"]["max_ton"]] .= d["GHP"]["max_ton"]
+                                heating_load_ton = thermal_load_ton .- cooling_load_ton
+                                # Make sure that the reduced heating load is not negative
+                                heating_load_ton[heating_load_ton .< 0] .= 0
+                                # If the updated peak thermal load is still more than user-defined GHP size, 
+                                # reduce cooling load as well
+                                updated_thermal_load_ton = heating_load_ton .+ cooling_load_ton
+                                updated_peak_thermal_load_ton = maximum(updated_thermal_load_ton)
+                                if updated_peak_thermal_load_ton > d["GHP"]["max_ton"]
+                                    updated_thermal_load_ton[updated_thermal_load_ton .>= d["GHP"]["max_ton"]] .= d["GHP"]["max_ton"]
+                                    cooling_load_ton = updated_thermal_load_ton .- heating_load_ton
+                                    ghpghx_inputs["cooling_thermal_load_ton"] = cooling_load_ton
+                                end
+                                heating_load_mmbtu = heating_load_ton .* TONNE_PER_MMBTU_HOUR
+                                ghpghx_inputs["heating_thermal_load_mmbtu_per_hr"] = heating_load_mmbtu
+                            end
+                        # if cooling load is not included, cut down heating load only and send to GhpGhx.jl
+                        else
+                            heating_load_mmbtu[heating_load_mmbtu .>= d["GHP"]["max_ton"] * TONNE_PER_MMBTU_HOUR] .= d["GHP"]["max_ton"] * TONNE_PER_MMBTU_HOUR
+                            ghpghx_inputs["heating_thermal_load_mmbtu_per_hr"] = heating_load_mmbtu
+                        end                         
+                    end
+                end
                 results, inputs_params = GhpGhx.ghp_model(ghpghx_inputs)
+                # If max_number_of_boreholes is specified, check if number of boreholes sized by GhpGhx.jl greater than user-specified max_number_of_boreholes,
+                # and if max_number_of_boreholes is less, reduce thermal load served by GHP until max_number_of_boreholes = number of boreholses sized by GhpGhx.jl                
+                if haskey(d["GHP"],"max_number_of_boreholes")
+                    optimal_number_of_boreholes = GhpGhx.get_results_for_reopt(results, inputs_params)["number_of_boreholes"]
+                    if optimal_number_of_boreholes > d["GHP"]["max_number_of_boreholes"]
+                        @info "Max number of boreholes specified is less than number of boreholes sized in GhpGhx.jl, reducing thermal load served by GHP further"
+                        max_iter = 10
+                        for iter = 1:max_iter
+                            borehole_ratio = d["GHP"]["max_number_of_boreholes"] / optimal_number_of_boreholes
+                            heating_load_mmbtu .*= borehole_ratio
+                            if haskey(ghpghx_inputs,"cooling_thermal_load_ton")
+                                cooling_load_ton .*= borehole_ratio
+                            # if cooling load is not included, cut down heating load only and send to GhpGhx.jl
+                            end
+                            ghpghx_inputs["heating_thermal_load_mmbtu_per_hr"] = heating_load_mmbtu                              
+                            ghpghx_inputs["cooling_thermal_load_ton"] = cooling_load_ton   
+
+                            # Rerun GhpGhx.jl
+                            results, inputs_params = GhpGhx.ghp_model(ghpghx_inputs)
+                            optimal_number_of_boreholes = GhpGhx.get_results_for_reopt(results, inputs_params)["number_of_boreholes"]
+                            # Solution is found if the new optimal number of boreholes sized by GhpGhx.jl = user-specified max number of boreholes,
+                            # Otherwise, continue solving until reaching max iteration
+                            if -0.5 < optimal_number_of_boreholes - d["GHP"]["max_number_of_boreholes"] < 0.5
+                                break
+                            else
+                                iter += 1
+                            end
+                        end
+                    end
+                end
+
                 # Create a dictionary of the results data needed for REopt
                 ghpghx_results = GhpGhx.get_results_for_reopt(results, inputs_params)
+                # Return results from GhpGhx.jl without load scaling if user does not provide GHP size or if user entered GHP size is greater than GHP size output
                 @info "GhpGhx.jl model solved" #with status $(results["status"])."
             catch e
                 @info e
@@ -644,9 +734,7 @@ function Scenario(d::Dict; flex_hvac_from_json=false)
             end                    
             append!(ghp_option_list, [GHP(ghpghx_response, ghp_inputs_removed_ghpghx_params)])
             # Print out ghpghx_response for loading into a future run without running GhpGhx.jl again
-            #open("scenarios/ghpghx_response.json","w") do f
-            #    JSON.print(f, ghpghx_response)
-            #end                
+            # open("scenarios/ghpghx_response.json","w") do f             
         end
     # If ghpghx_responses is included in inputs, do NOT run GhpGhx.jl model and use already-run ghpghx result as input to REopt
     elseif eval_ghp && get_ghpghx_from_input

src/results/existing_boiler.jl

@@ -5,7 +5,7 @@
 - `fuel_consumption_series_mmbtu_per_hour`  # Fuel consumption series [MMBtu/hr] 
 - `annual_fuel_consumption_mmbtu`  # Fuel consumed in a year [MMBtu]
 - `thermal_production_series_mmbtu_per_hour`  # Thermal energy production series [MMBtu/hr]
-- `annual_thermal_production_mmbtu`  # Thermal power production to TES (HotThermalStorage) series [MMBtu/hr]
+- `annual_thermal_production_mmbtu`  # Thermal power production in a year [MMBtu]
 - `thermal_to_storage_series_mmbtu_per_hour` # Thermal power production to TES (HotThermalStorage) series [MMBtu/hr]
 - `thermal_to_steamturbine_series_mmbtu_per_hour`  # Thermal power production to SteamTurbine series [MMBtu/hr]
 - `thermal_to_load_series_mmbtu_per_hour`  # Thermal power production to serve the heating load series [MMBtu/hr]

src/results/ghp.jl

@@ -15,6 +15,9 @@ GHP results:
 - `thermal_to_space_heating_load_series_mmbtu_per_hour`
 - `thermal_to_dhw_load_series_mmbtu_per_hour`
 - `thermal_to_load_series_ton`
+- `annual_thermal_production_mmbtu`  # GHP's heating thermal power production in a year [MMBtu]
+- `annual_thermal_production_tonhour`  # GHP's cooling thermal power production in a year [ton]
+
 """
 
 function add_ghp_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
@@ -25,6 +28,11 @@ function add_ghp_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
     # r["size_heat_pump_ton"] = 0.0
     # r["size_wwhp_heating_pump_ton"] = 0.0
     # r["size_wwhp_cooling_pump_ton"] = 0.0
+
+    # Set sizing factor = 1 if user inputs their own GHP size
+    if haskey(d, "GHP") && haskey(d["GHP"],"max_ton")
+        p.s.ghp_option_list[ghp_option_chosen].heatpump_capacity_sizing_factor_on_peak_load = 1.0
+    end
     if ghp_option_chosen >= 1
         r["ghpghx_chosen_outputs"] = p.s.ghp_option_list[ghp_option_chosen].ghpghx_response["outputs"]
 
@@ -42,13 +50,22 @@ function add_ghp_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
         r["space_heating_thermal_load_reduction_with_ghp_mmbtu_per_hour"] = round.(value.(HeatingThermalReductionWithGHP) ./ KWH_PER_MMBTU, digits=3)
         @expression(m, CoolingThermalReductionWithGHP[ts in p.time_steps],
 		    sum(p.cooling_thermal_load_reduction_with_ghp_kw[g,ts] * m[Symbol("binGHP"*_n)][g] for g in p.ghp_options))
+
+        @expression(m, HeatingThermalLoadServedWithGHP[ts in p.time_steps],
+		    sum(p.ghp_heating_thermal_load_served_kw[g,ts] * m[Symbol("binGHP"*_n)][g] for g in p.ghp_options))    
+        @expression(m, CoolingThermalLoadServedWithGHP[ts in p.time_steps],
+		    sum(p.ghp_cooling_thermal_load_served_kw[g,ts] * m[Symbol("binGHP"*_n)][g] for g in p.ghp_options)) 
+
         r["cooling_thermal_load_reduction_with_ghp_ton"] = round.(value.(CoolingThermalReductionWithGHP) ./ KWH_THERMAL_PER_TONHOUR, digits=3)
         r["ghx_residual_value_present_value"] = value(m[:ResidualGHXCapCost])
         r["avoided_capex_by_ghp_present_value"] = value(m[:AvoidedCapexByGHP])
-        r["thermal_to_space_heating_load_series_mmbtu_per_hour"] = d["HeatingLoad"]["space_heating_thermal_load_series_mmbtu_per_hour"] .- r["space_heating_thermal_load_reduction_with_ghp_mmbtu_per_hour"]
-        r["thermal_to_load_series_ton"] = d["CoolingLoad"]["load_series_ton"] .- r["cooling_thermal_load_reduction_with_ghp_ton"]
+        r["thermal_to_space_heating_load_series_mmbtu_per_hour"] = round.(value.(HeatingThermalLoadServedWithGHP) ./ KWH_PER_MMBTU, digits=3)
+        r["thermal_to_load_series_ton"] = round.(value.(CoolingThermalLoadServedWithGHP) ./ KWH_THERMAL_PER_TONHOUR, digits=3)
+        r["annual_thermal_production_mmbtu"] = sum(r["thermal_to_space_heating_load_series_mmbtu_per_hour"])
+        r["annual_thermal_production_tonhour"] = sum(r["thermal_to_load_series_ton"])
         if p.s.ghp_option_list[ghp_option_chosen].can_serve_dhw
             r["thermal_to_dhw_load_series_mmbtu_per_hour"] = d["HeatingLoad"]["dhw_thermal_load_series_mmbtu_per_hour"]
+            r["annual_thermal_production_mmbtu"] = r["annual_thermal_production_mmbtu"] + sum(r["thermal_to_dhw_load_series_mmbtu_per_hour"])
         else
             r["thermal_to_dhw_load_series_mmbtu_per_hour"] = zeros(length(p.time_steps))
         end
@@ -60,6 +77,8 @@ function add_ghp_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
         r["thermal_to_space_heating_load_series_mmbtu_per_hour"] = zeros(length(p.time_steps))
         r["thermal_to_load_series_ton"] = zeros(length(p.time_steps))
         r["thermal_to_dhw_load_series_mmbtu_per_hour"] = zeros(length(p.time_steps))
+        r["annual_thermal_production_mmbtu"] = 0.0
+        r["annual_thermal_production_tonhour"] = 0.0
     end
     d["GHP"] = r
     nothing