From 71e1c5adeab6165a8c3e6959582f5ba675f9f0bf Mon Sep 17 00:00:00 2001 From: wbecker Date: Fri, 28 Nov 2025 19:51:25 -0700 Subject: [PATCH 01/45] Add ramp rate limit parameter for CHP --- src/constraints/chp_constraints.jl | 21 +++++++++++++++++++++ src/core/chp.jl | 2 ++ 2 files changed, 23 insertions(+) diff --git a/src/constraints/chp_constraints.jl b/src/constraints/chp_constraints.jl index 0f25bf54f..b278c555a 100644 --- a/src/constraints/chp_constraints.jl +++ b/src/constraints/chp_constraints.jl @@ -135,6 +135,22 @@ function add_chp_rated_prod_constraint(m, p; _n="") end +function add_chp_ramp_rate_constraints(m, p; _n="") + # Ramp rate constraints limit how quickly CHP production can change between consecutive timesteps + # Ramp up constraint + @constraint(m, CHPRampUp[t in p.techs.chp, ts in p.time_steps[2:end]], + m[Symbol("dvRatedProduction"*_n)][t, ts] - m[Symbol("dvRatedProduction"*_n)][t, ts-1] <= + p.s.chp.ramp_rate_fraction_per_hour * m[Symbol("dvSize"*_n)][t] / p.s.settings.time_steps_per_hour + ) + + # Ramp down constraint + @constraint(m, CHPRampDown[t in p.techs.chp, ts in p.time_steps[2:end]], + m[Symbol("dvRatedProduction"*_n)][t, ts-1] - m[Symbol("dvRatedProduction"*_n)][t, ts] <= + p.s.chp.ramp_rate_fraction_per_hour * m[Symbol("dvSize"*_n)][t] / p.s.settings.time_steps_per_hour + ) +end + + """ add_chp_hourly_om_charges(m, p; _n="") @@ -197,6 +213,11 @@ function add_chp_constraints(m, p; _n="") add_chp_thermal_production_constraints(m, p; _n=_n) add_binCHPIsOnInTS_constraints(m, p; _n=_n) add_chp_rated_prod_constraint(m, p; _n=_n) + + # Add ramp rate constraints if ramp_rate_fraction_per_hour < 1.0 + if p.s.chp.ramp_rate_fraction_per_hour < 1.0 / p.s.settings.time_steps_per_hour + add_chp_ramp_rate_constraints(m, p; _n=_n) + end if p.s.chp.supplementary_firing_max_steam_ratio > 1.0 add_chp_supplementary_firing_constraints(m,p; _n=_n) diff --git a/src/core/chp.jl b/src/core/chp.jl index c59514a9f..121fb1aaa 100644 --- a/src/core/chp.jl +++ b/src/core/chp.jl @@ -30,6 +30,7 @@ conflict_res_min_allowable_fraction_of_max = 0.25 fuel_type::String = "natural_gas" # "restrict_to": ["natural_gas", "landfill_bio_gas", "propane", "diesel_oil"] om_cost_per_kw::Float64 = 0.0 # Annual CHP fixed operations and maintenance costs in \$/kw-yr om_cost_per_hr_per_kw_rated::Float64 = 0.0 # CHP non-fuel variable operations and maintenance costs in \$/hr/kw_rated + ramp_rate_fraction_per_hour::Float64 = 1.0 # Maximum rate of change in electric production per hour as a fraction of size_kw [kW/size_kw/hour]. supplementary_firing_capital_cost_per_kw::Float64 = 150.0 # Installed CHP supplementary firing system cost in \$/kW (based on rated electric power) supplementary_firing_max_steam_ratio::Float64 = 1.0 # Ratio of max fired steam to un-fired steam production. Relevant only for combustion_turbine prime_mover supplementary_firing_efficiency::Float64 = 0.92 # Thermal efficiency of the incremental steam production from supplementary firing. Relevant only for combustion_turbine prime_mover @@ -101,6 +102,7 @@ Base.@kwdef mutable struct CHP <: AbstractCHP fuel_type::String = "natural_gas" om_cost_per_kw::Float64 = 0.0 om_cost_per_hr_per_kw_rated::Float64 = 0.0 + ramp_rate_fraction_per_hour::Float64 = 1.0 electric_efficiency_half_load::Float64 = NaN # Assigned to electric_efficiency_full_load if not input thermal_efficiency_half_load::Float64 = NaN # Assigned to thermal_efficiency_full_load if not input supplementary_firing_capital_cost_per_kw::Float64 = 150.0 From e8c915331252e0d07bd9195ca14606a09c37b637 Mon Sep 17 00:00:00 2001 From: wbecker Date: Fri, 28 Nov 2025 19:52:13 -0700 Subject: [PATCH 02/45] Report curtailed CHP power (i.e. power produced but sent to a load bank) --- src/results/chp.jl | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/results/chp.jl b/src/results/chp.jl index 58f302804..b0480a858 100644 --- a/src/results/chp.jl +++ b/src/results/chp.jl @@ -67,9 +67,16 @@ function add_chp_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="") CHPtoBatt = zeros(length(p.time_steps)) end r["electric_to_storage_series_kw"] = round.(value.(CHPtoBatt), digits=3) + if p.s.chp.can_curtail + @expression(m, CHPtoCurtail[ts in p.time_steps], + sum(m[Symbol("dvCurtail"*_n)][t,ts] for t in p.techs.chp)) + else + CHPtoCurtail = zeros(length(p.time_steps)) + end + r["electric_curtailed_series_kw"] = round.(value.(CHPtoCurtail), digits=3) @expression(m, CHPtoLoad[ts in p.time_steps], sum(m[Symbol("dvRatedProduction"*_n)][t, ts] * p.production_factor[t, ts] * p.levelization_factor[t] - for t in p.techs.chp) - CHPtoBatt[ts] - CHPtoGrid[ts]) + for t in p.techs.chp) - CHPtoBatt[ts] - CHPtoGrid[ts] - CHPtoCurtail[ts]) r["electric_to_load_series_kw"] = round.(value.(CHPtoLoad), digits=3) # Thermal dispatch breakdown if !isempty(p.s.storage.types.hot) From 896d157fcd598d35582461c0f26e940157023c42 Mon Sep 17 00:00:00 2001 From: wbecker Date: Fri, 28 Nov 2025 19:58:06 -0700 Subject: [PATCH 03/45] Fix issue with dictkeys_tosymbols making booleans numeric --- src/core/utils.jl | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/core/utils.jl b/src/core/utils.jl index 95dc834a9..aa3e0e6a9 100644 --- a/src/core/utils.jl +++ b/src/core/utils.jl @@ -223,6 +223,17 @@ function dictkeys_tosymbols(d::Dict) end end end + # Convert numeric boolean values (0.0/1.0) to proper Bool type + if k in [ + "off_grid_flag", "add_soc_incentive", "include_climate_in_objective", + "include_health_in_objective", "include_export_cost_series_in_results" + ] && !isnothing(v) && !(typeof(v) <: Bool) + try + v = Bool(v) + catch + throw(@error("Unable to convert $k to Bool. Expected boolean or 0/1, got: $v")) + end + end if k in [ "fuel_cost_per_mmbtu", "wholesale_rate", "export_rate_beyond_net_metering_limit", # for ERP From 708250b6c7691a4b074f233a89f37c4d534e7613 Mon Sep 17 00:00:00 2001 From: wbecker Date: Fri, 28 Nov 2025 21:44:42 -0700 Subject: [PATCH 04/45] Add temp dev test files --- test/scenarios/nuclear_battery.json | 28 +++++ test/test_ramp.jl | 187 ++++++++++++++++++++++++++++ 2 files changed, 215 insertions(+) create mode 100644 test/scenarios/nuclear_battery.json create mode 100644 test/test_ramp.jl diff --git a/test/scenarios/nuclear_battery.json b/test/scenarios/nuclear_battery.json new file mode 100644 index 000000000..32204db44 --- /dev/null +++ b/test/scenarios/nuclear_battery.json @@ -0,0 +1,28 @@ +{ + "Site": { + "latitude": 37.78, + "longitude": -122.45 + }, + "ElectricLoad": {}, + "ElectricTariff": { + "blended_annual_energy_rate": 0.12, + "blended_annual_demand_rate": 10.0 + }, + "CHP": { + "fuel_cost_per_mmbtu": 0.1, + "prime_mover": "combustion_turbine", + "installed_cost_per_kw": 12000.0, + "om_cost_per_kw": 150.0, + "om_cost_per_kwh": 0.0, + "electric_efficiency_full_load": 0.25, + "electric_efficiency_half_load": 0.25, + "min_turn_down_fraction": 0.0, + "is_electric_only": true, + "can_curtail": false, + "unavailability_periods": [{"month":1, + "start_week_of_month":2, + "start_day_of_week":1, + "start_hour":1, + "duration_hours":0}] + } +} diff --git a/test/test_ramp.jl b/test/test_ramp.jl new file mode 100644 index 000000000..c4f558b89 --- /dev/null +++ b/test/test_ramp.jl @@ -0,0 +1,187 @@ +# using Revise +# using REopt +# using JSON +# using DelimitedFiles +# using PlotlyJS +# using Dates +# using Test +# using JuMP +# using HiGHS +# using DotEnv +# DotEnv.load!() + +# Notes to test for Off-Grid DC powered by Nuclear + Storage +# +# DONE: CHP electric only with ramp rate limit, no outage but force battery size to avoid utility serving the high ramp times +# Boiler + SteamTurbine (+ HighTempThermalStorage?) +# Add CHP and/or Boiler+ST to off-grid with SR requirements for nuclear +# Generator and Battery currently have unlimited SR **supply** +# PV has SR requirement - use that as proxy for nuclear and/or geothermal? +# Focus on analysis needs for project, but eventually consider GenericGenerator (+Heat?) with needed attributes + + +# Create a 15-minute interval load profile from hourly with sine wave sub-hourly variations +function upsample_load_variation(hourly_load; intra_hour_variation=0.2) + # Create 15-minute load profile with sine wave sub-hourly variations + fifteen_min_loads = Float64[] + for (hour_idx, hourly_load) in enumerate(hourly_load) + # Create 4 sub-hourly values per hour with sine wave variation + # Base load with ±10% sine wave variation within each hour + for quarter_hour in 1:4 + # Intra-hour phase angle for this quarter hour (0, π/2, π, 3π/2) + intra_hour_phase = (quarter_hour - 1) * π / 2 + # Sine wave variation: ±20% of hourly load + variation = intra_hour_variation * sin(intra_hour_phase + 2π * hour_idx / 24) # Daily cycle component + sub_hourly_load = hourly_load * (1.0 + variation) + push!(fifteen_min_loads, max(0.0, sub_hourly_load)) # Ensure non-negative + end + end + + return fifteen_min_loads +end + + +############### Nuclear as CHP with ramp rate and battery to support ################### + +# scenario = 1 # "Specify CHP+battery sizes 15-min load with variation, no outages" +# scenario = 2 # "Least-cost sizes CHP+battery for year long outage hourly load from CRB" +# scenario = 3 # "Force CHP/Nuclear to peak load size, but still size battery for ramp rate limit" +scenario = 3 +# Load scenario with CHP (electric only) cost and performance params similar to nuclear +input_data = JSON.parsefile("./scenarios/nuclear_battery.json") + +# New ramp rate input/constraint (fraction of capacity per hour) +input_data["CHP"]["ramp_rate_fraction_per_hour"] = 0.1 + +if scenario == 1 + # Create 15-minute load profile from hourly with sine wave intra-hour variation + sim_input = Dict( + "load_type" => "electric", + "doe_reference_name" => "Hospital", + "annual_kwh" => 4.0e6, + "latitude" => input_data["Site"]["latitude"], + "longitude" => input_data["Site"]["longitude"], + "year" => 2023) + hourly_loads_kw = simulated_load(sim_input)["loads_kw"] + fifteen_min_loads = upsample_load_variation(hourly_loads_kw; intra_hour_variation=0.05) + input_data["Settings"] = Dict("time_steps_per_hour" => 4) + input_data["ElectricLoad"] = Dict("loads_kw" => fifteen_min_loads, "year" => 2023) + + # You can adjust the fixed size here if desired + fixed_chp_size_kw = 800.0 + input_data["CHP"]["min_kw"] = fixed_chp_size_kw + input_data["CHP"]["max_kw"] = fixed_chp_size_kw + + # See if REopt sizes battery to support nuclear ramp rate limitation + input_data["ElectricStorage"] = Dict( + "min_kw" => 800.0, + "min_kwh" => 800.0, + "max_kw" => 800.0, + "max_kwh" => 800.0 + ) +elseif scenario == 2 || scenario == 3 + # Year-long outage scenario from CHP test scenarios + input_data["ElectricLoad"] = Dict("doe_reference_name" => "Hospital", "annual_kwh" => 4.0e6) + input_data["ElectricLoad"]["critical_load_fraction"] = 1.0 + input_data["ElectricUtility"] = Dict("outage_start_time_step" => 1, + "outage_end_time_step" => 8760) + input_data["ElectricStorage"] = Dict() + if scenario == 3 + # Force CHP to peak load size, but still size battery for ramp rate limit + # Estimate peak load from annual energy assuming capacity factor + estimated_peak_load_kw = 740.0 # Assuming 50% capacity factor + input_data["CHP"]["min_kw"] = estimated_peak_load_kw + input_data["CHP"]["max_kw"] = estimated_peak_load_kw + end +end + +# This was somehow being set to 0 in dictkeys_tosymbols in utils.py which was causing an error, +# but function was updated to avoid this +# input_data["Settings"]["off_grid_flag"] = false + +# Create scenario and inputs +s = Scenario(input_data) +inputs = REoptInputs(s) +ts_per_hour = s.settings.time_steps_per_hour + +# Run optimization with single model (no BAU comparison needed for fixed sizing) +m = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +results = run_reopt(m, inputs) + +# Print key results +println("\n===== Nuclear as CHP with ramp rate and battery to support =====") +println("\nCHP Results:") +println("CHP Size (kW): ", results["CHP"]["size_kw"]) +println("Annual CHP Production (kWh): ", round(results["CHP"]["annual_electric_production_kwh"], digits=0)) +println("CHP to Load (kWh): ", round(sum(results["CHP"]["electric_to_load_series_kw"])/ts_per_hour, digits=0)) +println("CHP to Battery (kWh): ", round(sum(results["CHP"]["electric_to_storage_series_kw"])/ts_per_hour, digits=0)) +println("CHP to Grid (kWh): ", round(sum(results["CHP"]["electric_to_grid_series_kw"])/ts_per_hour, digits=0)) +println("CHP Curtailed (kWh): ", round(sum(results["CHP"]["electric_curtailed_series_kw"])/ts_per_hour, digits=0)) +println("\nBattery Results:") +println("Battery Energy Size (kWh): ", results["ElectricStorage"]["size_kwh"]) +println("Battery Power Size (kW): ", results["ElectricStorage"]["size_kw"]) +println("Battery to Load (kWh): ", round(sum(results["ElectricStorage"]["storage_to_load_series_kw"])/ts_per_hour, digits=0)) +println("\nGrid Results:") +println("Utility to Load (kWh): ", round(sum(results["ElectricUtility"]["annual_energy_supplied_kwh"]), digits=0)) +# println("\nFinancial Results:") +# println("NPV (\$): ", round(results["Financial"]["npv"], digits=0)) +# println("Simple Payback (years): ", round(results["Financial"]["simple_payback_years"], digits=2)) +println("==========================================\n") + +# Create stacked area chart showing how CHP and battery serve the load +load_series = results["ElectricLoad"]["load_series_kw"] +chp_to_load = results["CHP"]["electric_to_load_series_kw"] +batt_to_load = results["ElectricStorage"]["storage_to_load_series_kw"] +grid_to_load = results["ElectricUtility"]["electric_to_load_series_kw"] + +# Create time steps array (in hours for x-axis) +time_hours = collect(1:length(load_series)) ./ ts_per_hour + +# Create stacked area chart +plt = plot( + [ + scatter( + x=time_hours, + y=grid_to_load, + mode="lines", + name="Grid to Load", + fill="tozeroy", + line=attr(width=0.5, color="rgb(128,128,128)"), + fillcolor="rgba(128,128,128,0.5)" + ), + scatter( + x=time_hours, + y=grid_to_load .+ chp_to_load, + mode="lines", + name="CHP to Load", + fill="tonexty", + line=attr(width=0.5, color="rgb(255,128,0)"), + fillcolor="rgba(255,128,0,0.6)" + ), + scatter( + x=time_hours, + y=grid_to_load .+ chp_to_load .+ batt_to_load, + mode="lines", + name="Battery to Load", + fill="tonexty", + line=attr(width=0.5, color="rgb(0,128,255)"), + fillcolor="rgba(0,128,255,0.6)" + ), + scatter( + x=time_hours, + y=load_series, + mode="lines", + name="Total Load", + line=attr(width=2, color="rgb(0,0,0)", dash="dot") + ) + ], + Layout( + title="Electric Load Service Breakdown - CHP with Ramp Rate and Battery Support", + xaxis=attr(title="Time (hours)", range=[0, 8760]), + yaxis=attr(title="Power (kW)"), + showlegend=true, + hovermode="x unified", + legend=attr(x=1.02, y=1, xanchor="left", yanchor="top") + ) +) +display(plt) \ No newline at end of file From 458f4d9ec7a1b629b20e3a5401fd1b07ef440a9c Mon Sep 17 00:00:00 2001 From: wbecker Date: Sat, 29 Nov 2025 11:14:14 -0700 Subject: [PATCH 05/45] Add CHP to off-grid analysis but with option to provide or require operating reserve --- .../operating_reserve_constraints.jl | 18 +++++++++++++++--- src/core/chp.jl | 16 +++++++++++++++- src/core/reopt_inputs.jl | 6 +++++- src/core/scenario.jl | 10 ++++++---- src/core/techs.jl | 7 +++++++ 5 files changed, 48 insertions(+), 9 deletions(-) diff --git a/src/constraints/operating_reserve_constraints.jl b/src/constraints/operating_reserve_constraints.jl index 0cf21ea8b..193c03293 100644 --- a/src/constraints/operating_reserve_constraints.jl +++ b/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) diff --git a/src/core/chp.jl b/src/core/chp.jl index c59514a9f..437f536ad 100644 --- a/src/core/chp.jl +++ b/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 diff --git a/src/core/reopt_inputs.jl b/src/core/reopt_inputs.jl index bbbd6d0f6..07916a750 100644 --- a/src/core/reopt_inputs.jl +++ b/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 diff --git a/src/core/scenario.jl b/src/core/scenario.jl index 30ea91200..e0c45bebc 100644 --- a/src/core/scenario.jl +++ b/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 diff --git a/src/core/techs.jl b/src/core/techs.jl index 10e6e966a..b8504e258 100644 --- a/src/core/techs.jl +++ b/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]) From b04c04bd19c7199ffaefb4bc758071be61245710 Mon Sep 17 00:00:00 2001 From: wbecker Date: Sat, 29 Nov 2025 11:19:57 -0700 Subject: [PATCH 06/45] Remove redundant CHP thermal result which was causing conflict with no heating load --- src/results/chp.jl | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/results/chp.jl b/src/results/chp.jl index 58f302804..242d05b42 100644 --- a/src/results/chp.jl +++ b/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], From 08eda3d50d1484883c8684422b8e9a3e8ac5d94f Mon Sep 17 00:00:00 2001 From: wbecker Date: Sat, 29 Nov 2025 11:21:08 -0700 Subject: [PATCH 07/45] Add temp dev test files for off-grid CHP --- test/scenarios/chp_offgrid.json | 34 +++++++++++ test/test_chp_offgrid.jl | 105 ++++++++++++++++++++++++++++++++ 2 files changed, 139 insertions(+) create mode 100644 test/scenarios/chp_offgrid.json create mode 100644 test/test_chp_offgrid.jl diff --git a/test/scenarios/chp_offgrid.json b/test/scenarios/chp_offgrid.json new file mode 100644 index 000000000..1d60c92f8 --- /dev/null +++ b/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 + } +} diff --git a/test/test_chp_offgrid.jl b/test/test_chp_offgrid.jl new file mode 100644 index 000000000..213de6b1a --- /dev/null +++ b/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") \ No newline at end of file From 1c4df386af315b0c411656d56e6011e0d1d0ae2b Mon Sep 17 00:00:00 2001 From: wbecker Date: Mon, 1 Dec 2025 22:02:54 -0700 Subject: [PATCH 08/45] Multiple CHPs inputs setup --- src/core/bau_inputs.jl | 6 +- src/core/bau_scenario.jl | 6 +- src/core/chp.jl | 7 ++ src/core/reopt_inputs.jl | 155 ++++++++++++++++++++------------------- src/core/scenario.jl | 55 ++++++++------ src/core/techs.jl | 24 +++--- 6 files changed, 141 insertions(+), 112 deletions(-) diff --git a/src/core/bau_inputs.jl b/src/core/bau_inputs.jl index 3464c39de..0f7311178 100644 --- a/src/core/bau_inputs.jl +++ b/src/core/bau_inputs.jl @@ -166,6 +166,9 @@ function BAUInputs(p::REoptInputs) heating_loads_served_by_tes = Dict{String,Array{String,1}}() unavailability = get_unavailability_by_tech(p.s, techs, p.time_steps) + # Initialize CHP-specific parameters as nested dictionary (empty for BAU) + chp_params = Dict{String, Dict{Symbol, Float64}}() + REoptInputs( bau_scenario, techs, @@ -234,7 +237,8 @@ function BAUInputs(p::REoptInputs) heating_loads_served_by_tes, unavailability, absorption_chillers_using_heating_load, - avoided_capex_by_ashp_present_value + avoided_capex_by_ashp_present_value, + chp_params ) end diff --git a/src/core/bau_scenario.jl b/src/core/bau_scenario.jl index a32b43880..264035f18 100644 --- a/src/core/bau_scenario.jl +++ b/src/core/bau_scenario.jl @@ -33,7 +33,8 @@ struct BAUScenario <: AbstractScenario cooling_load::CoolingLoad ghp_option_list::Array{Union{GHP, Nothing}, 1} # List of GHP objects (often just 1 element, but can be more) space_heating_thermal_load_reduction_with_ghp_kw::Union{Vector{Float64}, Nothing} - cooling_thermal_load_reduction_with_ghp_kw::Union{Vector{Float64}, Nothing} + cooling_thermal_load_reduction_with_ghp_kw::Union{Vector{Float64}, Nothing} + chps::Array{CHP, 1} # Empty array for BAU scenarios (no new CHP modeled) end @@ -150,6 +151,7 @@ function BAUScenario(s::Scenario) s.cooling_load, ghp_option_list, space_heating_thermal_load_reduction_with_ghp_kw, - cooling_thermal_load_reduction_with_ghp_kw + cooling_thermal_load_reduction_with_ghp_kw, + CHP[] # Empty array - no CHP in BAU scenario ) end \ No newline at end of file diff --git a/src/core/chp.jl b/src/core/chp.jl index c59514a9f..822322af5 100644 --- a/src/core/chp.jl +++ b/src/core/chp.jl @@ -80,6 +80,7 @@ conflict_res_min_allowable_fraction_of_max = 0.25 Base.@kwdef mutable struct CHP <: AbstractCHP # Required input fuel_cost_per_mmbtu::Union{<:Real, AbstractVector{<:Real}} = [] + name::String = "CHP" # for use with multiple CHPs # Inputs which defaults vary depending on prime_mover and size_class installed_cost_per_kw::Union{Float64, AbstractVector{Float64}} = Float64[] @@ -140,6 +141,7 @@ Base.@kwdef mutable struct CHP <: AbstractCHP emissions_factor_lb_NOx_per_mmbtu::Real = get(FUEL_DEFAULTS["emissions_factor_lb_NOx_per_mmbtu"],fuel_type,0) emissions_factor_lb_SO2_per_mmbtu::Real = get(FUEL_DEFAULTS["emissions_factor_lb_SO2_per_mmbtu"],fuel_type,0) emissions_factor_lb_PM25_per_mmbtu::Real = get(FUEL_DEFAULTS["emissions_factor_lb_PM25_per_mmbtu"],fuel_type,0) + fuel_cost_escalation_rate_fraction::Union{Nothing, Float64} = nothing end @@ -549,3 +551,8 @@ function get_size_class_from_size(chp_elec_size_heuristic_kw, class_bounds, n_cl end return size_class end + +# Get a specific CHP by name from an array of CHPs +function get_chp_by_name(name::String, chps::AbstractArray{CHP, 1}) + chps[findfirst(chp -> chp.name == name, chps)] +end diff --git a/src/core/reopt_inputs.jl b/src/core/reopt_inputs.jl index bbbd6d0f6..8eaf85f51 100644 --- a/src/core/reopt_inputs.jl +++ b/src/core/reopt_inputs.jl @@ -138,6 +138,8 @@ struct REoptInputs{ScenarioType <: AbstractScenario} <: AbstractInputs unavailability::Dict{String, Array{Float64,1}} # (techs.elec) absorption_chillers_using_heating_load::Dict{String,Array{String,1}} # ("AbsorptionChiller" or empty) avoided_capex_by_ashp_present_value::Dict{String, <:Real} # HVAC upgrade costs avoided (ASHP) + # CHP-specific parameters indexed by tech name, then parameter name + chp_params::Dict{String, Dict{Symbol, Float64}} # (techs.chp) end @@ -174,7 +176,8 @@ function REoptInputs(s::AbstractScenario) tech_renewable_energy_fraction, tech_emissions_factors_CO2, tech_emissions_factors_NOx, tech_emissions_factors_SO2, tech_emissions_factors_PM25, techs_operating_reserve_req_fraction, thermal_cop, fuel_cost_per_kwh, heating_cop, cooling_cop, heating_cf, cooling_cf, avoided_capex_by_ashp_present_value, - pbi_pwf, pbi_max_benefit, pbi_max_kw, pbi_benefit_per_kwh = setup_tech_inputs(s,time_steps) + pbi_pwf, pbi_max_benefit, pbi_max_kw, pbi_benefit_per_kwh, + chp_params = setup_tech_inputs(s,time_steps) months = 1:12 @@ -327,7 +330,8 @@ function REoptInputs(s::AbstractScenario) heating_loads_served_by_tes, unavailability, absorption_chillers_using_heating_load, - avoided_capex_by_ashp_present_value + avoided_capex_by_ashp_present_value, + chp_params ) end @@ -365,6 +369,9 @@ function setup_tech_inputs(s::AbstractScenario, time_steps) cooling_cf = Dict(t => zeros(length(time_steps)) for t in techs.cooling) cooling_cop = Dict(t => zeros(length(time_steps)) for t in techs.cooling) avoided_capex_by_ashp_present_value = Dict(t => 0.0 for t in techs.all) + + # Initialize CHP-specific parameters as nested dictionary + chp_params = Dict{String, Dict{Symbol, Float64}}() pbi_pwf = Dict{String, Any}() pbi_max_benefit = Dict{String, Any}() @@ -417,14 +424,15 @@ function setup_tech_inputs(s::AbstractScenario, time_steps) tech_renewable_energy_fraction, om_cost_per_kw, production_factor, fuel_cost_per_kwh, heating_cf) end - if "CHP" in techs.all + if !isempty(techs.chp) setup_chp_inputs(s, max_sizes, min_sizes, cap_cost_slope, om_cost_per_kw, production_factor, techs_by_exportbin, techs.segmented, n_segs_by_tech, seg_min_size, seg_max_size, seg_yint, techs, tech_renewable_energy_fraction, tech_emissions_factors_CO2, tech_emissions_factors_NOx, tech_emissions_factors_SO2, tech_emissions_factors_PM25, fuel_cost_per_kwh, - heating_cf, pbi_pwf, pbi_max_benefit, pbi_max_kw, pbi_benefit_per_kwh) - end - + heating_cf, pbi_pwf, pbi_max_benefit, pbi_max_kw, pbi_benefit_per_kwh, + chp_params) + end + if "ExistingChiller" in techs.all setup_existing_chiller_inputs(s, max_sizes, min_sizes, existing_sizes, cap_cost_slope, cooling_cop, cooling_cf) else @@ -496,7 +504,8 @@ function setup_tech_inputs(s::AbstractScenario, time_steps) tech_renewable_energy_fraction, tech_emissions_factors_CO2, tech_emissions_factors_NOx, tech_emissions_factors_SO2, tech_emissions_factors_PM25, techs_operating_reserve_req_fraction, thermal_cop, fuel_cost_per_kwh, heating_cop, cooling_cop, heating_cf, cooling_cf, avoided_capex_by_ashp_present_value, - pbi_pwf, pbi_max_benefit, pbi_max_kw, pbi_benefit_per_kwh + pbi_pwf, pbi_max_benefit, pbi_max_kw, pbi_benefit_per_kwh, + chp_params end @@ -638,43 +647,6 @@ function setup_pv_inputs(s::AbstractScenario, max_sizes, min_sizes, return nothing end -""" - function setup_chp_inputs(s::AbstractScenario, max_sizes, min_sizes, cap_cost_slope, om_cost_per_kw, - production_factor, techs_by_exportbin, segmented_techs, n_segs_by_tech, seg_min_size, seg_max_size, seg_yint, techs, - tech_renewable_energy_fraction, tech_emissions_factors_CO2, tech_emissions_factors_NOx, tech_emissions_factors_SO2, tech_emissions_factors_PM25, fuel_cost_per_kwh, - heating_cf - ) - -Update tech-indexed data arrays necessary to build the JuMP model with the values for CHP. -""" -function setup_chp_inputs(s::AbstractScenario, max_sizes, min_sizes, cap_cost_slope, om_cost_per_kw, - production_factor, techs_by_exportbin, segmented_techs, n_segs_by_tech, seg_min_size, seg_max_size, seg_yint, techs, - tech_renewable_energy_fraction, tech_emissions_factors_CO2, tech_emissions_factors_NOx, tech_emissions_factors_SO2, tech_emissions_factors_PM25, fuel_cost_per_kwh, - heating_cf - ) - max_sizes["CHP"] = s.chp.max_kw - min_sizes["CHP"] = s.chp.min_kw - update_cost_curve!(s.chp, "CHP", s.financial, - cap_cost_slope, segmented_techs, n_segs_by_tech, seg_min_size, seg_max_size, seg_yint - ) - om_cost_per_kw["CHP"] = s.chp.om_cost_per_kw - production_factor["CHP", :] = get_production_factor(s.chp, s.electric_load.year, s.electric_utility.outage_start_time_step, - s.electric_utility.outage_end_time_step, s.settings.time_steps_per_hour) - fillin_techs_by_exportbin(techs_by_exportbin, s.chp, "CHP") - if !s.chp.can_curtail - push!(techs.no_curtail, "CHP") - end - tech_renewable_energy_fraction["CHP"] = s.chp.fuel_renewable_energy_fraction - tech_emissions_factors_CO2["CHP"] = s.chp.emissions_factor_lb_CO2_per_mmbtu / KWH_PER_MMBTU # lb/mmtbu * mmtbu/kWh - tech_emissions_factors_NOx["CHP"] = s.chp.emissions_factor_lb_NOx_per_mmbtu / KWH_PER_MMBTU - tech_emissions_factors_SO2["CHP"] = s.chp.emissions_factor_lb_SO2_per_mmbtu / KWH_PER_MMBTU - tech_emissions_factors_PM25["CHP"] = s.chp.emissions_factor_lb_PM25_per_mmbtu / KWH_PER_MMBTU - chp_fuel_cost_per_kwh = s.chp.fuel_cost_per_mmbtu ./ KWH_PER_MMBTU - fuel_cost_per_kwh["CHP"] = per_hour_value_to_time_series(chp_fuel_cost_per_kwh, s.settings.time_steps_per_hour, "CHP") - heating_cf["CHP"] = ones(8760*s.settings.time_steps_per_hour) - return nothing -end - function setup_wind_inputs(s::AbstractScenario, max_sizes, min_sizes, existing_sizes, cap_cost_slope, om_cost_per_kw, production_factor, techs_by_exportbin, @@ -864,24 +836,32 @@ function setup_absorption_chiller_inputs(s::AbstractScenario, max_sizes, min_siz cooling_cop["AbsorptionChiller"] .= s.absorption_chiller.cop_electric cooling_cf["AbsorptionChiller"] .= 1.0 - if isnothing(s.chp) + if isempty(s.chps) thermal_factor = 1.0 - elseif s.chp.cooling_thermal_factor == 0.0 - throw(@error("The CHP cooling_thermal_factor is 0.0 which implies that CHP cannot serve AbsorptionChiller. If you - want to model CHP and AbsorptionChiller, you must specify a cooling_thermal_factor greater than 0.0")) else - thermal_factor = s.chp.cooling_thermal_factor + # Use the cooling_thermal_factor from the first CHP + if s.chps[1].cooling_thermal_factor == 0.0 + throw(@error("The CHP cooling_thermal_factor is 0.0 which implies that CHP cannot serve AbsorptionChiller. If you + want to model CHP and AbsorptionChiller, you must specify a cooling_thermal_factor greater than 0.0")) + else + thermal_factor = s.chps[1].cooling_thermal_factor + end end thermal_cop["AbsorptionChiller"] = s.absorption_chiller.cop_thermal * thermal_factor om_cost_per_kw["AbsorptionChiller"] = s.absorption_chiller.om_cost_per_kw return nothing end + """ function setup_chp_inputs(s::AbstractScenario, max_sizes, min_sizes, cap_cost_slope, om_cost_per_kw, production_factor, techs_by_exportbin, segmented_techs, n_segs_by_tech, seg_min_size, seg_max_size, seg_yint, techs, tech_renewable_energy_fraction, tech_emissions_factors_CO2, tech_emissions_factors_NOx, tech_emissions_factors_SO2, tech_emissions_factors_PM25, fuel_cost_per_kwh, - heating_cf, pbi_pwf, pbi_max_benefit, pbi_max_kw, pbi_benefit_per_kwh + heating_cf, pbi_pwf, pbi_max_benefit, pbi_max_kw, pbi_benefit_per_kwh, + chp_electric_efficiency_full_load, chp_electric_efficiency_half_load, + chp_thermal_efficiency_full_load, chp_thermal_efficiency_half_load, + chp_min_turn_down_fraction, chp_supplementary_firing_efficiency, chp_supplementary_firing_max_steam_ratio, + chp_om_cost_per_kwh ) Update tech-indexed data arrays necessary to build the JuMP model with the values for CHP. @@ -889,29 +869,44 @@ Update tech-indexed data arrays necessary to build the JuMP model with the value function setup_chp_inputs(s::AbstractScenario, max_sizes, min_sizes, cap_cost_slope, om_cost_per_kw, production_factor, techs_by_exportbin, segmented_techs, n_segs_by_tech, seg_min_size, seg_max_size, seg_yint, techs, tech_renewable_energy_fraction, tech_emissions_factors_CO2, tech_emissions_factors_NOx, tech_emissions_factors_SO2, tech_emissions_factors_PM25, fuel_cost_per_kwh, - heating_cf, pbi_pwf, pbi_max_benefit, pbi_max_kw, pbi_benefit_per_kwh + heating_cf, pbi_pwf, pbi_max_benefit, pbi_max_kw, pbi_benefit_per_kwh, + chp_params ) - max_sizes["CHP"] = s.chp.max_kw - min_sizes["CHP"] = s.chp.min_kw - update_cost_curve!(s.chp, "CHP", s.financial, - cap_cost_slope, segmented_techs, n_segs_by_tech, seg_min_size, seg_max_size, seg_yint - ) - om_cost_per_kw["CHP"] = s.chp.om_cost_per_kw - production_factor["CHP", :] = get_production_factor(s.chp, s.electric_load.year, s.electric_utility.outage_start_time_step, - s.electric_utility.outage_end_time_step, s.settings.time_steps_per_hour) - fillin_techs_by_exportbin(techs_by_exportbin, s.chp, "CHP") - if !s.chp.can_curtail - push!(techs.no_curtail, "CHP") - end - tech_renewable_energy_fraction["CHP"] = s.chp.fuel_renewable_energy_fraction - tech_emissions_factors_CO2["CHP"] = s.chp.emissions_factor_lb_CO2_per_mmbtu / KWH_PER_MMBTU # lb/mmtbu * mmtbu/kWh - tech_emissions_factors_NOx["CHP"] = s.chp.emissions_factor_lb_NOx_per_mmbtu / KWH_PER_MMBTU - tech_emissions_factors_SO2["CHP"] = s.chp.emissions_factor_lb_SO2_per_mmbtu / KWH_PER_MMBTU - tech_emissions_factors_PM25["CHP"] = s.chp.emissions_factor_lb_PM25_per_mmbtu / KWH_PER_MMBTU - chp_fuel_cost_per_kwh = s.chp.fuel_cost_per_mmbtu ./ KWH_PER_MMBTU - fuel_cost_per_kwh["CHP"] = per_hour_value_to_time_series(chp_fuel_cost_per_kwh, s.settings.time_steps_per_hour, "CHP") - heating_cf["CHP"] = ones(8760*s.settings.time_steps_per_hour) - setup_pbi_inputs!(techs, s.chp, "CHP", s.financial, pbi_pwf, pbi_max_benefit, pbi_max_kw, pbi_benefit_per_kwh) + for chp in s.chps + max_sizes[chp.name] = chp.max_kw + min_sizes[chp.name] = chp.min_kw + update_cost_curve!(chp, chp.name, s.financial, + cap_cost_slope, segmented_techs, n_segs_by_tech, seg_min_size, seg_max_size, seg_yint + ) + om_cost_per_kw[chp.name] = chp.om_cost_per_kw + production_factor[chp.name, :] = get_production_factor(chp, s.electric_load.year, s.electric_utility.outage_start_time_step, + s.electric_utility.outage_end_time_step, s.settings.time_steps_per_hour) + fillin_techs_by_exportbin(techs_by_exportbin, chp, chp.name) + if !chp.can_curtail + push!(techs.no_curtail, chp.name) + end + tech_renewable_energy_fraction[chp.name] = chp.fuel_renewable_energy_fraction + tech_emissions_factors_CO2[chp.name] = chp.emissions_factor_lb_CO2_per_mmbtu / KWH_PER_MMBTU # lb/mmtbu * mmtbu/kWh + tech_emissions_factors_NOx[chp.name] = chp.emissions_factor_lb_NOx_per_mmbtu / KWH_PER_MMBTU + tech_emissions_factors_SO2[chp.name] = chp.emissions_factor_lb_SO2_per_mmbtu / KWH_PER_MMBTU + tech_emissions_factors_PM25[chp.name] = chp.emissions_factor_lb_PM25_per_mmbtu / KWH_PER_MMBTU + chp_fuel_cost_per_kwh = chp.fuel_cost_per_mmbtu ./ KWH_PER_MMBTU + fuel_cost_per_kwh[chp.name] = per_hour_value_to_time_series(chp_fuel_cost_per_kwh, s.settings.time_steps_per_hour, chp.name) + heating_cf[chp.name] = ones(8760*s.settings.time_steps_per_hour) + setup_pbi_inputs!(techs, chp, chp.name, s.financial, pbi_pwf, pbi_max_benefit, pbi_max_kw, pbi_benefit_per_kwh) + + # Store CHP-specific efficiency and operational parameters in nested dict + chp_params[chp.name] = Dict{Symbol, Float64}( + :electric_efficiency_full_load => chp.electric_efficiency_full_load, + :electric_efficiency_half_load => chp.electric_efficiency_half_load, + :thermal_efficiency_full_load => chp.thermal_efficiency_full_load, + :thermal_efficiency_half_load => chp.thermal_efficiency_half_load, + :min_turn_down_fraction => chp.min_turn_down_fraction, + :supplementary_firing_efficiency => chp.supplementary_firing_efficiency, + :supplementary_firing_max_steam_ratio => chp.supplementary_firing_max_steam_ratio, + :om_cost_per_kwh => chp.om_cost_per_kwh + ) + end return nothing end @@ -1123,10 +1118,16 @@ function setup_present_worth_factors(s::AbstractScenario, techs::Techs) s.financial.offtaker_discount_rate_fraction ) end - if t == "CHP" - pwf_fuel["CHP"] = annuity( + if t in techs.chp + # Get the CHP object to check for custom fuel_cost_escalation_rate_fraction + chp = get_chp_by_name(t, s.chps) + # Use CHP-specific escalation rate if provided, otherwise use default from Financial + escalation_rate = isnothing(chp.fuel_cost_escalation_rate_fraction) ? + s.financial.chp_fuel_cost_escalation_rate_fraction : + chp.fuel_cost_escalation_rate_fraction + pwf_fuel[t] = annuity( s.financial.analysis_years, - s.financial.chp_fuel_cost_escalation_rate_fraction, + escalation_rate, s.financial.offtaker_discount_rate_fraction ) end @@ -1390,7 +1391,9 @@ function get_unavailability_by_tech(s::AbstractScenario, techs::Techs, time_step if !isempty(techs.elec) unavailability = Dict(tech => zeros(length(time_steps)) for tech in techs.elec) if !isempty(techs.chp) - unavailability["CHP"] = [s.chp.unavailability_hourly[i] for i in 1:8760 for _ in 1:s.settings.time_steps_per_hour] + for chp in s.chps + unavailability[chp.name] = [chp.unavailability_hourly[i] for i in 1:8760 for _ in 1:s.settings.time_steps_per_hour] + end end else unavailability = Dict(""=>Float64[]) diff --git a/src/core/scenario.jl b/src/core/scenario.jl index 30ea91200..570dd3338 100644 --- a/src/core/scenario.jl +++ b/src/core/scenario.jl @@ -16,7 +16,7 @@ struct Scenario <: AbstractScenario cooling_load::CoolingLoad existing_boiler::Union{ExistingBoiler, Nothing} boiler::Union{Boiler, Nothing} - chp::Union{CHP, Nothing} # use nothing for more items when they are not modeled? + chps::Array{CHP, 1} flexible_hvac::Union{FlexibleHVAC, Nothing} existing_chiller::Union{ExistingChiller, Nothing} absorption_chiller::Union{AbsorptionChiller, Nothing} @@ -399,28 +399,41 @@ function Scenario(d::Dict; flex_hvac_from_json=false) end - chp = nothing + chps = CHP[] chp_prime_mover = nothing if haskey(d, "CHP") - electric_only = get(d["CHP"], "is_electric_only", false) || get(d["CHP"], "thermal_efficiency_full_load", 0.5) == 0.0 - if !isnothing(existing_boiler) && !electric_only - total_fuel_heating_load_mmbtu_per_hour = (space_heating_load.loads_kw + dhw_load.loads_kw + process_heat_load.loads_kw) / existing_boiler.efficiency / KWH_PER_MMBTU - avg_boiler_fuel_load_mmbtu_per_hour = sum(total_fuel_heating_load_mmbtu_per_hour) / length(total_fuel_heating_load_mmbtu_per_hour) - chp = CHP(d["CHP"]; - avg_boiler_fuel_load_mmbtu_per_hour = avg_boiler_fuel_load_mmbtu_per_hour, - existing_boiler = existing_boiler, - electric_load_series_kw = electric_load.loads_kw, - year = electric_load.year, - sector = site.sector, - federal_procurement_type = site.federal_procurement_type) - else # Only if modeling CHP without heating_load and existing_boiler (for prime generator, electric-only) - 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) + chp_array = isa(d["CHP"], AbstractArray) ? d["CHP"] : [d["CHP"]] + for (i, chp_dict) in enumerate(chp_array) + electric_only = get(chp_dict, "is_electric_only", false) || get(chp_dict, "thermal_efficiency_full_load", 0.5) == 0.0 + + # Set default name if not provided + if !haskey(chp_dict, "name") + chp_dict["name"] = length(chp_array) > 1 ? "CHP$i" : "CHP" + end + + if !isnothing(existing_boiler) && !electric_only + total_fuel_heating_load_mmbtu_per_hour = (space_heating_load.loads_kw + dhw_load.loads_kw + process_heat_load.loads_kw) / existing_boiler.efficiency / KWH_PER_MMBTU + avg_boiler_fuel_load_mmbtu_per_hour = sum(total_fuel_heating_load_mmbtu_per_hour) / length(total_fuel_heating_load_mmbtu_per_hour) + chp = CHP(chp_dict; + avg_boiler_fuel_load_mmbtu_per_hour = avg_boiler_fuel_load_mmbtu_per_hour, + existing_boiler = existing_boiler, + electric_load_series_kw = electric_load.loads_kw, + year = electric_load.year, + sector = site.sector, + federal_procurement_type = site.federal_procurement_type) + else # Only if modeling CHP without heating_load and existing_boiler (for prime generator, electric-only) + chp = CHP(chp_dict, + electric_load_series_kw = electric_load.loads_kw, + year = electric_load.year, + sector = site.sector, + federal_procurement_type = site.federal_procurement_type) + end + push!(chps, chp) + end + # Store first CHP's prime mover for backward compatibility if needed + if !isempty(chps) + chp_prime_mover = chps[1].prime_mover end - chp_prime_mover = chp.prime_mover end max_cooling_demand_kw = 0 @@ -1022,7 +1035,7 @@ function Scenario(d::Dict; flex_hvac_from_json=false) cooling_load, existing_boiler, boiler, - chp, + chps, flexible_hvac, existing_chiller, absorption_chiller, diff --git a/src/core/techs.jl b/src/core/techs.jl index 10e6e966a..21e5c682a 100644 --- a/src/core/techs.jl +++ b/src/core/techs.jl @@ -186,21 +186,21 @@ function Techs(s::Scenario) end end - if !isnothing(s.chp) - push!(all_techs, "CHP") - push!(elec, "CHP") - push!(chp_techs, "CHP") - if s.chp.can_supply_steam_turbine - push!(techs_can_supply_steam_turbine, "CHP") + for chp in s.chps + push!(all_techs, chp.name) + push!(chp_techs, chp.name) + push!(elec, chp.name) + if chp.can_supply_steam_turbine + push!(techs_can_supply_steam_turbine, chp.name) end - if s.chp.can_serve_space_heating - push!(techs_can_serve_space_heating, "CHP") + if chp.can_serve_space_heating + push!(techs_can_serve_space_heating, chp.name) end - if s.chp.can_serve_dhw - push!(techs_can_serve_dhw, "CHP") + if chp.can_serve_dhw + push!(techs_can_serve_dhw, chp.name) end - if s.chp.can_serve_process_heat - push!(techs_can_serve_process_heat, "CHP") + if chp.can_serve_process_heat + push!(techs_can_serve_process_heat, chp.name) end end From 550f953146e238319e2893017d2bdccb26f45033 Mon Sep 17 00:00:00 2001 From: wbecker Date: Mon, 1 Dec 2025 22:03:43 -0700 Subject: [PATCH 09/45] Multiple CHP constraints indexing --- src/constraints/chp_constraints.jl | 278 ++++++++++-------- src/constraints/cost_curve_constraints.jl | 62 ++-- .../electric_utility_constraints.jl | 11 +- src/constraints/outage_constraints.jl | 4 +- 4 files changed, 201 insertions(+), 154 deletions(-) diff --git a/src/constraints/chp_constraints.jl b/src/constraints/chp_constraints.jl index 0f25bf54f..acacb4024 100644 --- a/src/constraints/chp_constraints.jl +++ b/src/constraints/chp_constraints.jl @@ -1,129 +1,157 @@ # REopt®, Copyright (c) Alliance for Sustainable Energy, LLC. See also https://github.com/NREL/REopt.jl/blob/master/LICENSE. function add_chp_fuel_burn_constraints(m, p; _n="") - # Fuel burn slope and intercept - fuel_burn_slope, fuel_burn_intercept = fuel_slope_and_intercept(; - electric_efficiency_full_load = p.s.chp.electric_efficiency_full_load, - electric_efficiency_half_load = p.s.chp.electric_efficiency_half_load, - fuel_higher_heating_value_kwh_per_unit=1 - ) + # Loop through each CHP and add constraints with tech-specific parameters + for t in p.techs.chp + # Fuel burn slope and intercept for this specific CHP + fuel_burn_slope, fuel_burn_intercept = fuel_slope_and_intercept(; + electric_efficiency_full_load = p.chp_params[t][:electric_efficiency_full_load], + electric_efficiency_half_load = p.chp_params[t][:electric_efficiency_half_load], + fuel_higher_heating_value_kwh_per_unit=1 + ) - # Fuel cost - m[:TotalCHPFuelCosts] = @expression(m, - sum(p.pwf_fuel[t] * m[:dvFuelUsage][t, ts] * p.fuel_cost_per_kwh[t][ts] for t in p.techs.chp, ts in p.time_steps) - ) - # Conditionally add dvFuelBurnYIntercept if coefficient p.FuelBurnYIntRate is greater than ~zero - if abs(fuel_burn_intercept) > 1.0E-7 - dv = "dvFuelBurnYIntercept"*_n - m[Symbol(dv)] = @variable(m, [p.techs.chp, p.time_steps], base_name=dv) - - #Constraint (1c1): Total Fuel burn for CHP **with** y-intercept fuel burn and supplementary firing - @constraint(m, CHPFuelBurnCon[t in p.techs.chp, ts in p.time_steps], - m[Symbol("dvFuelUsage"*_n)][t,ts] == p.hours_per_time_step * ( - m[Symbol("dvFuelBurnYIntercept"*_n)][t,ts] + - p.production_factor[t,ts] * fuel_burn_slope * m[Symbol("dvRatedProduction"*_n)][t,ts] + - m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] / p.s.chp.supplementary_firing_efficiency + # Conditionally add dvFuelBurnYIntercept if coefficient fuel_burn_intercept is greater than ~zero + if abs(fuel_burn_intercept) > 1.0E-7 + if !haskey(m, Symbol("dvFuelBurnYIntercept"*_n)) + dv = "dvFuelBurnYIntercept"*_n + m[Symbol(dv)] = @variable(m, [p.techs.chp, p.time_steps], base_name=dv) + end + + #Constraint (1c1): Total Fuel burn for CHP **with** y-intercept fuel burn and supplementary firing + @constraint(m, [ts in p.time_steps], + m[Symbol("dvFuelUsage"*_n)][t,ts] == p.hours_per_time_step * ( + m[Symbol("dvFuelBurnYIntercept"*_n)][t,ts] + + p.production_factor[t,ts] * fuel_burn_slope * m[Symbol("dvRatedProduction"*_n)][t,ts] + + m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] / p.chp_params[t][:supplementary_firing_efficiency] + ) ) - ) - #Constraint (1d): Y-intercept fuel burn for CHP - @constraint(m, CHPFuelBurnYIntCon[t in p.techs.chp, ts in p.time_steps], - fuel_burn_intercept * m[Symbol("dvSize"*_n)][t] - p.s.chp.max_kw * - (1-m[Symbol("binCHPIsOnInTS"*_n)][t,ts]) <= m[Symbol("dvFuelBurnYIntercept"*_n)][t,ts] - ) - else - #Constraint (1c2): Total Fuel burn for CHP **without** y-intercept fuel burn - @constraint(m, CHPFuelBurnConLinear[t in p.techs.chp, ts in p.time_steps], - m[Symbol("dvFuelUsage"*_n)][t,ts] == p.hours_per_time_step * ( - p.production_factor[t,ts] * fuel_burn_slope * m[Symbol("dvRatedProduction"*_n)][t,ts] + - m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] / p.s.chp.supplementary_firing_efficiency + #Constraint (1d): Y-intercept fuel burn for CHP + @constraint(m, [ts in p.time_steps], + fuel_burn_intercept * m[Symbol("dvSize"*_n)][t] - p.max_sizes[t] * + (1-m[Symbol("binCHPIsOnInTS"*_n)][t,ts]) <= m[Symbol("dvFuelBurnYIntercept"*_n)][t,ts] + ) + else + #Constraint (1c2): Total Fuel burn for CHP **without** y-intercept fuel burn + @constraint(m, [ts in p.time_steps], + m[Symbol("dvFuelUsage"*_n)][t,ts] == p.hours_per_time_step * ( + p.production_factor[t,ts] * fuel_burn_slope * m[Symbol("dvRatedProduction"*_n)][t,ts] + + m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] / p.chp_params[t][:supplementary_firing_efficiency] + ) ) - ) - end + end + end + + # Fuel cost (sum over all CHPs) + m[:TotalCHPFuelCosts] = @expression(m, + sum(p.pwf_fuel[t] * m[:dvFuelUsage][t, ts] * p.fuel_cost_per_kwh[t][ts] for t in p.techs.chp, ts in p.time_steps) + ) end function add_chp_thermal_production_constraints(m, p; _n="") - # Thermal production slope and intercept - thermal_prod_full_load = 1.0 / p.s.chp.electric_efficiency_full_load * p.s.chp.thermal_efficiency_full_load # [kWt/kWe] - thermal_prod_half_load = 0.5 / p.s.chp.electric_efficiency_half_load * p.s.chp.thermal_efficiency_half_load # [kWt/kWe] - thermal_prod_slope = (thermal_prod_full_load - thermal_prod_half_load) / (1.0 - 0.5) # [kWt/kWe] - thermal_prod_intercept = thermal_prod_full_load - thermal_prod_slope * 1.0 # [kWt/kWe_rated - - - # Conditionally add dvHeatingProductionYIntercept if coefficient p.s.chpThermalProdIntercept is greater than ~zero - if abs(thermal_prod_intercept) > 1.0E-7 - dv = "dvHeatingProductionYIntercept"*_n - m[Symbol(dv)] = @variable(m, [p.techs.chp, p.time_steps], base_name=dv) + # Loop through each CHP and add constraints with tech-specific thermal production parameters + for t in p.techs.chp + # Thermal production slope and intercept for this specific CHP + thermal_prod_full_load = 1.0 / p.chp_params[t][:electric_efficiency_full_load] * p.chp_params[t][:thermal_efficiency_full_load] # [kWt/kWe] + thermal_prod_half_load = 0.5 / p.chp_params[t][:electric_efficiency_half_load] * p.chp_params[t][:thermal_efficiency_half_load] # [kWt/kWe] + thermal_prod_slope = (thermal_prod_full_load - thermal_prod_half_load) / (1.0 - 0.5) # [kWt/kWe] + thermal_prod_intercept = thermal_prod_full_load - thermal_prod_slope * 1.0 # [kWt/kWe_rated] + + # Conditionally add dvHeatingProductionYIntercept if coefficient thermal_prod_intercept is greater than ~zero + if abs(thermal_prod_intercept) > 1.0E-7 + if !haskey(m, Symbol("dvHeatingProductionYIntercept"*_n)) + dv = "dvHeatingProductionYIntercept"*_n + m[Symbol(dv)] = @variable(m, [p.techs.chp, p.time_steps], base_name=dv) + end - #Constraint (2a-1): Upper Bounds on Thermal Production Y-Intercept - @constraint(m, CHPYInt2a1Con[t in p.techs.chp, ts in p.time_steps], - m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts] <= thermal_prod_intercept * m[Symbol("dvSize"*_n)][t] - ) - # Constraint (2a-2): Upper Bounds on Thermal Production Y-Intercept - @constraint(m, CHPYInt2a2Con[t in p.techs.chp, ts in p.time_steps], - m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts] <= thermal_prod_intercept * p.s.chp.max_kw - * m[Symbol("binCHPIsOnInTS"*_n)][t,ts] - ) - #Constraint (2b): Lower Bounds on Thermal Production Y-Intercept - @constraint(m, CHPYInt2bCon[t in p.techs.chp, ts in p.time_steps], - m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts] >= thermal_prod_intercept * m[Symbol("dvSize"*_n)][t] - - thermal_prod_intercept * p.s.chp.max_kw * (1 - m[Symbol("binCHPIsOnInTS"*_n)][t,ts]) - ) - # Constraint (2c): Thermal Production of CHP - # Note: p.HotWaterAmbientFactor[t,ts] * p.HotWaterThermalFactor[t,ts] removed from this but present in math - @constraint(m, CHPThermalProductionCon[t in p.techs.chp, ts in p.time_steps], - sum(m[Symbol("dvHeatingProduction"*_n)][t,q,ts] for q in p.heating_loads) == - thermal_prod_slope * p.production_factor[t,ts] * m[Symbol("dvRatedProduction"*_n)][t,ts] - + m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts] + - m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] - ) - else - @constraint(m, CHPThermalProductionConLinear[t in p.techs.chp, ts in p.time_steps], + #Constraint (2a-1): Upper Bounds on Thermal Production Y-Intercept + @constraint(m, [ts in p.time_steps], + m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts] <= thermal_prod_intercept * m[Symbol("dvSize"*_n)][t] + ) + # Constraint (2a-2): Upper Bounds on Thermal Production Y-Intercept + @constraint(m, [ts in p.time_steps], + m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts] <= thermal_prod_intercept * p.max_sizes[t] + * m[Symbol("binCHPIsOnInTS"*_n)][t,ts] + ) + #Constraint (2b): Lower Bounds on Thermal Production Y-Intercept + @constraint(m, [ts in p.time_steps], + m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts] >= thermal_prod_intercept * m[Symbol("dvSize"*_n)][t] + - thermal_prod_intercept * p.max_sizes[t] * (1 - m[Symbol("binCHPIsOnInTS"*_n)][t,ts]) + ) + # Constraint (2c): Thermal Production of CHP + @constraint(m, [ts in p.time_steps], sum(m[Symbol("dvHeatingProduction"*_n)][t,q,ts] for q in p.heating_loads) == - thermal_prod_slope * p.production_factor[t,ts] * m[Symbol("dvRatedProduction"*_n)][t,ts] + - m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] - ) + thermal_prod_slope * p.production_factor[t,ts] * m[Symbol("dvRatedProduction"*_n)][t,ts] + + m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts] + + m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] + ) + else + @constraint(m, [ts in p.time_steps], + sum(m[Symbol("dvHeatingProduction"*_n)][t,q,ts] for q in p.heating_loads) == + thermal_prod_slope * p.production_factor[t,ts] * m[Symbol("dvRatedProduction"*_n)][t,ts] + + m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] + ) + end end - end """ add_chp_supplementary_firing_constraints(m, p; _n="") Used by add_chp_constraints to add supplementary firing constraints if - p.s.chp.supplementary_firing_max_steam_ratio > 1.0 to add CHP supplementary firing operating constraints. + p.chp_params[t][:supplementary_firing_max_steam_ratio] > 1.0 to add CHP supplementary firing operating constraints. Else, the supplementary firing dispatch and size decision variables are set to zero. """ function add_chp_supplementary_firing_constraints(m, p; _n="") - thermal_prod_full_load = 1.0 / p.s.chp.electric_efficiency_full_load * p.s.chp.thermal_efficiency_full_load # [kWt/kWe] - thermal_prod_half_load = 0.5 / p.s.chp.electric_efficiency_half_load * p.s.chp.thermal_efficiency_half_load # [kWt/kWe] - thermal_prod_slope = (thermal_prod_full_load - thermal_prod_half_load) / (1.0 - 0.5) # [kWt/kWe] - - # Constrain upper limit of dvSupplementaryThermalProduction, using auxiliary variable for (size * useSupplementaryFiring) - @constraint(m, CHPSupplementaryFireCon[t in p.techs.chp, ts in p.time_steps], - m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] <= - (p.s.chp.supplementary_firing_max_steam_ratio - 1.0) * p.production_factor[t,ts] * (thermal_prod_slope * m[Symbol("dvSupplementaryFiringSize"*_n)][t] + m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts]) - ) - if solver_is_compatible_with_indicator_constraints(p.s.settings.solver_name) - # Constrain lower limit of 0 if CHP tech is off - @constraint(m, NoCHPSupplementaryFireOffCon[t in p.techs.chp, ts in p.time_steps], - !m[Symbol("binCHPIsOnInTS"*_n)][t,ts] => {m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] <= 0.0} - ) - else - #There's no upper bound specified for the CHP supplementary firing, so assume the entire heat load as a reasonable maximum that wouldn't be exceeded (but might not be the best possible value). - max_supplementary_firing_size = maximum(p.s.dhw_load.loads_kw .+ p.s.space_heating_load.loads_kw) - @constraint(m, NoCHPSupplementaryFireOffCon[t in p.techs.chp, ts in p.time_steps], - m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] <= (p.s.chp.supplementary_firing_max_steam_ratio - 1.0) * p.production_factor[t,ts] * (thermal_prod_slope * max_supplementary_firing_size + m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts]) - ) + # Check if the Y-intercept variable exists + has_y_intercept = haskey(m, Symbol("dvHeatingProductionYIntercept"*_n)) + + for t in p.techs.chp + thermal_prod_full_load = 1.0 / p.chp_params[t][:electric_efficiency_full_load] * p.chp_params[t][:thermal_efficiency_full_load] # [kWt/kWe] + thermal_prod_half_load = 0.5 / p.chp_params[t][:electric_efficiency_half_load] * p.chp_params[t][:thermal_efficiency_half_load] # [kWt/kWe] + thermal_prod_slope = (thermal_prod_full_load - thermal_prod_half_load) / (1.0 - 0.5) # [kWt/kWe] + thermal_prod_intercept = thermal_prod_full_load - thermal_prod_slope * 1.0 # [kWt/kWe_rated] + + # Constrain upper limit of dvSupplementaryThermalProduction + if has_y_intercept && abs(thermal_prod_intercept) > 1.0E-7 + @constraint(m, [ts in p.time_steps], + m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] <= + (p.chp_params[t][:supplementary_firing_max_steam_ratio] - 1.0) * p.production_factor[t,ts] * (thermal_prod_slope * m[Symbol("dvSupplementaryFiringSize"*_n)][t] + m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts]) + ) + else + @constraint(m, [ts in p.time_steps], + m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] <= + (p.chp_params[t][:supplementary_firing_max_steam_ratio] - 1.0) * p.production_factor[t,ts] * thermal_prod_slope * m[Symbol("dvSupplementaryFiringSize"*_n)][t] + ) + end + + if solver_is_compatible_with_indicator_constraints(p.s.settings.solver_name) + # Constrain lower limit of 0 if CHP tech is off + @constraint(m, [ts in p.time_steps], + !m[Symbol("binCHPIsOnInTS"*_n)][t,ts] => {m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] <= 0.0} + ) + else + #There's no upper bound specified for the CHP supplementary firing, so assume the entire heat load as a reasonable maximum that wouldn't be exceeded (but might not be the best possible value). + max_supplementary_firing_size = maximum(p.s.dhw_load.loads_kw .+ p.s.space_heating_load.loads_kw) + if has_y_intercept && abs(thermal_prod_intercept) > 1.0E-7 + @constraint(m, [ts in p.time_steps], + m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] <= (p.chp_params[t][:supplementary_firing_max_steam_ratio] - 1.0) * p.production_factor[t,ts] * (thermal_prod_slope * max_supplementary_firing_size + m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts]) + ) + else + @constraint(m, [ts in p.time_steps], + m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] <= (p.chp_params[t][:supplementary_firing_max_steam_ratio] - 1.0) * p.production_factor[t,ts] * thermal_prod_slope * max_supplementary_firing_size + ) + end + end end end function add_binCHPIsOnInTS_constraints(m, p; _n="") # Note, min_turn_down_fraction for CHP is only enforced in p.time_steps_with_grid @constraint(m, [t in p.techs.chp, ts in p.time_steps_with_grid], - m[Symbol("dvRatedProduction"*_n)][t, ts] <= p.s.chp.max_kw * m[Symbol("binCHPIsOnInTS"*_n)][t, ts] + m[Symbol("dvRatedProduction"*_n)][t, ts] <= p.max_sizes[t] * m[Symbol("binCHPIsOnInTS"*_n)][t, ts] ) @constraint(m, [t in p.techs.chp, ts in p.time_steps_with_grid], - p.s.chp.min_turn_down_fraction * m[Symbol("dvSize"*_n)][t] - m[Symbol("dvRatedProduction"*_n)][t, ts] <= - p.s.chp.max_kw * (1 - m[Symbol("binCHPIsOnInTS"*_n)][t, ts]) + p.chp_params[t][:min_turn_down_fraction] * m[Symbol("dvSize"*_n)][t] - m[Symbol("dvRatedProduction"*_n)][t, ts] <= + p.max_sizes[t] * (1 - m[Symbol("binCHPIsOnInTS"*_n)][t, ts]) ) end @@ -145,22 +173,28 @@ function add_chp_hourly_om_charges(m, p; _n="") dv = "dvOMByHourBySizeCHP"*_n m[Symbol(dv)] = @variable(m, [p.techs.chp, p.time_steps], base_name=dv, lower_bound=0) - #Constraint CHP-hourly-om-a: om per hour, per time step >= per_unit_size_cost * size for when on, >= zero when off - @constraint(m, CHPHourlyOMBySizeA[t in p.techs.chp, ts in p.time_steps], - p.s.chp.om_cost_per_hr_per_kw_rated * m[Symbol("dvSize"*_n)][t] - - p.s.chp.max_kw * p.s.chp.om_cost_per_hr_per_kw_rated * (1-m[Symbol("binCHPIsOnInTS"*_n)][t,ts]) - <= m[Symbol("dvOMByHourBySizeCHP"*_n)][t, ts] - ) - #Constraint CHP-hourly-om-b: om per hour, per time step <= per_unit_size_cost * size for each hour - @constraint(m, CHPHourlyOMBySizeB[t in p.techs.chp, ts in p.time_steps], - p.s.chp.om_cost_per_hr_per_kw_rated * m[Symbol("dvSize"*_n)][t] - >= m[Symbol("dvOMByHourBySizeCHP"*_n)][t, ts] - ) - #Constraint CHP-hourly-om-c: om per hour, per time step <= zero when off, <= per_unit_size_cost*max_size - @constraint(m, CHPHourlyOMBySizeC[t in p.techs.chp, ts in p.time_steps], - p.s.chp.max_kw * p.s.chp.om_cost_per_hr_per_kw_rated * m[Symbol("binCHPIsOnInTS"*_n)][t,ts] - >= m[Symbol("dvOMByHourBySizeCHP"*_n)][t, ts] - ) + for t in p.techs.chp + # Find the CHP object for this tech to get om_cost_per_hr_per_kw_rated + chp_idx = findfirst(chp -> chp.name == t, p.s.chps) + om_cost_per_hr_per_kw = p.s.chps[chp_idx].om_cost_per_hr_per_kw_rated + + #Constraint CHP-hourly-om-a: om per hour, per time step >= per_unit_size_cost * size for when on, >= zero when off + @constraint(m, [ts in p.time_steps], + om_cost_per_hr_per_kw * m[Symbol("dvSize"*_n)][t] - + p.max_sizes[t] * om_cost_per_hr_per_kw * (1-m[Symbol("binCHPIsOnInTS"*_n)][t,ts]) + <= m[Symbol("dvOMByHourBySizeCHP"*_n)][t, ts] + ) + #Constraint CHP-hourly-om-b: om per hour, per time step <= per_unit_size_cost * size for each hour + @constraint(m, [ts in p.time_steps], + om_cost_per_hr_per_kw * m[Symbol("dvSize"*_n)][t] + >= m[Symbol("dvOMByHourBySizeCHP"*_n)][t, ts] + ) + #Constraint CHP-hourly-om-c: om per hour, per time step <= zero when off, <= per_unit_size_cost*max_size + @constraint(m, [ts in p.time_steps], + p.max_sizes[t] * om_cost_per_hr_per_kw * m[Symbol("binCHPIsOnInTS"*_n)][t,ts] + >= m[Symbol("dvOMByHourBySizeCHP"*_n)][t, ts] + ) + end m[:TotalHourlyCHPOMCosts] = @expression(m, p.third_party_factor * p.pwf_om * sum(m[Symbol(dv)][t, ts] * p.hours_per_time_step for t in p.techs.chp, ts in p.time_steps)) @@ -184,12 +218,14 @@ function add_chp_constraints(m, p; _n="") m[:TotalHourlyCHPOMCosts] = 0 m[:TotalCHPFuelCosts] = 0 + # Sum om_cost_per_kwh for each CHP tech m[:TotalCHPPerUnitProdOMCosts] = @expression(m, p.third_party_factor * p.pwf_om * - sum(p.s.chp.om_cost_per_kwh * p.hours_per_time_step * + sum(p.chp_params[t][:om_cost_per_kwh] * p.hours_per_time_step * m[:dvRatedProduction][t, ts] for t in p.techs.chp, ts in p.time_steps) ) - if p.s.chp.om_cost_per_hr_per_kw_rated > 1.0E-7 + # Check if any CHP has hourly O&M charges + if any(chp.om_cost_per_hr_per_kw_rated > 1.0E-7 for chp in p.s.chps) add_chp_hourly_om_charges(m, p) end @@ -198,10 +234,12 @@ function add_chp_constraints(m, p; _n="") add_binCHPIsOnInTS_constraints(m, p; _n=_n) add_chp_rated_prod_constraint(m, p; _n=_n) - if p.s.chp.supplementary_firing_max_steam_ratio > 1.0 - add_chp_supplementary_firing_constraints(m,p; _n=_n) - else - for t in p.techs.chp + # Add supplementary firing constraints - function handles per-tech logic + add_chp_supplementary_firing_constraints(m,p; _n=_n) + + # Fix supplementary firing variables to zero for CHPs without supplementary firing + for t in p.techs.chp + if p.chp_params[t][:supplementary_firing_max_steam_ratio] <= 1.0 fix(m[Symbol("dvSupplementaryFiringSize"*_n)][t], 0.0, force=true) for ts in p.time_steps fix(m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts], 0.0, force=true) diff --git a/src/constraints/cost_curve_constraints.jl b/src/constraints/cost_curve_constraints.jl index da1633aa1..bcb49d626 100644 --- a/src/constraints/cost_curve_constraints.jl +++ b/src/constraints/cost_curve_constraints.jl @@ -159,38 +159,44 @@ function initial_capex_no_incentives(m::JuMP.AbstractModel, p::REoptInputs; _n=" ) end - if "CHP" in p.techs.all + if !isempty(p.techs.chp) m[:CHPCapexNoIncentives] = JuMP.GenericAffExpr{Float64, JuMP.VariableRef}() - cost_list = p.s.chp.installed_cost_per_kw - size_list = p.s.chp.tech_sizes_for_cost_curve + + # Loop through each CHP and apply its specific cost curve + for chp in p.s.chps + t = chp.name + cost_list = chp.installed_cost_per_kw + size_list = chp.tech_sizes_for_cost_curve - t="CHP" - if t in p.techs.segmented && !isempty(size_list) - # Use "no incentives" version of p.cap_cost_slope and p.seg_yint - cost_slope_no_inc = [cost_list[1]] - seg_yint_no_inc = [0.0] - for s in range(2, stop=length(size_list)) - tmp_slope = round((cost_list[s] * size_list[s] - cost_list[s-1] * size_list[s-1]) / - (size_list[s] - size_list[s-1]), digits=0) - tmp_y_int = round(cost_list[s-1] * size_list[s-1] - tmp_slope * size_list[s-1], digits=0) - append!(cost_slope_no_inc, tmp_slope) - append!(seg_yint_no_inc, tmp_y_int) - end - append!(cost_slope_no_inc, cost_list[end]) - append!(seg_yint_no_inc, 0.0) + if t in p.techs.segmented && !isempty(size_list) + # Use "no incentives" version of p.cap_cost_slope and p.seg_yint + cost_slope_no_inc = [cost_list[1]] + seg_yint_no_inc = [0.0] + for s in range(2, stop=length(size_list)) + tmp_slope = round((cost_list[s] * size_list[s] - cost_list[s-1] * size_list[s-1]) / + (size_list[s] - size_list[s-1]), digits=0) + tmp_y_int = round(cost_list[s-1] * size_list[s-1] - tmp_slope * size_list[s-1], digits=0) + append!(cost_slope_no_inc, tmp_slope) + append!(seg_yint_no_inc, tmp_y_int) + end + append!(cost_slope_no_inc, cost_list[end]) + append!(seg_yint_no_inc, 0.0) - add_to_expression!(m[:CHPCapexNoIncentives], - sum(cost_slope_no_inc[s] * m[Symbol("dvSegmentSystemSize"*t)][s] + - seg_yint_no_inc[s] * m[Symbol("binSegment"*t)][s] for s in eachindex(cost_slope_no_inc)) - ) - else - add_to_expression!(m[:CHPCapexNoIncentives], cost_list * m[Symbol("dvPurchaseSize"*_n)]["CHP"]) - end - if p.s.chp.supplementary_firing_capital_cost_per_kw > 0 - add_to_expression!(m[:CHPCapexNoIncentives], - p.s.chp.supplementary_firing_capital_cost_per_kw * m[Symbol("dvSupplementaryFiringSize"*_n)]["CHP"] - ) + add_to_expression!(m[:CHPCapexNoIncentives], + sum(cost_slope_no_inc[s] * m[Symbol("dvSegmentSystemSize"*t)][s] + + seg_yint_no_inc[s] * m[Symbol("binSegment"*t)][s] for s in eachindex(cost_slope_no_inc)) + ) + else + add_to_expression!(m[:CHPCapexNoIncentives], cost_list * m[Symbol("dvPurchaseSize"*_n)][t]) + end + + if chp.supplementary_firing_capital_cost_per_kw > 0 + add_to_expression!(m[:CHPCapexNoIncentives], + chp.supplementary_firing_capital_cost_per_kw * m[Symbol("dvSupplementaryFiringSize"*_n)][t] + ) + end end + add_to_expression!(m[:InitialCapexNoIncentives], m[:CHPCapexNoIncentives]) end diff --git a/src/constraints/electric_utility_constraints.jl b/src/constraints/electric_utility_constraints.jl index 07dbbab76..27e4bd645 100644 --- a/src/constraints/electric_utility_constraints.jl +++ b/src/constraints/electric_utility_constraints.jl @@ -189,13 +189,16 @@ If the monthly demand rate is tiered than also adds binMonthlyDemandTier and con function add_monthly_peak_constraint(m, p; _n="") ## Constraint (11d): Monthly peak demand is >= demand at each hour in the month - if (!isempty(p.techs.chp)) && !(p.s.chp.reduces_demand_charges) + # Get CHPs that do NOT reduce demand charges + chps_not_reducing_demand = String[chp.name for chp in p.s.chps if !chp.reduces_demand_charges] + + if !isempty(chps_not_reducing_demand) @constraint(m, [mth in p.months, ts in p.s.electric_tariff.time_steps_monthly[mth]], sum(m[Symbol("dvPeakDemandMonth"*_n)][mth, t] for t in 1:p.s.electric_tariff.n_monthly_demand_tiers) >= sum(m[Symbol("dvGridPurchase"*_n)][ts, tier] for tier in 1:p.s.electric_tariff.n_energy_tiers) + - sum(p.production_factor[t, ts] * p.levelization_factor[t] * m[Symbol("dvRatedProduction"*_n)][t, ts] for t in p.techs.chp) - - sum(sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for b in p.s.storage.types.elec) for t in p.techs.chp) - - sum(sum(m[Symbol("dvProductionToGrid")][t,u,ts] for u in p.export_bins_by_tech[t]) for t in p.techs.chp) + sum(p.production_factor[t, ts] * p.levelization_factor[t] * m[Symbol("dvRatedProduction"*_n)][t, ts] for t in chps_not_reducing_demand) - + sum(sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for b in p.s.storage.types.elec) for t in chps_not_reducing_demand) - + sum(sum(m[Symbol("dvProductionToGrid")][t,u,ts] for u in p.export_bins_by_tech[t]) for t in chps_not_reducing_demand) ) else diff --git a/src/constraints/outage_constraints.jl b/src/constraints/outage_constraints.jl index d7fc82a44..f14a9be91 100644 --- a/src/constraints/outage_constraints.jl +++ b/src/constraints/outage_constraints.jl @@ -176,8 +176,8 @@ end function add_MG_CHP_fuel_burn_constraints(m, p; _n="") # Fuel burn slope and intercept fuel_burn_slope, fuel_burn_intercept = fuel_slope_and_intercept(; - electric_efficiency_full_load = p.s.chp.electric_efficiency_full_load, - electric_efficiency_half_load = p.s.chp.electric_efficiency_half_load, + electric_efficiency_full_load = p.s.chps[1].electric_efficiency_full_load, + electric_efficiency_half_load = p.s.chps[1].electric_efficiency_half_load, fuel_higher_heating_value_kwh_per_unit=1 ) From 4f4cc6ea4d64822609b2afdf561c8b1a49aaad03 Mon Sep 17 00:00:00 2001 From: wbecker Date: Mon, 1 Dec 2025 22:04:12 -0700 Subject: [PATCH 10/45] Multiple CHPs reopt.jl handling --- src/core/reopt.jl | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/src/core/reopt.jl b/src/core/reopt.jl index 6af9ef3c1..a955561a3 100644 --- a/src/core/reopt.jl +++ b/src/core/reopt.jl @@ -153,6 +153,9 @@ function run_reopt(ms::AbstractArray{T, 1}, p::REoptInputs) where T <: JuMP.Abst if !isempty(p.techs.pv) organize_multiple_pv_results(p, results_dict) end + if !isempty(p.techs.chp) + organize_multiple_chp_results(p, results_dict) + end return results_dict else throw(@error("REopt scenarios solved either with errors or non-optimal solutions.")) @@ -284,11 +287,17 @@ function build_reopt!(m::JuMP.AbstractModel, p::REoptInputs) m[:TotalFuelCosts] += m[:TotalCHPFuelCosts] m[:TotalPerUnitHourOMCosts] += m[:TotalHourlyCHPOMCosts] - if p.s.chp.standby_rate_per_kw_per_month > 1.0e-7 - m[:TotalCHPStandbyCharges] += sum(p.pwf_e * 12 * p.s.chp.standby_rate_per_kw_per_month * m[:dvSize][t] for t in p.techs.chp) + # Add standby charges for each CHP + for chp in p.s.chps + if chp.standby_rate_per_kw_per_month > 1.0e-7 + m[:TotalCHPStandbyCharges] += p.pwf_e * 12 * chp.standby_rate_per_kw_per_month * m[:dvSize][chp.name] + end end - m[:TotalTechCapCosts] += sum(p.s.chp.supplementary_firing_capital_cost_per_kw * m[:dvSupplementaryFiringSize][t] for t in p.techs.chp) + # Add supplementary firing capital costs for each CHP + for chp in p.s.chps + m[:TotalTechCapCosts] += chp.supplementary_firing_capital_cost_per_kw * m[:dvSupplementaryFiringSize][chp.name] + end end if !isempty(setdiff(p.techs.heating, p.techs.elec)) @@ -620,6 +629,9 @@ function run_reopt(m::JuMP.AbstractModel, p::REoptInputs; organize_pvs=true) if organize_pvs && !isempty(p.techs.pv) # do not want to organize_pvs when running BAU case in parallel b/c then proform code fails organize_multiple_pv_results(p, results) end + if organize_pvs && !isempty(p.techs.chp) # same logic as PV + organize_multiple_chp_results(p, results) + end # add error messages (if any) and warnings to results dict results["Messages"] = logger_to_dict() From 3da54c780f5c1c4547294383bc08f7a556f1288d Mon Sep 17 00:00:00 2001 From: wbecker Date: Mon, 1 Dec 2025 22:05:14 -0700 Subject: [PATCH 11/45] Multiple CHPs results and outagesim --- src/outagesim/backup_reliability.jl | 4 +- src/results/chp.jl | 161 ++++++++++++++++------------ src/results/proforma.jl | 10 +- src/results/results.jl | 2 +- 4 files changed, 101 insertions(+), 76 deletions(-) diff --git a/src/outagesim/backup_reliability.jl b/src/outagesim/backup_reliability.jl index 9692da169..6841f5ca7 100644 --- a/src/outagesim/backup_reliability.jl +++ b/src/outagesim/backup_reliability.jl @@ -798,8 +798,8 @@ function backup_reliability_reopt_inputs(;d::Dict, p::REoptInputs, r::Dict = Dic end if prime_kw > 0 fuel_slope, fuel_intercept = fuel_slope_and_intercept( - electric_efficiency_full_load=p.s.chp.electric_efficiency_full_load, - electric_efficiency_half_load=p.s.chp.electric_efficiency_half_load, + electric_efficiency_full_load=p.s.chps[1].electric_efficiency_full_load, + electric_efficiency_half_load=p.s.chps[1].electric_efficiency_half_load, fuel_higher_heating_value_kwh_per_unit=1 ) r2[:generator_fuel_burn_rate_per_kwh] = [fuel_slope] diff --git a/src/results/chp.jl b/src/results/chp.jl index 58f302804..5f3aede16 100644 --- a/src/results/chp.jl +++ b/src/results/chp.jl @@ -31,108 +31,107 @@ function add_chp_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="") # Adds the `CHP` results to the dictionary passed back from `run_reopt` using the solved model `m` and the `REoptInputs` for node `_n`. # Note: the node number is an empty string if evaluating a single `Site`. + + # Add each CHP's results using its name as the key + for chp_name in p.techs.chp + r = get_chp_results_for_tech(m, p, chp_name, _n) + d[chp_name] = r + end + nothing +end + +function get_chp_results_for_tech(m::JuMP.AbstractModel, p::REoptInputs, chp_name::String, _n::String) r = Dict{String, Any}() - r["size_kw"] = value(sum(m[Symbol("dvSize"*_n)][t] for t in p.techs.chp)) - r["size_supplemental_firing_kw"] = value(sum(m[Symbol("dvSupplementaryFiringSize"*_n)][t] for t in p.techs.chp)) - @expression(m, CHPFuelUsedKWH, sum(m[Symbol("dvFuelUsage"*_n)][t, ts] for t in p.techs.chp, ts in p.time_steps)) - r["annual_fuel_consumption_mmbtu"] = round(value(CHPFuelUsedKWH) / KWH_PER_MMBTU, digits=3) - @expression(m, Year1CHPElecProd, - p.hours_per_time_step * sum(m[Symbol("dvRatedProduction"*_n)][t,ts] * p.production_factor[t, ts] - for t in p.techs.chp, ts in p.time_steps)) - r["annual_electric_production_kwh"] = round(value(Year1CHPElecProd), digits=3) + + # Find the CHP object for this name + chp_idx = findfirst(chp -> chp.name == chp_name, p.s.chps) + if isnothing(chp_idx) + @warn "CHP named $chp_name not found in scenario" + return r + end + chp = p.s.chps[chp_idx] + + r["size_kw"] = value(m[Symbol("dvSize"*_n)][chp_name]) + r["size_supplemental_firing_kw"] = value(m[Symbol("dvSupplementaryFiringSize"*_n)][chp_name]) + CHPFuelUsedKWH = sum(value(m[Symbol("dvFuelUsage"*_n)][chp_name, ts]) for ts in p.time_steps) + r["annual_fuel_consumption_mmbtu"] = round(CHPFuelUsedKWH / KWH_PER_MMBTU, digits=3) + Year1CHPElecProd = p.hours_per_time_step * sum(value(m[Symbol("dvRatedProduction"*_n)][chp_name,ts]) * p.production_factor[chp_name, ts] + for ts in p.time_steps) + r["annual_electric_production_kwh"] = round(Year1CHPElecProd, digits=3) - @expression(m, CHPThermalProdKW[ts in p.time_steps], - sum(sum(m[Symbol("dvHeatingProduction"*_n)][t,q,ts] - m[Symbol("dvProductionToWaste"*_n)][t,q,ts] for q in p.heating_loads) + - m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] for t in p.techs.chp)) + CHPThermalProdKW = [sum(value(m[Symbol("dvHeatingProduction"*_n)][chp_name,q,ts]) - value(m[Symbol("dvProductionToWaste"*_n)][chp_name,q,ts]) for q in p.heating_loads) + + value(m[Symbol("dvSupplementaryThermalProduction"*_n)][chp_name,ts]) for ts in p.time_steps] - r["thermal_production_series_mmbtu_per_hour"] = round.(value.(CHPThermalProdKW) / KWH_PER_MMBTU, digits=5) + r["thermal_production_series_mmbtu_per_hour"] = round.(CHPThermalProdKW / KWH_PER_MMBTU, digits=5) r["annual_thermal_production_mmbtu"] = round(p.hours_per_time_step * sum(r["thermal_production_series_mmbtu_per_hour"]), digits=3) - @expression(m, CHPElecProdTotal[ts in p.time_steps], - sum(m[Symbol("dvRatedProduction"*_n)][t,ts] * p.production_factor[t, ts] for t in p.techs.chp)) - r["electric_production_series_kw"] = round.(value.(CHPElecProdTotal), digits=3) + CHPElecProdTotal = [value(m[Symbol("dvRatedProduction"*_n)][chp_name,ts]) * p.production_factor[chp_name, ts] for ts in p.time_steps] + r["electric_production_series_kw"] = round.(CHPElecProdTotal, digits=3) # Electric dispatch breakdown if !isempty(p.s.electric_tariff.export_bins) - @expression(m, CHPtoGrid[ts in p.time_steps], sum(m[Symbol("dvProductionToGrid"*_n)][t,u,ts] - for t in p.techs.chp, u in p.export_bins_by_tech[t])) + CHPtoGrid = [sum(value(m[Symbol("dvProductionToGrid"*_n)][chp_name,u,ts]) + for u in p.export_bins_by_tech[chp_name]) for ts in p.time_steps] else CHPtoGrid = zeros(length(p.time_steps)) end - r["electric_to_grid_series_kw"] = round.(value.(CHPtoGrid), digits=3) + r["electric_to_grid_series_kw"] = round.(CHPtoGrid, digits=3) if !isempty(p.s.storage.types.elec) - @expression(m, CHPtoBatt[ts in p.time_steps], - sum(m[Symbol("dvProductionToStorage"*_n)]["ElectricStorage",t,ts] for t in p.techs.chp)) + CHPtoBatt = [value(m[Symbol("dvProductionToStorage"*_n)]["ElectricStorage",chp_name,ts]) for ts in p.time_steps] else CHPtoBatt = zeros(length(p.time_steps)) end - r["electric_to_storage_series_kw"] = round.(value.(CHPtoBatt), digits=3) - @expression(m, CHPtoLoad[ts in p.time_steps], - sum(m[Symbol("dvRatedProduction"*_n)][t, ts] * p.production_factor[t, ts] * p.levelization_factor[t] - for t in p.techs.chp) - CHPtoBatt[ts] - CHPtoGrid[ts]) - r["electric_to_load_series_kw"] = round.(value.(CHPtoLoad), digits=3) + r["electric_to_storage_series_kw"] = round.(CHPtoBatt, digits=3) + CHPtoLoad = [value(m[Symbol("dvRatedProduction"*_n)][chp_name, ts]) * p.production_factor[chp_name, ts] * p.levelization_factor[chp_name] - CHPtoBatt[ts] - CHPtoGrid[ts] for ts in p.time_steps] + r["electric_to_load_series_kw"] = round.(CHPtoLoad, digits=3) # Thermal dispatch breakdown if !isempty(p.s.storage.types.hot) - @expression(m, CHPToHotTES[ts in p.time_steps], - sum(m[Symbol("dvHeatToStorage"*_n)][b, t, q, ts] for b in p.s.storage.types.hot, t in p.techs.chp, q in p.heating_loads)) - @expression(m, CHPToHotTESByQuality[q in p.heating_loads, ts in p.time_steps], - sum(m[Symbol("dvHeatToStorage"*_n)][b, t, q, ts] for b in p.s.storage.types.hot, t in p.techs.chp)) + CHPToHotTES = [sum(value(m[Symbol("dvHeatToStorage"*_n)][b, chp_name, q, ts]) for b in p.s.storage.types.hot, q in p.heating_loads) for ts in p.time_steps] + CHPToHotTESByQuality = Dict(q => [sum(value(m[Symbol("dvHeatToStorage"*_n)][b, chp_name, q, ts]) for b in p.s.storage.types.hot) for ts in p.time_steps] for q in p.heating_loads) else - @expression(m, CHPToHotTES[ts in p.time_steps], 0.0) - @expression(m, CHPToHotTESByQuality[q in p.heating_loads, ts in p.time_steps], 0.0) + CHPToHotTES = zeros(length(p.time_steps)) + CHPToHotTESByQuality = Dict(q => zeros(length(p.time_steps)) for q in p.heating_loads) end - r["thermal_to_storage_series_mmbtu_per_hour"] = round.(value.(CHPToHotTES / KWH_PER_MMBTU), digits=5) - @expression(m, CHPThermalToWasteKW[ts in p.time_steps], - sum(m[Symbol("dvProductionToWaste"*_n)][t,q,ts] for q in p.heating_loads, t in p.techs.chp)) - @expression(m, CHPThermalToWasteByQualityKW[q in p.heating_loads, ts in p.time_steps], - sum(m[Symbol("dvProductionToWaste"*_n)][t,q,ts] for t in p.techs.chp)) - r["thermal_curtailed_series_mmbtu_per_hour"] = round.(value.(CHPThermalToWasteKW) / KWH_PER_MMBTU, digits=5) - if !isempty(p.techs.steam_turbine) && p.s.chp.can_supply_steam_turbine - @expression(m, CHPToSteamTurbineKW[ts in p.time_steps], sum(m[Symbol("dvThermalToSteamTurbine"*_n)][t,q,ts] for t in p.techs.chp, q in p.heating_loads)) - @expression(m, CHPToSteamTurbineByQualityKW[q in p.heating_loads, ts in p.time_steps], sum(m[Symbol("dvThermalToSteamTurbine"*_n)][t,q,ts] for t in p.techs.chp)) + r["thermal_to_storage_series_mmbtu_per_hour"] = round.(CHPToHotTES / KWH_PER_MMBTU, digits=5) + CHPThermalToWasteKW = [sum(value(m[Symbol("dvProductionToWaste"*_n)][chp_name,q,ts]) for q in p.heating_loads) for ts in p.time_steps] + CHPThermalToWasteByQualityKW = Dict(q => [value(m[Symbol("dvProductionToWaste"*_n)][chp_name,q,ts]) for ts in p.time_steps] for q in p.heating_loads) + r["thermal_curtailed_series_mmbtu_per_hour"] = round.(CHPThermalToWasteKW / KWH_PER_MMBTU, digits=5) + if !isempty(p.techs.steam_turbine) && chp.can_supply_steam_turbine + CHPToSteamTurbineKW = [sum(value(m[Symbol("dvThermalToSteamTurbine"*_n)][chp_name,q,ts]) for q in p.heating_loads) for ts in p.time_steps] + CHPToSteamTurbineByQualityKW = Dict(q => [value(m[Symbol("dvThermalToSteamTurbine"*_n)][chp_name,q,ts]) for ts in p.time_steps] for q in p.heating_loads) else CHPToSteamTurbineKW = zeros(length(p.time_steps)) - @expression(m, CHPToSteamTurbineByQualityKW[q in p.heating_loads, ts in p.time_steps], 0.0) + CHPToSteamTurbineByQualityKW = Dict(q => zeros(length(p.time_steps)) for q in p.heating_loads) end - r["thermal_to_steamturbine_series_mmbtu_per_hour"] = round.(value.(CHPToSteamTurbineKW) / KWH_PER_MMBTU, digits=5) - @expression(m, CHPThermalToLoadKW[ts in p.time_steps], - 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) + r["thermal_to_steamturbine_series_mmbtu_per_hour"] = round.(CHPToSteamTurbineKW / KWH_PER_MMBTU, digits=5) + CHPThermalToLoadKW = [sum(value(m[Symbol("dvHeatingProduction"*_n)][chp_name,q,ts]) for q in p.heating_loads) + value(m[Symbol("dvSupplementaryThermalProduction"*_n)][chp_name,ts]) - CHPToHotTES[ts] - CHPToSteamTurbineKW[ts] - CHPThermalToWasteKW[ts] for ts in p.time_steps] + r["thermal_to_load_series_mmbtu_per_hour"] = round.(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], - m[:dvHeatingProduction]["CHP","DomesticHotWater",ts] - CHPToHotTESByQuality["DomesticHotWater",ts] - CHPToSteamTurbineByQualityKW["DomesticHotWater",ts] - CHPThermalToWasteByQualityKW["DomesticHotWater",ts] - ) + if "DomesticHotWater" in p.heating_loads && chp.can_serve_dhw + CHPToDHWKW = [value(m[:dvHeatingProduction][chp_name,"DomesticHotWater",ts]) - CHPToHotTESByQuality["DomesticHotWater"][ts] - CHPToSteamTurbineByQualityKW["DomesticHotWater"][ts] - CHPThermalToWasteByQualityKW["DomesticHotWater"][ts] + for ts in p.time_steps] else - @expression(m, CHPToDHWKW[ts in p.time_steps], 0.0) + CHPToDHWKW = zeros(length(p.time_steps)) end - r["thermal_to_dhw_load_series_mmbtu_per_hour"] = round.(value.(CHPToDHWKW ./ KWH_PER_MMBTU), digits=5) + r["thermal_to_dhw_load_series_mmbtu_per_hour"] = round.(CHPToDHWKW ./ KWH_PER_MMBTU, digits=5) - if "SpaceHeating" in p.heating_loads && p.s.chp.can_serve_space_heating - @expression(m, CHPToSpaceHeatingKW[ts in p.time_steps], - m[:dvHeatingProduction]["CHP","SpaceHeating",ts] - CHPToHotTESByQuality["SpaceHeating",ts] - CHPToSteamTurbineByQualityKW["SpaceHeating",ts] - CHPThermalToWasteByQualityKW["SpaceHeating",ts] - ) + if "SpaceHeating" in p.heating_loads && chp.can_serve_space_heating + CHPToSpaceHeatingKW = [value(m[:dvHeatingProduction][chp_name,"SpaceHeating",ts]) - CHPToHotTESByQuality["SpaceHeating"][ts] - CHPToSteamTurbineByQualityKW["SpaceHeating"][ts] - CHPThermalToWasteByQualityKW["SpaceHeating"][ts] + for ts in p.time_steps] else - @expression(m, CHPToSpaceHeatingKW[ts in p.time_steps], 0.0) + CHPToSpaceHeatingKW = zeros(length(p.time_steps)) end - r["thermal_to_space_heating_load_series_mmbtu_per_hour"] = round.(value.(CHPToSpaceHeatingKW ./ KWH_PER_MMBTU), digits=5) + r["thermal_to_space_heating_load_series_mmbtu_per_hour"] = round.(CHPToSpaceHeatingKW ./ KWH_PER_MMBTU, digits=5) - if "ProcessHeat" in p.heating_loads && p.s.chp.can_serve_process_heat - @expression(m, CHPToProcessHeatKW[ts in p.time_steps], - m[:dvHeatingProduction]["CHP","ProcessHeat",ts] - CHPToHotTESByQuality["ProcessHeat",ts] - CHPToSteamTurbineByQualityKW["ProcessHeat",ts] - CHPThermalToWasteByQualityKW["ProcessHeat",ts] - ) + if "ProcessHeat" in p.heating_loads && chp.can_serve_process_heat + CHPToProcessHeatKW = [value(m[:dvHeatingProduction][chp_name,"ProcessHeat",ts]) - CHPToHotTESByQuality["ProcessHeat"][ts] - CHPToSteamTurbineByQualityKW["ProcessHeat"][ts] - CHPThermalToWasteByQualityKW["ProcessHeat"][ts] + for ts in p.time_steps] else - @expression(m, CHPToProcessHeatKW[ts in p.time_steps], 0.0) + CHPToProcessHeatKW = zeros(length(p.time_steps)) end - r["thermal_to_process_heat_load_series_mmbtu_per_hour"] = round.(value.(CHPToProcessHeatKW ./ KWH_PER_MMBTU), digits=5) + r["thermal_to_process_heat_load_series_mmbtu_per_hour"] = round.(CHPToProcessHeatKW ./ KWH_PER_MMBTU, digits=5) - r["year_one_fuel_cost_before_tax"] = round(value(m[:TotalCHPFuelCosts] / p.pwf_fuel["CHP"]), digits=3) + r["year_one_fuel_cost_before_tax"] = round(value(m[:TotalCHPFuelCosts] / p.pwf_fuel[chp_name]), digits=3) r["year_one_fuel_cost_after_tax"] = r["year_one_fuel_cost_before_tax"] * (1 - p.s.financial.offtaker_tax_rate_fraction) r["lifecycle_fuel_cost_after_tax"] = round(value(m[:TotalCHPFuelCosts]) * (1- p.s.financial.offtaker_tax_rate_fraction), digits=3) #Standby charges and hourly O&M @@ -141,6 +140,26 @@ function add_chp_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="") r["lifecycle_standby_cost_after_tax"] = round(value(m[Symbol("TotalCHPStandbyCharges")]) * (1 - p.s.financial.offtaker_tax_rate_fraction), digits=0) r["initial_capital_costs"] = round(value(m[Symbol("CHPCapexNoIncentives")]), digits=2) - d["CHP"] = r + return r +end + + +""" + organize_multiple_chp_results(p::REoptInputs, d::Dict) + +The last step in results processing: if more than one CHP was modeled then move their results from the top +level keys (that use each CHP.name) to an array of results with "CHP" as the top key in the results dict `d`. +""" +function organize_multiple_chp_results(p::REoptInputs, d::Dict) + if length(p.techs.chp) == 1 && p.techs.chp[1] == "CHP" + return nothing + end + chps = Dict[] + for chpname in p.techs.chp + d[chpname]["name"] = chpname # add name to results dict to distinguish each CHP + push!(chps, d[chpname]) + delete!(d, chpname) + end + d["CHP"] = chps nothing end diff --git a/src/results/proforma.jl b/src/results/proforma.jl index a6cfc007f..8fb1f40fa 100644 --- a/src/results/proforma.jl +++ b/src/results/proforma.jl @@ -133,8 +133,14 @@ function proforma_results(p::REoptInputs, d::Dict) end # calculate CHP o+m costs, incentives, and depreciation - if "CHP" in keys(d) && d["CHP"]["size_kw"] > 0 - update_metrics(m, p, p.s.chp, "CHP", d, third_party) + for chp_name in p.techs.chp + if chp_name in keys(d) && get(d[chp_name], "size_kw", 0) > 0 + # Find the CHP object for this name + chp_idx = findfirst(chp -> chp.name == chp_name, p.s.chps) + if !isnothing(chp_idx) + update_metrics(m, p, p.s.chps[chp_idx], chp_name, d, third_party) + end + end end # calculate ExistingBoiler o+m costs (just fuel, no non-fuel operating costs currently) diff --git a/src/results/results.jl b/src/results/results.jl index 4847830bc..ffe634106 100644 --- a/src/results/results.jl +++ b/src/results/results.jl @@ -44,7 +44,7 @@ function reopt_results(m::JuMP.AbstractModel, p::REoptInputs; _n="") add_wind_results(m, p, d; _n) end - if "CHP" in p.techs.all + if !isempty(p.techs.chp) add_chp_results(m, p, d; _n) end From cf725e79b23be2261a56345003aac9cd1afd098b Mon Sep 17 00:00:00 2001 From: wbecker Date: Mon, 1 Dec 2025 22:06:28 -0700 Subject: [PATCH 12/45] Update references to s.chp(s) in runtests.jl --- test/runtests.jl | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 1da1eeb2d..081c97a95 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -37,7 +37,7 @@ else # run HiGHS tests s = Scenario(input_data) @test s.financial.owner_tax_rate_fraction == 0.26 @test s.financial.elec_cost_escalation_rate_fraction == 0.0166 - for tech_struct in (s.pvs[1], s.wind, s.chp, s.steam_turbine) + for tech_struct in (s.pvs[1], s.wind, s.chps[1], s.steam_turbine) for incentive_input_name in (:macrs_option_years, :macrs_bonus_fraction) @test getfield(tech_struct, incentive_input_name) != 0 end @@ -58,7 +58,7 @@ else # run HiGHS tests s = Scenario(input_data) @test s.financial.owner_tax_rate_fraction == 0.26 @test s.financial.chp_fuel_cost_escalation_rate_fraction == 0.00581 #national avg - for tech_struct in (s.pvs[1], s.wind, s.chp, s.steam_turbine) + for tech_struct in (s.pvs[1], s.wind, s.chps[1], s.steam_turbine) for incentive_input_name in (:macrs_option_years, :macrs_bonus_fraction) @test getfield(tech_struct, incentive_input_name) != 0 end @@ -79,7 +79,7 @@ else # run HiGHS tests s = Scenario(input_data) @test s.financial.owner_tax_rate_fraction == 0.0 @test s.financial.elec_cost_escalation_rate_fraction == -0.00088 - for tech_struct in (s.pvs[1], s.wind, s.chp, s.ghp_option_list[1], s.steam_turbine) + for tech_struct in (s.pvs[1], s.wind, s.chps[1], s.ghp_option_list[1], s.steam_turbine) for incentive_input_name in (:macrs_option_years, :macrs_bonus_fraction, :federal_itc_fraction) default = 0 try @@ -1990,8 +1990,8 @@ else # run HiGHS tests @test round(total_chiller_electric_consumption, digits=0) ≈ 320544.0 atol=1.0 # loads_kw is **electric**, loads_kw_thermal is **thermal** #Test CHP defaults use average fuel load, size class 2 for recip_engine - @test inputs.s.chp.min_allowable_kw ≈ 50.0 atol=0.01 - @test inputs.s.chp.om_cost_per_kwh ≈ 0.0235 atol=0.0001 + @test inputs.s.chps[1].min_allowable_kw ≈ 50.0 atol=0.01 + @test inputs.s.chps[1].om_cost_per_kwh ≈ 0.0235 atol=0.0001 delete!(input_data, "SpaceHeatingLoad") delete!(input_data, "DomesticHotWaterLoad") @@ -2007,7 +2007,7 @@ else # run HiGHS tests @test round(total_chiller_electric_consumption, digits=0) ≈ 3876410 atol=1.0 # Check that without heating load or max_kw input, CHP.max_kw gets set based on peak electric load - @test inputs.s.chp.max_kw ≈ maximum(inputs.s.electric_load.loads_kw) atol=0.01 + @test inputs.s.chps[1].max_kw ≈ maximum(inputs.s.electric_load.loads_kw) atol=0.01 input_data["SpaceHeatingLoad"] = Dict{Any, Any}("monthly_mmbtu" => repeat([1000.0], 12)) input_data["DomesticHotWaterLoad"] = Dict{Any, Any}("monthly_mmbtu" => repeat([1000.0], 12)) @@ -2017,8 +2017,8 @@ else # run HiGHS tests inputs = REoptInputs(s) #Test CHP defaults use average fuel load, size class changes to 3 - @test inputs.s.chp.min_allowable_kw ≈ 125.0 atol=0.1 - @test inputs.s.chp.om_cost_per_kwh ≈ 0.021 atol=0.0001 + @test inputs.s.chps[1].min_allowable_kw ≈ 125.0 atol=0.1 + @test inputs.s.chps[1].om_cost_per_kwh ≈ 0.021 atol=0.0001 #Update CHP prime_mover and test new defaults input_data["CHP"]["prime_mover"] = "combustion_turbine" input_data["CHP"]["size_class"] = 1 @@ -2028,8 +2028,8 @@ else # run HiGHS tests s = Scenario(input_data) inputs = REoptInputs(s) - @test inputs.s.chp.min_allowable_kw ≈ 2000.0 atol=0.1 - @test inputs.s.chp.om_cost_per_kwh ≈ 0.014499999999999999 atol=0.0001 + @test inputs.s.chps[1].min_allowable_kw ≈ 2000.0 atol=0.1 + @test inputs.s.chps[1].om_cost_per_kwh ≈ 0.014499999999999999 atol=0.0001 total_heating_fuel_load_mmbtu = (sum(inputs.s.space_heating_load.loads_kw) + sum(inputs.s.dhw_load.loads_kw)) / input_data["ExistingBoiler"]["efficiency"] / REopt.KWH_PER_MMBTU @@ -2078,14 +2078,14 @@ else # run HiGHS tests s = Scenario(input_data) inputs = REoptInputs(s) # Costs are 75% of CHP - @test inputs.s.chp.installed_cost_per_kw ≈ (0.75*installed_cost_chp) atol=1.0 - @test inputs.s.chp.om_cost_per_kwh ≈ (0.75*0.0145) atol=0.0001 - @test inputs.s.chp.federal_itc_fraction ≈ 0.0 atol=0.0001 + @test inputs.s.chps[1].installed_cost_per_kw ≈ (0.75*installed_cost_chp) atol=1.0 + @test inputs.s.chps[1].om_cost_per_kwh ≈ (0.75*0.0145) atol=0.0001 + @test inputs.s.chps[1].federal_itc_fraction ≈ 0.0 atol=0.0001 # Thermal efficiency set to zero - @test inputs.s.chp.thermal_efficiency_full_load == 0 - @test inputs.s.chp.thermal_efficiency_half_load == 0 + @test inputs.s.chps[1].thermal_efficiency_full_load == 0 + @test inputs.s.chps[1].thermal_efficiency_half_load == 0 # Max size based on electric load, not heating load - @test inputs.s.chp.max_kw ≈ maximum(inputs.s.electric_load.loads_kw) atol=0.001 + @test inputs.s.chps[1].max_kw ≈ maximum(inputs.s.electric_load.loads_kw) atol=0.001 end @testset "Hybrid/blended heating and cooling loads" begin From c8d60b49ce9b8eabd00f6af1b9baf33fd46a8083 Mon Sep 17 00:00:00 2001 From: wbecker Date: Tue, 2 Dec 2025 14:42:00 -0700 Subject: [PATCH 13/45] Add chps to MPCScenario struct --- src/mpc/scenario.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/mpc/scenario.jl b/src/mpc/scenario.jl index db6c6d096..8ad33426f 100644 --- a/src/mpc/scenario.jl +++ b/src/mpc/scenario.jl @@ -11,6 +11,7 @@ struct MPCScenario <: AbstractScenario cooling_load::MPCCoolingLoad limits::MPCLimits node::Int + chps::Array{CHP, 1} # Empty array for MPC scenarios (no CHP modeled) end @@ -31,6 +32,7 @@ Method for creating the MPCScenario struct: cooling_load::MPCCoolingLoad limits::MPCLimits node::Int + chps::Array{CHP, 1} end ``` @@ -125,6 +127,7 @@ function MPCScenario(d::Dict) generator, cooling_load, limits, - node + node, + CHP[] # Empty array - no CHP in MPC scenarios ) end From c8b5ba1f5d1dee22d19ff5c4fb10a34b0dba07bd Mon Sep 17 00:00:00 2001 From: wbecker Date: Tue, 2 Dec 2025 14:43:03 -0700 Subject: [PATCH 14/45] Handle empty export bins for CHP --- src/results/chp.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/results/chp.jl b/src/results/chp.jl index 5f3aede16..0da0976e7 100644 --- a/src/results/chp.jl +++ b/src/results/chp.jl @@ -69,7 +69,7 @@ function get_chp_results_for_tech(m::JuMP.AbstractModel, p::REoptInputs, chp_nam CHPElecProdTotal = [value(m[Symbol("dvRatedProduction"*_n)][chp_name,ts]) * p.production_factor[chp_name, ts] for ts in p.time_steps] r["electric_production_series_kw"] = round.(CHPElecProdTotal, digits=3) # Electric dispatch breakdown - if !isempty(p.s.electric_tariff.export_bins) + if !isempty(p.s.electric_tariff.export_bins) && haskey(p.export_bins_by_tech, chp_name) && !isempty(p.export_bins_by_tech[chp_name]) CHPtoGrid = [sum(value(m[Symbol("dvProductionToGrid"*_n)][chp_name,u,ts]) for u in p.export_bins_by_tech[chp_name]) for ts in p.time_steps] else From 4ad72f84b622ac2cf8b5230c1629b2e7933b1e00 Mon Sep 17 00:00:00 2001 From: wbecker Date: Tue, 2 Dec 2025 14:43:33 -0700 Subject: [PATCH 15/45] Add tolerance on payback test --- test/runtests.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 081c97a95..ed53950cf 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3915,7 +3915,7 @@ else # run HiGHS tests capital_costs_after_non_discounted_incentives = results["Financial"]["capital_costs_after_non_discounted_incentives"] # Calculated payback from above-two metrics payback = capital_costs_after_non_discounted_incentives / savings - @test round(results["Financial"]["simple_payback_years"], digits=2) ≈ round(payback, digits=2) + @test round(results["Financial"]["simple_payback_years"], digits=2) ≈ round(payback, digits=2) rtol=0.01 finalize(backend(m1)) empty!(m1) finalize(backend(m2)) @@ -3932,7 +3932,7 @@ else # run HiGHS tests m2 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) results = run_reopt([m1,m2], inputs) payback = results["Financial"]["capital_costs_after_non_discounted_incentives"] / results["Financial"]["year_one_total_operating_cost_savings_after_tax"] - @test round(results["Financial"]["simple_payback_years"], digits=2) ≈ round(payback, digits=2) + @test round(results["Financial"]["simple_payback_years"], digits=2) ≈ round(payback, digits=2) rtol=0.01 finalize(backend(m1)) empty!(m1) finalize(backend(m2)) From 5e1a5660ecde67d300f2162c3ce7dd665c7f4978 Mon Sep 17 00:00:00 2001 From: wbecker Date: Tue, 2 Dec 2025 15:01:44 -0700 Subject: [PATCH 16/45] Updates with multiple CHPs in test_with_xpress (obsolete anyway) --- test/test_with_xpress.jl | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/test/test_with_xpress.jl b/test/test_with_xpress.jl index 47cfcb312..d0adff595 100644 --- a/test/test_with_xpress.jl +++ b/test/test_with_xpress.jl @@ -260,6 +260,7 @@ end empty!(m2) GC.gc() end +end @testset "FlexibleHVAC" begin @@ -956,8 +957,8 @@ end @test round(total_chiller_electric_consumption, digits=0) ≈ 320544.0 atol=1.0 # loads_kw is **electric**, loads_kw_thermal is **thermal** #Test CHP defaults use average fuel load, size class 2 for recip_engine - @test inputs.s.chp.min_allowable_kw ≈ 50.0 atol=0.01 - @test inputs.s.chp.om_cost_per_kwh ≈ 0.0235 atol=0.0001 + @test inputs.s.chps[1].min_allowable_kw ≈ 50.0 atol=0.01 + @test inputs.s.chps[1].om_cost_per_kwh ≈ 0.0235 atol=0.0001 delete!(input_data, "SpaceHeatingLoad") delete!(input_data, "DomesticHotWaterLoad") @@ -973,7 +974,7 @@ end @test round(total_chiller_electric_consumption, digits=0) ≈ 3876410 atol=1.0 # Check that without heating load or max_kw input, CHP.max_kw gets set based on peak electric load - @test inputs.s.chp.max_kw ≈ maximum(inputs.s.electric_load.loads_kw) atol=0.01 + @test inputs.s.chps[1].max_kw ≈ maximum(inputs.s.electric_load.loads_kw) atol=0.01 input_data["SpaceHeatingLoad"] = Dict{Any, Any}("monthly_mmbtu" => repeat([1000.0], 12)) input_data["DomesticHotWaterLoad"] = Dict{Any, Any}("monthly_mmbtu" => repeat([1000.0], 12)) @@ -983,8 +984,8 @@ end inputs = REoptInputs(s) #Test CHP defaults use average fuel load, size class changes to 3 - @test inputs.s.chp.min_allowable_kw ≈ 125.0 atol=0.1 - @test inputs.s.chp.om_cost_per_kwh ≈ 0.021 atol=0.0001 + @test inputs.s.chps[1].min_allowable_kw ≈ 125.0 atol=0.1 + @test inputs.s.chps[1].om_cost_per_kwh ≈ 0.021 atol=0.0001 #Update CHP prime_mover and test new defaults input_data["CHP"]["prime_mover"] = "combustion_turbine" input_data["CHP"]["size_class"] = 1 @@ -994,8 +995,8 @@ end s = Scenario(input_data) inputs = REoptInputs(s) - @test inputs.s.chp.min_allowable_kw ≈ 2000.0 atol=0.1 - @test inputs.s.chp.om_cost_per_kwh ≈ 0.014499999999999999 atol=0.0001 + @test inputs.s.chps[1].min_allowable_kw ≈ 2000.0 atol=0.1 + @test inputs.s.chps[1].om_cost_per_kwh ≈ 0.014499999999999999 atol=0.0001 total_heating_fuel_load_mmbtu = (sum(inputs.s.space_heating_load.loads_kw) + sum(inputs.s.dhw_load.loads_kw)) / input_data["ExistingBoiler"]["efficiency"] / REopt.KWH_PER_MMBTU @@ -1044,14 +1045,14 @@ end s = Scenario(input_data) inputs = REoptInputs(s) # Costs are 75% of CHP - @test inputs.s.chp.installed_cost_per_kw ≈ (0.75*installed_cost_chp) atol=1.0 - @test inputs.s.chp.om_cost_per_kwh ≈ (0.75*0.0145) atol=0.0001 - @test inputs.s.chp.federal_itc_fraction ≈ 0.0 atol=0.0001 + @test inputs.s.chps[1].installed_cost_per_kw ≈ (0.75*installed_cost_chp) atol=1.0 + @test inputs.s.chps[1].om_cost_per_kwh ≈ (0.75*0.0145) atol=0.0001 + @test inputs.s.chps[1].federal_itc_fraction ≈ 0.0 atol=0.0001 # Thermal efficiency set to zero - @test inputs.s.chp.thermal_efficiency_full_load == 0 - @test inputs.s.chp.thermal_efficiency_half_load == 0 + @test inputs.s.chps[1].thermal_efficiency_full_load == 0 + @test inputs.s.chps[1].thermal_efficiency_half_load == 0 # Max size based on electric load, not heating load - @test inputs.s.chp.max_kw ≈ maximum(inputs.s.electric_load.loads_kw) atol=0.001 + @test inputs.s.chps[1].max_kw ≈ maximum(inputs.s.electric_load.loads_kw) atol=0.001 end @testset "Hybrid/blended heating and cooling loads" begin From edaed8e2d2e9b323679154ae543f03910ad00e9d Mon Sep 17 00:00:00 2001 From: wbecker Date: Tue, 2 Dec 2025 15:08:43 -0700 Subject: [PATCH 17/45] Add temp dev test files for multiple CHPs --- test/scenarios/multiple_chps.json | 55 +++++++++++++++++++++++++++++++ test/test_multiple_chps.jl | 34 +++++++++++++++++++ 2 files changed, 89 insertions(+) create mode 100644 test/scenarios/multiple_chps.json create mode 100644 test/test_multiple_chps.jl diff --git a/test/scenarios/multiple_chps.json b/test/scenarios/multiple_chps.json new file mode 100644 index 000000000..6d07508b7 --- /dev/null +++ b/test/scenarios/multiple_chps.json @@ -0,0 +1,55 @@ +{ + "Site": { + "latitude": 34.5794343, + "longitude": -118.1164613 + }, + "Financial": { + "om_cost_escalation_rate_fraction": 0.025, + "elec_cost_escalation_rate_fraction": 0.023, + "offtaker_tax_rate_fraction": 0.26, + "offtaker_discount_rate_fraction": 0.083, + "analysis_years": 25 + }, + "ElectricLoad": { + "doe_reference_name": "Hospital", + "annual_kwh": 10000000.0, + "year": 2017 + }, + "ElectricTariff": { + "blended_annual_energy_rate": 0.12, + "blended_annual_demand_rate": 15.0 + }, + "DomesticHotWaterLoad": { + "doe_reference_name": "Hospital", + "annual_mmbtu": 500000.0 + }, + "SpaceHeatingLoad": { + "doe_reference_name": "Hospital", + "annual_mmbtu": 1000000.0 + }, + "ExistingBoiler": { + "fuel_cost_per_mmbtu": 12.0 + }, + "CHP": [ + { + "name": "CHP_recip_engine", + "prime_mover": "recip_engine", + "size_class": 1, + "min_kw": 100.0, + "max_kw": 500.0, + "fuel_cost_per_mmbtu": 8.0, + "can_serve_dhw": true, + "can_serve_space_heating": true + }, + { + "name": "CHP_micro_turbine", + "prime_mover": "micro_turbine", + "size_class": 2, + "min_kw": 50.0, + "max_kw": 300.0, + "fuel_cost_per_mmbtu": 8.0, + "can_serve_dhw": true, + "can_serve_space_heating": false + } + ] +} diff --git a/test/test_multiple_chps.jl b/test/test_multiple_chps.jl new file mode 100644 index 000000000..cb6fe2977 --- /dev/null +++ b/test/test_multiple_chps.jl @@ -0,0 +1,34 @@ +using Revise +using REopt +using JSON +using DelimitedFiles +using PlotlyJS +using Dates +using Test +using JuMP +using HiGHS +using DotEnv +DotEnv.load!() + +############### Multiple CHPs Test ################### +input_data = JSON.parsefile("./scenarios/multiple_chps.json") +s = Scenario(input_data) +inputs = REoptInputs(s) + +m1 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +results = run_reopt(m1, inputs) +# m2 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +# results = run_reopt([m1,m2], inputs) + +# Check that both CHPs are in results +@test length(results["CHP"]) == 2 +CHP_recip_engine = results["CHP"][findfirst(chp -> chp["name"] == "CHP_recip_engine", results["CHP"])] +CHP_micro_turbine = results["CHP"][findfirst(chp -> chp["name"] == "CHP_micro_turbine", results["CHP"])] + +# Check that each CHP has sizing results +@test CHP_recip_engine["size_kw"] > 10.0 +@test CHP_micro_turbine["size_kw"] > 10.0 + +# Check that results include electric production for each CHP +@test sum(CHP_recip_engine["electric_production_series_kw"]) > 5000.0 +@test sum(CHP_micro_turbine["electric_production_series_kw"]) > 5000.0 \ No newline at end of file From 3543e377339ba89ac3da0ad7a6b985183342ea33 Mon Sep 17 00:00:00 2001 From: wbecker Date: Tue, 2 Dec 2025 16:39:26 -0700 Subject: [PATCH 18/45] Fix test for chps reference --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index ed53950cf..7179d1d2d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2070,7 +2070,7 @@ else # run HiGHS tests "blended_annual_demand_rate" => 0.0 ) s_chp = Scenario(input_data) inputs_chp = REoptInputs(s) - installed_cost_chp = s_chp.chp.installed_cost_per_kw + installed_cost_chp = s_chp.chps[1].installed_cost_per_kw # Now get prime generator (electric only) input_data["CHP"]["is_electric_only"] = true From da3322a901ea0c1867a954bcc5c83d0f17344a0c Mon Sep 17 00:00:00 2001 From: wbecker Date: Sun, 14 Dec 2025 19:46:32 -0700 Subject: [PATCH 19/45] Fix Boiler results for non-hourly timesteps and unit conversions --- src/results/boiler.jl | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/results/boiler.jl b/src/results/boiler.jl index 3b0beaee0..bb5ee7c83 100644 --- a/src/results/boiler.jl +++ b/src/results/boiler.jl @@ -24,13 +24,12 @@ function add_boiler_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n=" r["size_mmbtu_per_hour"] = round(value(m[Symbol("dvSize"*_n)]["Boiler"]) / KWH_PER_MMBTU, digits=3) r["fuel_consumption_series_mmbtu_per_hour"] = round.(value.(m[:dvFuelUsage]["Boiler", ts] for ts in p.time_steps) / KWH_PER_MMBTU, digits=3) - r["annual_fuel_consumption_mmbtu"] = round(sum(r["fuel_consumption_series_mmbtu_per_hour"]), digits=3) + r["annual_fuel_consumption_mmbtu"] = round(sum(r["fuel_consumption_series_mmbtu_per_hour"]) / p.s.settings.time_steps_per_hour, digits=3) r["thermal_production_series_mmbtu_per_hour"] = round.(sum(value.(m[:dvHeatingProduction]["Boiler", q, ts] for ts in p.time_steps) for q in p.heating_loads) ./ KWH_PER_MMBTU, digits=5) - r["annual_thermal_production_mmbtu"] = round(sum(r["thermal_production_series_mmbtu_per_hour"]), digits=3) - - if !isempty(p.s.storage.types.hot) + r["annual_thermal_production_mmbtu"] = round(sum(r["thermal_production_series_mmbtu_per_hour"]) / p.s.settings.time_steps_per_hour, digits=3) + if !isempty(p.s.storage.types.hot) @expression(m, NewBoilerToHotTESKW[ts in p.time_steps], sum(m[:dvHeatToStorage][b,"Boiler",q,ts] for b in p.s.storage.types.hot, q in p.heating_loads) ) @@ -48,7 +47,7 @@ function add_boiler_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n=" NewBoilerToSteamTurbine = zeros(length(p.time_steps)) @expression(m, NewBoilerToSteamTurbineByQuality[q in p.heating_loads, ts in p.time_steps], 0.0) end - r["thermal_to_steamturbine_series_mmbtu_per_hour"] = round.(value.(NewBoilerToSteamTurbine), digits=3) + r["thermal_to_steamturbine_series_mmbtu_per_hour"] = round.(value.(NewBoilerToSteamTurbine) ./ KWH_PER_MMBTU, digits=3) BoilerToLoad = @expression(m, [ts in p.time_steps], sum(value.(m[:dvHeatingProduction]["Boiler", q, ts]) for q in p.heating_loads) - NewBoilerToHotTESKW[ts] - NewBoilerToSteamTurbine[ts] From c8d92668f6ec90408a58d800d8256ec2aba8c7dd Mon Sep 17 00:00:00 2001 From: wbecker Date: Sun, 14 Dec 2025 19:49:46 -0700 Subject: [PATCH 20/45] Add temp dev Boiler + SteamTurbine test files --- test/scenarios/boiler_steamturbine.json | 51 +++++++++++++++++++++++++ test/test_boiler_steamturbine.jl | 47 +++++++++++++++++++++++ 2 files changed, 98 insertions(+) create mode 100644 test/scenarios/boiler_steamturbine.json create mode 100644 test/test_boiler_steamturbine.jl diff --git a/test/scenarios/boiler_steamturbine.json b/test/scenarios/boiler_steamturbine.json new file mode 100644 index 000000000..9893ae6b3 --- /dev/null +++ b/test/scenarios/boiler_steamturbine.json @@ -0,0 +1,51 @@ +{ + "Site": { + "latitude": 37.78, + "longitude": -122.45 + }, + "ElectricLoad": { + "doe_reference_name": "Hospital", + "annual_kwh": 5000000.0 + }, + "SpaceHeatingLoad": { + "doe_reference_name": "Hospital", + "annual_mmbtu": 10000.0 + }, + "DomesticHotWaterLoad": { + "doe_reference_name": "Hospital", + "annual_mmbtu": 2000.0 + }, + "ElectricTariff": { + "blended_annual_energy_rate": 0.10, + "blended_annual_demand_rate": 10.0 + }, + "ExistingBoiler":{ + "fuel_cost_per_mmbtu": 8.0, + "efficiency": 0.75, + "can_supply_steam_turbine": false + }, + "Boiler": { + "min_mmbtu_per_hour": 0.0, + "max_mmbtu_per_hour": 100.0, + "efficiency": 0.8, + "fuel_type": "natural_gas", + "fuel_cost_per_mmbtu": 8.0, + "installed_cost_per_mmbtu_per_hour": 50000.0, + "om_cost_per_mmbtu_per_hour": 500.0, + "om_cost_per_mmbtu": 0.5, + "can_supply_steam_turbine": true + }, + "SteamTurbine": { + "min_kw": 0.0, + "max_kw": 5000.0, + "electric_produced_to_thermal_consumed_ratio": 0.30, + "thermal_produced_to_thermal_consumed_ratio": 0.50, + "installed_cost_per_kw": 1000.0, + "om_cost_per_kw": 10.0, + "om_cost_per_kwh": 0.01, + "can_wholesale": false, + "can_curtail": false, + "macrs_bonus_fraction": 0.0, + "macrs_option_years": 5 + } +} diff --git a/test/test_boiler_steamturbine.jl b/test/test_boiler_steamturbine.jl new file mode 100644 index 000000000..c09b73b63 --- /dev/null +++ b/test/test_boiler_steamturbine.jl @@ -0,0 +1,47 @@ +using Revise +using REopt +using JSON +using DelimitedFiles +using PlotlyJS +using Dates +using Test +using JuMP +using HiGHS +using DotEnv +DotEnv.load!() + +############### Boiler + SteamTurbine with simplified parameters ################### +# This test uses the simpler SteamTurbine setup with electric_produced_to_thermal_consumed_ratio +# and thermal_produced_to_thermal_consumed_ratio instead of detailed steam parameters + +input_data = JSON.parsefile("./scenarios/boiler_steamturbine.json") +s = Scenario(input_data) +inputs = REoptInputs(s) + +m1 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.001, "output_flag" => false, "log_to_console" => false)) +m2 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.001, "output_flag" => false, "log_to_console" => false)) +results = run_reopt([m1,m2], inputs) + +println("\n=== Boiler + SteamTurbine Test Results ===") +println("Financial NPV: ", round(results["Financial"]["npv"], digits=2)) +println("Financial LCC: ", round(results["Financial"]["lcc"], digits=2)) +println("\nBoiler:") +println(" Size (MMBtu/hr): ", round(results["Boiler"]["size_mmbtu_per_hour"], digits=2)) +println(" Annual thermal production (MMBtu): ", round(results["Boiler"]["annual_thermal_production_mmbtu"], digits=2)) +println(" Annual fuel consumption (MMBtu): ", round(results["Boiler"]["annual_fuel_consumption_mmbtu"], digits=2)) +println(" Thermal to SteamTurbine (MMBtu): ", round(sum(results["Boiler"]["thermal_to_steamturbine_series_mmbtu_per_hour"]), digits=2)) +println("\nSteamTurbine:") +println(" Size (kW): ", round(results["SteamTurbine"]["size_kw"], digits=2)) +println(" Annual electric production (kWh): ", round(results["SteamTurbine"]["annual_electric_production_kwh"], digits=2)) +println(" Annual thermal consumption (MMBtu): ", round(results["SteamTurbine"]["annual_thermal_consumption_mmbtu"], digits=2)) +println(" Annual thermal production (MMBtu): ", round(results["SteamTurbine"]["annual_thermal_production_mmbtu"], digits=2)) + +# Verify energy balance: Boiler thermal to ST should match ST thermal consumption +boiler_to_st = sum(results["Boiler"]["thermal_to_steamturbine_series_mmbtu_per_hour"]) +st_thermal_in = results["SteamTurbine"]["annual_thermal_consumption_mmbtu"] +println("\nEnergy Balance Check:") +println(" Boiler->SteamTurbine: ", round(boiler_to_st, digits=2), " MMBtu") +println(" SteamTurbine thermal in: ", round(st_thermal_in, digits=2), " MMBtu") +println(" Difference: ", round(abs(boiler_to_st - st_thermal_in), digits=2), " MMBtu") + +println("=============================================\n") From a52f4f17836db4b701b8f855e1c8f28070810a81 Mon Sep 17 00:00:00 2001 From: wbecker Date: Mon, 15 Dec 2025 15:18:01 -0700 Subject: [PATCH 21/45] Avoid binCHPIsOnInTS binary dV array if not needed --- src/constraints/chp_constraints.jl | 94 +++++++++++++++++------------- 1 file changed, 54 insertions(+), 40 deletions(-) diff --git a/src/constraints/chp_constraints.jl b/src/constraints/chp_constraints.jl index 0f25bf54f..67af76526 100644 --- a/src/constraints/chp_constraints.jl +++ b/src/constraints/chp_constraints.jl @@ -1,17 +1,11 @@ # REopt®, Copyright (c) Alliance for Sustainable Energy, LLC. See also https://github.com/NREL/REopt.jl/blob/master/LICENSE. -function add_chp_fuel_burn_constraints(m, p; _n="") - # Fuel burn slope and intercept - fuel_burn_slope, fuel_burn_intercept = fuel_slope_and_intercept(; - electric_efficiency_full_load = p.s.chp.electric_efficiency_full_load, - electric_efficiency_half_load = p.s.chp.electric_efficiency_half_load, - fuel_higher_heating_value_kwh_per_unit=1 - ) - +function add_chp_fuel_burn_constraints(m, p; _n="", binary_created=true, fuel_burn_slope, fuel_burn_intercept) # Fuel cost m[:TotalCHPFuelCosts] = @expression(m, sum(p.pwf_fuel[t] * m[:dvFuelUsage][t, ts] * p.fuel_cost_per_kwh[t][ts] for t in p.techs.chp, ts in p.time_steps) ) # Conditionally add dvFuelBurnYIntercept if coefficient p.FuelBurnYIntRate is greater than ~zero + # Note: if intercept exists, binary_created will always be true (checked in add_chp_constraints) if abs(fuel_burn_intercept) > 1.0E-7 dv = "dvFuelBurnYIntercept"*_n m[Symbol(dv)] = @variable(m, [p.techs.chp, p.time_steps], base_name=dv) @@ -40,15 +34,9 @@ function add_chp_fuel_burn_constraints(m, p; _n="") end end -function add_chp_thermal_production_constraints(m, p; _n="") - # Thermal production slope and intercept - thermal_prod_full_load = 1.0 / p.s.chp.electric_efficiency_full_load * p.s.chp.thermal_efficiency_full_load # [kWt/kWe] - thermal_prod_half_load = 0.5 / p.s.chp.electric_efficiency_half_load * p.s.chp.thermal_efficiency_half_load # [kWt/kWe] - thermal_prod_slope = (thermal_prod_full_load - thermal_prod_half_load) / (1.0 - 0.5) # [kWt/kWe] - thermal_prod_intercept = thermal_prod_full_load - thermal_prod_slope * 1.0 # [kWt/kWe_rated - - +function add_chp_thermal_production_constraints(m, p; _n="", binary_created=true, thermal_prod_slope, thermal_prod_intercept) # Conditionally add dvHeatingProductionYIntercept if coefficient p.s.chpThermalProdIntercept is greater than ~zero + # Note: if intercept exists, binary_created will always be true (checked in add_chp_constraints) if abs(thermal_prod_intercept) > 1.0E-7 dv = "dvHeatingProductionYIntercept"*_n m[Symbol(dv)] = @variable(m, [p.techs.chp, p.time_steps], base_name=dv) @@ -92,16 +80,20 @@ Used by add_chp_constraints to add supplementary firing constraints if p.s.chp.supplementary_firing_max_steam_ratio > 1.0 to add CHP supplementary firing operating constraints. Else, the supplementary firing dispatch and size decision variables are set to zero. """ -function add_chp_supplementary_firing_constraints(m, p; _n="") - thermal_prod_full_load = 1.0 / p.s.chp.electric_efficiency_full_load * p.s.chp.thermal_efficiency_full_load # [kWt/kWe] - thermal_prod_half_load = 0.5 / p.s.chp.electric_efficiency_half_load * p.s.chp.thermal_efficiency_half_load # [kWt/kWe] - thermal_prod_slope = (thermal_prod_full_load - thermal_prod_half_load) / (1.0 - 0.5) # [kWt/kWe] - +function add_chp_supplementary_firing_constraints(m, p; _n="", thermal_prod_slope, thermal_prod_intercept) # Constrain upper limit of dvSupplementaryThermalProduction, using auxiliary variable for (size * useSupplementaryFiring) - @constraint(m, CHPSupplementaryFireCon[t in p.techs.chp, ts in p.time_steps], - m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] <= - (p.s.chp.supplementary_firing_max_steam_ratio - 1.0) * p.production_factor[t,ts] * (thermal_prod_slope * m[Symbol("dvSupplementaryFiringSize"*_n)][t] + m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts]) - ) + # Note: dvHeatingProductionYIntercept only exists when thermal_prod_intercept > 0 + if abs(thermal_prod_intercept) > 1.0E-7 + @constraint(m, CHPSupplementaryFireCon[t in p.techs.chp, ts in p.time_steps], + m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] <= + (p.s.chp.supplementary_firing_max_steam_ratio - 1.0) * p.production_factor[t,ts] * (thermal_prod_slope * m[Symbol("dvSupplementaryFiringSize"*_n)][t] + m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts]) + ) + else + @constraint(m, CHPSupplementaryFireCon[t in p.techs.chp, ts in p.time_steps], + m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] <= + (p.s.chp.supplementary_firing_max_steam_ratio - 1.0) * p.production_factor[t,ts] * thermal_prod_slope * m[Symbol("dvSupplementaryFiringSize"*_n)][t] + ) + end if solver_is_compatible_with_indicator_constraints(p.s.settings.solver_name) # Constrain lower limit of 0 if CHP tech is off @constraint(m, NoCHPSupplementaryFireOffCon[t in p.techs.chp, ts in p.time_steps], @@ -174,13 +166,32 @@ end Used in src/reopt.jl to add_chp_constraints if !isempty(p.techs.chp) to add CHP operating constraints and cost expressions. """ -function add_chp_constraints(m, p; _n="") - # TODO if chp.min_turn_down_fraction is 0.0, and there is no fuel burn or thermal y-intercept, we don't need the binary below - @warn """Adding binary variable to model CHP. - Some solvers are very slow with integer variables""" - @variables m begin - binCHPIsOnInTS[p.techs.chp, p.time_steps], Bin # 1 If technology t is operating in time step; 0 otherwise - end +function add_chp_constraints(m, p; _n="") + # Calculate fuel burn and thermal production slopes and intercepts which are passed into some constraint functions + fuel_burn_slope, fuel_burn_intercept = fuel_slope_and_intercept(; + electric_efficiency_full_load = p.s.chp.electric_efficiency_full_load, + electric_efficiency_half_load = p.s.chp.electric_efficiency_half_load, + fuel_higher_heating_value_kwh_per_unit=1 + ) + + thermal_prod_full_load = 1.0 / p.s.chp.electric_efficiency_full_load * p.s.chp.thermal_efficiency_full_load + thermal_prod_half_load = 0.5 / p.s.chp.electric_efficiency_half_load * p.s.chp.thermal_efficiency_half_load + thermal_prod_slope = (thermal_prod_full_load - thermal_prod_half_load) / (1.0 - 0.5) + thermal_prod_intercept = thermal_prod_full_load - thermal_prod_slope * 1.0 + + # Check if binary variable is needed which creates binCHPIsOnInTS binary decision variable and activates other constraints + binary_needed = (abs(fuel_burn_intercept) > 1.0E-7) || + (abs(thermal_prod_intercept) > 1.0E-7) || + (p.s.chp.min_turn_down_fraction > 1.0E-7) || + (p.s.chp.om_cost_per_hr_per_kw_rated > 1.0E-7) + + if binary_needed + @warn """Adding binary variable binCHPIsOnInTS to model CHP. + Some solvers are very slow with integer variables""" + @variables m begin + binCHPIsOnInTS[p.techs.chp, p.time_steps], Bin # 1 If technology t is operating in time step; 0 otherwise + end + end m[:TotalHourlyCHPOMCosts] = 0 m[:TotalCHPFuelCosts] = 0 @@ -189,17 +200,20 @@ function add_chp_constraints(m, p; _n="") m[:dvRatedProduction][t, ts] for t in p.techs.chp, ts in p.time_steps) ) - if p.s.chp.om_cost_per_hr_per_kw_rated > 1.0E-7 - add_chp_hourly_om_charges(m, p) - end - - add_chp_fuel_burn_constraints(m, p; _n=_n) - add_chp_thermal_production_constraints(m, p; _n=_n) - add_binCHPIsOnInTS_constraints(m, p; _n=_n) + # These constraints are always needed + add_chp_fuel_burn_constraints(m, p; _n=_n, binary_created=binary_needed, fuel_burn_slope=fuel_burn_slope, fuel_burn_intercept=fuel_burn_intercept) + add_chp_thermal_production_constraints(m, p; _n=_n, binary_created=binary_needed, thermal_prod_slope=thermal_prod_slope, thermal_prod_intercept=thermal_prod_intercept) add_chp_rated_prod_constraint(m, p; _n=_n) + # These constraints and decision variables created within the functions are only needed if binary_needed or more specific parameters is non-zero + if binary_needed + add_binCHPIsOnInTS_constraints(m, p; _n=_n) + end + if p.s.chp.om_cost_per_hr_per_kw_rated > 1.0E-7 + add_chp_hourly_om_charges(m, p; _n=_n) + end if p.s.chp.supplementary_firing_max_steam_ratio > 1.0 - add_chp_supplementary_firing_constraints(m,p; _n=_n) + add_chp_supplementary_firing_constraints(m, p; _n=_n, thermal_prod_slope=thermal_prod_slope, thermal_prod_intercept=thermal_prod_intercept) else for t in p.techs.chp fix(m[Symbol("dvSupplementaryFiringSize"*_n)][t], 0.0, force=true) From b86fcf94eaf3d35ad49259aed4cf689d8d8bc1fb Mon Sep 17 00:00:00 2001 From: wbecker Date: Mon, 15 Dec 2025 15:35:34 -0700 Subject: [PATCH 22/45] Add temp dev test file for avoiding CHP binary --- test/test_avoid_chp_binary.jl | 205 ++++++++++++++++++++++++++++++++++ 1 file changed, 205 insertions(+) create mode 100644 test/test_avoid_chp_binary.jl diff --git a/test/test_avoid_chp_binary.jl b/test/test_avoid_chp_binary.jl new file mode 100644 index 000000000..199683cd8 --- /dev/null +++ b/test/test_avoid_chp_binary.jl @@ -0,0 +1,205 @@ + +using Revise +using REopt +using JSON +using DelimitedFiles +using PlotlyJS +using Dates +using Test +using JuMP +using HiGHS +using DotEnv +DotEnv.load!() + +############### CHP Binary Creation Tests ################### +# Tests to verify that binCHPIsOnInTS is only created when needed + +@testset "CHP Binary Creation Tests" begin + + @testset "CHP with min_turn_down_fraction > 0 (binary SHOULD be created)" begin + # Test that binCHPIsOnInTS gets created when min_turn_down_fraction > 0 + data = JSON.parsefile("./scenarios/chp_unavailability_outage.json") + + # Set min_turn_down_fraction to a non-zero value + data["CHP"]["min_turn_down_fraction"] = 0.5 + + # Remove outage to simplify the test + delete!(data, "ElectricUtility") + data["ElectricUtility"] = Dict("net_metering_limit_kw" => 0.0, "co2_from_avert" => true) + + s = Scenario(data) + inputs = REoptInputs(s) + m = Model(optimizer_with_attributes(HiGHS.Optimizer, + "output_flag" => false, + "log_to_console" => false, + "mip_rel_gap" => 0.01)) + results = run_reopt(m, inputs) + + # Verify the binary variable was created + @test haskey(m.obj_dict, :binCHPIsOnInTS) + binary_var = m[:binCHPIsOnInTS] + + # Verify model solved successfully + @test results["CHP"]["size_kw"] > 0.0 + + println("✓ Test 1 passed: Binary created with min_turn_down_fraction = 0.5") + finalize(backend(m)) + empty!(m) + end + + @testset "CHP with min_turn_down_fraction = 0 (binary should NOT be created)" begin + # Test that binCHPIsOnInTS is NOT created when all conditions are zero/equal + data = JSON.parsefile("./scenarios/chp_unavailability_outage.json") + + # Set min_turn_down_fraction to zero + data["CHP"]["min_turn_down_fraction"] = 0.0 + + # Ensure no intercepts by making efficiencies equal (no part-load curve) + # When electric_efficiency_full_load == electric_efficiency_half_load, intercept = 0 + data["CHP"]["electric_efficiency_full_load"] = 0.35 + data["CHP"]["electric_efficiency_half_load"] = 0.35 + data["CHP"]["thermal_efficiency_full_load"] = 0.45 + data["CHP"]["thermal_efficiency_half_load"] = 0.45 + + # Ensure no hourly O&M cost + data["CHP"]["om_cost_per_hr_per_kw_rated"] = 0.0 + + # Remove outage to simplify + delete!(data, "ElectricUtility") + data["ElectricUtility"] = Dict("net_metering_limit_kw" => 0.0, "co2_from_avert" => true) + + s = Scenario(data) + inputs = REoptInputs(s) + m = Model(optimizer_with_attributes(HiGHS.Optimizer, + "output_flag" => false, + "log_to_console" => false, + "mip_rel_gap" => 0.01)) + results = run_reopt(m, inputs) + + # Verify the binary variable was NOT created + @test !haskey(m.obj_dict, :binCHPIsOnInTS) + + # Verify model still solved successfully + @test results["CHP"]["size_kw"] > 0.0 + + println("✓ Test 2 passed: Binary NOT created with min_turn_down_fraction = 0 and no intercepts") + finalize(backend(m)) + empty!(m) + end + + @testset "CHP with fuel_burn_intercept > 0 (binary SHOULD be created)" begin + # Test that binary is created when fuel burn has non-zero intercept + data = JSON.parsefile("./scenarios/chp_unavailability_outage.json") + + # Set min_turn_down_fraction to zero + data["CHP"]["min_turn_down_fraction"] = 0.0 + + # Create non-zero fuel burn intercept by having different efficiencies + data["CHP"]["electric_efficiency_full_load"] = 0.40 + data["CHP"]["electric_efficiency_half_load"] = 0.30 # Different from full load + data["CHP"]["thermal_efficiency_full_load"] = 0.45 + data["CHP"]["thermal_efficiency_half_load"] = 0.45 # Keep thermal same to isolate fuel burn + + # Ensure no hourly O&M cost + data["CHP"]["om_cost_per_hr_per_kw_rated"] = 0.0 + + delete!(data, "ElectricUtility") + data["ElectricUtility"] = Dict("net_metering_limit_kw" => 0.0, "co2_from_avert" => true) + + s = Scenario(data) + inputs = REoptInputs(s) + m = Model(optimizer_with_attributes(HiGHS.Optimizer, + "output_flag" => false, + "log_to_console" => false, + "mip_rel_gap" => 0.01)) + results = run_reopt(m, inputs) + + # Verify the binary variable was created + @test haskey(m.obj_dict, :binCHPIsOnInTS) + + # verify model solved successfully + @test results["Financial"]["lcc"] != 0.0 + + println("✓ Test 3 passed: Binary created with non-zero fuel_burn_intercept") + finalize(backend(m)) + empty!(m) + end + + @testset "CHP with thermal_prod_intercept > 0 (binary SHOULD be created)" begin + # Test that binary is created when thermal production has non-zero intercept + data = JSON.parsefile("./scenarios/chp_unavailability_outage.json") + + # Set min_turn_down_fraction to zero + data["CHP"]["min_turn_down_fraction"] = 0.0 + + # Create non-zero thermal prod intercept by having different thermal efficiencies + data["CHP"]["electric_efficiency_full_load"] = 0.35 + data["CHP"]["electric_efficiency_half_load"] = 0.35 # Keep electric same + data["CHP"]["thermal_efficiency_full_load"] = 0.40 + data["CHP"]["thermal_efficiency_half_load"] = 0.50 # Different from full load + + # Ensure no hourly O&M cost + data["CHP"]["om_cost_per_hr_per_kw_rated"] = 0.0 + + delete!(data, "ElectricUtility") + data["ElectricUtility"] = Dict("net_metering_limit_kw" => 0.0, "co2_from_avert" => true) + + s = Scenario(data) + inputs = REoptInputs(s) + m = Model(optimizer_with_attributes(HiGHS.Optimizer, + "output_flag" => false, + "log_to_console" => false, + "mip_rel_gap" => 0.01)) + results = run_reopt(m, inputs) + + # Verify the binary variable was created + @test haskey(m.obj_dict, :binCHPIsOnInTS) + + # Verify model solved successfully + @test results["Financial"]["lcc"] != 0.0 + + println("✓ Test 4 passed: Binary created with non-zero thermal_prod_intercept") + finalize(backend(m)) + empty!(m) + end + + @testset "CHP with om_cost_per_hr_per_kw_rated > 0 (binary SHOULD be created)" begin + # Test that binary is created when hourly O&M cost is non-zero + data = JSON.parsefile("./scenarios/chp_unavailability_outage.json") + + # Set min_turn_down_fraction to zero + data["CHP"]["min_turn_down_fraction"] = 0.0 + + # Make efficiencies equal (no intercepts) + data["CHP"]["electric_efficiency_full_load"] = 0.35 + data["CHP"]["electric_efficiency_half_load"] = 0.35 + data["CHP"]["thermal_efficiency_full_load"] = 0.45 + data["CHP"]["thermal_efficiency_half_load"] = 0.45 + + # Set non-zero hourly O&M cost + data["CHP"]["om_cost_per_hr_per_kw_rated"] = 0.01 + + delete!(data, "ElectricUtility") + data["ElectricUtility"] = Dict("net_metering_limit_kw" => 0.0, "co2_from_avert" => true) + + s = Scenario(data) + inputs = REoptInputs(s) + m = Model(optimizer_with_attributes(HiGHS.Optimizer, + "output_flag" => false, + "log_to_console" => false, + "mip_rel_gap" => 0.01)) + results = run_reopt(m, inputs) + + # Verify the binary variable was created + @test haskey(m.obj_dict, :binCHPIsOnInTS) + + # Verify model solved successfully + @test results["Financial"]["lcc"] != 0.0 + + println("✓ Test 5 passed: Binary created with om_cost_per_hr_per_kw_rated = 0.01") + finalize(backend(m)) + empty!(m) + end + + println("\n===== All CHP Binary Creation Tests Passed! =====\n") +end \ No newline at end of file From 290efa2f1ab29bd4d33af3e9b8be6c9b5a7b4994 Mon Sep 17 00:00:00 2001 From: wbecker Date: Mon, 15 Dec 2025 23:02:41 -0700 Subject: [PATCH 23/45] Add optional production_factor_series input for CHP to limit electric production for each timestep Could represent run-of-river hydro or OA temp-dependent turbine power --- src/core/chp.jl | 2 ++ src/core/production_factor.jl | 6 ++++++ 2 files changed, 8 insertions(+) diff --git a/src/core/chp.jl b/src/core/chp.jl index c59514a9f..77a979371 100644 --- a/src/core/chp.jl +++ b/src/core/chp.jl @@ -22,6 +22,7 @@ conflict_res_min_allowable_fraction_of_max = 0.25 cooling_thermal_factor::Float64 = NaN # only needed with cooling load unavailability_periods::AbstractVector{Dict} = Dict[] # CHP unavailability periods for scheduled and unscheduled maintenance, list of dictionaries with keys of "['month', 'start_week_of_month', 'start_day_of_week', 'start_hour', 'duration_hours'] all values are one-indexed and start_day_of_week uses 1 for Monday, 7 for Sunday unavailability_hourly::AbstractVector{Int64} = Int64[] # Hourly 8760 profile of unavailability (1) and availability (0) + production_factor_series::Union{Nothing, Vector{<:Real}} = nothing # Optional user-provided production factor time-series (length of time_steps_per_hour * 8760). If provided, this will override the production factor calculated from unavailability. # Optional inputs: size_class::Union{Int, Nothing} = nothing # CHP size class for using appropriate default inputs, with size_class=0 using an average of all other size class data @@ -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 + production_factor_series::Union{Nothing, Vector{<:Real}} = nothing macrs_option_years::Int = 5 macrs_bonus_fraction::Float64 = 1.0 diff --git a/src/core/production_factor.jl b/src/core/production_factor.jl index 6fc19e337..07749ccbf 100644 --- a/src/core/production_factor.jl +++ b/src/core/production_factor.jl @@ -278,9 +278,15 @@ end production_factor for CHP accounts for unavailability (`unavailability_periods`) of CHP due to scheduled (mostly off-peak) and "unscheduled" (on-peak) maintenance. Note: this same prod_factor should be applied to electric and thermal production + +If the user provides their own production_factor_series, that will be returned instead of generating from unavailability. """ function get_production_factor(chp::AbstractCHP, year::Int=2017, outage_start_time_step::Int=0, outage_end_time_step::Int=0, ts_per_hour::Int=1) + if !(isnothing(chp.production_factor_series)) + return chp.production_factor_series + end + prod_factor = [1.0 - chp.unavailability_hourly[i] for i in 1:8760 for _ in 1:ts_per_hour] # Ignore unavailability in time_step if it intersects with an outage interval(s) From 7c460e48ae98f5bc8f46da6373a7e4f2e53c0fc4 Mon Sep 17 00:00:00 2001 From: wbecker Date: Mon, 15 Dec 2025 23:03:02 -0700 Subject: [PATCH 24/45] Add temp dev test file for CHP production_factor_series input --- test/test_chp_prodfactor.jl | 59 +++++++++++++++++++++++++++++++++++++ 1 file changed, 59 insertions(+) create mode 100644 test/test_chp_prodfactor.jl diff --git a/test/test_chp_prodfactor.jl b/test/test_chp_prodfactor.jl new file mode 100644 index 000000000..80e833811 --- /dev/null +++ b/test/test_chp_prodfactor.jl @@ -0,0 +1,59 @@ +using Revise +using REopt +using JSON +using DelimitedFiles +using PlotlyJS +using Dates +using Test +using JuMP +using HiGHS +using DotEnv +DotEnv.load!() + + +########### CHP production_factor_series test ############# +@testset "CHP with production_factor_series" begin + # Create a simple scenario with electric-only CHP and a custom production factor + input_data = Dict( + "Site" => Dict( + "latitude" => 39.7, + "longitude" => -104.9 + ), + "ElectricLoad" => Dict( + "loads_kw" => repeat([100.0], 8760), + "year" => 2025 + ), + "ElectricTariff" => Dict( + "urdb_label" => "5ed6c1a15457a3367add15ae" + ), + "CHP" => Dict( + "is_electric_only" => true, + "fuel_cost_per_mmbtu" => 4.0, + "max_kw" => 100.0, + "min_kw" => 100.0, + "production_factor_series" => vcat(repeat([0.5], 4380), repeat([1.0], 4380)) # 50% for first half year, 100% for second half + ) + ) + + s = Scenario(input_data) + p = REoptInputs(s) + + # Verify the production factor series was properly assigned + @test !isnothing(s.chp.production_factor_series) + @test length(s.chp.production_factor_series) == 8760 + @test s.chp.production_factor_series[1] ≈ 0.5 atol=0.001 + @test s.chp.production_factor_series[8760] ≈ 1.0 atol=0.001 + + # Run the optimization + m = Model(optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false, "log_to_console" => false, "mip_rel_gap" => 0.01)) + results = run_reopt(m, p) + + # Check that production matches the pattern (lower in first half, higher in second half) + first_half_avg = sum(results["CHP"]["electric_production_series_kw"][1:4380]) / 4380 + second_half_avg = sum(results["CHP"]["electric_production_series_kw"][4381:8760]) / 4380 + + println("CHP production_factor_series test passed!") + println("CHP size: ", results["CHP"]["size_kw"], " kW") + println("First half avg production: ", round(first_half_avg, digits=2), " kW") + println("Second half avg production: ", round(second_half_avg, digits=2), " kW") +end \ No newline at end of file From f2a483b58cfca8e7cd6d074bfe08f7345ddbc46d Mon Sep 17 00:00:00 2001 From: wbecker Date: Fri, 30 Jan 2026 12:00:18 -0700 Subject: [PATCH 25/45] Add organize_chps argument to run_reopt() --- src/core/reopt.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/core/reopt.jl b/src/core/reopt.jl index a955561a3..f42ed18ed 100644 --- a/src/core/reopt.jl +++ b/src/core/reopt.jl @@ -88,8 +88,8 @@ end Method for use with Threads when running BAU in parallel with optimal scenario. """ function run_reopt(t::Tuple{JuMP.AbstractModel, AbstractInputs}) - run_reopt(t[1], t[2]; organize_pvs=false) - # must organize_pvs after adding proforma results + run_reopt(t[1], t[2]; organize_pvs=false, organize_chps=false) +# must organize_pvs/chps after adding proforma results end @@ -598,7 +598,7 @@ function build_reopt!(m::JuMP.AbstractModel, p::REoptInputs) end -function run_reopt(m::JuMP.AbstractModel, p::REoptInputs; organize_pvs=true) +function run_reopt(m::JuMP.AbstractModel, p::REoptInputs; organize_pvs=true, organize_chps=true) try build_reopt!(m, p) @@ -629,7 +629,7 @@ function run_reopt(m::JuMP.AbstractModel, p::REoptInputs; organize_pvs=true) if organize_pvs && !isempty(p.techs.pv) # do not want to organize_pvs when running BAU case in parallel b/c then proform code fails organize_multiple_pv_results(p, results) end - if organize_pvs && !isempty(p.techs.chp) # same logic as PV + if organize_chps && !isempty(p.techs.chp) # same logic as PV organize_multiple_chp_results(p, results) end From 88bb5acc96700d53dd242f142eed70347e4669f7 Mon Sep 17 00:00:00 2001 From: wbecker Date: Fri, 30 Jan 2026 17:18:57 -0700 Subject: [PATCH 26/45] Fix proforma_results() for multiple CHPs --- src/results/proforma.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/results/proforma.jl b/src/results/proforma.jl index 8fb1f40fa..b67398a0a 100644 --- a/src/results/proforma.jl +++ b/src/results/proforma.jl @@ -372,8 +372,8 @@ function update_metrics(m::Metrics, p::REoptInputs, tech::AbstractTech, tech_nam total_kw = results[tech_name]["size_kw"] existing_kw = :existing_kw in fieldnames(typeof(tech)) ? tech.existing_kw : 0 new_kw = total_kw - existing_kw - if tech_name == "CHP" - capital_cost = results["CHP"]["initial_capital_costs"] + if tech_name in [chp.name for chp in p.s.chps] + capital_cost = results[tech_name]["initial_capital_costs"] elseif tech_name in [pv.name for pv in p.s.pvs] # Check if it's a PV technology capital_cost = get_pv_initial_capex(p, tech, new_kw) else @@ -382,9 +382,9 @@ function update_metrics(m::Metrics, p::REoptInputs, tech::AbstractTech, tech_nam # owner is responsible for only new technologies operating and maintenance cost in optimal case # CHP doesn't have existing CHP, and it has different O&M cost parameters - if tech_name == "CHP" - hours_operating = sum(results["CHP"]["electric_production_series_kw"] .> 0.0) / (8760 * p.s.settings.time_steps_per_hour) - annual_om = -1 * (results["CHP"]["annual_electric_production_kwh"] * tech.om_cost_per_kwh + + if tech_name in [chp.name for chp in p.s.chps] + hours_operating = sum(results[tech_name]["electric_production_series_kw"] .> 0.0) / (8760 * p.s.settings.time_steps_per_hour) + annual_om = -1 * (results[tech_name]["annual_electric_production_kwh"] * tech.om_cost_per_kwh + new_kw * tech.om_cost_per_kw + new_kw * tech.om_cost_per_hr_per_kw_rated * hours_operating) elseif third_party @@ -397,9 +397,9 @@ function update_metrics(m::Metrics, p::REoptInputs, tech::AbstractTech, tech_nam m.om_series += escalate_om(annual_om) m.om_series_bau += escalate_om(-1 * existing_kw * tech.om_cost_per_kw) - if tech_name == "CHP" + if tech_name in [chp.name for chp in p.s.chps] escalate_fuel(val, esc_rate) = [val * (1 + esc_rate)^yr for yr in 1:years] - fuel_cost = results["CHP"]["year_one_fuel_cost_before_tax"] + fuel_cost = results[tech_name]["year_one_fuel_cost_before_tax"] m.fuel_cost_series += escalate_fuel(-1 * fuel_cost, p.s.financial.chp_fuel_cost_escalation_rate_fraction) end @@ -417,7 +417,7 @@ function update_metrics(m::Metrics, p::REoptInputs, tech::AbstractTech, tech_nam pbi_series = Float64[] pbi_series_bau = Float64[] existing_energy_bau = third_party ? get(results[tech_name], "year_one_energy_produced_kwh_bau", 0) : 0 - if tech_name == "CHP" + if tech_name in [chp.name for chp in p.s.chps] year_one_energy = results[tech_name]["annual_electric_production_kwh"] else year_one_energy = "year_one_energy_produced_kwh" in keys(results[tech_name]) ? results[tech_name]["year_one_energy_produced_kwh"] : results[tech_name]["annual_energy_produced_kwh"] From ef251c82d9b85dbdbffd6ddd9b4dd5fa411edf0a Mon Sep 17 00:00:00 2001 From: wbecker Date: Sat, 31 Jan 2026 09:28:43 -0700 Subject: [PATCH 27/45] Generalize function to get initial CapEx with a cost curve To be used for PV, CHP, and any future techs using a cost curve --- src/results/financial.jl | 34 ++++++++++++++++++++-------------- src/results/proforma.jl | 2 +- src/results/pv.jl | 2 +- 3 files changed, 22 insertions(+), 16 deletions(-) diff --git a/src/results/financial.jl b/src/results/financial.jl index 537d605f9..cd44a961e 100644 --- a/src/results/financial.jl +++ b/src/results/financial.jl @@ -252,8 +252,8 @@ function calculate_lcoe(p::REoptInputs, tech_results::Dict, tech::AbstractTech) federal_tax_rate_fraction = p.s.financial.offtaker_tax_rate_fraction end capital_costs = if typeof(tech) == PV && :tech_sizes_for_cost_curve in fieldnames(typeof(tech)) - # Use PV-specific cost curve calculation for PV tech - get_pv_initial_capex(p, tech, new_kw) + # Use cost curve calculation for PV tech + get_tech_initial_capex(tech, new_kw) else # Use simple calculation for other techs like Wind new_kw * tech.installed_cost_per_kw @@ -361,29 +361,35 @@ function get_depreciation_schedule(p::REoptInputs, tech::Union{AbstractTech,Abst return depreciation_schedule end -function get_pv_initial_capex(p::REoptInputs, pv::AbstractTech, size_kw::Float64) - cost_list = pv.installed_cost_per_kw - size_list = pv.tech_sizes_for_cost_curve - pv_size = size_kw +""" + get_tech_initial_capex(tech::AbstractTech, size_kw::Float64) + +Calculate initial capital costs for a technology with either a simple cost per kW or a segmented cost curve. +Works for PV, CHP, and other technologies with cost curves defined by tech_sizes_for_cost_curve. +""" +function get_tech_initial_capex(tech::AbstractTech, size_kw::Float64) + cost_list = tech.installed_cost_per_kw + size_list = tech.tech_sizes_for_cost_curve initial_capex = 0.0 if typeof(cost_list) == Vector{Float64} - if pv_size <= size_list[1] - initial_capex = pv_size * cost_list[1] - elseif pv_size > size_list[end] - initial_capex = pv_size * cost_list[end] + if size_kw <= size_list[1] + initial_capex = size_kw * cost_list[1] + elseif size_kw > size_list[end] + initial_capex = size_kw * cost_list[end] else for s in 2:length(size_list) - if (pv_size > size_list[s-1]) && (pv_size <= size_list[s]) + if (size_kw > size_list[s-1]) && (size_kw <= size_list[s]) slope = (cost_list[s] * size_list[s] - cost_list[s-1] * size_list[s-1]) / (size_list[s] - size_list[s-1]) - initial_capex = cost_list[s-1] * size_list[s-1] + (pv_size - size_list[s-1]) * slope + initial_capex = cost_list[s-1] * size_list[s-1] + (size_kw - size_list[s-1]) * slope end end end else - initial_capex = cost_list * pv_size + initial_capex = cost_list * size_kw end return initial_capex -end \ No newline at end of file +end + diff --git a/src/results/proforma.jl b/src/results/proforma.jl index b67398a0a..896624e56 100644 --- a/src/results/proforma.jl +++ b/src/results/proforma.jl @@ -375,7 +375,7 @@ function update_metrics(m::Metrics, p::REoptInputs, tech::AbstractTech, tech_nam if tech_name in [chp.name for chp in p.s.chps] capital_cost = results[tech_name]["initial_capital_costs"] elseif tech_name in [pv.name for pv in p.s.pvs] # Check if it's a PV technology - capital_cost = get_pv_initial_capex(p, tech, new_kw) + capital_cost = get_tech_initial_capex(tech, new_kw) else capital_cost = new_kw * tech.installed_cost_per_kw end diff --git a/src/results/pv.jl b/src/results/pv.jl index b76996cf4..c5db17af6 100644 --- a/src/results/pv.jl +++ b/src/results/pv.jl @@ -38,7 +38,7 @@ function add_pv_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="") r["installed_cost_per_kw"] = pv_tech.installed_cost_per_kw else # Get cost from cost curve - r["installed_cost_per_kw"] = get_pv_initial_capex(p, pv_tech, r["size_kw"]) / r["size_kw"] + r["installed_cost_per_kw"] = get_tech_initial_capex(pv_tech, r["size_kw"]) / r["size_kw"] end if !isnothing(pv_tech.size_class) && !isempty(pv_tech.tech_sizes_for_cost_curve) From 6be57bb1fbe787329e0c37dba4ad94a1db37b5f7 Mon Sep 17 00:00:00 2001 From: wbecker Date: Sat, 31 Jan 2026 09:39:37 -0700 Subject: [PATCH 28/45] Fix CHP om_cost_per_hr_per_kw_rated accounting in proforma.jl --- src/results/proforma.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/results/proforma.jl b/src/results/proforma.jl index 896624e56..eb0534916 100644 --- a/src/results/proforma.jl +++ b/src/results/proforma.jl @@ -383,7 +383,7 @@ function update_metrics(m::Metrics, p::REoptInputs, tech::AbstractTech, tech_nam # owner is responsible for only new technologies operating and maintenance cost in optimal case # CHP doesn't have existing CHP, and it has different O&M cost parameters if tech_name in [chp.name for chp in p.s.chps] - hours_operating = sum(results[tech_name]["electric_production_series_kw"] .> 0.0) / (8760 * p.s.settings.time_steps_per_hour) + hours_operating = sum(results[tech_name]["electric_production_series_kw"] .> 0.0) * p.hours_per_time_step annual_om = -1 * (results[tech_name]["annual_electric_production_kwh"] * tech.om_cost_per_kwh + new_kw * tech.om_cost_per_kw + new_kw * tech.om_cost_per_hr_per_kw_rated * hours_operating) From 4f3aa7e34ff1aae628c527a01da8ac995d60e21e Mon Sep 17 00:00:00 2001 From: wbecker Date: Sat, 31 Jan 2026 09:53:25 -0700 Subject: [PATCH 29/45] Fix individual CHP reporting when multiple CHPs --- src/results/chp.jl | 24 +++++++++++++++++++----- 1 file changed, 19 insertions(+), 5 deletions(-) diff --git a/src/results/chp.jl b/src/results/chp.jl index 0da0976e7..c33def301 100644 --- a/src/results/chp.jl +++ b/src/results/chp.jl @@ -131,14 +131,28 @@ function get_chp_results_for_tech(m::JuMP.AbstractModel, p::REoptInputs, chp_nam end r["thermal_to_process_heat_load_series_mmbtu_per_hour"] = round.(CHPToProcessHeatKW ./ KWH_PER_MMBTU, digits=5) - r["year_one_fuel_cost_before_tax"] = round(value(m[:TotalCHPFuelCosts] / p.pwf_fuel[chp_name]), digits=3) + # Calculate individual CHP fuel cost (not the total across all CHPs) + chp_fuel_cost_lifecycle = sum(p.pwf_fuel[chp_name] * value(m[:dvFuelUsage][chp_name, ts]) * p.fuel_cost_per_kwh[chp_name][ts] for ts in p.time_steps) + r["year_one_fuel_cost_before_tax"] = round(chp_fuel_cost_lifecycle / p.pwf_fuel[chp_name], digits=3) r["year_one_fuel_cost_after_tax"] = r["year_one_fuel_cost_before_tax"] * (1 - p.s.financial.offtaker_tax_rate_fraction) - r["lifecycle_fuel_cost_after_tax"] = round(value(m[:TotalCHPFuelCosts]) * (1- p.s.financial.offtaker_tax_rate_fraction), digits=3) + r["lifecycle_fuel_cost_after_tax"] = round(chp_fuel_cost_lifecycle * (1- p.s.financial.offtaker_tax_rate_fraction), digits=3) + #Standby charges and hourly O&M - r["year_one_standby_cost_before_tax"] = round(value(m[Symbol("TotalCHPStandbyCharges")]) / p.pwf_e, digits=0) + r["year_one_standby_cost_before_tax"] = round(chp.standby_rate_per_kw_per_month * 12 * value(m[Symbol("dvSize"*_n)][chp_name]), digits=0) r["year_one_standby_cost_after_tax"] = r["year_one_standby_cost_before_tax"] * (1 - p.s.financial.offtaker_tax_rate_fraction) - r["lifecycle_standby_cost_after_tax"] = round(value(m[Symbol("TotalCHPStandbyCharges")]) * (1 - p.s.financial.offtaker_tax_rate_fraction), digits=0) - r["initial_capital_costs"] = round(value(m[Symbol("CHPCapexNoIncentives")]), digits=2) + r["lifecycle_standby_cost_after_tax"] = round(r["year_one_standby_cost_before_tax"] * p.pwf_e * (1 - p.s.financial.offtaker_tax_rate_fraction), digits=0) + + # Calculate individual CHP capital costs + size_kw = value(m[Symbol("dvSize"*_n)][chp_name]) + suppl_firing_size = value(m[Symbol("dvSupplementaryFiringSize"*_n)][chp_name]) + + # Base CHP capital cost using generalized function + capital_cost = get_tech_initial_capex(chp, size_kw) + + # Add supplementary firing capital cost + capital_cost += chp.supplementary_firing_capital_cost_per_kw * suppl_firing_size + + r["initial_capital_costs"] = round(capital_cost, digits=2) return r end From bfe9f3b1c793ffc514a294f616ccb08f9eac9fa2 Mon Sep 17 00:00:00 2001 From: wbecker Date: Sat, 31 Jan 2026 09:57:08 -0700 Subject: [PATCH 30/45] Update test_multiple_chps after proforma alignment --- test/scenarios/multiple_chps.json | 39 ++++++++++++++++------- test/test_multiple_chps.jl | 52 ++++++++++++++++++++----------- 2 files changed, 62 insertions(+), 29 deletions(-) diff --git a/test/scenarios/multiple_chps.json b/test/scenarios/multiple_chps.json index 6d07508b7..1d67a77b6 100644 --- a/test/scenarios/multiple_chps.json +++ b/test/scenarios/multiple_chps.json @@ -16,16 +16,16 @@ "year": 2017 }, "ElectricTariff": { - "blended_annual_energy_rate": 0.12, + "blended_annual_energy_rate": 0.15, "blended_annual_demand_rate": 15.0 }, "DomesticHotWaterLoad": { "doe_reference_name": "Hospital", - "annual_mmbtu": 500000.0 + "annual_mmbtu": 5000.0 }, "SpaceHeatingLoad": { "doe_reference_name": "Hospital", - "annual_mmbtu": 1000000.0 + "annual_mmbtu": 10000.0 }, "ExistingBoiler": { "fuel_cost_per_mmbtu": 12.0 @@ -34,22 +34,39 @@ { "name": "CHP_recip_engine", "prime_mover": "recip_engine", - "size_class": 1, - "min_kw": 100.0, - "max_kw": 500.0, + "size_class": 3, + "installed_cost_per_kw": 3000.0, + "macrs_bonus_fraction": 0.0, + "macrs_option_years": 0, + "min_turn_down_fraction": 0.0, + "min_kw": 150.0, + "max_kw": 150.0, "fuel_cost_per_mmbtu": 8.0, "can_serve_dhw": true, - "can_serve_space_heating": true + "can_serve_space_heating": true, + "standby_rate_per_kw_per_month": 5.0, + "om_cost_per_kwh": 0.05, + "om_cost_per_kw": 10.0, + "om_cost_per_hr_per_kw_rated": 0.05 }, { "name": "CHP_micro_turbine", "prime_mover": "micro_turbine", - "size_class": 2, - "min_kw": 50.0, - "max_kw": 300.0, + "size_class": 3, + "installed_cost_per_kw": 3500.0, + "macrs_bonus_fraction": 0.0, + "macrs_option_years": 0, + "min_turn_down_fraction": 0.0, + "min_kw": 100.0, + "max_kw": 100.0, "fuel_cost_per_mmbtu": 8.0, "can_serve_dhw": true, - "can_serve_space_heating": false + "can_serve_space_heating": false, + "standby_rate_per_kw_per_month": 5.0, + "om_cost_per_kw": 10.0, + "om_cost_per_kwh": 0.05, + "om_cost_per_hr_per_kw_rated": 0.05, + "supplementary_firing_max_steam_ratio": 1.2 } ] } diff --git a/test/test_multiple_chps.jl b/test/test_multiple_chps.jl index cb6fe2977..94fb398c8 100644 --- a/test/test_multiple_chps.jl +++ b/test/test_multiple_chps.jl @@ -1,26 +1,26 @@ -using Revise -using REopt -using JSON -using DelimitedFiles -using PlotlyJS -using Dates -using Test -using JuMP -using HiGHS -using DotEnv -DotEnv.load!() +# using Revise +# using REopt +# using JSON +# using DelimitedFiles +# using PlotlyJS +# using Dates +# using Test +# using JuMP +# using HiGHS +# using DotEnv +# DotEnv.load!() ############### Multiple CHPs Test ################### + +# With BAU Scenario input_data = JSON.parsefile("./scenarios/multiple_chps.json") s = Scenario(input_data) inputs = REoptInputs(s) -m1 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) -results = run_reopt(m1, inputs) -# m2 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) -# results = run_reopt([m1,m2], inputs) +m1 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.05, "output_flag" => false, "log_to_console" => false)) +m2 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.05, "output_flag" => false, "log_to_console" => false)) +results = run_reopt([m1,m2], inputs) -# Check that both CHPs are in results @test length(results["CHP"]) == 2 CHP_recip_engine = results["CHP"][findfirst(chp -> chp["name"] == "CHP_recip_engine", results["CHP"])] CHP_micro_turbine = results["CHP"][findfirst(chp -> chp["name"] == "CHP_micro_turbine", results["CHP"])] @@ -30,5 +30,21 @@ CHP_micro_turbine = results["CHP"][findfirst(chp -> chp["name"] == "CHP_micro_tu @test CHP_micro_turbine["size_kw"] > 10.0 # Check that results include electric production for each CHP -@test sum(CHP_recip_engine["electric_production_series_kw"]) > 5000.0 -@test sum(CHP_micro_turbine["electric_production_series_kw"]) > 5000.0 \ No newline at end of file +@test sum(CHP_recip_engine["electric_production_series_kw"]) > 500.0 +@test sum(CHP_micro_turbine["electric_production_series_kw"]) > 500.0 + +# Test that LCC matches sum of discounted cashflows for optimal case (within 0.1% tolerance) +# Note, cashflows are negative where lcc is positive, so negating cashflows in comparison +optimal_lcc = results["Financial"]["lcc"] +optimal_cashflow_sum = -1*sum(results["Financial"]["offtaker_discounted_annual_free_cashflows"]) +@test isapprox(optimal_lcc, optimal_cashflow_sum, rtol=0.001) + +# Test that LCC matches sum of discounted cashflows for BAU case (within 0.1% tolerance) +bau_lcc = results["Financial"]["lcc_bau"] +bau_cashflow_sum = -1*sum(results["Financial"]["offtaker_discounted_annual_free_cashflows_bau"]) +@test isapprox(bau_lcc, bau_cashflow_sum, rtol=0.001) + +# Add a test to compare npv with difference in optimal and bau cashflows +npv_calculated = bau_cashflow_sum - optimal_cashflow_sum +npv_reported = results["Financial"]["npv"] +@test isapprox(npv_calculated, npv_reported, rtol=0.001) From b4e0c388c49f06041ff04ebb27801fbb6c000f41 Mon Sep 17 00:00:00 2001 From: wbecker Date: Sat, 31 Jan 2026 09:57:49 -0700 Subject: [PATCH 31/45] Add test_multiple_chps.jl to runtests.jl for CI Actions --- test/runtests.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/test/runtests.jl b/test/runtests.jl index 7179d1d2d..171ca1d1f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1144,6 +1144,10 @@ else # run HiGHS tests empty!(m2) GC.gc() end + + @testset "Multiple CHPs" begin + include("test_multiple_chps.jl") + end end @testset verbose=true "FlexibleHVAC" begin From f92edcb187d75c7da4b0023014e68438be1a7498 Mon Sep 17 00:00:00 2001 From: wbecker Date: Sat, 31 Jan 2026 10:45:55 -0700 Subject: [PATCH 32/45] Add test_avoid_chp_binaries.jl to runtests.jl --- test/runtests.jl | 4 ++++ test/test_avoid_chp_binary.jl | 22 +++++++++++----------- 2 files changed, 15 insertions(+), 11 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 171ca1d1f..7f692e917 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1148,6 +1148,10 @@ else # run HiGHS tests @testset "Multiple CHPs" begin include("test_multiple_chps.jl") end + + @testset "Avoid CHP Binary" begin + include("test_avoid_chp_binary.jl") + end end @testset verbose=true "FlexibleHVAC" begin diff --git a/test/test_avoid_chp_binary.jl b/test/test_avoid_chp_binary.jl index 199683cd8..0179b89f1 100644 --- a/test/test_avoid_chp_binary.jl +++ b/test/test_avoid_chp_binary.jl @@ -1,15 +1,15 @@ -using Revise -using REopt -using JSON -using DelimitedFiles -using PlotlyJS -using Dates -using Test -using JuMP -using HiGHS -using DotEnv -DotEnv.load!() +# using Revise +# using REopt +# using JSON +# using DelimitedFiles +# using PlotlyJS +# using Dates +# using Test +# using JuMP +# using HiGHS +# using DotEnv +# DotEnv.load!() ############### CHP Binary Creation Tests ################### # Tests to verify that binCHPIsOnInTS is only created when needed From 1cebf4aa664feb7482f77f05d22e17ce479188a2 Mon Sep 17 00:00:00 2001 From: wbecker Date: Sat, 31 Jan 2026 13:42:28 -0700 Subject: [PATCH 33/45] Update expected solar dataset query dataset --- test/runtests.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 7f692e917..5c03e8cec 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -145,18 +145,18 @@ else # run HiGHS tests dataset, distance, datasource = REopt.call_solar_dataset_api(latitude, longitude, radius) @test dataset == "nsrdb" - # 3. Oulu, Findland + # 3. Oulu, Finland latitude, longitude = 65.0102196310875, 25.465387094897675 radius = 0 dataset, distance, datasource = REopt.call_solar_dataset_api(latitude, longitude, radius) - @test dataset == "intl" + @test dataset == "nsrdb" # 4. Fairbanks, AK site = "Fairbanks" latitude, longitude = 64.84112047064114, -147.71570239058084 radius = 20 dataset, distance, datasource = REopt.call_solar_dataset_api(latitude, longitude, radius) - @test dataset == "tmy3" + @test dataset == "nsrdb" end @testset "ASHP min allowable size and COP, CF Profiles" begin From b263720992736a0a3365e56fcb2715e9d919a416 Mon Sep 17 00:00:00 2001 From: wbecker Date: Sun, 1 Feb 2026 22:51:00 -0700 Subject: [PATCH 34/45] Fix issues after merging offgrid-chp --- .../operating_reserve_constraints.jl | 8 +++----- src/core/reopt_inputs.jl | 18 +++++++++++------- src/core/techs.jl | 6 +++--- 3 files changed, 17 insertions(+), 15 deletions(-) diff --git a/src/constraints/operating_reserve_constraints.jl b/src/constraints/operating_reserve_constraints.jl index 193c03293..78bc2e455 100644 --- a/src/constraints/operating_reserve_constraints.jl +++ b/src/constraints/operating_reserve_constraints.jl @@ -42,11 +42,9 @@ function add_operating_reserve_constraints(m, p; _n="") 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 + @constraint(m, [t in intersect(p.techs.chp, p.techs.providing_oper_res), ts in p.time_steps_without_grid], + m[Symbol("dvOpResFromTechs"*_n)][t,ts] <= m[:binCHPIsOnInTS][t, ts] * p.max_sizes[t] + ) 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) diff --git a/src/core/reopt_inputs.jl b/src/core/reopt_inputs.jl index 6816aef87..151b048d1 100644 --- a/src/core/reopt_inputs.jl +++ b/src/core/reopt_inputs.jl @@ -839,12 +839,16 @@ function setup_absorption_chiller_inputs(s::AbstractScenario, max_sizes, min_siz if isempty(s.chps) thermal_factor = 1.0 else - # Use the cooling_thermal_factor from the first CHP - if s.chps[1].cooling_thermal_factor == 0.0 - throw(@error("The CHP cooling_thermal_factor is 0.0 which implies that CHP cannot serve AbsorptionChiller. If you - want to model CHP and AbsorptionChiller, you must specify a cooling_thermal_factor greater than 0.0")) + # Use the maximum cooling_thermal_factor from all CHPs that can supply the absorption chiller + chps_serving_absorption_chiller = filter(chp -> chp.can_serve_cooling, s.chps) + if isempty(chps_serving_absorption_chiller) + thermal_factor = 1.0 else - thermal_factor = s.chps[1].cooling_thermal_factor + thermal_factor = maximum(chp.cooling_thermal_factor for chp in chps_serving_absorption_chiller) + if thermal_factor == 0.0 + throw(@error("All CHPs have cooling_thermal_factor of 0.0 which implies that CHP cannot serve AbsorptionChiller. If you + want to model CHP and AbsorptionChiller, you must specify a cooling_thermal_factor greater than 0.0 for at least one CHP")) + end end end thermal_cop["AbsorptionChiller"] = s.absorption_chiller.cop_thermal * thermal_factor @@ -1362,8 +1366,8 @@ function setup_operating_reserve_fraction(s::AbstractScenario, techs_operating_r 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 + for chp in s.chps + techs_operating_reserve_req_fraction[chp.name] = chp.operating_reserve_required_fraction end return nothing diff --git a/src/core/techs.jl b/src/core/techs.jl index 3c69f0cd7..839f4facf 100644 --- a/src/core/techs.jl +++ b/src/core/techs.jl @@ -203,10 +203,10 @@ function Techs(s::Scenario) push!(techs_can_serve_process_heat, chp.name) end if s.settings.off_grid_flag - if s.chp.operating_reserve_required_fraction > 0.0 - push!(requiring_oper_res, "CHP") + if chp.operating_reserve_required_fraction > 0.0 + push!(requiring_oper_res, chp.name) else - push!(providing_oper_res, "CHP") + push!(providing_oper_res, chp.name) end end end From 986156176c87fabf638f159fc5a3c24389031dc6 Mon Sep 17 00:00:00 2001 From: wbecker Date: Sun, 1 Feb 2026 23:26:43 -0700 Subject: [PATCH 35/45] Add offgrid CHP tests to runtests --- test/runtests.jl | 4 ++++ test/test_chp_offgrid.jl | 26 +++++++++++++------------- 2 files changed, 17 insertions(+), 13 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 5c03e8cec..94d23cf62 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1152,6 +1152,10 @@ else # run HiGHS tests @testset "Avoid CHP Binary" begin include("test_avoid_chp_binary.jl") end + + @testset "Off-Grid CHP" begin + include("test_chp_offgrid.jl") + end end @testset verbose=true "FlexibleHVAC" begin diff --git a/test/test_chp_offgrid.jl b/test/test_chp_offgrid.jl index 213de6b1a..bfc8fa7dc 100644 --- a/test/test_chp_offgrid.jl +++ b/test/test_chp_offgrid.jl @@ -1,14 +1,14 @@ -using Revise -using REopt -using JSON -using DelimitedFiles -using PlotlyJS -using Dates -using Test -using JuMP -using HiGHS -using DotEnv -DotEnv.load!() +# 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 ################### @@ -30,7 +30,7 @@ s = Scenario(input_data) inputs = REoptInputs(s) println("\nScenario Setup:") -println(" - CHP operating_reserve_required_fraction: ", s.chp.operating_reserve_required_fraction) +println(" - CHP operating_reserve_required_fraction: ", s.chps[1].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) @@ -49,7 +49,7 @@ println(" - Load met fraction: ", round(results["ElectricLoad"]["offgrid_load_m # 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) +chp_or_required = sum(chp_production_to_load .* s.chps[1].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 From 090c09ac33e2ad4d0d2167cca7497a0eed5b48af Mon Sep 17 00:00:00 2001 From: wbecker Date: Mon, 2 Feb 2026 09:05:47 -0700 Subject: [PATCH 36/45] Fix issue with CHP -> AbsChl fieldname --- src/core/reopt_inputs.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/core/reopt_inputs.jl b/src/core/reopt_inputs.jl index 151b048d1..520c99bb1 100644 --- a/src/core/reopt_inputs.jl +++ b/src/core/reopt_inputs.jl @@ -839,12 +839,12 @@ function setup_absorption_chiller_inputs(s::AbstractScenario, max_sizes, min_siz if isempty(s.chps) thermal_factor = 1.0 else - # Use the maximum cooling_thermal_factor from all CHPs that can supply the absorption chiller - chps_serving_absorption_chiller = filter(chp -> chp.can_serve_cooling, s.chps) - if isempty(chps_serving_absorption_chiller) + # Use the maximum cooling_thermal_factor from all CHPs that have thermal output (not electric-only) + chps_with_thermal = filter(chp -> !chp.is_electric_only, s.chps) + if isempty(chps_with_thermal) thermal_factor = 1.0 else - thermal_factor = maximum(chp.cooling_thermal_factor for chp in chps_serving_absorption_chiller) + thermal_factor = maximum(chp.cooling_thermal_factor for chp in chps_with_thermal) if thermal_factor == 0.0 throw(@error("All CHPs have cooling_thermal_factor of 0.0 which implies that CHP cannot serve AbsorptionChiller. If you want to model CHP and AbsorptionChiller, you must specify a cooling_thermal_factor greater than 0.0 for at least one CHP")) From ddb906b3c4e35570cb08bf56fd4b8c82f6e9482b Mon Sep 17 00:00:00 2001 From: wbecker Date: Mon, 2 Feb 2026 11:28:01 -0700 Subject: [PATCH 37/45] Fix merge issues after merging ramp-chp --- src/constraints/chp_constraints.jl | 29 ++++++++++++++++------------- src/core/reopt_inputs.jl | 3 ++- 2 files changed, 18 insertions(+), 14 deletions(-) diff --git a/src/constraints/chp_constraints.jl b/src/constraints/chp_constraints.jl index aab8fac1d..f66c3af74 100644 --- a/src/constraints/chp_constraints.jl +++ b/src/constraints/chp_constraints.jl @@ -165,17 +165,20 @@ end function add_chp_ramp_rate_constraints(m, p; _n="") # Ramp rate constraints limit how quickly CHP production can change between consecutive timesteps - # Ramp up constraint - @constraint(m, CHPRampUp[t in p.techs.chp, ts in p.time_steps[2:end]], - m[Symbol("dvRatedProduction"*_n)][t, ts] - m[Symbol("dvRatedProduction"*_n)][t, ts-1] <= - p.s.chp.ramp_rate_fraction_per_hour * m[Symbol("dvSize"*_n)][t] / p.s.settings.time_steps_per_hour - ) - - # Ramp down constraint - @constraint(m, CHPRampDown[t in p.techs.chp, ts in p.time_steps[2:end]], - m[Symbol("dvRatedProduction"*_n)][t, ts-1] - m[Symbol("dvRatedProduction"*_n)][t, ts] <= - p.s.chp.ramp_rate_fraction_per_hour * m[Symbol("dvSize"*_n)][t] / p.s.settings.time_steps_per_hour - ) + # Loop through each CHP and add constraints with tech-specific ramp rate parameters + for t in p.techs.chp + # Ramp up constraint + @constraint(m, [ts in p.time_steps[2:end]], + m[Symbol("dvRatedProduction"*_n)][t, ts] - m[Symbol("dvRatedProduction"*_n)][t, ts-1] <= + p.chp_params[t][:ramp_rate_fraction_per_hour] * m[Symbol("dvSize"*_n)][t] / p.s.settings.time_steps_per_hour + ) + + # Ramp down constraint + @constraint(m, [ts in p.time_steps[2:end]], + m[Symbol("dvRatedProduction"*_n)][t, ts-1] - m[Symbol("dvRatedProduction"*_n)][t, ts] <= + p.chp_params[t][:ramp_rate_fraction_per_hour] * m[Symbol("dvSize"*_n)][t] / p.s.settings.time_steps_per_hour + ) + end end @@ -275,8 +278,8 @@ function add_chp_constraints(m, p; _n="") add_chp_thermal_production_constraints(m, p; _n=_n) add_chp_rated_prod_constraint(m, p; _n=_n) - # Add ramp rate constraints if ramp_rate_fraction_per_hour < 1.0 - if p.s.chp.ramp_rate_fraction_per_hour < 1.0 / p.s.settings.time_steps_per_hour + # Add ramp rate constraints if any CHP has ramp_rate_fraction_per_hour < 1.0 + if any(p.chp_params[t][:ramp_rate_fraction_per_hour] < 1.0 / p.s.settings.time_steps_per_hour for t in p.techs.chp) add_chp_ramp_rate_constraints(m, p; _n=_n) end diff --git a/src/core/reopt_inputs.jl b/src/core/reopt_inputs.jl index 520c99bb1..ad47ae614 100644 --- a/src/core/reopt_inputs.jl +++ b/src/core/reopt_inputs.jl @@ -908,7 +908,8 @@ function setup_chp_inputs(s::AbstractScenario, max_sizes, min_sizes, cap_cost_sl :min_turn_down_fraction => chp.min_turn_down_fraction, :supplementary_firing_efficiency => chp.supplementary_firing_efficiency, :supplementary_firing_max_steam_ratio => chp.supplementary_firing_max_steam_ratio, - :om_cost_per_kwh => chp.om_cost_per_kwh + :om_cost_per_kwh => chp.om_cost_per_kwh, + :ramp_rate_fraction_per_hour => chp.ramp_rate_fraction_per_hour ) end return nothing From c089e14caff3b4b74b08c2ff2dff9115f8ca1080 Mon Sep 17 00:00:00 2001 From: wbecker Date: Mon, 2 Feb 2026 11:30:17 -0700 Subject: [PATCH 38/45] Add test for CHP ramp to runtests.jl --- test/runtests.jl | 4 ++ test/test_ramp.jl | 114 +++++++++++++++++++++++----------------------- 2 files changed, 61 insertions(+), 57 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 94d23cf62..e9add94b3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1156,6 +1156,10 @@ else # run HiGHS tests @testset "Off-Grid CHP" begin include("test_chp_offgrid.jl") end + + @testset "CHP Ramp Rate" begin + include("test_ramp.jl") + end end @testset verbose=true "FlexibleHVAC" begin diff --git a/test/test_ramp.jl b/test/test_ramp.jl index c4f558b89..f1f96b1b7 100644 --- a/test/test_ramp.jl +++ b/test/test_ramp.jl @@ -128,60 +128,60 @@ println("Utility to Load (kWh): ", round(sum(results["ElectricUtility"]["annual_ # println("Simple Payback (years): ", round(results["Financial"]["simple_payback_years"], digits=2)) println("==========================================\n") -# Create stacked area chart showing how CHP and battery serve the load -load_series = results["ElectricLoad"]["load_series_kw"] -chp_to_load = results["CHP"]["electric_to_load_series_kw"] -batt_to_load = results["ElectricStorage"]["storage_to_load_series_kw"] -grid_to_load = results["ElectricUtility"]["electric_to_load_series_kw"] - -# Create time steps array (in hours for x-axis) -time_hours = collect(1:length(load_series)) ./ ts_per_hour - -# Create stacked area chart -plt = plot( - [ - scatter( - x=time_hours, - y=grid_to_load, - mode="lines", - name="Grid to Load", - fill="tozeroy", - line=attr(width=0.5, color="rgb(128,128,128)"), - fillcolor="rgba(128,128,128,0.5)" - ), - scatter( - x=time_hours, - y=grid_to_load .+ chp_to_load, - mode="lines", - name="CHP to Load", - fill="tonexty", - line=attr(width=0.5, color="rgb(255,128,0)"), - fillcolor="rgba(255,128,0,0.6)" - ), - scatter( - x=time_hours, - y=grid_to_load .+ chp_to_load .+ batt_to_load, - mode="lines", - name="Battery to Load", - fill="tonexty", - line=attr(width=0.5, color="rgb(0,128,255)"), - fillcolor="rgba(0,128,255,0.6)" - ), - scatter( - x=time_hours, - y=load_series, - mode="lines", - name="Total Load", - line=attr(width=2, color="rgb(0,0,0)", dash="dot") - ) - ], - Layout( - title="Electric Load Service Breakdown - CHP with Ramp Rate and Battery Support", - xaxis=attr(title="Time (hours)", range=[0, 8760]), - yaxis=attr(title="Power (kW)"), - showlegend=true, - hovermode="x unified", - legend=attr(x=1.02, y=1, xanchor="left", yanchor="top") - ) -) -display(plt) \ No newline at end of file +# # Create stacked area chart showing how CHP and battery serve the load +# load_series = results["ElectricLoad"]["load_series_kw"] +# chp_to_load = results["CHP"]["electric_to_load_series_kw"] +# batt_to_load = results["ElectricStorage"]["storage_to_load_series_kw"] +# grid_to_load = results["ElectricUtility"]["electric_to_load_series_kw"] + +# # Create time steps array (in hours for x-axis) +# time_hours = collect(1:length(load_series)) ./ ts_per_hour + +# # Create stacked area chart +# plt = plot( +# [ +# scatter( +# x=time_hours, +# y=grid_to_load, +# mode="lines", +# name="Grid to Load", +# fill="tozeroy", +# line=attr(width=0.5, color="rgb(128,128,128)"), +# fillcolor="rgba(128,128,128,0.5)" +# ), +# scatter( +# x=time_hours, +# y=grid_to_load .+ chp_to_load, +# mode="lines", +# name="CHP to Load", +# fill="tonexty", +# line=attr(width=0.5, color="rgb(255,128,0)"), +# fillcolor="rgba(255,128,0,0.6)" +# ), +# scatter( +# x=time_hours, +# y=grid_to_load .+ chp_to_load .+ batt_to_load, +# mode="lines", +# name="Battery to Load", +# fill="tonexty", +# line=attr(width=0.5, color="rgb(0,128,255)"), +# fillcolor="rgba(0,128,255,0.6)" +# ), +# scatter( +# x=time_hours, +# y=load_series, +# mode="lines", +# name="Total Load", +# line=attr(width=2, color="rgb(0,0,0)", dash="dot") +# ) +# ], +# Layout( +# title="Electric Load Service Breakdown - CHP with Ramp Rate and Battery Support", +# xaxis=attr(title="Time (hours)", range=[0, 8760]), +# yaxis=attr(title="Power (kW)"), +# showlegend=true, +# hovermode="x unified", +# legend=attr(x=1.02, y=1, xanchor="left", yanchor="top") +# ) +# ) +# display(plt) \ No newline at end of file From db47435191cd675c0c6a7f420156c720d357b132 Mon Sep 17 00:00:00 2001 From: wbecker Date: Mon, 2 Feb 2026 16:59:11 -0700 Subject: [PATCH 39/45] Add CHP production factor test to runtests.jl --- test/runtests.jl | 5 +++++ test/test_chp_prodfactor.jl | 30 +++++++++++++++--------------- 2 files changed, 20 insertions(+), 15 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index e9add94b3..0abaa40fe 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1160,6 +1160,11 @@ else # run HiGHS tests @testset "CHP Ramp Rate" begin include("test_ramp.jl") end + + @testset "CHP Production Factor" begin + include("test_chp_prodfactor.jl") + end + end @testset verbose=true "FlexibleHVAC" begin diff --git a/test/test_chp_prodfactor.jl b/test/test_chp_prodfactor.jl index 80e833811..5209eec34 100644 --- a/test/test_chp_prodfactor.jl +++ b/test/test_chp_prodfactor.jl @@ -1,14 +1,14 @@ -using Revise -using REopt -using JSON -using DelimitedFiles -using PlotlyJS -using Dates -using Test -using JuMP -using HiGHS -using DotEnv -DotEnv.load!() +# using Revise +# using REopt +# using JSON +# using DelimitedFiles +# using PlotlyJS +# using Dates +# using Test +# using JuMP +# using HiGHS +# using DotEnv +# DotEnv.load!() ########### CHP production_factor_series test ############# @@ -39,10 +39,10 @@ DotEnv.load!() p = REoptInputs(s) # Verify the production factor series was properly assigned - @test !isnothing(s.chp.production_factor_series) - @test length(s.chp.production_factor_series) == 8760 - @test s.chp.production_factor_series[1] ≈ 0.5 atol=0.001 - @test s.chp.production_factor_series[8760] ≈ 1.0 atol=0.001 + @test !isnothing(s.chps[1].production_factor_series) + @test length(s.chps[1].production_factor_series) == 8760 + @test s.chps[1].production_factor_series[1] ≈ 0.5 atol=0.001 + @test s.chps[1].production_factor_series[8760] ≈ 1.0 atol=0.001 # Run the optimization m = Model(optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false, "log_to_console" => false, "mip_rel_gap" => 0.01)) From 4ef6458b313f91375d90e4560bdd3d764652b9b3 Mon Sep 17 00:00:00 2001 From: wbecker Date: Sun, 3 May 2026 13:59:31 -0600 Subject: [PATCH 40/45] Fix issues after merging develop 1 Co-authored-by: Copilot --- src/constraints/chp_constraints.jl | 8 +- src/core/scenario.jl | 139 ++++++++++++----------------- test/runtests.jl | 2 +- 3 files changed, 63 insertions(+), 86 deletions(-) diff --git a/src/constraints/chp_constraints.jl b/src/constraints/chp_constraints.jl index f1f5847c3..833adbc90 100644 --- a/src/constraints/chp_constraints.jl +++ b/src/constraints/chp_constraints.jl @@ -228,7 +228,9 @@ load or send to waste in dispatch. """ function add_chp_to_absorption_chiller_only_constraints(m, p; _n="") monthly_timesteps = get_monthly_time_steps(p.s.electric_load.year; time_steps_per_hour=p.s.settings.time_steps_per_hour) - for mth in p.s.chp.months_serving_absorption_chiller_only + chp_with_ac_only = filter(chp -> chp.serve_absorption_chiller_only, p.s.chps) + months_to_restrict = unique(vcat([chp.months_serving_absorption_chiller_only for chp in chp_with_ac_only]...)) + for mth in months_to_restrict @constraint(m, [t in p.techs.chp, q in [p.s.absorption_chiller.heating_load_input], ts in monthly_timesteps[mth]], m[Symbol("dvProductionToWaste"*_n)]["CHP",q,ts] + m[Symbol("dvHeatToAbsorptionChiller"*_n)]["CHP",q,ts] == m[Symbol("dvHeatingProduction"*_n)]["CHP",q,ts] ) @@ -368,11 +370,11 @@ function add_chp_constraints(m, p; _n="") end end - if !isempty(p.techs.absorption_chiller) && p.s.chp.serve_absorption_chiller_only + if !isempty(p.techs.absorption_chiller) && any(chp.serve_absorption_chiller_only for chp in p.s.chps) add_chp_to_absorption_chiller_only_constraints(m, p; _n=_n) end - if p.s.chp.follow_electrical_load + if any(chp.follow_electrical_load for chp in p.s.chps) add_chp_electrical_load_following_constraints(m, p; _n=_n) end end diff --git a/src/core/scenario.jl b/src/core/scenario.jl index 21ec94c13..7adc5cb15 100644 --- a/src/core/scenario.jl +++ b/src/core/scenario.jl @@ -408,45 +408,6 @@ function Scenario(d::Dict; flex_hvac_from_json=false) end - chps = CHP[] - chp_prime_mover = nothing - if haskey(d, "CHP") - chp_array = isa(d["CHP"], AbstractArray) ? d["CHP"] : [d["CHP"]] - for (i, chp_dict) in enumerate(chp_array) - electric_only = get(chp_dict, "is_electric_only", false) || get(chp_dict, "thermal_efficiency_full_load", 0.5) == 0.0 - - # Set default name if not provided - if !haskey(chp_dict, "name") - chp_dict["name"] = length(chp_array) > 1 ? "CHP$i" : "CHP" - end - - if !isnothing(existing_boiler) && !electric_only - total_fuel_heating_load_mmbtu_per_hour = (space_heating_load.loads_kw + dhw_load.loads_kw + process_heat_load.loads_kw) / existing_boiler.efficiency / KWH_PER_MMBTU - avg_boiler_fuel_load_mmbtu_per_hour = sum(total_fuel_heating_load_mmbtu_per_hour) / length(total_fuel_heating_load_mmbtu_per_hour) - chp = CHP(chp_dict; - avg_boiler_fuel_load_mmbtu_per_hour = avg_boiler_fuel_load_mmbtu_per_hour, - existing_boiler = existing_boiler, - electric_load_series_kw = electric_load.loads_kw, - year = electric_load.year, - sector = site.sector, - 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(chp_dict; - electric_load_series_kw = electric_load.loads_kw, - year = electric_load.year, - sector = site.sector, - federal_procurement_type = site.federal_procurement_type, - off_grid_flag = settings.off_grid_flag) - end - push!(chps, chp) - end - # Store first CHP's prime mover for backward compatibility if needed - if !isempty(chps) - chp_prime_mover = chps[1].prime_mover - end - end - max_cooling_demand_kw = 0 if haskey(d, "CoolingLoad") && !haskey(d, "FlexibleHVAC") d["CoolingLoad"] = convert(Dict{String, Any}, d["CoolingLoad"]) @@ -504,51 +465,65 @@ function Scenario(d::Dict; flex_hvac_from_json=false) end - chp = nothing + chps = CHP[] + chp_prime_mover = nothing if haskey(d, "CHP") - electric_only = get(d["CHP"], "is_electric_only", false) || get(d["CHP"], "thermal_efficiency_full_load", 0.5) == 0.0 - - # If AbsorptionChiller is evaluated, need cooling load -> absorption chiller heating consumption estimate for CHP sizing heuristic - avg_cooling_load_kw = nothing - absorption_chiller_cop = nothing - # User can override by explicitly setting include_cooling_in_chp_size = false - if "include_cooling_in_chp_size" in keys(d["CHP"]) - include_cooling_in_size = pop!(d["CHP"], "include_cooling_in_chp_size") - else - include_cooling_in_size = haskey(d, "AbsorptionChiller") - end - - if max_cooling_demand_kw > 0 && include_cooling_in_size - # Use already-processed cooling_load object - avg_cooling_load_kw = sum(cooling_load.loads_kw_thermal) / length(cooling_load.loads_kw_thermal) - # Get absorption chiller COP if specified, otherwise will use default - if haskey(d, "AbsorptionChiller") && haskey(d["AbsorptionChiller"], "cop_thermal") - absorption_chiller_cop = d["AbsorptionChiller"]["cop_thermal"] + chp_array = isa(d["CHP"], AbstractArray) ? d["CHP"] : [d["CHP"]] + for (i, chp_dict) in enumerate(chp_array) + electric_only = get(chp_dict, "is_electric_only", false) || get(chp_dict, "thermal_efficiency_full_load", 0.5) == 0.0 + + # Set default name if not provided + if !haskey(chp_dict, "name") + chp_dict["name"] = length(chp_array) > 1 ? "CHP$i" : "CHP" + end + + # If AbsorptionChiller is evaluated, include cooling load in CHP sizing heuristic + # User can override by explicitly setting include_cooling_in_chp_size = false in the CHP dict + avg_cooling_load_kw = nothing + absorption_chiller_cop = nothing + if "include_cooling_in_chp_size" in keys(chp_dict) + include_cooling_in_size = pop!(chp_dict, "include_cooling_in_chp_size") + else + include_cooling_in_size = haskey(d, "AbsorptionChiller") + end + + if max_cooling_demand_kw > 0 && include_cooling_in_size + avg_cooling_load_kw = sum(cooling_load.loads_kw_thermal) / length(cooling_load.loads_kw_thermal) + if haskey(d, "AbsorptionChiller") && haskey(d["AbsorptionChiller"], "cop_thermal") + absorption_chiller_cop = d["AbsorptionChiller"]["cop_thermal"] + end + end + + if !isnothing(existing_boiler) && !electric_only + total_fuel_heating_load_mmbtu_per_hour = (space_heating_load.loads_kw + dhw_load.loads_kw + process_heat_load.loads_kw) / existing_boiler.efficiency / KWH_PER_MMBTU + avg_boiler_fuel_load_mmbtu_per_hour = sum(total_fuel_heating_load_mmbtu_per_hour) / length(total_fuel_heating_load_mmbtu_per_hour) + chp = CHP(chp_dict; + avg_boiler_fuel_load_mmbtu_per_hour = avg_boiler_fuel_load_mmbtu_per_hour, + existing_boiler = existing_boiler, + electric_load_series_kw = electric_load.loads_kw, + avg_cooling_load_kw = avg_cooling_load_kw, + absorption_chiller_cop = absorption_chiller_cop, + include_cooling_in_size = include_cooling_in_size, + year = electric_load.year, + sector = site.sector, + 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(chp_dict; + electric_load_series_kw = electric_load.loads_kw, + avg_cooling_load_kw = avg_cooling_load_kw, + absorption_chiller_cop = absorption_chiller_cop, + include_cooling_in_size = include_cooling_in_size, + year = electric_load.year, + sector = site.sector, + federal_procurement_type = site.federal_procurement_type, + off_grid_flag = settings.off_grid_flag) end + push!(chps, chp) end - - if !isnothing(existing_boiler) && !electric_only - total_fuel_heating_load_mmbtu_per_hour = (space_heating_load.loads_kw + dhw_load.loads_kw + process_heat_load.loads_kw) / existing_boiler.efficiency / KWH_PER_MMBTU - avg_boiler_fuel_load_mmbtu_per_hour = sum(total_fuel_heating_load_mmbtu_per_hour) / length(total_fuel_heating_load_mmbtu_per_hour) - chp = CHP(d["CHP"]; - avg_boiler_fuel_load_mmbtu_per_hour = avg_boiler_fuel_load_mmbtu_per_hour, - existing_boiler = existing_boiler, - electric_load_series_kw = electric_load.loads_kw, - avg_cooling_load_kw = avg_cooling_load_kw, - absorption_chiller_cop = absorption_chiller_cop, - include_cooling_in_size = include_cooling_in_size, - year = electric_load.year, - sector = site.sector, - federal_procurement_type = site.federal_procurement_type) - else # Only if modeling CHP without heating_load and existing_boiler (for prime generator, electric-only) - chp = CHP(d["CHP"]; - electric_load_series_kw = electric_load.loads_kw, - avg_cooling_load_kw = avg_cooling_load_kw, - absorption_chiller_cop = absorption_chiller_cop, - include_cooling_in_size = include_cooling_in_size, - year = electric_load.year, - sector = site.sector, - federal_procurement_type = site.federal_procurement_type) + # Store first CHP's prime mover for AbsorptionChiller construction below + if !isempty(chps) + chp_prime_mover = chps[1].prime_mover end end @@ -569,7 +544,7 @@ function Scenario(d::Dict; flex_hvac_from_json=false) if haskey(d, "AbsorptionChiller") absorption_chiller = AbsorptionChiller(d["AbsorptionChiller"]; existing_boiler = existing_boiler, - chp_prime_mover = isnothing(chp) ? nothing : chp.prime_mover, + chp_prime_mover = chp_prime_mover, cooling_load = cooling_load) end end diff --git a/test/runtests.jl b/test/runtests.jl index dae7c3bb6..2ea849578 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -317,7 +317,7 @@ else # run HiGHS tests ) ) s = Scenario(d) - @test s.chp.max_kw ≈ case3_max_size_kw atol=1e-3 + @test s.chps[1].max_kw ≈ case3_max_size_kw atol=1e-3 #Case 5: with higher electric load, still sized according to thermal input avg_electric_load_kw = 2000.0 From 33a84a9b27963d044dab1e070fccaad6ae7d4c68 Mon Sep 17 00:00:00 2001 From: Byron Pullutasig <115118857+bpulluta@users.noreply.github.com> Date: Wed, 3 Jun 2026 14:28:12 -0600 Subject: [PATCH 41/45] added ability to add existing_chp --- src/core/chp.jl | 6 +++-- src/core/reopt_inputs.jl | 9 ++++---- src/results/chp.jl | 4 ++-- test/runtests.jl | 4 ++++ test/test_chp_existing_kw.jl | 44 ++++++++++++++++++++++++++++++++++++ 5 files changed, 59 insertions(+), 8 deletions(-) create mode 100644 test/test_chp_existing_kw.jl diff --git a/src/core/chp.jl b/src/core/chp.jl index 042689bc5..6ea42065f 100644 --- a/src/core/chp.jl +++ b/src/core/chp.jl @@ -25,8 +25,9 @@ conflict_res_min_allowable_fraction_of_max = 0.25 # Optional inputs: size_class::Union{Int, Nothing} = nothing # CHP size class for using appropriate default inputs, with size_class=0 using an average of all other size class data - min_kw::Float64 = 0.0 # Minimum CHP size (based on electric) constraint for optimization - max_kw::Float64 = NaN # Maximum CHP size (based on electric) constraint for optimization. Determined by heuristic sizing based on heating load or electric load. + existing_kw::Float64 = 0.0 # Existing CHP electric capacity (based on rated electric power) + min_kw::Float64 = 0.0 # Minimum NEW CHP size (based on electric) to add beyond existing_kw + max_kw::Float64 = NaN # Maximum NEW CHP size (based on electric) to add beyond existing_kw. Determined by heuristic sizing based on heating load or electric load. fuel_type::String = "natural_gas" # "restrict_to": ["natural_gas", "landfill_bio_gas", "propane", "diesel_oil"] om_cost_per_kw::Float64 = 0.0 # Annual CHP fixed operations and maintenance costs in \$/kw-yr om_cost_per_hr_per_kw_rated::Float64 = 0.0 # CHP non-fuel variable operations and maintenance costs in \$/hr/kw_rated @@ -103,6 +104,7 @@ Base.@kwdef mutable struct CHP <: AbstractCHP # Optional inputs: prime_mover::Union{String, Nothing} = nothing size_class::Union{Int, Nothing} = nothing + existing_kw::Float64 = 0.0 min_kw::Float64 = 0.0 max_kw::Float64 = NaN fuel_type::String = "natural_gas" diff --git a/src/core/reopt_inputs.jl b/src/core/reopt_inputs.jl index 878e29fa6..29e47202c 100644 --- a/src/core/reopt_inputs.jl +++ b/src/core/reopt_inputs.jl @@ -435,7 +435,7 @@ function setup_tech_inputs(s::AbstractScenario, time_steps) end if !isempty(techs.chp) - setup_chp_inputs(s, max_sizes, min_sizes, cap_cost_slope, om_cost_per_kw, + setup_chp_inputs(s, max_sizes, min_sizes, existing_sizes, cap_cost_slope, om_cost_per_kw, production_factor, techs_by_exportbin, techs.segmented, n_segs_by_tech, seg_min_size, seg_max_size, seg_yint, techs, tech_renewable_energy_fraction, tech_emissions_factors_CO2, tech_emissions_factors_NOx, tech_emissions_factors_SO2, tech_emissions_factors_PM25, fuel_cost_per_kwh, @@ -885,15 +885,16 @@ end Update tech-indexed data arrays necessary to build the JuMP model with the values for CHP. """ -function setup_chp_inputs(s::AbstractScenario, max_sizes, min_sizes, cap_cost_slope, om_cost_per_kw, +function setup_chp_inputs(s::AbstractScenario, max_sizes, min_sizes, existing_sizes, cap_cost_slope, om_cost_per_kw, production_factor, techs_by_exportbin, segmented_techs, n_segs_by_tech, seg_min_size, seg_max_size, seg_yint, techs, tech_renewable_energy_fraction, tech_emissions_factors_CO2, tech_emissions_factors_NOx, tech_emissions_factors_SO2, tech_emissions_factors_PM25, fuel_cost_per_kwh, heating_cf, pbi_pwf, pbi_max_benefit, pbi_max_kw, pbi_benefit_per_kwh, chp_params ) for chp in s.chps - max_sizes[chp.name] = chp.max_kw - min_sizes[chp.name] = chp.min_kw + max_sizes[chp.name] = chp.existing_kw + chp.max_kw + min_sizes[chp.name] = chp.existing_kw + chp.min_kw + existing_sizes[chp.name] = chp.existing_kw update_cost_curve!(chp, chp.name, s.financial, cap_cost_slope, segmented_techs, n_segs_by_tech, seg_min_size, seg_max_size, seg_yint ) diff --git a/src/results/chp.jl b/src/results/chp.jl index 340448edc..cb83a743e 100644 --- a/src/results/chp.jl +++ b/src/results/chp.jl @@ -163,8 +163,8 @@ function get_chp_results_for_tech(m::JuMP.AbstractModel, p::REoptInputs, chp_nam r["year_one_standby_cost_after_tax"] = r["year_one_standby_cost_before_tax"] * (1 - p.s.financial.offtaker_tax_rate_fraction) r["lifecycle_standby_cost_after_tax"] = round(r["year_one_standby_cost_before_tax"] * p.pwf_e * (1 - p.s.financial.offtaker_tax_rate_fraction), digits=0) - # Calculate individual CHP capital costs - size_kw = value(m[Symbol("dvSize"*_n)][chp_name]) + # Calculate individual CHP capital costs using purchased CHP capacity (excluding existing_kw) + size_kw = value(m[Symbol("dvPurchaseSize"*_n)][chp_name]) suppl_firing_size = value(m[Symbol("dvSupplementaryFiringSize"*_n)][chp_name]) # Base CHP capital cost using generalized function diff --git a/test/runtests.jl b/test/runtests.jl index 2ea849578..dcab15147 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1416,6 +1416,10 @@ else # run HiGHS tests @testset "CHP Production Factor" begin include("test_chp_prodfactor.jl") end + + @testset "CHP Existing KW" begin + include("test_chp_existing_kw.jl") + end end diff --git a/test/test_chp_existing_kw.jl b/test/test_chp_existing_kw.jl new file mode 100644 index 000000000..47c038416 --- /dev/null +++ b/test/test_chp_existing_kw.jl @@ -0,0 +1,44 @@ +# REopt®, Copyright (c) Alliance for Energy Innovation, LLC. See also https://github.com/NatLabRockies/REopt.jl/blob/master/LICENSE. + +input_base = JSON.parsefile("./scenarios/chp_sizing.json") +input_base["CHP"]["min_kw"] = 0.0 +input_base["CHP"]["max_kw"] = 200.0 + +input_existing = deepcopy(input_base) +input_existing["CHP"]["existing_kw"] = 100.0 + +s_base = Scenario(input_base) +p_base = REoptInputs(s_base) + +s_existing = Scenario(input_existing) +p_existing = REoptInputs(s_existing) + +chp_name = s_existing.chps[1].name + +# Bounds semantics are aligned with Generator/PV: min/max are NEW capacity above existing. +@test p_existing.existing_sizes[chp_name] == 100.0 +@test p_existing.min_sizes[chp_name] == 100.0 +@test p_existing.max_sizes[chp_name] == 300.0 + +m_base = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +r_base = run_reopt(m_base, p_base) + +m_existing = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +r_existing = run_reopt(m_existing, p_existing) + +size_total = r_existing["CHP"]["size_kw"] +new_purchase = size_total - input_existing["CHP"]["existing_kw"] + +# Technical feasibility checks +@test new_purchase <= input_existing["CHP"]["max_kw"] + 1e-3 +@test size_total <= input_existing["CHP"]["existing_kw"] + input_existing["CHP"]["max_kw"] + 1e-3 + +# Financial consistency check: existing CHP should not be charged as new CHP capex +@test r_existing["CHP"]["initial_capital_costs"] < r_base["CHP"]["initial_capital_costs"] +@test r_existing["Financial"]["initial_capital_costs"] < r_base["Financial"]["initial_capital_costs"] + +finalize(backend(m_base)) +empty!(m_base) +finalize(backend(m_existing)) +empty!(m_existing) +GC.gc() From 65aed807a17b695f6f153a9d15ae7aecfcf4002f Mon Sep 17 00:00:00 2001 From: Byron Pullutasig <115118857+bpulluta@users.noreply.github.com> Date: Thu, 4 Jun 2026 13:40:08 -0600 Subject: [PATCH 42/45] Strengthen existing CHP regression with exact capex assertion --- test/test_chp_existing_kw.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/test_chp_existing_kw.jl b/test/test_chp_existing_kw.jl index 47c038416..74e3840a6 100644 --- a/test/test_chp_existing_kw.jl +++ b/test/test_chp_existing_kw.jl @@ -28,12 +28,15 @@ r_existing = run_reopt(m_existing, p_existing) size_total = r_existing["CHP"]["size_kw"] new_purchase = size_total - input_existing["CHP"]["existing_kw"] +expected_existing_chp_capex = REopt.get_tech_initial_capex(s_existing.chps[1], new_purchase) + + s_existing.chps[1].supplementary_firing_capital_cost_per_kw * r_existing["CHP"]["size_supplemental_firing_kw"] # Technical feasibility checks @test new_purchase <= input_existing["CHP"]["max_kw"] + 1e-3 @test size_total <= input_existing["CHP"]["existing_kw"] + input_existing["CHP"]["max_kw"] + 1e-3 # Financial consistency check: existing CHP should not be charged as new CHP capex +@test r_existing["CHP"]["initial_capital_costs"] ≈ expected_existing_chp_capex atol=1e-2 @test r_existing["CHP"]["initial_capital_costs"] < r_base["CHP"]["initial_capital_costs"] @test r_existing["Financial"]["initial_capital_costs"] < r_base["Financial"]["initial_capital_costs"] From 354882c4b865ab5dc1cc876ac054a55a09e2215b Mon Sep 17 00:00:00 2001 From: wbecker Date: Wed, 1 Jul 2026 15:34:07 -0600 Subject: [PATCH 43/45] Fix CHP cooling-size keyword after merge Use include_cooling_in_chp_size when constructing CHP objects in Scenario for multi-CHP inputs. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com> --- src/core/scenario.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/core/scenario.jl b/src/core/scenario.jl index 55110ce07..07db4a2f6 100644 --- a/src/core/scenario.jl +++ b/src/core/scenario.jl @@ -481,12 +481,12 @@ function Scenario(d::Dict; flex_hvac_from_json=false) avg_cooling_load_kw = nothing absorption_chiller_cop = nothing if "include_cooling_in_chp_size" in keys(chp_dict) - include_cooling_in_size = pop!(chp_dict, "include_cooling_in_chp_size") + include_cooling_in_chp_size = pop!(chp_dict, "include_cooling_in_chp_size") else - include_cooling_in_size = haskey(d, "AbsorptionChiller") + include_cooling_in_chp_size = haskey(d, "AbsorptionChiller") end - if max_cooling_demand_kw > 0 && include_cooling_in_size + if max_cooling_demand_kw > 0 && include_cooling_in_chp_size avg_cooling_load_kw = sum(cooling_load.loads_kw_thermal) / length(cooling_load.loads_kw_thermal) if haskey(d, "AbsorptionChiller") && haskey(d["AbsorptionChiller"], "cop_thermal") absorption_chiller_cop = d["AbsorptionChiller"]["cop_thermal"] @@ -502,7 +502,7 @@ function Scenario(d::Dict; flex_hvac_from_json=false) electric_load_series_kw = electric_load.loads_kw, avg_cooling_load_kw = avg_cooling_load_kw, absorption_chiller_cop = absorption_chiller_cop, - include_cooling_in_size = include_cooling_in_size, + include_cooling_in_chp_size = include_cooling_in_chp_size, year = electric_load.year, sector = site.sector, federal_procurement_type = site.federal_procurement_type, @@ -512,7 +512,7 @@ function Scenario(d::Dict; flex_hvac_from_json=false) electric_load_series_kw = electric_load.loads_kw, avg_cooling_load_kw = avg_cooling_load_kw, absorption_chiller_cop = absorption_chiller_cop, - include_cooling_in_size = include_cooling_in_size, + include_cooling_in_chp_size = include_cooling_in_chp_size, year = electric_load.year, sector = site.sector, federal_procurement_type = site.federal_procurement_type, From f7100fe3bc2aee49bca0202ff6f9c77f881c5fb6 Mon Sep 17 00:00:00 2001 From: Byron Pullutasig <115118857+bpulluta@users.noreply.github.com> Date: Wed, 8 Jul 2026 18:34:06 -0400 Subject: [PATCH 44/45] added tests and cleaned existing chp naming convention --- src/core/bau_inputs.jl | 40 +++++++++- src/core/bau_scenario.jl | 17 +++- src/core/techs.jl | 21 +++++ test/test_chp_existing_kw.jl | 150 +++++++++++++++++++++++++++++------ 4 files changed, 197 insertions(+), 31 deletions(-) diff --git a/src/core/bau_inputs.jl b/src/core/bau_inputs.jl index de08b0148..5d222ca98 100644 --- a/src/core/bau_inputs.jl +++ b/src/core/bau_inputs.jl @@ -28,6 +28,7 @@ function BAUInputs(p::REoptInputs) thermal_cop = Dict{String, Float64}() heating_cop = Dict{String, Array{Float64,1}}() cooling_cop = Dict{String, Array{Float64,1}}() + chp_params = Dict{String, Dict{Symbol, Float64}}() heating_cf = Dict{String, Array{Float64,1}}() cooling_cf = Dict{String, Array{Float64,1}}() avoided_capex_by_ashp_present_value = Dict(t => 0.0 for t in techs.all) @@ -92,6 +93,42 @@ function BAUInputs(p::REoptInputs) fuel_cost_per_kwh["Generator"] = p.fuel_cost_per_kwh["Generator"] end + for chp in bau_scenario.chps + chp_name = chp.name + max_sizes[chp_name] = chp.existing_kw + min_sizes[chp_name] = chp.existing_kw + existing_sizes[chp_name] = chp.existing_kw + cap_cost_slope[chp_name] = 0.0 + om_cost_per_kw[chp_name] = chp.om_cost_per_kw + production_factor[chp_name, :] = p.production_factor[chp_name, :] + tech_renewable_energy_fraction[chp_name] = chp.fuel_renewable_energy_fraction + tech_emissions_factors_CO2[chp_name] = chp.emissions_factor_lb_CO2_per_mmbtu / KWH_PER_MMBTU + tech_emissions_factors_NOx[chp_name] = chp.emissions_factor_lb_NOx_per_mmbtu / KWH_PER_MMBTU + tech_emissions_factors_SO2[chp_name] = chp.emissions_factor_lb_SO2_per_mmbtu / KWH_PER_MMBTU + tech_emissions_factors_PM25[chp_name] = chp.emissions_factor_lb_PM25_per_mmbtu / KWH_PER_MMBTU + fillin_techs_by_exportbin(techs_by_exportbin, chp, chp_name) + if chp_name in p.techs.pbi + push!(techs.pbi, chp_name) + end + if !chp.can_curtail + push!(techs.no_curtail, chp_name) + end + chp_fuel_cost_per_kwh = chp.fuel_cost_per_mmbtu ./ KWH_PER_MMBTU + fuel_cost_per_kwh[chp_name] = per_hour_value_to_time_series(chp_fuel_cost_per_kwh, p.s.settings.time_steps_per_hour, chp_name) + heating_cf[chp_name] = ones(8760 * p.s.settings.time_steps_per_hour) + chp_params[chp_name] = Dict{Symbol, Float64}( + :electric_efficiency_full_load => chp.electric_efficiency_full_load, + :electric_efficiency_half_load => chp.electric_efficiency_half_load, + :thermal_efficiency_full_load => chp.thermal_efficiency_full_load, + :thermal_efficiency_half_load => chp.thermal_efficiency_half_load, + :min_turn_down_fraction => chp.min_turn_down_fraction, + :supplementary_firing_efficiency => chp.supplementary_firing_efficiency, + :supplementary_firing_max_steam_ratio => chp.supplementary_firing_max_steam_ratio, + :om_cost_per_kwh => chp.om_cost_per_kwh, + :ramp_rate_fraction_per_hour => chp.ramp_rate_fraction_per_hour + ) + end + if "ExistingBoiler" in techs.all setup_existing_boiler_inputs(bau_scenario, max_sizes, min_sizes, existing_sizes, cap_cost_slope, boiler_efficiency, tech_renewable_energy_fraction, tech_emissions_factors_CO2, tech_emissions_factors_NOx, tech_emissions_factors_SO2, tech_emissions_factors_PM25, fuel_cost_per_kwh, @@ -166,9 +203,6 @@ function BAUInputs(p::REoptInputs) heating_loads_served_by_tes = Dict{String,Array{String,1}}() unavailability = get_unavailability_by_tech(p.s, techs, p.time_steps) - # Initialize CHP-specific parameters as nested dictionary (empty for BAU) - chp_params = Dict{String, Dict{Symbol, Float64}}() - REoptInputs( bau_scenario, techs, diff --git a/src/core/bau_scenario.jl b/src/core/bau_scenario.jl index bd4254c5f..482a35504 100644 --- a/src/core/bau_scenario.jl +++ b/src/core/bau_scenario.jl @@ -34,7 +34,7 @@ struct BAUScenario <: AbstractScenario ghp_option_list::Array{Union{GHP, Nothing}, 1} # List of GHP objects (often just 1 element, but can be more) space_heating_thermal_load_reduction_with_ghp_kw::Union{Vector{Float64}, Nothing} cooling_thermal_load_reduction_with_ghp_kw::Union{Vector{Float64}, Nothing} - chps::Array{CHP, 1} # Empty array for BAU scenarios (no new CHP modeled) + chps::Array{CHP, 1} end @@ -70,8 +70,9 @@ end Constructs the BAUScenario (used to create the Business-as-usual inputs) based on the Scenario for the optimized case. The following assumptions are made for the BAU scenario: -- sets the `PV` and `Generator` min_kw and max_kw values to the existing_kw values +- sets the `PV`, `Generator`, and `CHP` min_kw and max_kw values to the existing_kw values - sets wind and storage max_kw values to zero (existing wind and storage cannot be modeled) +- only includes `PV`, `Generator`, and `CHP` systems that have existing_kw > 0 """ function BAUScenario(s::Scenario) @@ -130,6 +131,16 @@ function BAUScenario(s::Scenario) =# site = bau_site(s.site) + # set CHP.max_kw to existing_kw + chps = CHP[] + for chp in s.chps + if chp.existing_kw > 0 + bau_chp = deepcopy(chp) + bau_chp.supplementary_firing_capital_cost_per_kw = 0.0 + push!(chps, bau_chp) + end + end + return BAUScenario( s.settings, site, @@ -152,6 +163,6 @@ function BAUScenario(s::Scenario) ghp_option_list, space_heating_thermal_load_reduction_with_ghp_kw, cooling_thermal_load_reduction_with_ghp_kw, - CHP[] # Empty array - no CHP in BAU scenario + chps ) end \ No newline at end of file diff --git a/src/core/techs.jl b/src/core/techs.jl index 9026e1787..adf676d9a 100644 --- a/src/core/techs.jl +++ b/src/core/techs.jl @@ -58,6 +58,27 @@ function Techs(p::REoptInputs, s::BAUScenario) push!(electric_chillers, "ExistingChiller") end + for chp in s.chps + push!(all_techs, chp.name) + push!(chp_techs, chp.name) + push!(elec, chp.name) + if !chp.is_electric_only + push!(heating_techs, chp.name) + end + if chp.can_serve_space_heating + push!(techs_can_serve_space_heating, chp.name) + end + if chp.can_serve_dhw + push!(techs_can_serve_dhw, chp.name) + end + if chp.can_serve_process_heat + push!(techs_can_serve_process_heat, chp.name) + end + if chp.can_supply_steam_turbine + push!(techs_can_supply_steam_turbine, chp.name) + end + end + cooling_techs = union(electric_chillers, absorption_chillers) fuel_burning_techs = union(gentechs, boiler_techs, chp_techs) thermal_techs = union(heating_techs, boiler_techs, cooling_techs) diff --git a/test/test_chp_existing_kw.jl b/test/test_chp_existing_kw.jl index 74e3840a6..9dfa6bc7d 100644 --- a/test/test_chp_existing_kw.jl +++ b/test/test_chp_existing_kw.jl @@ -1,4 +1,11 @@ -# REopt®, Copyright (c) Alliance for Energy Innovation, LLC. See also https://github.com/NatLabRockies/REopt.jl/blob/master/LICENSE. +# using Revise +# using REopt +# using JSON +# using Test +# using JuMP +# using HiGHS + +############### Existing CHP Sizing ################### input_base = JSON.parsefile("./scenarios/chp_sizing.json") input_base["CHP"]["min_kw"] = 0.0 @@ -7,41 +14,134 @@ input_base["CHP"]["max_kw"] = 200.0 input_existing = deepcopy(input_base) input_existing["CHP"]["existing_kw"] = 100.0 -s_base = Scenario(input_base) -p_base = REoptInputs(s_base) - s_existing = Scenario(input_existing) p_existing = REoptInputs(s_existing) - chp_name = s_existing.chps[1].name -# Bounds semantics are aligned with Generator/PV: min/max are NEW capacity above existing. +# Bounds: max_kw is NEW capacity; internal max_sizes = existing + max_kw @test p_existing.existing_sizes[chp_name] == 100.0 @test p_existing.min_sizes[chp_name] == 100.0 @test p_existing.max_sizes[chp_name] == 300.0 -m_base = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) -r_base = run_reopt(m_base, p_base) +m1 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +r_base = run_reopt(m1, REoptInputs(Scenario(input_base))) -m_existing = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) -r_existing = run_reopt(m_existing, p_existing) +m2 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +r_existing = run_reopt(m2, p_existing) -size_total = r_existing["CHP"]["size_kw"] -new_purchase = size_total - input_existing["CHP"]["existing_kw"] -expected_existing_chp_capex = REopt.get_tech_initial_capex(s_existing.chps[1], new_purchase) + - s_existing.chps[1].supplementary_firing_capital_cost_per_kw * r_existing["CHP"]["size_supplemental_firing_kw"] +new_purchase = r_existing["CHP"]["size_kw"] - 100.0 +expected_capex = REopt.get_tech_initial_capex(s_existing.chps[1], new_purchase) + + s_existing.chps[1].supplementary_firing_capital_cost_per_kw * r_existing["CHP"]["size_supplemental_firing_kw"] -# Technical feasibility checks -@test new_purchase <= input_existing["CHP"]["max_kw"] + 1e-3 -@test size_total <= input_existing["CHP"]["existing_kw"] + input_existing["CHP"]["max_kw"] + 1e-3 - -# Financial consistency check: existing CHP should not be charged as new CHP capex -@test r_existing["CHP"]["initial_capital_costs"] ≈ expected_existing_chp_capex atol=1e-2 +@test r_existing["CHP"]["initial_capital_costs"] ≈ expected_capex atol=1e-2 @test r_existing["CHP"]["initial_capital_costs"] < r_base["CHP"]["initial_capital_costs"] @test r_existing["Financial"]["initial_capital_costs"] < r_base["Financial"]["initial_capital_costs"] -finalize(backend(m_base)) -empty!(m_base) -finalize(backend(m_existing)) -empty!(m_existing) -GC.gc() +############### BAU Scenario ################### + +bau_s = REopt.BAUScenario(s_existing) +@test length(bau_s.chps) == 1 +@test bau_s.chps[1].existing_kw == 100.0 +@test bau_s.chps[1].supplementary_firing_capital_cost_per_kw == 0.0 + +bau_inputs = REopt.BAUInputs(p_existing) +@test bau_inputs.max_sizes[chp_name] == 100.0 +@test bau_inputs.min_sizes[chp_name] == 100.0 +@test bau_inputs.cap_cost_slope[chp_name] == 0.0 + +# With BAU: existing CHP reduces BAU LCC and yields smaller incremental NPV +m3 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +m4 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +r_full_existing = run_reopt([m3, m4], p_existing) + +m5 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +m6 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +r_full_base = run_reopt([m5, m6], REoptInputs(Scenario(input_base))) + +@test r_full_existing["Financial"]["lcc_bau"] < r_full_base["Financial"]["lcc_bau"] +@test r_full_existing["Financial"]["npv"] < r_full_base["Financial"]["npv"] +@test r_full_existing["Financial"]["npv"] > 0 + +############### Existing Only (No additional NEW CHP) ################### + +input_noexp = JSON.parsefile("./scenarios/chp_sizing.json") +input_noexp["CHP"]["existing_kw"] = 150.0 +input_noexp["CHP"]["min_kw"] = 0.0 +input_noexp["CHP"]["max_kw"] = 0.0 + +m7 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +m8 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +r_noexp = run_reopt([m7, m8], REoptInputs(Scenario(input_noexp))) + +@test r_noexp["CHP"]["size_kw"] ≈ 150.0 atol=0.1 +@test r_noexp["CHP"]["initial_capital_costs"] == 0.0 +@test r_noexp["CHP"]["annual_electric_production_kwh"] > 0 +@test r_noexp["CHP"]["annual_thermal_production_mmbtu"] > 0 +@test r_noexp["CHP"]["year_one_fuel_cost_before_tax"] > 0 +@test abs(r_noexp["Financial"]["npv"]) < 1000.0 + +############### Outage Resilience ################### + +input_outage = JSON.parsefile("./scenarios/chp_sizing.json") +input_outage["CHP"]["existing_kw"] = 150.0 +input_outage["CHP"]["min_kw"] = 0.0 +input_outage["CHP"]["max_kw"] = 0.0 +input_outage["Generator"] = Dict( + "existing_kw" => 100.0, "min_kw" => 0.0, "max_kw" => 0.0, + "installed_cost_per_kw" => 500.0, "om_cost_per_kw" => 10.0, + "fuel_cost_per_gallon" => 3.0, + "electric_efficiency_full_load" => 0.32, "electric_efficiency_half_load" => 0.32, + "fuel_avail_gal" => 1000.0, + "only_runs_during_grid_outage" => true, "sells_energy_back_to_grid" => false +) +input_outage["ElectricUtility"] = Dict( + "outage_start_time_step" => 5000, "outage_end_time_step" => 5048, + "co2_from_avert" => true +) +input_outage["ElectricLoad"] = Dict( + "doe_reference_name" => "Hospital", "annual_kwh" => 1000000.0, + "critical_load_fraction" => 0.5 +) + +m9 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +m10 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +r_outage = run_reopt([m9, m10], REoptInputs(Scenario(input_outage))) + +@test r_outage["CHP"]["size_kw"] ≈ 150.0 atol=0.1 +@test r_outage["Generator"]["size_kw"] ≈ 100.0 atol=0.1 +@test r_outage["CHP"]["initial_capital_costs"] == 0.0 +@test r_outage["CHP"]["annual_electric_production_kwh"] > 0 +@test r_outage["ElectricLoad"]["bau_critical_load_met"] == true +@test abs(r_outage["Financial"]["npv"]) < 5000.0 + +############### Multiple CHPs with Existing ################### + +input_multi = JSON.parsefile("./scenarios/multiple_chps.json") +input_multi["CHP"][1]["existing_kw"] = 150.0 +input_multi["CHP"][1]["min_kw"] = 0.0 +input_multi["CHP"][1]["max_kw"] = 0.0 +input_multi["CHP"][2]["existing_kw"] = 0.0 +input_multi["CHP"][2]["min_kw"] = 100.0 +input_multi["CHP"][2]["max_kw"] = 100.0 + +s_multi = Scenario(input_multi) +p_multi = REoptInputs(s_multi) + +bau_multi = REopt.BAUScenario(s_multi) +@test length(bau_multi.chps) == 1 +@test bau_multi.chps[1].name == "CHP_recip_engine" + +m11 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +m12 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +r_multi = run_reopt([m11, m12], p_multi) + +chp_results = r_multi["CHP"] +recip = first(filter(c -> c["name"] == "CHP_recip_engine", chp_results)) +micro = first(filter(c -> c["name"] == "CHP_micro_turbine", chp_results)) + +@test recip["size_kw"] ≈ 150.0 atol=0.1 +@test recip["initial_capital_costs"] == 0.0 +@test micro["size_kw"] ≈ 100.0 atol=0.1 +@test micro["initial_capital_costs"] > 0.0 +@test recip["annual_electric_production_kwh"] > 0 +@test micro["annual_electric_production_kwh"] > 0 From 73e771516d12b7afa134e3cb89eb5aa92b16f064 Mon Sep 17 00:00:00 2001 From: Byron Pullutasig <115118857+bpulluta@users.noreply.github.com> Date: Thu, 9 Jul 2026 11:40:15 -0400 Subject: [PATCH 45/45] fixed proforma by adding validity/sign change before calling fzero --- src/results/proforma.jl | 21 ++++++++++++++++----- 1 file changed, 16 insertions(+), 5 deletions(-) diff --git a/src/results/proforma.jl b/src/results/proforma.jl index e4d995126..8353a3ae3 100644 --- a/src/results/proforma.jl +++ b/src/results/proforma.jl @@ -522,10 +522,21 @@ function irr(cashflows::AbstractArray{<:Real, 1}) return 0.0 end f(r) = npv(cashflows, r) - rate = 0.0 - try - rate = fzero(f, [0.0, 0.99]) - finally - return round(rate, digits=3) + a, b = 0.0, 0.99 + fa = f(a) + fb = f(b) + + # Only use bracketed solve when endpoint signs differ. + # Otherwise return 0.0 instead of throwing an ArgumentError from Roots. + if !(isfinite(fa) && isfinite(fb)) || fa * fb >= 0 + return 0.0 + end + + rate = try + fzero(f, [a, b]) + catch + return 0.0 end + + return round(rate, digits=3) end