Off-grid CHP

src/constraints/operating_reserve_constraints.jl

@@ -8,9 +8,15 @@ function add_operating_reserve_constraints(m, p; _n="")
         sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for b in p.s.storage.types.elec) -
         m[Symbol("dvCurtail"*_n)][t, ts]
     )
+    # 1b. Production going to load from requiring_oper_res (separate calculation for techs that require reserves)
+    m[:ProductionToLoadRequiringOR] = @expression(m, [t in p.techs.requiring_oper_res, ts in p.time_steps_without_grid],
+        p.production_factor[t, ts] * p.levelization_factor[t] * m[Symbol("dvRatedProduction"*_n)][t,ts] -
+        sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for b in p.s.storage.types.elec) -
+        m[Symbol("dvCurtail"*_n)][t, ts]
+    )
     # 2. Total OR required by requiring_oper_res & Load 
     m[:OpResRequired] = @expression(m, [ts in p.time_steps_without_grid],
-        sum(m[:ProductionToLoadOR][t,ts] * p.techs_operating_reserve_req_fraction[t] for t in p.techs.requiring_oper_res)
+        sum(m[:ProductionToLoadRequiringOR][t,ts] * p.techs_operating_reserve_req_fraction[t] for t in p.techs.requiring_oper_res)
         + p.s.electric_load.critical_loads_kw[ts] * m[Symbol("dvOffgridLoadServedFraction"*_n)][ts] * p.s.electric_load.operating_reserve_required_fraction
     )
     # 3. Operating reserve provided - battery  
@@ -28,13 +34,19 @@ function add_operating_reserve_constraints(m, p; _n="")
     )
 
     # 5a. Upper bound on dvOpResFromTechs (for generator techs).  Note: will need to add new constraints for each new tech that can provide operating reserves
