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
     )
   
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.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()
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
 
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
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..0da0976e7 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]))
+    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
         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
 	
diff --git a/test/runtests.jl b/test/runtests.jl
index 1da1eeb2d..7179d1d2d 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
@@ -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
@@ -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
@@ -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))
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
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