Add existing CHP capacity support

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/operating_reserve_constraints.jl

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

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

@@ -28,6 +28,7 @@ function BAUInputs(p::REoptInputs)
     thermal_cop = Dict{String, Float64}()
     heating_cop = Dict{String, Array{Float64,1}}()
     cooling_cop = Dict{String, Array{Float64,1}}()
+    chp_params = Dict{String, Dict{Symbol, Float64}}()
     heating_cf = Dict{String, Array{Float64,1}}() 
     cooling_cf = Dict{String, Array{Float64,1}}() 
     avoided_capex_by_ashp_present_value = Dict(t => 0.0 for t in techs.all)
@@ -92,6 +93,42 @@ function BAUInputs(p::REoptInputs)
         fuel_cost_per_kwh["Generator"] = p.fuel_cost_per_kwh["Generator"]        
     end
 
+    for chp in bau_scenario.chps
+        chp_name = chp.name
+        max_sizes[chp_name] = chp.existing_kw
+        min_sizes[chp_name] = chp.existing_kw
+        existing_sizes[chp_name] = chp.existing_kw
+        cap_cost_slope[chp_name] = 0.0
+        om_cost_per_kw[chp_name] = chp.om_cost_per_kw
+        production_factor[chp_name, :] = p.production_factor[chp_name, :]
+        tech_renewable_energy_fraction[chp_name] = chp.fuel_renewable_energy_fraction
+        tech_emissions_factors_CO2[chp_name] = chp.emissions_factor_lb_CO2_per_mmbtu / KWH_PER_MMBTU
+        tech_emissions_factors_NOx[chp_name] = chp.emissions_factor_lb_NOx_per_mmbtu / KWH_PER_MMBTU
+        tech_emissions_factors_SO2[chp_name] = chp.emissions_factor_lb_SO2_per_mmbtu / KWH_PER_MMBTU
+        tech_emissions_factors_PM25[chp_name] = chp.emissions_factor_lb_PM25_per_mmbtu / KWH_PER_MMBTU
+        fillin_techs_by_exportbin(techs_by_exportbin, chp, chp_name)
+        if chp_name in p.techs.pbi
+            push!(techs.pbi, chp_name)
+        end
+        if !chp.can_curtail
+            push!(techs.no_curtail, chp_name)
+        end
+        chp_fuel_cost_per_kwh = chp.fuel_cost_per_mmbtu ./ KWH_PER_MMBTU
+        fuel_cost_per_kwh[chp_name] = per_hour_value_to_time_series(chp_fuel_cost_per_kwh, p.s.settings.time_steps_per_hour, chp_name)
+        heating_cf[chp_name] = ones(8760 * p.s.settings.time_steps_per_hour)
+        chp_params[chp_name] = Dict{Symbol, Float64}(
+            :electric_efficiency_full_load => chp.electric_efficiency_full_load,
+            :electric_efficiency_half_load => chp.electric_efficiency_half_load,
+            :thermal_efficiency_full_load => chp.thermal_efficiency_full_load,
+            :thermal_efficiency_half_load => chp.thermal_efficiency_half_load,
+            :min_turn_down_fraction => chp.min_turn_down_fraction,
+            :supplementary_firing_efficiency => chp.supplementary_firing_efficiency,
+            :supplementary_firing_max_steam_ratio => chp.supplementary_firing_max_steam_ratio,
+            :om_cost_per_kwh => chp.om_cost_per_kwh,
+            :ramp_rate_fraction_per_hour => chp.ramp_rate_fraction_per_hour
+        )
+    end
+
     if "ExistingBoiler" in techs.all
         setup_existing_boiler_inputs(bau_scenario, max_sizes, min_sizes, existing_sizes, cap_cost_slope, boiler_efficiency,
             tech_renewable_energy_fraction, tech_emissions_factors_CO2, tech_emissions_factors_NOx, tech_emissions_factors_SO2, tech_emissions_factors_PM25, fuel_cost_per_kwh,
@@ -234,7 +271,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}
 end
 
 
@@ -69,8 +70,9 @@ end
 Constructs the BAUScenario (used to create the Business-as-usual inputs) based on the Scenario for the optimized case.
 
 The following assumptions are made for the BAU scenario: 
-- sets the `PV` and `Generator` min_kw and max_kw values to the existing_kw values
+- sets the `PV`, `Generator`, and `CHP` min_kw and max_kw values to the existing_kw values
 - sets wind and storage max_kw values to zero (existing wind and storage cannot be modeled)
+- only includes `PV`, `Generator`, and `CHP` systems that have existing_kw > 0
 """
 function BAUScenario(s::Scenario)
 
@@ -129,6 +131,16 @@ function BAUScenario(s::Scenario)
     =#
     site = bau_site(s.site)
 
+    # set CHP.max_kw to existing_kw
+    chps = CHP[]
+    for chp in s.chps
+        if chp.existing_kw > 0
+            bau_chp = deepcopy(chp)
+            bau_chp.supplementary_firing_capital_cost_per_kw = 0.0
+            push!(chps, bau_chp)
+        end
+    end
+
     return BAUScenario(
         s.settings,
         site, 
@@ -150,6 +162,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,
+        chps
     )
 end