-    @constraint(m, [t in p.techs.gen, ts in p.time_steps_without_grid],
+    @constraint(m, [t in intersect(p.techs.gen, p.techs.providing_oper_res), ts in p.time_steps_without_grid],
         m[Symbol("dvOpResFromTechs"*_n)][t,ts] <= m[:binGenIsOnInTS][t, ts] * p.max_sizes[t] 
     )
     # 5b. Upper bound on dvOpResFromTechs (for pv techs)
-    @constraint(m, [t in p.techs.pv, ts in p.time_steps_without_grid],
+    @constraint(m, [t in intersect(p.techs.pv, p.techs.providing_oper_res), ts in p.time_steps_without_grid],
         m[Symbol("dvOpResFromTechs"*_n)][t,ts] <= p.max_sizes[t] 
     )
+    # 5c. Upper bound on dvOpResFromTechs (for CHP techs)
+    if "CHP" in p.techs.providing_oper_res
+        @constraint(m, [ts in p.time_steps_without_grid],
+            m[Symbol("dvOpResFromTechs"*_n)]["CHP",ts] <= m[:binCHPIsOnInTS]["CHP", ts] * p.max_sizes["CHP"] 
+        )
+    end
 
     m[:OpResProvided] = @expression(m, [ts in p.time_steps_without_grid],
         sum(m[Symbol("dvOpResFromTechs"*_n)][t,ts] for t in p.techs.providing_oper_res)

src/core/chp.jl

@@ -40,6 +40,7 @@ conflict_res_min_allowable_fraction_of_max = 0.25
     can_serve_space_heating::Bool = true # If CHP can supply heat to the space heating load
     can_serve_process_heat::Bool = true # If CHP can supply heat to the process heating load
     is_electric_only::Bool = false # If CHP is a prime generator that does not supply heat
+    operating_reserve_required_fraction::Real = off_grid_flag ? 0.0 : 0.0 # If off grid, 10%, else 0%. Applied to each time_step as a % of CHP generation. When positive, CHP requires reserves from Generator or Battery; when 0, CHP can provide reserves.
 
     macrs_option_years::Int = 5 # Notes: this value cannot be 0 if aiming to apply 100% bonus depreciation; default may change if Site.sector is not "commercial/industrial"
     macrs_bonus_fraction::Float64 = 1.0 #Note: default may change if Site.sector is not "commercial/industrial"
@@ -113,6 +114,7 @@ Base.@kwdef mutable struct CHP <: AbstractCHP
     can_serve_space_heating::Bool = true
     can_serve_process_heat::Bool = true
     is_electric_only::Bool = false
+    operating_reserve_required_fraction::Real = 0.0
 
     macrs_option_years::Int = 5
     macrs_bonus_fraction::Float64 = 1.0
@@ -149,7 +151,8 @@ function CHP(d::Dict;
             electric_load_series_kw::Array{<:Real,1}=Real[],
             year::Int64=2017,
             sector::String,
-            federal_procurement_type::String)
+            federal_procurement_type::String,
+            off_grid_flag::Bool=false)
     # If array inputs are coming from Julia JSON.parsefile (reader), they have type Vector{Any}; convert to expected type here
     for (k,v) in d
         if typeof(v) <: AbstractVector{Any} && k != "unavailability_periods"
@@ -274,6 +277,17 @@ function CHP(d::Dict;
         setproperty!(chp, :thermal_efficiency_half_load, 0.0)
     end
 
+    # Set operating_reserve_required_fraction based on off_grid_flag
+    if haskey(d, "operating_reserve_required_fraction")
+        chp.operating_reserve_required_fraction = d["operating_reserve_required_fraction"]
+    end
+
+    # Validate operating_reserve_required_fraction for on-grid scenarios
+    if !off_grid_flag && !(chp.operating_reserve_required_fraction == 0.0)
+        @warn "CHP operating_reserve_required_fraction applies only when off_grid_flag is true. Setting operating_reserve_required_fraction to 0.0 for this on-grid analysis."
+        chp.operating_reserve_required_fraction = 0.0
+    end
+
     return chp
 end
 

src/core/reopt_inputs.jl

@@ -1354,13 +1354,17 @@ function setup_ghp_inputs(s::AbstractScenario, time_steps, time_steps_without_gr
 end
 
 function setup_operating_reserve_fraction(s::AbstractScenario, techs_operating_reserve_req_fraction)
-    # currently only PV and Wind require operating reserves
+    # Currently PV, Wind, and CHP can require or provide operating reserves in off-grid scenarios
     for pv in s.pvs 
         techs_operating_reserve_req_fraction[pv.name] = pv.operating_reserve_required_fraction
     end
 
     techs_operating_reserve_req_fraction["Wind"] = s.wind.operating_reserve_required_fraction
 
+    if !isnothing(s.chp)
+        techs_operating_reserve_req_fraction["CHP"] = s.chp.operating_reserve_required_fraction
+    end
+
     return nothing
 end
 

src/core/scenario.jl

@@ -83,7 +83,7 @@ function Scenario(d::Dict; flex_hvac_from_json=false)
 
     # Check that only PV, electric storage, and generator are modeled for off-grid
     if settings.off_grid_flag
-        offgrid_allowed_keys = ["PV", "Wind", "ElectricStorage", "Generator", "Settings", "Site", "Financial", "ElectricLoad", "ElectricTariff", "ElectricUtility"]
+        offgrid_allowed_keys = ["PV", "Wind", "ElectricStorage", "Generator", "CHP", "Settings", "Site", "Financial", "ElectricLoad", "ElectricTariff", "ElectricUtility"]
         unallowed_keys = setdiff(keys(d), offgrid_allowed_keys) 
         if !isempty(unallowed_keys)
             throw(@error("The following key(s) are not permitted when `off_grid_flag` is true: $unallowed_keys."))
@@ -412,13 +412,15 @@ function Scenario(d::Dict; flex_hvac_from_json=false)
                     electric_load_series_kw = electric_load.loads_kw,
                     year = electric_load.year,
                     sector = site.sector,
-                    federal_procurement_type = site.federal_procurement_type)
+                    federal_procurement_type = site.federal_procurement_type,
+                    off_grid_flag = settings.off_grid_flag)
         else # Only if modeling CHP without heating_load and existing_boiler (for prime generator, electric-only)
-            chp = CHP(d["CHP"],
+            chp = CHP(d["CHP"];
                     electric_load_series_kw = electric_load.loads_kw,
                     year = electric_load.year,
                     sector = site.sector,
-                    federal_procurement_type = site.federal_procurement_type)
+                    federal_procurement_type = site.federal_procurement_type,
+                    off_grid_flag = settings.off_grid_flag)
         end
         chp_prime_mover = chp.prime_mover
     end

src/core/techs.jl

@@ -202,6 +202,13 @@ function Techs(s::Scenario)
         if s.chp.can_serve_process_heat
             push!(techs_can_serve_process_heat, "CHP")
         end
+        if s.settings.off_grid_flag
+            if s.chp.operating_reserve_required_fraction > 0.0
+                push!(requiring_oper_res, "CHP")
+            else
+                push!(providing_oper_res, "CHP")
+            end
+        end
     end
 
     if !isempty(s.ghp_option_list) && !isnothing(s.ghp_option_list[1])

src/results/chp.jl

@@ -99,11 +99,6 @@ function add_chp_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
         sum(sum(m[Symbol("dvHeatingProduction"*_n)][t,q,ts] for q in p.heating_loads) + m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts]
             for t in p.techs.chp) - CHPToHotTES[ts] - CHPToSteamTurbineKW[ts] - CHPThermalToWasteKW[ts])
     r["thermal_to_load_series_mmbtu_per_hour"] = round.(value.(CHPThermalToLoadKW ./ KWH_PER_MMBTU), digits=5)
