Multiple CHPs

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
 

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

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
     )
 

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
 

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

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

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
 
 
@@ -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))
@@ -589,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)
@@ -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_chps && !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()