-
-	CHPToLoadKW = @expression(m, [ts in p.time_steps],
-		sum(value.(m[:dvHeatingProduction]["CHP",q,ts] for q in p.heating_loads)) - CHPToHotTES[ts] - CHPToSteamTurbineKW[ts]
-    )
-	r["thermal_to_load_series_mmbtu_per_hour"] = round.(value.(CHPThermalToLoadKW ./ KWH_PER_MMBTU), digits=5)
 
     if "DomesticHotWater" in p.heating_loads && p.s.chp.can_serve_dhw
         @expression(m, CHPToDHWKW[ts in p.time_steps], 

test/scenarios/chp_offgrid.json

@@ -0,0 +1,34 @@
+{
+    "Settings": {
+        "time_steps_per_hour": 1,
+        "off_grid_flag": true
+    },
+    "Site": {
+        "latitude": 37.78,
+        "longitude": -122.45
+    },
+    "ElectricLoad": {
+        "doe_reference_name": "Hospital",
+        "annual_kwh": 1000000.0
+    },
+    "CHP": {
+        "prime_mover": "recip_engine",
+        "size_class": 2,
+        "min_kw": 50.0,
+        "max_kw": 500.0,
+        "installed_cost_per_kw": 2300.0,
+        "om_cost_per_kwh": 0.0,
+        "om_cost_per_kw": 50.0,
+        "fuel_cost_per_mmbtu": 8.0,
+        "is_electric_only": true
+    },
+    "ElectricStorage": {
+        "min_kw": 0.0,
+        "max_kw": 1000.0,
+        "min_kwh": 0.0,
+        "max_kwh": 5000.0,
+        "soc_min_fraction": 0.2,
+        "soc_init_fraction": 1.0,
+        "can_grid_charge": false
+    }
+}

test/test_chp_offgrid.jl

@@ -0,0 +1,105 @@
+using Revise
+using REopt
+using JSON
+using DelimitedFiles
+using PlotlyJS
+using Dates
+using Test
+using JuMP
+using HiGHS
+using DotEnv
+DotEnv.load!()
+
+###############   Off-Grid CHP with Operating Reserves Test    ###################
+
+# Test off-grid scenario with CHP and ElectricStorage where CHP requires 10% operating reserves
+# This tests that ElectricStorage provides the required reserves for CHP generation
+
+input_data = JSON.parsefile("./scenarios/chp_offgrid.json")
+input_data["ElectricLoad"]["min_load_met_annual_fraction"] = 0.999
+input_data["ElectricLoad"]["operating_reserve_required_fraction"] = 0.0
+# TODO see if no battery shows up when below is set to 0 so CHP can provide SR
+input_data["CHP"]["operating_reserve_required_fraction"] = 0.2
+
+println("\n" * "="^80)
+println("Testing Off-Grid CHP with Operating Reserves")
+println("="^80)
+
+# Create scenario and run optimization
+s = Scenario(input_data)
+inputs = REoptInputs(s)
+
+println("\nScenario Setup:")
+println("  - CHP operating_reserve_required_fraction: ", s.chp.operating_reserve_required_fraction)
+println("  - ElectricLoad operating_reserve_required_fraction: ", s.electric_load.operating_reserve_required_fraction)
+println("  - Off-grid flag: ", s.settings.off_grid_flag)
+
+m = Model(optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false, "log_to_console" => false, "mip_rel_gap" => 0.01))
+results = run_reopt(m, inputs)
+
+# Check that the model solved successfully
+@test termination_status(m) == MOI.OPTIMAL || termination_status(m) == MOI.LOCALLY_SOLVED
+
+println("\nOptimization Results:")
+println("  - Termination status: ", termination_status(m))
+println("  - CHP size: ", round(results["CHP"]["size_kw"], digits=1), " kW")
+println("  - Battery power: ", round(results["ElectricStorage"]["size_kw"], digits=1), " kW")
+println("  - Battery energy: ", round(results["ElectricStorage"]["size_kwh"], digits=1), " kWh")
+println("  - Load met fraction: ", round(results["ElectricLoad"]["offgrid_load_met_fraction"], digits=4))
+
+# Calculate operating reserve requirements
+chp_production_to_load = results["CHP"]["electric_to_load_series_kw"]
+chp_or_required = sum(chp_production_to_load .* s.chp.operating_reserve_required_fraction)
+load_or_required = sum(s.electric_load.critical_loads_kw .* s.electric_load.operating_reserve_required_fraction)
+total_or_required = chp_or_required + load_or_required
+
+# Get operating reserve provided
+or_provided = sum(results["ElectricLoad"]["offgrid_annual_oper_res_provided_series_kwh"])
+
+println("\nOperating Reserve Analysis:")
+println("  - CHP OR required: ", round(chp_or_required, digits=1), " kWh")
+println("  - Load OR required: ", round(load_or_required, digits=1), " kWh")
+println("  - Total OR required: ", round(total_or_required, digits=1), " kWh")
+println("  - Total OR provided: ", round(or_provided, digits=1), " kWh")
+println("  - OR margin: ", round(or_provided - total_or_required, digits=1), " kWh")
+
+# Test 1: Operating reserves provided >= operating reserves required (with tolerance for solver gap)
+@test or_provided >= total_or_required * (1 - 0.02)  # Allow 2% gap due to mip_rel_gap and rounding
+println("\n✓ Test 1 PASSED: Operating reserves provided (", round(or_provided, digits=1), 
+        " kWh) >= required (", round(total_or_required, digits=1), " kWh) within tolerance")
+
+# Test 2: CHP produces electricity in off-grid scenario
+chp_annual_production = sum(results["CHP"]["electric_to_load_series_kw"]) + 
+                        sum(results["CHP"]["electric_to_storage_series_kw"])
+@test chp_annual_production > 0
+println("✓ Test 2 PASSED: CHP produces electricity (", round(chp_annual_production, digits=1), " kWh)")
+
+# Test 3: Battery is sized to provide operating reserves
+battery_sized = results["ElectricStorage"]["size_kw"] > 0 || results["ElectricStorage"]["size_kwh"] > 0
+@test battery_sized
+println("✓ Test 3 PASSED: Battery is sized (", round(results["ElectricStorage"]["size_kw"], digits=1), 
+        " kW, ", round(results["ElectricStorage"]["size_kwh"], digits=1), " kWh)")
+
+# Test 4: Load is met according to minimum fraction
+@test results["ElectricLoad"]["offgrid_load_met_fraction"] >= s.electric_load.min_load_met_annual_fraction
+println("✓ Test 4 PASSED: Load met fraction (", round(results["ElectricLoad"]["offgrid_load_met_fraction"], digits=4), 
+        ") >= minimum (", s.electric_load.min_load_met_annual_fraction, ")")
+
+# Test 5: No grid interaction in off-grid scenario
+@test results["ElectricUtility"]["annual_energy_supplied_kwh"] ≈ 0.0 atol=0.01
+println("✓ Test 5 PASSED: No grid interaction (", results["ElectricUtility"]["annual_energy_supplied_kwh"], " kWh)")
+
+# Test 6: Financial calculation includes off-grid costs
+f = results["Financial"]
+lcc_check = f["lifecycle_generation_tech_capital_costs"] + f["lifecycle_storage_capital_costs"] + 
+            f["lifecycle_om_costs_after_tax"] + f["lifecycle_fuel_costs_after_tax"] + 
+            f["lifecycle_chp_standby_cost_after_tax"] + f["lifecycle_elecbill_after_tax"] + 
+            f["lifecycle_offgrid_other_annual_costs_after_tax"] + f["lifecycle_offgrid_other_capital_costs"] + 
+            f["lifecycle_outage_cost"] + f["lifecycle_MG_upgrade_and_fuel_cost"] - 
+            f["lifecycle_production_incentive_after_tax"]
+@test lcc_check ≈ f["lcc"] atol=1.0
+println("✓ Test 6 PASSED: Financial calculations consistent (LCC: \$", round(f["lcc"], digits=0), ")")
+
+println("\n" * "="^80)
+println("All tests PASSED for Off-Grid CHP with Operating Reserves!")
+println("="^80 * "\n")