diff --git a/src/constraints/chp_constraints.jl b/src/constraints/chp_constraints.jl index 7a94efcb4..833adbc90 100644 --- a/src/constraints/chp_constraints.jl +++ b/src/constraints/chp_constraints.jl @@ -1,129 +1,157 @@ # REopt®, Copyright (c) Alliance for Energy Innovation, LLC. See also https://github.com/NatLabRockies/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 - ) - # 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) + ) + + # 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 + ) + + # 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, 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 + #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 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 + # 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 - # 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) - - #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] + # 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, 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]) - ) + # 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 @@ -135,6 +163,25 @@ function add_chp_rated_prod_constraint(m, p; _n="") end +function add_chp_ramp_rate_constraints(m, p; _n="") + # Ramp rate constraints limit how quickly CHP production can change between consecutive timesteps + # Loop through each CHP and add constraints with tech-specific ramp rate parameters + for t in p.techs.chp + # Ramp up constraint + @constraint(m, [ts in p.time_steps[2:end]], + m[Symbol("dvRatedProduction"*_n)][t, ts] - m[Symbol("dvRatedProduction"*_n)][t, ts-1] <= + p.chp_params[t][:ramp_rate_fraction_per_hour] * m[Symbol("dvSize"*_n)][t] / p.s.settings.time_steps_per_hour + ) + + # Ramp down constraint + @constraint(m, [ts in p.time_steps[2:end]], + m[Symbol("dvRatedProduction"*_n)][t, ts-1] - m[Symbol("dvRatedProduction"*_n)][t, ts] <= + p.chp_params[t][:ramp_rate_fraction_per_hour] * m[Symbol("dvSize"*_n)][t] / p.s.settings.time_steps_per_hour + ) + end +end + + """ add_chp_hourly_om_charges(m, p; _n="") @@ -145,22 +192,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)) @@ -175,7 +228,9 @@ load or send to waste in dispatch. """ function add_chp_to_absorption_chiller_only_constraints(m, p; _n="") monthly_timesteps = get_monthly_time_steps(p.s.electric_load.year; time_steps_per_hour=p.s.settings.time_steps_per_hour) - for mth in p.s.chp.months_serving_absorption_chiller_only + chp_with_ac_only = filter(chp -> chp.serve_absorption_chiller_only, p.s.chps) + months_to_restrict = unique(vcat([chp.months_serving_absorption_chiller_only for chp in chp_with_ac_only]...)) + for mth in months_to_restrict @constraint(m, [t in p.techs.chp, q in [p.s.absorption_chiller.heating_load_input], ts in monthly_timesteps[mth]], m[Symbol("dvProductionToWaste"*_n)]["CHP",q,ts] + m[Symbol("dvHeatToAbsorptionChiller"*_n)]["CHP",q,ts] == m[Symbol("dvHeatingProduction"*_n)]["CHP",q,ts] ) @@ -236,34 +291,78 @@ end Used in src/reopt.jl to add_chp_constraints if !isempty(p.techs.chp) to add CHP operating constraints and cost expressions. """ -function add_chp_constraints(m, p; _n="") - # TODO if chp.min_turn_down_fraction is 0.0, and there is no fuel burn or thermal y-intercept, we don't need the binary below - @warn """Adding binary variable to model CHP. - Some solvers are very slow with integer variables""" - @variables m begin - binCHPIsOnInTS[p.techs.chp, p.time_steps], Bin # 1 If technology t is operating in time step; 0 otherwise - end +function add_chp_constraints(m, p; _n="") + # Check if binary variable is needed by evaluating any CHP's parameters + # Binary is needed if any CHP has non-zero intercepts, min_turn_down, or hourly O&M + binary_needed = false + for t in p.techs.chp + # Calculate fuel burn slopes/intercepts for this 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 + ) + + # Calculate thermal production slopes/intercepts for this CHP + thermal_prod_full_load = 1.0 / p.chp_params[t][:electric_efficiency_full_load] * p.chp_params[t][:thermal_efficiency_full_load] + thermal_prod_half_load = 0.5 / p.chp_params[t][:electric_efficiency_half_load] * p.chp_params[t][:thermal_efficiency_half_load] + thermal_prod_slope = (thermal_prod_full_load - thermal_prod_half_load) / (1.0 - 0.5) + thermal_prod_intercept = thermal_prod_full_load - thermal_prod_slope * 1.0 + + # Check if this CHP needs binary variables + chp_idx = findfirst(chp -> chp.name == t, p.s.chps) + if (abs(fuel_burn_intercept) > 1.0E-7) || + (abs(thermal_prod_intercept) > 1.0E-7) || + (p.chp_params[t][:min_turn_down_fraction] > 1.0E-7) || + (p.s.chps[chp_idx].om_cost_per_hr_per_kw_rated > 1.0E-7) + binary_needed = true + break + end + end + + # Create binary variable if needed + if binary_needed + @warn """Adding binary variable binCHPIsOnInTS to model CHP. + Some solvers are very slow with integer variables""" + @variables m begin + binCHPIsOnInTS[p.techs.chp, p.time_steps], Bin # 1 If technology t is operating in time step; 0 otherwise + end + end m[:TotalHourlyCHPOMCosts] = 0 m[:TotalCHPFuelCosts] = 0 + # 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 - add_chp_hourly_om_charges(m, p) - end - + # These constraints are always needed add_chp_fuel_burn_constraints(m, p; _n=_n) add_chp_thermal_production_constraints(m, p; _n=_n) - add_binCHPIsOnInTS_constraints(m, p; _n=_n) add_chp_rated_prod_constraint(m, p; _n=_n) + + # Add ramp rate constraints if any CHP has ramp_rate_fraction_per_hour < 1.0 + if any(p.chp_params[t][:ramp_rate_fraction_per_hour] < 1.0 / p.s.settings.time_steps_per_hour for t in p.techs.chp) + add_chp_ramp_rate_constraints(m, p; _n=_n) + end - 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 + # These constraints are only needed if binary was created + if binary_needed + add_binCHPIsOnInTS_constraints(m, p; _n=_n) + + # 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; _n=_n) + end + end + + # 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) @@ -271,11 +370,11 @@ function add_chp_constraints(m, p; _n="") end end - if !isempty(p.techs.absorption_chiller) && p.s.chp.serve_absorption_chiller_only + if !isempty(p.techs.absorption_chiller) && any(chp.serve_absorption_chiller_only for chp in p.s.chps) add_chp_to_absorption_chiller_only_constraints(m, p; _n=_n) end - if p.s.chp.follow_electrical_load + if any(chp.follow_electrical_load for chp in p.s.chps) add_chp_electrical_load_following_constraints(m, p; _n=_n) end end diff --git a/src/constraints/cost_curve_constraints.jl b/src/constraints/cost_curve_constraints.jl index 913213357..6ff5ebd11 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 bd04118d5..7c2443fb0 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/operating_reserve_constraints.jl b/src/constraints/operating_reserve_constraints.jl index 568d3c340..e01af7753 100644 --- a/src/constraints/operating_reserve_constraints.jl +++ b/src/constraints/operating_reserve_constraints.jl @@ -8,9 +8,15 @@ function add_operating_reserve_constraints(m, p; _n="") sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for b in p.s.storage.types.elec) - m[Symbol("dvCurtail"*_n)][t, ts] ) + # 1b. Production going to load from requiring_oper_res (separate calculation for techs that require reserves) + m[:ProductionToLoadRequiringOR] = @expression(m, [t in p.techs.requiring_oper_res, ts in p.time_steps_without_grid], + p.production_factor[t, ts] * p.levelization_factor[t] * m[Symbol("dvRatedProduction"*_n)][t,ts] - + sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for b in p.s.storage.types.elec) - + m[Symbol("dvCurtail"*_n)][t, ts] + ) # 2. Total OR required by requiring_oper_res & Load m[:OpResRequired] = @expression(m, [ts in p.time_steps_without_grid], - sum(m[:ProductionToLoadOR][t,ts] * p.techs_operating_reserve_req_fraction[t] for t in p.techs.requiring_oper_res) + sum(m[:ProductionToLoadRequiringOR][t,ts] * p.techs_operating_reserve_req_fraction[t] for t in p.techs.requiring_oper_res) + p.s.electric_load.critical_loads_kw[ts] * m[Symbol("dvOffgridLoadServedFraction"*_n)][ts] * p.s.electric_load.operating_reserve_required_fraction ) # 3. Operating reserve provided - battery @@ -28,13 +34,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) diff --git a/src/constraints/outage_constraints.jl b/src/constraints/outage_constraints.jl index 8d5cb66de..e7c50bd7d 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 a94a3a91d..5d222ca98 100644 --- a/src/core/bau_inputs.jl +++ b/src/core/bau_inputs.jl @@ -28,6 +28,7 @@ function BAUInputs(p::REoptInputs) thermal_cop = Dict{String, Float64}() heating_cop = Dict{String, Array{Float64,1}}() cooling_cop = Dict{String, Array{Float64,1}}() + chp_params = Dict{String, Dict{Symbol, Float64}}() heating_cf = Dict{String, Array{Float64,1}}() cooling_cf = Dict{String, Array{Float64,1}}() avoided_capex_by_ashp_present_value = Dict(t => 0.0 for t in techs.all) @@ -92,6 +93,42 @@ function BAUInputs(p::REoptInputs) fuel_cost_per_kwh["Generator"] = p.fuel_cost_per_kwh["Generator"] end + for chp in bau_scenario.chps + chp_name = chp.name + max_sizes[chp_name] = chp.existing_kw + min_sizes[chp_name] = chp.existing_kw + existing_sizes[chp_name] = chp.existing_kw + cap_cost_slope[chp_name] = 0.0 + om_cost_per_kw[chp_name] = chp.om_cost_per_kw + production_factor[chp_name, :] = p.production_factor[chp_name, :] + tech_renewable_energy_fraction[chp_name] = chp.fuel_renewable_energy_fraction + tech_emissions_factors_CO2[chp_name] = chp.emissions_factor_lb_CO2_per_mmbtu / KWH_PER_MMBTU + tech_emissions_factors_NOx[chp_name] = chp.emissions_factor_lb_NOx_per_mmbtu / KWH_PER_MMBTU + tech_emissions_factors_SO2[chp_name] = chp.emissions_factor_lb_SO2_per_mmbtu / KWH_PER_MMBTU + tech_emissions_factors_PM25[chp_name] = chp.emissions_factor_lb_PM25_per_mmbtu / KWH_PER_MMBTU + fillin_techs_by_exportbin(techs_by_exportbin, chp, chp_name) + if chp_name in p.techs.pbi + push!(techs.pbi, chp_name) + end + if !chp.can_curtail + push!(techs.no_curtail, chp_name) + end + chp_fuel_cost_per_kwh = chp.fuel_cost_per_mmbtu ./ KWH_PER_MMBTU + fuel_cost_per_kwh[chp_name] = per_hour_value_to_time_series(chp_fuel_cost_per_kwh, p.s.settings.time_steps_per_hour, chp_name) + heating_cf[chp_name] = ones(8760 * p.s.settings.time_steps_per_hour) + chp_params[chp_name] = Dict{Symbol, Float64}( + :electric_efficiency_full_load => chp.electric_efficiency_full_load, + :electric_efficiency_half_load => chp.electric_efficiency_half_load, + :thermal_efficiency_full_load => chp.thermal_efficiency_full_load, + :thermal_efficiency_half_load => chp.thermal_efficiency_half_load, + :min_turn_down_fraction => chp.min_turn_down_fraction, + :supplementary_firing_efficiency => chp.supplementary_firing_efficiency, + :supplementary_firing_max_steam_ratio => chp.supplementary_firing_max_steam_ratio, + :om_cost_per_kwh => chp.om_cost_per_kwh, + :ramp_rate_fraction_per_hour => chp.ramp_rate_fraction_per_hour + ) + end + if "ExistingBoiler" in techs.all setup_existing_boiler_inputs(bau_scenario, max_sizes, min_sizes, existing_sizes, cap_cost_slope, boiler_efficiency, tech_renewable_energy_fraction, tech_emissions_factors_CO2, tech_emissions_factors_NOx, tech_emissions_factors_SO2, tech_emissions_factors_PM25, fuel_cost_per_kwh, @@ -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 diff --git a/src/core/bau_scenario.jl b/src/core/bau_scenario.jl index aa8b07035..482a35504 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} 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 \ No newline at end of file diff --git a/src/core/chp.jl b/src/core/chp.jl index 07ff763a4..99805e285 100644 --- a/src/core/chp.jl +++ b/src/core/chp.jl @@ -25,11 +25,13 @@ conflict_res_min_allowable_fraction_of_max = 0.25 # Optional inputs: size_class::Union{Int, Nothing} = nothing # CHP size class for using appropriate default inputs, with size_class=0 using an average of all other size class data - min_kw::Float64 = 0.0 # Minimum CHP size (based on electric) constraint for optimization - max_kw::Float64 = NaN # Maximum CHP size (based on electric) constraint for optimization. Determined by heuristic sizing based on heating load or electric load. + existing_kw::Float64 = 0.0 # Existing CHP electric capacity (based on rated electric power) + min_kw::Float64 = 0.0 # Minimum NEW CHP size (based on electric) to add beyond existing_kw + max_kw::Float64 = NaN # Maximum NEW CHP size (based on electric) to add beyond existing_kw. Determined by heuristic sizing based on heating load or electric load. fuel_type::String = "natural_gas" # "restrict_to": ["natural_gas", "landfill_bio_gas", "propane", "diesel_oil"] om_cost_per_kw::Float64 = 0.0 # Annual CHP fixed operations and maintenance costs in \$/kw-yr om_cost_per_hr_per_kw_rated::Float64 = 0.0 # CHP non-fuel variable operations and maintenance costs in \$/hr/kw_rated + ramp_rate_fraction_per_hour::Float64 = 1.0 # Maximum rate of change in electric production per hour as a fraction of size_kw [kW/size_kw/hour]. supplementary_firing_capital_cost_per_kw::Float64 = 150.0 # Installed CHP supplementary firing system cost in \$/kW (based on rated electric power) supplementary_firing_max_steam_ratio::Float64 = 1.0 # Ratio of max fired steam to un-fired steam production. Relevant only for combustion_turbine prime_mover supplementary_firing_efficiency::Float64 = 0.92 # Thermal efficiency of the incremental steam production from supplementary firing. Relevant only for combustion_turbine prime_mover @@ -40,6 +42,8 @@ conflict_res_min_allowable_fraction_of_max = 0.25 can_serve_space_heating::Bool = true # If CHP can supply heat to the space heating load can_serve_process_heat::Bool = true # If CHP can supply heat to the process heating load is_electric_only::Bool = false # If CHP is a prime generator that does not supply heat + operating_reserve_required_fraction::Real = off_grid_flag ? 0.0 : 0.0 # If off grid, 10%, else 0%. Applied to each time_step as a % of CHP generation. When positive, CHP requires reserves from Generator or Battery; when 0, CHP can provide reserves. + production_factor_series::Union{Nothing, Vector{<:Real}} = nothing # Optional user-provided production factor time-series (length of time_steps_per_hour * 8760). If provided, this will override the production factor calculated from unavailability. serve_absorption_chiller_only::Bool = false # If CHP produced heat either serves absorption chiller or sends it to waste; only applies to the months specified in months_serving_absorption_chiller_only if true months_serving_absorption_chiller_only::AbstractVector{Int64} = Int64[] # months in which CHP only sevres the absorption chiller, with 1=January and 12=December; only applied when serve_absorption_chiller_only = true follow_electrical_load::Bool = false # If CHP follows the electrical load by running at capacity or meeting the load only. @@ -83,6 +87,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[] @@ -99,11 +104,13 @@ Base.@kwdef mutable struct CHP <: AbstractCHP # Optional inputs: prime_mover::Union{String, Nothing} = nothing size_class::Union{Int, Nothing} = nothing + existing_kw::Float64 = 0.0 min_kw::Float64 = 0.0 max_kw::Float64 = NaN fuel_type::String = "natural_gas" om_cost_per_kw::Float64 = 0.0 om_cost_per_hr_per_kw_rated::Float64 = 0.0 + ramp_rate_fraction_per_hour::Float64 = 1.0 electric_efficiency_half_load::Float64 = NaN # Assigned to electric_efficiency_full_load if not input thermal_efficiency_half_load::Float64 = NaN # Assigned to thermal_efficiency_full_load if not input supplementary_firing_capital_cost_per_kw::Float64 = 150.0 @@ -116,6 +123,8 @@ Base.@kwdef mutable struct CHP <: AbstractCHP can_serve_space_heating::Bool = true can_serve_process_heat::Bool = true is_electric_only::Bool = false + operating_reserve_required_fraction::Real = 0.0 + production_factor_series::Union{Nothing, Vector{<:Real}} = nothing serve_absorption_chiller_only::Bool = false months_serving_absorption_chiller_only::AbstractVector{Int64} = Int64[] follow_electrical_load::Bool = false @@ -146,6 +155,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 @@ -158,7 +168,8 @@ function CHP(d::Dict; include_cooling_in_chp_size::Bool=false, year::Int64=2017, sector::String, - federal_procurement_type::String) + federal_procurement_type::String, + off_grid_flag::Bool=false) # If array inputs are coming from Julia JSON.parsefile (reader), they have type Vector{Any}; convert to expected type here for (k,v) in d if typeof(v) <: AbstractVector{Any} && k != "unavailability_periods" @@ -287,6 +298,17 @@ function CHP(d::Dict; setproperty!(chp, :thermal_efficiency_half_load, 0.0) end + # Set operating_reserve_required_fraction based on off_grid_flag + if haskey(d, "operating_reserve_required_fraction") + chp.operating_reserve_required_fraction = d["operating_reserve_required_fraction"] + end + + # Validate operating_reserve_required_fraction for on-grid scenarios + if !off_grid_flag && !(chp.operating_reserve_required_fraction == 0.0) + @warn "CHP operating_reserve_required_fraction applies only when off_grid_flag is true. Setting operating_reserve_required_fraction to 0.0 for this on-grid analysis." + chp.operating_reserve_required_fraction = 0.0 + end + if chp.serve_absorption_chiller_only && isempty(chp.months_serving_absorption_chiller_only) @warn "CHP.serve_absorption_chiller_only is set to true, but no months are specified. All months will be enforced." chp.months_serving_absorption_chiller_only = [1,2,3,4,5,6,7,8,9,10,11,12] @@ -309,7 +331,7 @@ function CHP(d::Dict; end end end - + return chp end @@ -620,3 +642,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/production_factor.jl b/src/core/production_factor.jl index 8bce6219f..6da9f0252 100644 --- a/src/core/production_factor.jl +++ b/src/core/production_factor.jl @@ -278,9 +278,15 @@ end production_factor for CHP accounts for unavailability (`unavailability_periods`) of CHP due to scheduled (mostly off-peak) and "unscheduled" (on-peak) maintenance. Note: this same prod_factor should be applied to electric and thermal production + +If the user provides their own production_factor_series, that will be returned instead of generating from unavailability. """ function get_production_factor(chp::AbstractCHP, year::Int=2017, outage_start_time_step::Int=0, outage_end_time_step::Int=0, ts_per_hour::Int=1) + if !(isnothing(chp.production_factor_series)) + return chp.production_factor_series + end + prod_factor = [1.0 - chp.unavailability_hourly[i] for i in 1:8760 for _ in 1:ts_per_hour] # Ignore unavailability in time_step if it intersects with an outage interval(s) diff --git a/src/core/reopt.jl b/src/core/reopt.jl index 83a1cce32..258cadab0 100644 --- a/src/core/reopt.jl +++ b/src/core/reopt.jl @@ -87,8 +87,8 @@ end Method for use with Threads when running BAU in parallel with optimal scenario. """ function run_reopt(t::Tuple{JuMP.AbstractModel, AbstractInputs}) - run_reopt(t[1], t[2]; organize_pvs=false) - # must organize_pvs after adding proforma results + run_reopt(t[1], t[2]; organize_pvs=false, organize_chps=false) +# must organize_pvs/chps after adding proforma results end @@ -151,6 +151,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.")) @@ -283,11 +286,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)) @@ -588,7 +597,7 @@ function build_reopt!(m::JuMP.AbstractModel, p::REoptInputs) end -function run_reopt(m::JuMP.AbstractModel, p::REoptInputs; organize_pvs=true) +function run_reopt(m::JuMP.AbstractModel, p::REoptInputs; organize_pvs=true, organize_chps=true) try build_reopt!(m, p) @@ -619,6 +628,9 @@ function run_reopt(m::JuMP.AbstractModel, p::REoptInputs; organize_pvs=true) if organize_pvs && !isempty(p.techs.pv) # do not want to organize_pvs when running BAU case in parallel b/c then proform code fails organize_multiple_pv_results(p, results) end + if organize_chps && !isempty(p.techs.chp) # same logic as PV + organize_multiple_chp_results(p, results) + end # add error messages (if any) and warnings to results dict results["Messages"] = logger_to_dict() diff --git a/src/core/reopt_inputs.jl b/src/core/reopt_inputs.jl index d29961ee4..29e47202c 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 @@ -336,7 +339,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 @@ -374,6 +378,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}() @@ -427,14 +434,15 @@ function setup_tech_inputs(s::AbstractScenario, time_steps) om_cost_per_kw, production_factor, fuel_cost_per_kwh, heating_cf) end - if "CHP" in techs.all - setup_chp_inputs(s, max_sizes, min_sizes, cap_cost_slope, om_cost_per_kw, + if !isempty(techs.chp) + setup_chp_inputs(s, max_sizes, min_sizes, existing_sizes, cap_cost_slope, om_cost_per_kw, production_factor, techs_by_exportbin, techs.segmented, n_segs_by_tech, seg_min_size, seg_max_size, seg_yint, techs, tech_renewable_energy_fraction, tech_emissions_factors_CO2, tech_emissions_factors_NOx, tech_emissions_factors_SO2, tech_emissions_factors_PM25, fuel_cost_per_kwh, - 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 @@ -506,7 +514,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 @@ -648,43 +657,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, @@ -879,54 +851,83 @@ 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 maximum cooling_thermal_factor from all CHPs that have thermal output (not electric-only) + chps_with_thermal = filter(chp -> !chp.is_electric_only, s.chps) + if isempty(chps_with_thermal) + thermal_factor = 1.0 + else + thermal_factor = maximum(chp.cooling_thermal_factor for chp in chps_with_thermal) + if thermal_factor == 0.0 + throw(@error("All CHPs have cooling_thermal_factor of 0.0 which implies that CHP cannot serve AbsorptionChiller. If you + want to model CHP and AbsorptionChiller, you must specify a cooling_thermal_factor greater than 0.0 for at least one CHP")) + end + 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. """ -function setup_chp_inputs(s::AbstractScenario, max_sizes, min_sizes, cap_cost_slope, om_cost_per_kw, +function setup_chp_inputs(s::AbstractScenario, max_sizes, min_sizes, existing_sizes, cap_cost_slope, om_cost_per_kw, production_factor, techs_by_exportbin, segmented_techs, n_segs_by_tech, seg_min_size, seg_max_size, seg_yint, techs, tech_renewable_energy_fraction, tech_emissions_factors_CO2, tech_emissions_factors_NOx, tech_emissions_factors_SO2, tech_emissions_factors_PM25, fuel_cost_per_kwh, - heating_cf, pbi_pwf, pbi_max_benefit, pbi_max_kw, pbi_benefit_per_kwh - ) - 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 + heating_cf, pbi_pwf, pbi_max_benefit, pbi_max_kw, pbi_benefit_per_kwh, + chp_params ) - 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.existing_kw + chp.max_kw + min_sizes[chp.name] = chp.existing_kw + chp.min_kw + existing_sizes[chp.name] = chp.existing_kw + update_cost_curve!(chp, chp.name, s.financial, + cap_cost_slope, segmented_techs, n_segs_by_tech, seg_min_size, seg_max_size, seg_yint + ) + 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, + :ramp_rate_fraction_per_hour => chp.ramp_rate_fraction_per_hour + ) + end return nothing end @@ -1138,10 +1139,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 @@ -1369,13 +1376,17 @@ function setup_ghp_inputs(s::AbstractScenario, time_steps, time_steps_without_gr end function setup_operating_reserve_fraction(s::AbstractScenario, techs_operating_reserve_req_fraction) - # currently only PV and Wind require operating reserves + # Currently PV, Wind, and CHP can require or provide operating reserves in off-grid scenarios for pv in s.pvs techs_operating_reserve_req_fraction[pv.name] = pv.operating_reserve_required_fraction end techs_operating_reserve_req_fraction["Wind"] = s.wind.operating_reserve_required_fraction + for chp in s.chps + techs_operating_reserve_req_fraction[chp.name] = chp.operating_reserve_required_fraction + end + return nothing end @@ -1405,7 +1416,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 609caf69a..07db4a2f6 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} @@ -82,7 +82,7 @@ function Scenario(d::Dict; flex_hvac_from_json=false) # Check that only PV, electric storage, and generator are modeled for off-grid if settings.off_grid_flag - offgrid_allowed_keys = ["PV", "Wind", "ElectricStorage", "Generator", "Settings", "Site", "Financial", "ElectricLoad", "ElectricTariff", "ElectricUtility"] + offgrid_allowed_keys = ["PV", "Wind", "ElectricStorage", "Generator", "CHP", "Settings", "Site", "Financial", "ElectricLoad", "ElectricTariff", "ElectricUtility"] unallowed_keys = setdiff(keys(d), offgrid_allowed_keys) if !isempty(unallowed_keys) throw(@error("The following key(s) are not permitted when `off_grid_flag` is true: $unallowed_keys.")) @@ -464,51 +464,65 @@ function Scenario(d::Dict; flex_hvac_from_json=false) end - chp = nothing + chps = CHP[] + chp_prime_mover = nothing if haskey(d, "CHP") - electric_only = get(d["CHP"], "is_electric_only", false) || get(d["CHP"], "thermal_efficiency_full_load", 0.5) == 0.0 - - # If AbsorptionChiller is evaluated, need cooling load -> absorption chiller heating consumption estimate for CHP sizing heuristic - avg_cooling_load_kw = nothing - absorption_chiller_cop = nothing - # User can override by explicitly setting include_cooling_in_chp_size = false - if "include_cooling_in_chp_size" in keys(d["CHP"]) - include_cooling_in_chp_size = pop!(d["CHP"], "include_cooling_in_chp_size") - else - include_cooling_in_chp_size = haskey(d, "AbsorptionChiller") - end - - if max_cooling_demand_kw > 0 && include_cooling_in_chp_size - # Use already-processed cooling_load object - avg_cooling_load_kw = sum(cooling_load.loads_kw_thermal) / length(cooling_load.loads_kw_thermal) - # Get absorption chiller COP if specified, otherwise will use default - if haskey(d, "AbsorptionChiller") && haskey(d["AbsorptionChiller"], "cop_thermal") - absorption_chiller_cop = d["AbsorptionChiller"]["cop_thermal"] + chp_array = isa(d["CHP"], AbstractArray) ? d["CHP"] : [d["CHP"]] + for (i, chp_dict) in enumerate(chp_array) + electric_only = get(chp_dict, "is_electric_only", false) || get(chp_dict, "thermal_efficiency_full_load", 0.5) == 0.0 + + # Set default name if not provided + if !haskey(chp_dict, "name") + chp_dict["name"] = length(chp_array) > 1 ? "CHP$i" : "CHP" + end + + # If AbsorptionChiller is evaluated, include cooling load in CHP sizing heuristic + # User can override by explicitly setting include_cooling_in_chp_size = false in the CHP dict + avg_cooling_load_kw = nothing + absorption_chiller_cop = nothing + if "include_cooling_in_chp_size" in keys(chp_dict) + include_cooling_in_chp_size = pop!(chp_dict, "include_cooling_in_chp_size") + else + include_cooling_in_chp_size = haskey(d, "AbsorptionChiller") + end + + if max_cooling_demand_kw > 0 && include_cooling_in_chp_size + avg_cooling_load_kw = sum(cooling_load.loads_kw_thermal) / length(cooling_load.loads_kw_thermal) + if haskey(d, "AbsorptionChiller") && haskey(d["AbsorptionChiller"], "cop_thermal") + absorption_chiller_cop = d["AbsorptionChiller"]["cop_thermal"] + end end + + if !isnothing(existing_boiler) && !electric_only + total_fuel_heating_load_mmbtu_per_hour = (space_heating_load.loads_kw + dhw_load.loads_kw + process_heat_load.loads_kw) / existing_boiler.efficiency / KWH_PER_MMBTU + avg_boiler_fuel_load_mmbtu_per_hour = sum(total_fuel_heating_load_mmbtu_per_hour) / length(total_fuel_heating_load_mmbtu_per_hour) + chp = CHP(chp_dict; + avg_boiler_fuel_load_mmbtu_per_hour = avg_boiler_fuel_load_mmbtu_per_hour, + existing_boiler = existing_boiler, + electric_load_series_kw = electric_load.loads_kw, + avg_cooling_load_kw = avg_cooling_load_kw, + absorption_chiller_cop = absorption_chiller_cop, + include_cooling_in_chp_size = include_cooling_in_chp_size, + year = electric_load.year, + sector = site.sector, + federal_procurement_type = site.federal_procurement_type, + off_grid_flag = settings.off_grid_flag) + else # Only if modeling CHP without heating_load and existing_boiler (for prime generator, electric-only) + chp = CHP(chp_dict; + electric_load_series_kw = electric_load.loads_kw, + avg_cooling_load_kw = avg_cooling_load_kw, + absorption_chiller_cop = absorption_chiller_cop, + include_cooling_in_chp_size = include_cooling_in_chp_size, + year = electric_load.year, + sector = site.sector, + federal_procurement_type = site.federal_procurement_type, + off_grid_flag = settings.off_grid_flag) + end + push!(chps, chp) end - - if !isnothing(existing_boiler) && !electric_only - total_fuel_heating_load_mmbtu_per_hour = (space_heating_load.loads_kw + dhw_load.loads_kw + process_heat_load.loads_kw) / existing_boiler.efficiency / KWH_PER_MMBTU - avg_boiler_fuel_load_mmbtu_per_hour = sum(total_fuel_heating_load_mmbtu_per_hour) / length(total_fuel_heating_load_mmbtu_per_hour) - chp = CHP(d["CHP"]; - avg_boiler_fuel_load_mmbtu_per_hour = avg_boiler_fuel_load_mmbtu_per_hour, - existing_boiler = existing_boiler, - electric_load_series_kw = electric_load.loads_kw, - avg_cooling_load_kw = avg_cooling_load_kw, - absorption_chiller_cop = absorption_chiller_cop, - include_cooling_in_chp_size = include_cooling_in_chp_size, - year = electric_load.year, - sector = site.sector, - federal_procurement_type = site.federal_procurement_type) - else # Only if modeling CHP without heating_load and existing_boiler (for prime generator, electric-only, or only sending heat to absorption chiller) - chp = CHP(d["CHP"]; - electric_load_series_kw = electric_load.loads_kw, - avg_cooling_load_kw = avg_cooling_load_kw, - absorption_chiller_cop = absorption_chiller_cop, - include_cooling_in_chp_size = include_cooling_in_chp_size, - year = electric_load.year, - sector = site.sector, - federal_procurement_type = site.federal_procurement_type) + # Store first CHP's prime mover for AbsorptionChiller construction below + if !isempty(chps) + chp_prime_mover = chps[1].prime_mover end end @@ -529,7 +543,7 @@ function Scenario(d::Dict; flex_hvac_from_json=false) if haskey(d, "AbsorptionChiller") absorption_chiller = AbsorptionChiller(d["AbsorptionChiller"]; existing_boiler = existing_boiler, - chp_prime_mover = isnothing(chp) ? nothing : chp.prime_mover, + chp_prime_mover = chp_prime_mover, cooling_load = cooling_load) end end @@ -1055,7 +1069,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 6cd303aee..adf676d9a 100644 --- a/src/core/techs.jl +++ b/src/core/techs.jl @@ -58,6 +58,27 @@ function Techs(p::REoptInputs, s::BAUScenario) push!(electric_chillers, "ExistingChiller") end + for chp in s.chps + push!(all_techs, chp.name) + push!(chp_techs, chp.name) + push!(elec, chp.name) + if !chp.is_electric_only + push!(heating_techs, chp.name) + end + if chp.can_serve_space_heating + push!(techs_can_serve_space_heating, chp.name) + end + if chp.can_serve_dhw + push!(techs_can_serve_dhw, chp.name) + end + if chp.can_serve_process_heat + push!(techs_can_serve_process_heat, chp.name) + end + if chp.can_supply_steam_turbine + push!(techs_can_supply_steam_turbine, chp.name) + end + end + cooling_techs = union(electric_chillers, absorption_chillers) fuel_burning_techs = union(gentechs, boiler_techs, chp_techs) thermal_techs = union(heating_techs, boiler_techs, cooling_techs) @@ -186,21 +207,28 @@ 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 + if s.settings.off_grid_flag + if chp.operating_reserve_required_fraction > 0.0 + push!(requiring_oper_res, chp.name) + else + push!(providing_oper_res, chp.name) + end end end diff --git a/src/core/utils.jl b/src/core/utils.jl index 7ffb793a2..5477c19b2 100644 --- a/src/core/utils.jl +++ b/src/core/utils.jl @@ -222,6 +222,17 @@ function dictkeys_tosymbols(d::Dict) end end end + # Convert numeric boolean values (0.0/1.0) to proper Bool type + if k in [ + "off_grid_flag", "add_soc_incentive", "include_climate_in_objective", + "include_health_in_objective", "include_export_cost_series_in_results" + ] && !isnothing(v) && !(typeof(v) <: Bool) + try + v = Bool(v) + catch + throw(@error("Unable to convert $k to Bool. Expected boolean or 0/1, got: $v")) + end + end if k in [ "fuel_cost_per_mmbtu", "wholesale_rate", "export_rate_beyond_net_metering_limit", # for ERP diff --git a/src/mpc/scenario.jl b/src/mpc/scenario.jl index fbf95a504..f84f8f541 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 15ef48e26..fd5d6f5ab 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/boiler.jl b/src/results/boiler.jl index 84beea330..394cbeca2 100644 --- a/src/results/boiler.jl +++ b/src/results/boiler.jl @@ -28,8 +28,7 @@ function add_boiler_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n=" r["size_mmbtu_per_hour"] = round(value(m[Symbol("dvSize"*_n)]["Boiler"]) / KWH_PER_MMBTU, digits=3) r["fuel_consumption_series_mmbtu_per_hour"] = round.(value.(m[:dvFuelUsage]["Boiler", ts] for ts in p.time_steps) / KWH_PER_MMBTU, digits=3) - r["annual_fuel_consumption_mmbtu"] = round(sum(r["fuel_consumption_series_mmbtu_per_hour"]), digits=3) - + r["annual_fuel_consumption_mmbtu"] = round(sum(r["fuel_consumption_series_mmbtu_per_hour"]) / p.s.settings.time_steps_per_hour, digits=3) if !isempty(p.s.storage.types.hot) @expression(m, NewBoilerToHotTESKW[ts in p.time_steps], @@ -49,7 +48,7 @@ function add_boiler_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n=" NewBoilerToSteamTurbine = zeros(length(p.time_steps)) @expression(m, NewBoilerToSteamTurbineByQuality[q in p.heating_loads, ts in p.time_steps], 0.0) end - r["thermal_to_steamturbine_series_mmbtu_per_hour"] = round.(value.(NewBoilerToSteamTurbine), digits=3) + r["thermal_to_steamturbine_series_mmbtu_per_hour"] = round.(value.(NewBoilerToSteamTurbine) ./ KWH_PER_MMBTU, digits=3) if "AbsorptionChiller" in p.techs.cooling @expression(m, NewBoilertoAbsorptionChillerKW[ts in p.time_steps], sum(m[:dvHeatToAbsorptionChiller]["Boiler",q,ts] for q in p.heating_loads)) diff --git a/src/results/chp.jl b/src/results/chp.jl index 40c1dff7e..89ce82790 100644 --- a/src/results/chp.jl +++ b/src/results/chp.jl @@ -35,114 +35,166 @@ 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, 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) + 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.(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) + + 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) + + # Check if this specific CHP can curtail + if chp.can_curtail + CHPtoCurtail = [value(m[Symbol("dvCurtail"*_n)][chp_name,ts]) for ts in p.time_steps] + else + CHPtoCurtail = zeros(length(p.time_steps)) + end + r["electric_curtailed_series_kw"] = round.(CHPtoCurtail, 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] - CHPtoCurtail[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) + r["thermal_to_steamturbine_series_mmbtu_per_hour"] = round.(CHPToSteamTurbineKW / KWH_PER_MMBTU, digits=5) if "AbsorptionChiller" in p.techs.cooling - @expression(m, CHPtoAbsorptionChillerKW[ts in p.time_steps], sum(m[:dvHeatToAbsorptionChiller][t,q,ts] for t in p.techs.chp, q in p.heating_loads)) - @expression(m, CHPtoAbsorptionChillerByQualityKW[q in p.heating_loads, ts in p.time_steps], sum(m[:dvHeatToAbsorptionChiller][t,q,ts] for t in p.techs.chp)) + CHPtoAbsorptionChillerKW = [sum(value(m[Symbol("dvHeatToAbsorptionChiller"*_n)][chp_name,q,ts]) for q in p.heating_loads) for ts in p.time_steps] + CHPtoAbsorptionChillerByQualityKW = Dict(q => [value(m[Symbol("dvHeatToAbsorptionChiller"*_n)][chp_name,q,ts]) for ts in p.time_steps] for q in p.heating_loads) else - @expression(m, CHPtoAbsorptionChillerKW[ts in p.time_steps], 0.0) - @expression(m, CHPtoAbsorptionChillerByQualityKW[q in p.heating_loads, ts in p.time_steps], 0.0) + CHPtoAbsorptionChillerKW = zeros(length(p.time_steps)) + CHPtoAbsorptionChillerByQualityKW = Dict(q => zeros(length(p.time_steps)) for q in p.heating_loads) end - r["thermal_to_absorption_chiller_series_mmbtu_per_hour"] = round.(value.(CHPtoAbsorptionChillerKW) / KWH_PER_MMBTU, digits=5) - - if "DomesticHotWater" in p.heating_loads && p.s.chp.can_serve_dhw && sum(p.heating_loads_kw["DomesticHotWater"]) > 0.0 - @expression(m, CHPToDHWKW[ts in p.time_steps], - m[:dvHeatingProduction]["CHP","DomesticHotWater",ts] - CHPToHotTESByQuality["DomesticHotWater",ts] - CHPToSteamTurbineByQualityKW["DomesticHotWater",ts] - - CHPThermalToWasteByQualityKW["DomesticHotWater",ts] - CHPtoAbsorptionChillerByQualityKW["DomesticHotWater",ts] - ) + r["thermal_to_absorption_chiller_series_mmbtu_per_hour"] = round.(CHPtoAbsorptionChillerKW / 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] - CHPtoAbsorptionChillerKW[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 && chp.can_serve_dhw && sum(p.heating_loads_kw["DomesticHotWater"]) > 0.0 + CHPToDHWKW = [value(m[:dvHeatingProduction][chp_name,"DomesticHotWater",ts]) - CHPToHotTESByQuality["DomesticHotWater"][ts] - CHPToSteamTurbineByQualityKW["DomesticHotWater"][ts] - CHPThermalToWasteByQualityKW["DomesticHotWater"][ts] - CHPtoAbsorptionChillerByQualityKW["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.(max.(0.0, value.(CHPToDHWKW ./ KWH_PER_MMBTU)), digits=5) + r["thermal_to_dhw_load_series_mmbtu_per_hour"] = round.(max.(0.0, CHPToDHWKW ./ KWH_PER_MMBTU), digits=5) - if "SpaceHeating" in p.heating_loads && p.s.chp.can_serve_space_heating && sum(p.heating_loads_kw["SpaceHeating"]) > 0.0 - @expression(m, CHPToSpaceHeatingKW[ts in p.time_steps], - m[:dvHeatingProduction]["CHP","SpaceHeating",ts] - CHPToHotTESByQuality["SpaceHeating",ts] - CHPToSteamTurbineByQualityKW["SpaceHeating",ts] - - CHPThermalToWasteByQualityKW["SpaceHeating",ts] - CHPtoAbsorptionChillerByQualityKW["SpaceHeating",ts] - ) + if "SpaceHeating" in p.heating_loads && chp.can_serve_space_heating && sum(p.heating_loads_kw["SpaceHeating"]) > 0.0 + CHPToSpaceHeatingKW = [value(m[:dvHeatingProduction][chp_name,"SpaceHeating",ts]) - CHPToHotTESByQuality["SpaceHeating"][ts] - CHPToSteamTurbineByQualityKW["SpaceHeating"][ts] - CHPThermalToWasteByQualityKW["SpaceHeating"][ts] - CHPtoAbsorptionChillerByQualityKW["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.(max.(0.0, value.(CHPToSpaceHeatingKW ./ KWH_PER_MMBTU)), digits=5) + r["thermal_to_space_heating_load_series_mmbtu_per_hour"] = round.(max.(0.0, CHPToSpaceHeatingKW ./ KWH_PER_MMBTU), digits=5) - if "ProcessHeat" in p.heating_loads && p.s.chp.can_serve_process_heat && sum(p.heating_loads_kw["ProcessHeat"]) > 0.0 - @expression(m, CHPToProcessHeatKW[ts in p.time_steps], - m[:dvHeatingProduction]["CHP","ProcessHeat",ts] - CHPToHotTESByQuality["ProcessHeat",ts] - CHPToSteamTurbineByQualityKW["ProcessHeat",ts] - - CHPThermalToWasteByQualityKW["ProcessHeat",ts] - CHPtoAbsorptionChillerByQualityKW["ProcessHeat",ts] - ) + if "ProcessHeat" in p.heating_loads && chp.can_serve_process_heat && sum(p.heating_loads_kw["ProcessHeat"]) > 0.0 + CHPToProcessHeatKW = [value(m[:dvHeatingProduction][chp_name,"ProcessHeat",ts]) - CHPToHotTESByQuality["ProcessHeat"][ts] - CHPToSteamTurbineByQualityKW["ProcessHeat"][ts] - CHPThermalToWasteByQualityKW["ProcessHeat"][ts] - CHPtoAbsorptionChillerByQualityKW["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.(max.(0.0, value.(CHPToProcessHeatKW ./ KWH_PER_MMBTU)), digits=5) - - r["thermal_to_load_series_mmbtu_per_hour"] = r["thermal_to_dhw_load_series_mmbtu_per_hour"] .+ r["thermal_to_space_heating_load_series_mmbtu_per_hour"] .+ r["thermal_to_process_heat_load_series_mmbtu_per_hour"] - r["thermal_production_series_mmbtu_per_hour"] = r["thermal_to_load_series_mmbtu_per_hour"] .+ r["thermal_to_absorption_chiller_series_mmbtu_per_hour"] .+ r["thermal_to_storage_series_mmbtu_per_hour"] .+ r["thermal_to_steamturbine_series_mmbtu_per_hour"] - r["annual_thermal_production_mmbtu"] = round(p.hours_per_time_step * sum(r["thermal_production_series_mmbtu_per_hour"]), digits=3) + r["thermal_to_process_heat_load_series_mmbtu_per_hour"] = round.(max.(0.0, CHPToProcessHeatKW ./ KWH_PER_MMBTU), digits=5) - r["year_one_fuel_cost_before_tax"] = round(value(m[:TotalCHPFuelCosts] / p.pwf_fuel["CHP"]), digits=3) + # Calculate individual CHP fuel cost (not the total across all CHPs) + chp_fuel_cost_lifecycle = sum(p.pwf_fuel[chp_name] * value(m[:dvFuelUsage][chp_name, ts]) * p.fuel_cost_per_kwh[chp_name][ts] for ts in p.time_steps) + r["year_one_fuel_cost_before_tax"] = round(chp_fuel_cost_lifecycle / p.pwf_fuel[chp_name], digits=3) r["year_one_fuel_cost_after_tax"] = r["year_one_fuel_cost_before_tax"] * (1 - p.s.financial.offtaker_tax_rate_fraction) - r["lifecycle_fuel_cost_after_tax"] = round(value(m[:TotalCHPFuelCosts]) * (1- p.s.financial.offtaker_tax_rate_fraction), digits=3) + r["lifecycle_fuel_cost_after_tax"] = round(chp_fuel_cost_lifecycle * (1- p.s.financial.offtaker_tax_rate_fraction), digits=3) + #Standby charges and hourly O&M - r["year_one_standby_cost_before_tax"] = round(value(m[Symbol("TotalCHPStandbyCharges")]) / p.pwf_e, digits=0) + r["year_one_standby_cost_before_tax"] = round(chp.standby_rate_per_kw_per_month * 12 * value(m[Symbol("dvSize"*_n)][chp_name]), digits=0) r["year_one_standby_cost_after_tax"] = r["year_one_standby_cost_before_tax"] * (1 - p.s.financial.offtaker_tax_rate_fraction) - r["lifecycle_standby_cost_after_tax"] = round(value(m[Symbol("TotalCHPStandbyCharges")]) * (1 - p.s.financial.offtaker_tax_rate_fraction), digits=0) - r["initial_capital_costs"] = round(value(m[Symbol("CHPCapexNoIncentives")]), digits=2) + r["lifecycle_standby_cost_after_tax"] = round(r["year_one_standby_cost_before_tax"] * p.pwf_e * (1 - p.s.financial.offtaker_tax_rate_fraction), digits=0) + + # Calculate individual CHP capital costs using purchased CHP capacity (excluding existing_kw) + size_kw = value(m[Symbol("dvPurchaseSize"*_n)][chp_name]) + suppl_firing_size = value(m[Symbol("dvSupplementaryFiringSize"*_n)][chp_name]) + + # Base CHP capital cost using generalized function + capital_cost = get_tech_initial_capex(chp, size_kw) + + # Add supplementary firing capital cost + capital_cost += chp.supplementary_firing_capital_cost_per_kw * suppl_firing_size + + r["initial_capital_costs"] = round(capital_cost, digits=2) - 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/financial.jl b/src/results/financial.jl index 724b9900f..8c6f98bf1 100644 --- a/src/results/financial.jl +++ b/src/results/financial.jl @@ -252,8 +252,8 @@ function calculate_lcoe(p::REoptInputs, tech_results::Dict, tech::AbstractTech) federal_tax_rate_fraction = p.s.financial.offtaker_tax_rate_fraction end capital_costs = if typeof(tech) == PV && :tech_sizes_for_cost_curve in fieldnames(typeof(tech)) - # Use PV-specific cost curve calculation for PV tech - get_pv_initial_capex(p, tech, new_kw) + # Use cost curve calculation for PV tech + get_tech_initial_capex(tech, new_kw) else # Use simple calculation for other techs like Wind new_kw * tech.installed_cost_per_kw @@ -361,29 +361,35 @@ function get_depreciation_schedule(p::REoptInputs, tech::Union{AbstractTech,Abst return depreciation_schedule end -function get_pv_initial_capex(p::REoptInputs, pv::AbstractTech, size_kw::Float64) - cost_list = pv.installed_cost_per_kw - size_list = pv.tech_sizes_for_cost_curve - pv_size = size_kw +""" + get_tech_initial_capex(tech::AbstractTech, size_kw::Float64) + +Calculate initial capital costs for a technology with either a simple cost per kW or a segmented cost curve. +Works for PV, CHP, and other technologies with cost curves defined by tech_sizes_for_cost_curve. +""" +function get_tech_initial_capex(tech::AbstractTech, size_kw::Float64) + cost_list = tech.installed_cost_per_kw + size_list = tech.tech_sizes_for_cost_curve initial_capex = 0.0 if typeof(cost_list) == Vector{Float64} - if pv_size <= size_list[1] - initial_capex = pv_size * cost_list[1] - elseif pv_size > size_list[end] - initial_capex = pv_size * cost_list[end] + if size_kw <= size_list[1] + initial_capex = size_kw * cost_list[1] + elseif size_kw > size_list[end] + initial_capex = size_kw * cost_list[end] else for s in 2:length(size_list) - if (pv_size > size_list[s-1]) && (pv_size <= size_list[s]) + if (size_kw > size_list[s-1]) && (size_kw <= size_list[s]) slope = (cost_list[s] * size_list[s] - cost_list[s-1] * size_list[s-1]) / (size_list[s] - size_list[s-1]) - initial_capex = cost_list[s-1] * size_list[s-1] + (pv_size - size_list[s-1]) * slope + initial_capex = cost_list[s-1] * size_list[s-1] + (size_kw - size_list[s-1]) * slope end end end else - initial_capex = cost_list * pv_size + initial_capex = cost_list * size_kw end return initial_capex -end \ No newline at end of file +end + diff --git a/src/results/proforma.jl b/src/results/proforma.jl index 450d53759..8353a3ae3 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) @@ -366,19 +372,19 @@ function update_metrics(m::Metrics, p::REoptInputs, tech::AbstractTech, tech_nam total_kw = results[tech_name]["size_kw"] existing_kw = :existing_kw in fieldnames(typeof(tech)) ? tech.existing_kw : 0 new_kw = total_kw - existing_kw - if tech_name == "CHP" - capital_cost = results["CHP"]["initial_capital_costs"] + if tech_name in [chp.name for chp in p.s.chps] + capital_cost = results[tech_name]["initial_capital_costs"] elseif tech_name in [pv.name for pv in p.s.pvs] # Check if it's a PV technology - capital_cost = get_pv_initial_capex(p, tech, new_kw) + capital_cost = get_tech_initial_capex(tech, new_kw) else capital_cost = new_kw * tech.installed_cost_per_kw end # owner is responsible for only new technologies operating and maintenance cost in optimal case # CHP doesn't have existing CHP, and it has different O&M cost parameters - if tech_name == "CHP" - hours_operating = sum(results["CHP"]["electric_production_series_kw"] .> 0.0) / (8760 * p.s.settings.time_steps_per_hour) - annual_om = -1 * (results["CHP"]["annual_electric_production_kwh"] * tech.om_cost_per_kwh + + if tech_name in [chp.name for chp in p.s.chps] + hours_operating = sum(results[tech_name]["electric_production_series_kw"] .> 0.0) * p.hours_per_time_step + annual_om = -1 * (results[tech_name]["annual_electric_production_kwh"] * tech.om_cost_per_kwh + new_kw * tech.om_cost_per_kw + new_kw * tech.om_cost_per_hr_per_kw_rated * hours_operating) elseif third_party @@ -391,9 +397,9 @@ function update_metrics(m::Metrics, p::REoptInputs, tech::AbstractTech, tech_nam m.om_series += escalate_om(annual_om) m.om_series_bau += escalate_om(-1 * existing_kw * tech.om_cost_per_kw) - if tech_name == "CHP" + if tech_name in [chp.name for chp in p.s.chps] escalate_fuel(val, esc_rate) = [val * (1 + esc_rate)^yr for yr in 1:years] - fuel_cost = results["CHP"]["year_one_fuel_cost_before_tax"] + fuel_cost = results[tech_name]["year_one_fuel_cost_before_tax"] m.fuel_cost_series += escalate_fuel(-1 * fuel_cost, p.s.financial.chp_fuel_cost_escalation_rate_fraction) end @@ -411,7 +417,7 @@ function update_metrics(m::Metrics, p::REoptInputs, tech::AbstractTech, tech_nam pbi_series = Float64[] pbi_series_bau = Float64[] existing_energy_bau = third_party ? get(results[tech_name], "year_one_energy_produced_kwh_bau", 0) : 0 - if tech_name == "CHP" + if tech_name in [chp.name for chp in p.s.chps] year_one_energy = results[tech_name]["annual_electric_production_kwh"] else year_one_energy = "year_one_energy_produced_kwh" in keys(results[tech_name]) ? results[tech_name]["year_one_energy_produced_kwh"] : results[tech_name]["annual_energy_produced_kwh"] @@ -516,10 +522,21 @@ function irr(cashflows::AbstractArray{<:Real, 1}) return 0.0 end f(r) = npv(cashflows, r) - rate = 0.0 - try - rate = fzero(f, [0.0, 0.99]) - finally - return round(rate, digits=3) + a, b = 0.0, 0.99 + fa = f(a) + fb = f(b) + + # Only use bracketed solve when endpoint signs differ. + # Otherwise return 0.0 instead of throwing an ArgumentError from Roots. + if !(isfinite(fa) && isfinite(fb)) || fa * fb >= 0 + return 0.0 end + + rate = try + fzero(f, [a, b]) + catch + return 0.0 + end + + return round(rate, digits=3) end diff --git a/src/results/pv.jl b/src/results/pv.jl index 034cb5e05..58d8a14fe 100644 --- a/src/results/pv.jl +++ b/src/results/pv.jl @@ -38,7 +38,7 @@ function add_pv_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="") r["installed_cost_per_kw"] = pv_tech.installed_cost_per_kw else # Get cost from cost curve - r["installed_cost_per_kw"] = get_pv_initial_capex(p, pv_tech, r["size_kw"]) / r["size_kw"] + r["installed_cost_per_kw"] = get_tech_initial_capex(pv_tech, r["size_kw"]) / r["size_kw"] end if !isnothing(pv_tech.size_class) && !isempty(pv_tech.tech_sizes_for_cost_curve) diff --git a/src/results/results.jl b/src/results/results.jl index f0ccf8532..48c083643 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 a6d3da742..7ed3aa633 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -41,7 +41,7 @@ else # run HiGHS tests @test s.financial.elec_cost_escalation_rate_fraction == 0.0166 @test s.financial.om_cost_escalation_rate_fraction == 0.025 @test s.financial.offtaker_discount_rate_fraction == 0.0624 - 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 @@ -66,7 +66,7 @@ else # run HiGHS tests @test s.financial.elec_cost_escalation_rate_fraction == 0.0074 #national avg @test s.financial.om_cost_escalation_rate_fraction == 0.015 @test s.financial.offtaker_discount_rate_fraction == 0.045 - 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 @@ -91,7 +91,7 @@ else # run HiGHS tests @test s.financial.elec_cost_escalation_rate_fraction == -0.0062 @test s.financial.om_cost_escalation_rate_fraction == 0.015 @test s.financial.offtaker_discount_rate_fraction == 0.045 - 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 @@ -157,7 +157,7 @@ else # run HiGHS tests dataset, distance, datasource = REopt.call_solar_dataset_api(latitude, longitude, radius) @test dataset == "nsrdb" - # 3. Oulu, Findland + # 3. Oulu, Finland latitude, longitude = 65.0102196310875, 25.465387094897675 radius = 0 dataset, distance, datasource = REopt.call_solar_dataset_api(latitude, longitude, radius) @@ -317,7 +317,7 @@ else # run HiGHS tests ) ) s = Scenario(d) - @test s.chp.max_kw ≈ case3_max_size_kw atol=1e-3 + @test s.chps[1].max_kw ≈ case3_max_size_kw atol=1e-3 #Case 5: with higher electric load, still sized according to thermal input avg_electric_load_kw = 2000.0 @@ -1397,6 +1397,31 @@ else # run HiGHS tests empty!(m2) GC.gc() end + + @testset "Multiple CHPs" begin + include("test_multiple_chps.jl") + end + + @testset "Avoid CHP Binary" begin + include("test_avoid_chp_binary.jl") + end + + @testset "Off-Grid CHP" begin + include("test_chp_offgrid.jl") + end + + @testset "CHP Ramp Rate" begin + include("test_ramp.jl") + end + + @testset "CHP Production Factor" begin + include("test_chp_prodfactor.jl") + end + + @testset "CHP Existing KW" begin + include("test_chp_existing_kw.jl") + end + end @testset verbose=true "FlexibleHVAC" begin @@ -2248,8 +2273,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.027 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.027 atol=0.0001 delete!(input_data, "SpaceHeatingLoad") delete!(input_data, "DomesticHotWaterLoad") @@ -2265,7 +2290,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)) @@ -2275,8 +2300,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.023 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.023 atol=0.0001 #Update CHP prime_mover and test new defaults input_data["CHP"]["prime_mover"] = "combustion_turbine" input_data["CHP"]["size_class"] = 1 @@ -2286,8 +2311,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.014999999999999999 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.014999999999999999 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 @@ -2328,7 +2353,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 @@ -2336,14 +2361,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.015) 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.015) 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 @@ -4182,7 +4207,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)) @@ -4199,7 +4224,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/boiler_steamturbine.json b/test/scenarios/boiler_steamturbine.json new file mode 100644 index 000000000..9893ae6b3 --- /dev/null +++ b/test/scenarios/boiler_steamturbine.json @@ -0,0 +1,51 @@ +{ + "Site": { + "latitude": 37.78, + "longitude": -122.45 + }, + "ElectricLoad": { + "doe_reference_name": "Hospital", + "annual_kwh": 5000000.0 + }, + "SpaceHeatingLoad": { + "doe_reference_name": "Hospital", + "annual_mmbtu": 10000.0 + }, + "DomesticHotWaterLoad": { + "doe_reference_name": "Hospital", + "annual_mmbtu": 2000.0 + }, + "ElectricTariff": { + "blended_annual_energy_rate": 0.10, + "blended_annual_demand_rate": 10.0 + }, + "ExistingBoiler":{ + "fuel_cost_per_mmbtu": 8.0, + "efficiency": 0.75, + "can_supply_steam_turbine": false + }, + "Boiler": { + "min_mmbtu_per_hour": 0.0, + "max_mmbtu_per_hour": 100.0, + "efficiency": 0.8, + "fuel_type": "natural_gas", + "fuel_cost_per_mmbtu": 8.0, + "installed_cost_per_mmbtu_per_hour": 50000.0, + "om_cost_per_mmbtu_per_hour": 500.0, + "om_cost_per_mmbtu": 0.5, + "can_supply_steam_turbine": true + }, + "SteamTurbine": { + "min_kw": 0.0, + "max_kw": 5000.0, + "electric_produced_to_thermal_consumed_ratio": 0.30, + "thermal_produced_to_thermal_consumed_ratio": 0.50, + "installed_cost_per_kw": 1000.0, + "om_cost_per_kw": 10.0, + "om_cost_per_kwh": 0.01, + "can_wholesale": false, + "can_curtail": false, + "macrs_bonus_fraction": 0.0, + "macrs_option_years": 5 + } +} diff --git a/test/scenarios/chp_offgrid.json b/test/scenarios/chp_offgrid.json new file mode 100644 index 000000000..1d60c92f8 --- /dev/null +++ b/test/scenarios/chp_offgrid.json @@ -0,0 +1,34 @@ +{ + "Settings": { + "time_steps_per_hour": 1, + "off_grid_flag": true + }, + "Site": { + "latitude": 37.78, + "longitude": -122.45 + }, + "ElectricLoad": { + "doe_reference_name": "Hospital", + "annual_kwh": 1000000.0 + }, + "CHP": { + "prime_mover": "recip_engine", + "size_class": 2, + "min_kw": 50.0, + "max_kw": 500.0, + "installed_cost_per_kw": 2300.0, + "om_cost_per_kwh": 0.0, + "om_cost_per_kw": 50.0, + "fuel_cost_per_mmbtu": 8.0, + "is_electric_only": true + }, + "ElectricStorage": { + "min_kw": 0.0, + "max_kw": 1000.0, + "min_kwh": 0.0, + "max_kwh": 5000.0, + "soc_min_fraction": 0.2, + "soc_init_fraction": 1.0, + "can_grid_charge": false + } +} diff --git a/test/scenarios/multiple_chps.json b/test/scenarios/multiple_chps.json new file mode 100644 index 000000000..1d67a77b6 --- /dev/null +++ b/test/scenarios/multiple_chps.json @@ -0,0 +1,72 @@ +{ + "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.15, + "blended_annual_demand_rate": 15.0 + }, + "DomesticHotWaterLoad": { + "doe_reference_name": "Hospital", + "annual_mmbtu": 5000.0 + }, + "SpaceHeatingLoad": { + "doe_reference_name": "Hospital", + "annual_mmbtu": 10000.0 + }, + "ExistingBoiler": { + "fuel_cost_per_mmbtu": 12.0 + }, + "CHP": [ + { + "name": "CHP_recip_engine", + "prime_mover": "recip_engine", + "size_class": 3, + "installed_cost_per_kw": 3000.0, + "macrs_bonus_fraction": 0.0, + "macrs_option_years": 0, + "min_turn_down_fraction": 0.0, + "min_kw": 150.0, + "max_kw": 150.0, + "fuel_cost_per_mmbtu": 8.0, + "can_serve_dhw": true, + "can_serve_space_heating": true, + "standby_rate_per_kw_per_month": 5.0, + "om_cost_per_kwh": 0.05, + "om_cost_per_kw": 10.0, + "om_cost_per_hr_per_kw_rated": 0.05 + }, + { + "name": "CHP_micro_turbine", + "prime_mover": "micro_turbine", + "size_class": 3, + "installed_cost_per_kw": 3500.0, + "macrs_bonus_fraction": 0.0, + "macrs_option_years": 0, + "min_turn_down_fraction": 0.0, + "min_kw": 100.0, + "max_kw": 100.0, + "fuel_cost_per_mmbtu": 8.0, + "can_serve_dhw": true, + "can_serve_space_heating": false, + "standby_rate_per_kw_per_month": 5.0, + "om_cost_per_kw": 10.0, + "om_cost_per_kwh": 0.05, + "om_cost_per_hr_per_kw_rated": 0.05, + "supplementary_firing_max_steam_ratio": 1.2 + } + ] +} diff --git a/test/scenarios/nuclear_battery.json b/test/scenarios/nuclear_battery.json new file mode 100644 index 000000000..32204db44 --- /dev/null +++ b/test/scenarios/nuclear_battery.json @@ -0,0 +1,28 @@ +{ + "Site": { + "latitude": 37.78, + "longitude": -122.45 + }, + "ElectricLoad": {}, + "ElectricTariff": { + "blended_annual_energy_rate": 0.12, + "blended_annual_demand_rate": 10.0 + }, + "CHP": { + "fuel_cost_per_mmbtu": 0.1, + "prime_mover": "combustion_turbine", + "installed_cost_per_kw": 12000.0, + "om_cost_per_kw": 150.0, + "om_cost_per_kwh": 0.0, + "electric_efficiency_full_load": 0.25, + "electric_efficiency_half_load": 0.25, + "min_turn_down_fraction": 0.0, + "is_electric_only": true, + "can_curtail": false, + "unavailability_periods": [{"month":1, + "start_week_of_month":2, + "start_day_of_week":1, + "start_hour":1, + "duration_hours":0}] + } +} diff --git a/test/test_avoid_chp_binary.jl b/test/test_avoid_chp_binary.jl new file mode 100644 index 000000000..0179b89f1 --- /dev/null +++ b/test/test_avoid_chp_binary.jl @@ -0,0 +1,205 @@ + +# using Revise +# using REopt +# using JSON +# using DelimitedFiles +# using PlotlyJS +# using Dates +# using Test +# using JuMP +# using HiGHS +# using DotEnv +# DotEnv.load!() + +############### CHP Binary Creation Tests ################### +# Tests to verify that binCHPIsOnInTS is only created when needed + +@testset "CHP Binary Creation Tests" begin + + @testset "CHP with min_turn_down_fraction > 0 (binary SHOULD be created)" begin + # Test that binCHPIsOnInTS gets created when min_turn_down_fraction > 0 + data = JSON.parsefile("./scenarios/chp_unavailability_outage.json") + + # Set min_turn_down_fraction to a non-zero value + data["CHP"]["min_turn_down_fraction"] = 0.5 + + # Remove outage to simplify the test + delete!(data, "ElectricUtility") + data["ElectricUtility"] = Dict("net_metering_limit_kw" => 0.0, "co2_from_avert" => true) + + s = Scenario(data) + inputs = REoptInputs(s) + m = Model(optimizer_with_attributes(HiGHS.Optimizer, + "output_flag" => false, + "log_to_console" => false, + "mip_rel_gap" => 0.01)) + results = run_reopt(m, inputs) + + # Verify the binary variable was created + @test haskey(m.obj_dict, :binCHPIsOnInTS) + binary_var = m[:binCHPIsOnInTS] + + # Verify model solved successfully + @test results["CHP"]["size_kw"] > 0.0 + + println("✓ Test 1 passed: Binary created with min_turn_down_fraction = 0.5") + finalize(backend(m)) + empty!(m) + end + + @testset "CHP with min_turn_down_fraction = 0 (binary should NOT be created)" begin + # Test that binCHPIsOnInTS is NOT created when all conditions are zero/equal + data = JSON.parsefile("./scenarios/chp_unavailability_outage.json") + + # Set min_turn_down_fraction to zero + data["CHP"]["min_turn_down_fraction"] = 0.0 + + # Ensure no intercepts by making efficiencies equal (no part-load curve) + # When electric_efficiency_full_load == electric_efficiency_half_load, intercept = 0 + data["CHP"]["electric_efficiency_full_load"] = 0.35 + data["CHP"]["electric_efficiency_half_load"] = 0.35 + data["CHP"]["thermal_efficiency_full_load"] = 0.45 + data["CHP"]["thermal_efficiency_half_load"] = 0.45 + + # Ensure no hourly O&M cost + data["CHP"]["om_cost_per_hr_per_kw_rated"] = 0.0 + + # Remove outage to simplify + delete!(data, "ElectricUtility") + data["ElectricUtility"] = Dict("net_metering_limit_kw" => 0.0, "co2_from_avert" => true) + + s = Scenario(data) + inputs = REoptInputs(s) + m = Model(optimizer_with_attributes(HiGHS.Optimizer, + "output_flag" => false, + "log_to_console" => false, + "mip_rel_gap" => 0.01)) + results = run_reopt(m, inputs) + + # Verify the binary variable was NOT created + @test !haskey(m.obj_dict, :binCHPIsOnInTS) + + # Verify model still solved successfully + @test results["CHP"]["size_kw"] > 0.0 + + println("✓ Test 2 passed: Binary NOT created with min_turn_down_fraction = 0 and no intercepts") + finalize(backend(m)) + empty!(m) + end + + @testset "CHP with fuel_burn_intercept > 0 (binary SHOULD be created)" begin + # Test that binary is created when fuel burn has non-zero intercept + data = JSON.parsefile("./scenarios/chp_unavailability_outage.json") + + # Set min_turn_down_fraction to zero + data["CHP"]["min_turn_down_fraction"] = 0.0 + + # Create non-zero fuel burn intercept by having different efficiencies + data["CHP"]["electric_efficiency_full_load"] = 0.40 + data["CHP"]["electric_efficiency_half_load"] = 0.30 # Different from full load + data["CHP"]["thermal_efficiency_full_load"] = 0.45 + data["CHP"]["thermal_efficiency_half_load"] = 0.45 # Keep thermal same to isolate fuel burn + + # Ensure no hourly O&M cost + data["CHP"]["om_cost_per_hr_per_kw_rated"] = 0.0 + + delete!(data, "ElectricUtility") + data["ElectricUtility"] = Dict("net_metering_limit_kw" => 0.0, "co2_from_avert" => true) + + s = Scenario(data) + inputs = REoptInputs(s) + m = Model(optimizer_with_attributes(HiGHS.Optimizer, + "output_flag" => false, + "log_to_console" => false, + "mip_rel_gap" => 0.01)) + results = run_reopt(m, inputs) + + # Verify the binary variable was created + @test haskey(m.obj_dict, :binCHPIsOnInTS) + + # verify model solved successfully + @test results["Financial"]["lcc"] != 0.0 + + println("✓ Test 3 passed: Binary created with non-zero fuel_burn_intercept") + finalize(backend(m)) + empty!(m) + end + + @testset "CHP with thermal_prod_intercept > 0 (binary SHOULD be created)" begin + # Test that binary is created when thermal production has non-zero intercept + data = JSON.parsefile("./scenarios/chp_unavailability_outage.json") + + # Set min_turn_down_fraction to zero + data["CHP"]["min_turn_down_fraction"] = 0.0 + + # Create non-zero thermal prod intercept by having different thermal efficiencies + data["CHP"]["electric_efficiency_full_load"] = 0.35 + data["CHP"]["electric_efficiency_half_load"] = 0.35 # Keep electric same + data["CHP"]["thermal_efficiency_full_load"] = 0.40 + data["CHP"]["thermal_efficiency_half_load"] = 0.50 # Different from full load + + # Ensure no hourly O&M cost + data["CHP"]["om_cost_per_hr_per_kw_rated"] = 0.0 + + delete!(data, "ElectricUtility") + data["ElectricUtility"] = Dict("net_metering_limit_kw" => 0.0, "co2_from_avert" => true) + + s = Scenario(data) + inputs = REoptInputs(s) + m = Model(optimizer_with_attributes(HiGHS.Optimizer, + "output_flag" => false, + "log_to_console" => false, + "mip_rel_gap" => 0.01)) + results = run_reopt(m, inputs) + + # Verify the binary variable was created + @test haskey(m.obj_dict, :binCHPIsOnInTS) + + # Verify model solved successfully + @test results["Financial"]["lcc"] != 0.0 + + println("✓ Test 4 passed: Binary created with non-zero thermal_prod_intercept") + finalize(backend(m)) + empty!(m) + end + + @testset "CHP with om_cost_per_hr_per_kw_rated > 0 (binary SHOULD be created)" begin + # Test that binary is created when hourly O&M cost is non-zero + data = JSON.parsefile("./scenarios/chp_unavailability_outage.json") + + # Set min_turn_down_fraction to zero + data["CHP"]["min_turn_down_fraction"] = 0.0 + + # Make efficiencies equal (no intercepts) + data["CHP"]["electric_efficiency_full_load"] = 0.35 + data["CHP"]["electric_efficiency_half_load"] = 0.35 + data["CHP"]["thermal_efficiency_full_load"] = 0.45 + data["CHP"]["thermal_efficiency_half_load"] = 0.45 + + # Set non-zero hourly O&M cost + data["CHP"]["om_cost_per_hr_per_kw_rated"] = 0.01 + + delete!(data, "ElectricUtility") + data["ElectricUtility"] = Dict("net_metering_limit_kw" => 0.0, "co2_from_avert" => true) + + s = Scenario(data) + inputs = REoptInputs(s) + m = Model(optimizer_with_attributes(HiGHS.Optimizer, + "output_flag" => false, + "log_to_console" => false, + "mip_rel_gap" => 0.01)) + results = run_reopt(m, inputs) + + # Verify the binary variable was created + @test haskey(m.obj_dict, :binCHPIsOnInTS) + + # Verify model solved successfully + @test results["Financial"]["lcc"] != 0.0 + + println("✓ Test 5 passed: Binary created with om_cost_per_hr_per_kw_rated = 0.01") + finalize(backend(m)) + empty!(m) + end + + println("\n===== All CHP Binary Creation Tests Passed! =====\n") +end \ No newline at end of file diff --git a/test/test_boiler_steamturbine.jl b/test/test_boiler_steamturbine.jl new file mode 100644 index 000000000..c09b73b63 --- /dev/null +++ b/test/test_boiler_steamturbine.jl @@ -0,0 +1,47 @@ +using Revise +using REopt +using JSON +using DelimitedFiles +using PlotlyJS +using Dates +using Test +using JuMP +using HiGHS +using DotEnv +DotEnv.load!() + +############### Boiler + SteamTurbine with simplified parameters ################### +# This test uses the simpler SteamTurbine setup with electric_produced_to_thermal_consumed_ratio +# and thermal_produced_to_thermal_consumed_ratio instead of detailed steam parameters + +input_data = JSON.parsefile("./scenarios/boiler_steamturbine.json") +s = Scenario(input_data) +inputs = REoptInputs(s) + +m1 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.001, "output_flag" => false, "log_to_console" => false)) +m2 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.001, "output_flag" => false, "log_to_console" => false)) +results = run_reopt([m1,m2], inputs) + +println("\n=== Boiler + SteamTurbine Test Results ===") +println("Financial NPV: ", round(results["Financial"]["npv"], digits=2)) +println("Financial LCC: ", round(results["Financial"]["lcc"], digits=2)) +println("\nBoiler:") +println(" Size (MMBtu/hr): ", round(results["Boiler"]["size_mmbtu_per_hour"], digits=2)) +println(" Annual thermal production (MMBtu): ", round(results["Boiler"]["annual_thermal_production_mmbtu"], digits=2)) +println(" Annual fuel consumption (MMBtu): ", round(results["Boiler"]["annual_fuel_consumption_mmbtu"], digits=2)) +println(" Thermal to SteamTurbine (MMBtu): ", round(sum(results["Boiler"]["thermal_to_steamturbine_series_mmbtu_per_hour"]), digits=2)) +println("\nSteamTurbine:") +println(" Size (kW): ", round(results["SteamTurbine"]["size_kw"], digits=2)) +println(" Annual electric production (kWh): ", round(results["SteamTurbine"]["annual_electric_production_kwh"], digits=2)) +println(" Annual thermal consumption (MMBtu): ", round(results["SteamTurbine"]["annual_thermal_consumption_mmbtu"], digits=2)) +println(" Annual thermal production (MMBtu): ", round(results["SteamTurbine"]["annual_thermal_production_mmbtu"], digits=2)) + +# Verify energy balance: Boiler thermal to ST should match ST thermal consumption +boiler_to_st = sum(results["Boiler"]["thermal_to_steamturbine_series_mmbtu_per_hour"]) +st_thermal_in = results["SteamTurbine"]["annual_thermal_consumption_mmbtu"] +println("\nEnergy Balance Check:") +println(" Boiler->SteamTurbine: ", round(boiler_to_st, digits=2), " MMBtu") +println(" SteamTurbine thermal in: ", round(st_thermal_in, digits=2), " MMBtu") +println(" Difference: ", round(abs(boiler_to_st - st_thermal_in), digits=2), " MMBtu") + +println("=============================================\n") diff --git a/test/test_chp_existing_kw.jl b/test/test_chp_existing_kw.jl new file mode 100644 index 000000000..9dfa6bc7d --- /dev/null +++ b/test/test_chp_existing_kw.jl @@ -0,0 +1,147 @@ +# using Revise +# using REopt +# using JSON +# using Test +# using JuMP +# using HiGHS + +############### Existing CHP Sizing ################### + +input_base = JSON.parsefile("./scenarios/chp_sizing.json") +input_base["CHP"]["min_kw"] = 0.0 +input_base["CHP"]["max_kw"] = 200.0 + +input_existing = deepcopy(input_base) +input_existing["CHP"]["existing_kw"] = 100.0 + +s_existing = Scenario(input_existing) +p_existing = REoptInputs(s_existing) +chp_name = s_existing.chps[1].name + +# Bounds: max_kw is NEW capacity; internal max_sizes = existing + max_kw +@test p_existing.existing_sizes[chp_name] == 100.0 +@test p_existing.min_sizes[chp_name] == 100.0 +@test p_existing.max_sizes[chp_name] == 300.0 + +m1 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +r_base = run_reopt(m1, REoptInputs(Scenario(input_base))) + +m2 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +r_existing = run_reopt(m2, p_existing) + +new_purchase = r_existing["CHP"]["size_kw"] - 100.0 +expected_capex = REopt.get_tech_initial_capex(s_existing.chps[1], new_purchase) + + s_existing.chps[1].supplementary_firing_capital_cost_per_kw * r_existing["CHP"]["size_supplemental_firing_kw"] + +@test r_existing["CHP"]["initial_capital_costs"] ≈ expected_capex atol=1e-2 +@test r_existing["CHP"]["initial_capital_costs"] < r_base["CHP"]["initial_capital_costs"] +@test r_existing["Financial"]["initial_capital_costs"] < r_base["Financial"]["initial_capital_costs"] + +############### BAU Scenario ################### + +bau_s = REopt.BAUScenario(s_existing) +@test length(bau_s.chps) == 1 +@test bau_s.chps[1].existing_kw == 100.0 +@test bau_s.chps[1].supplementary_firing_capital_cost_per_kw == 0.0 + +bau_inputs = REopt.BAUInputs(p_existing) +@test bau_inputs.max_sizes[chp_name] == 100.0 +@test bau_inputs.min_sizes[chp_name] == 100.0 +@test bau_inputs.cap_cost_slope[chp_name] == 0.0 + +# With BAU: existing CHP reduces BAU LCC and yields smaller incremental NPV +m3 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +m4 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +r_full_existing = run_reopt([m3, m4], p_existing) + +m5 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +m6 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +r_full_base = run_reopt([m5, m6], REoptInputs(Scenario(input_base))) + +@test r_full_existing["Financial"]["lcc_bau"] < r_full_base["Financial"]["lcc_bau"] +@test r_full_existing["Financial"]["npv"] < r_full_base["Financial"]["npv"] +@test r_full_existing["Financial"]["npv"] > 0 + +############### Existing Only (No additional NEW CHP) ################### + +input_noexp = JSON.parsefile("./scenarios/chp_sizing.json") +input_noexp["CHP"]["existing_kw"] = 150.0 +input_noexp["CHP"]["min_kw"] = 0.0 +input_noexp["CHP"]["max_kw"] = 0.0 + +m7 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +m8 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +r_noexp = run_reopt([m7, m8], REoptInputs(Scenario(input_noexp))) + +@test r_noexp["CHP"]["size_kw"] ≈ 150.0 atol=0.1 +@test r_noexp["CHP"]["initial_capital_costs"] == 0.0 +@test r_noexp["CHP"]["annual_electric_production_kwh"] > 0 +@test r_noexp["CHP"]["annual_thermal_production_mmbtu"] > 0 +@test r_noexp["CHP"]["year_one_fuel_cost_before_tax"] > 0 +@test abs(r_noexp["Financial"]["npv"]) < 1000.0 + +############### Outage Resilience ################### + +input_outage = JSON.parsefile("./scenarios/chp_sizing.json") +input_outage["CHP"]["existing_kw"] = 150.0 +input_outage["CHP"]["min_kw"] = 0.0 +input_outage["CHP"]["max_kw"] = 0.0 +input_outage["Generator"] = Dict( + "existing_kw" => 100.0, "min_kw" => 0.0, "max_kw" => 0.0, + "installed_cost_per_kw" => 500.0, "om_cost_per_kw" => 10.0, + "fuel_cost_per_gallon" => 3.0, + "electric_efficiency_full_load" => 0.32, "electric_efficiency_half_load" => 0.32, + "fuel_avail_gal" => 1000.0, + "only_runs_during_grid_outage" => true, "sells_energy_back_to_grid" => false +) +input_outage["ElectricUtility"] = Dict( + "outage_start_time_step" => 5000, "outage_end_time_step" => 5048, + "co2_from_avert" => true +) +input_outage["ElectricLoad"] = Dict( + "doe_reference_name" => "Hospital", "annual_kwh" => 1000000.0, + "critical_load_fraction" => 0.5 +) + +m9 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +m10 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +r_outage = run_reopt([m9, m10], REoptInputs(Scenario(input_outage))) + +@test r_outage["CHP"]["size_kw"] ≈ 150.0 atol=0.1 +@test r_outage["Generator"]["size_kw"] ≈ 100.0 atol=0.1 +@test r_outage["CHP"]["initial_capital_costs"] == 0.0 +@test r_outage["CHP"]["annual_electric_production_kwh"] > 0 +@test r_outage["ElectricLoad"]["bau_critical_load_met"] == true +@test abs(r_outage["Financial"]["npv"]) < 5000.0 + +############### Multiple CHPs with Existing ################### + +input_multi = JSON.parsefile("./scenarios/multiple_chps.json") +input_multi["CHP"][1]["existing_kw"] = 150.0 +input_multi["CHP"][1]["min_kw"] = 0.0 +input_multi["CHP"][1]["max_kw"] = 0.0 +input_multi["CHP"][2]["existing_kw"] = 0.0 +input_multi["CHP"][2]["min_kw"] = 100.0 +input_multi["CHP"][2]["max_kw"] = 100.0 + +s_multi = Scenario(input_multi) +p_multi = REoptInputs(s_multi) + +bau_multi = REopt.BAUScenario(s_multi) +@test length(bau_multi.chps) == 1 +@test bau_multi.chps[1].name == "CHP_recip_engine" + +m11 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +m12 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +r_multi = run_reopt([m11, m12], p_multi) + +chp_results = r_multi["CHP"] +recip = first(filter(c -> c["name"] == "CHP_recip_engine", chp_results)) +micro = first(filter(c -> c["name"] == "CHP_micro_turbine", chp_results)) + +@test recip["size_kw"] ≈ 150.0 atol=0.1 +@test recip["initial_capital_costs"] == 0.0 +@test micro["size_kw"] ≈ 100.0 atol=0.1 +@test micro["initial_capital_costs"] > 0.0 +@test recip["annual_electric_production_kwh"] > 0 +@test micro["annual_electric_production_kwh"] > 0 diff --git a/test/test_chp_offgrid.jl b/test/test_chp_offgrid.jl new file mode 100644 index 000000000..bfc8fa7dc --- /dev/null +++ b/test/test_chp_offgrid.jl @@ -0,0 +1,105 @@ +# using Revise +# using REopt +# using JSON +# using DelimitedFiles +# using PlotlyJS +# using Dates +# using Test +# using JuMP +# using HiGHS +# using DotEnv +# DotEnv.load!() + +############### Off-Grid CHP with Operating Reserves Test ################### + +# Test off-grid scenario with CHP and ElectricStorage where CHP requires 10% operating reserves +# This tests that ElectricStorage provides the required reserves for CHP generation + +input_data = JSON.parsefile("./scenarios/chp_offgrid.json") +input_data["ElectricLoad"]["min_load_met_annual_fraction"] = 0.999 +input_data["ElectricLoad"]["operating_reserve_required_fraction"] = 0.0 +# TODO see if no battery shows up when below is set to 0 so CHP can provide SR +input_data["CHP"]["operating_reserve_required_fraction"] = 0.2 + +println("\n" * "="^80) +println("Testing Off-Grid CHP with Operating Reserves") +println("="^80) + +# Create scenario and run optimization +s = Scenario(input_data) +inputs = REoptInputs(s) + +println("\nScenario Setup:") +println(" - CHP operating_reserve_required_fraction: ", s.chps[1].operating_reserve_required_fraction) +println(" - ElectricLoad operating_reserve_required_fraction: ", s.electric_load.operating_reserve_required_fraction) +println(" - Off-grid flag: ", s.settings.off_grid_flag) + +m = Model(optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false, "log_to_console" => false, "mip_rel_gap" => 0.01)) +results = run_reopt(m, inputs) + +# Check that the model solved successfully +@test termination_status(m) == MOI.OPTIMAL || termination_status(m) == MOI.LOCALLY_SOLVED + +println("\nOptimization Results:") +println(" - Termination status: ", termination_status(m)) +println(" - CHP size: ", round(results["CHP"]["size_kw"], digits=1), " kW") +println(" - Battery power: ", round(results["ElectricStorage"]["size_kw"], digits=1), " kW") +println(" - Battery energy: ", round(results["ElectricStorage"]["size_kwh"], digits=1), " kWh") +println(" - Load met fraction: ", round(results["ElectricLoad"]["offgrid_load_met_fraction"], digits=4)) + +# Calculate operating reserve requirements +chp_production_to_load = results["CHP"]["electric_to_load_series_kw"] +chp_or_required = sum(chp_production_to_load .* s.chps[1].operating_reserve_required_fraction) +load_or_required = sum(s.electric_load.critical_loads_kw .* s.electric_load.operating_reserve_required_fraction) +total_or_required = chp_or_required + load_or_required + +# Get operating reserve provided +or_provided = sum(results["ElectricLoad"]["offgrid_annual_oper_res_provided_series_kwh"]) + +println("\nOperating Reserve Analysis:") +println(" - CHP OR required: ", round(chp_or_required, digits=1), " kWh") +println(" - Load OR required: ", round(load_or_required, digits=1), " kWh") +println(" - Total OR required: ", round(total_or_required, digits=1), " kWh") +println(" - Total OR provided: ", round(or_provided, digits=1), " kWh") +println(" - OR margin: ", round(or_provided - total_or_required, digits=1), " kWh") + +# Test 1: Operating reserves provided >= operating reserves required (with tolerance for solver gap) +@test or_provided >= total_or_required * (1 - 0.02) # Allow 2% gap due to mip_rel_gap and rounding +println("\n✓ Test 1 PASSED: Operating reserves provided (", round(or_provided, digits=1), + " kWh) >= required (", round(total_or_required, digits=1), " kWh) within tolerance") + +# Test 2: CHP produces electricity in off-grid scenario +chp_annual_production = sum(results["CHP"]["electric_to_load_series_kw"]) + + sum(results["CHP"]["electric_to_storage_series_kw"]) +@test chp_annual_production > 0 +println("✓ Test 2 PASSED: CHP produces electricity (", round(chp_annual_production, digits=1), " kWh)") + +# Test 3: Battery is sized to provide operating reserves +battery_sized = results["ElectricStorage"]["size_kw"] > 0 || results["ElectricStorage"]["size_kwh"] > 0 +@test battery_sized +println("✓ Test 3 PASSED: Battery is sized (", round(results["ElectricStorage"]["size_kw"], digits=1), + " kW, ", round(results["ElectricStorage"]["size_kwh"], digits=1), " kWh)") + +# Test 4: Load is met according to minimum fraction +@test results["ElectricLoad"]["offgrid_load_met_fraction"] >= s.electric_load.min_load_met_annual_fraction +println("✓ Test 4 PASSED: Load met fraction (", round(results["ElectricLoad"]["offgrid_load_met_fraction"], digits=4), + ") >= minimum (", s.electric_load.min_load_met_annual_fraction, ")") + +# Test 5: No grid interaction in off-grid scenario +@test results["ElectricUtility"]["annual_energy_supplied_kwh"] ≈ 0.0 atol=0.01 +println("✓ Test 5 PASSED: No grid interaction (", results["ElectricUtility"]["annual_energy_supplied_kwh"], " kWh)") + +# Test 6: Financial calculation includes off-grid costs +f = results["Financial"] +lcc_check = f["lifecycle_generation_tech_capital_costs"] + f["lifecycle_storage_capital_costs"] + + f["lifecycle_om_costs_after_tax"] + f["lifecycle_fuel_costs_after_tax"] + + f["lifecycle_chp_standby_cost_after_tax"] + f["lifecycle_elecbill_after_tax"] + + f["lifecycle_offgrid_other_annual_costs_after_tax"] + f["lifecycle_offgrid_other_capital_costs"] + + f["lifecycle_outage_cost"] + f["lifecycle_MG_upgrade_and_fuel_cost"] - + f["lifecycle_production_incentive_after_tax"] +@test lcc_check ≈ f["lcc"] atol=1.0 +println("✓ Test 6 PASSED: Financial calculations consistent (LCC: \$", round(f["lcc"], digits=0), ")") + +println("\n" * "="^80) +println("All tests PASSED for Off-Grid CHP with Operating Reserves!") +println("="^80 * "\n") \ No newline at end of file diff --git a/test/test_chp_prodfactor.jl b/test/test_chp_prodfactor.jl new file mode 100644 index 000000000..5209eec34 --- /dev/null +++ b/test/test_chp_prodfactor.jl @@ -0,0 +1,59 @@ +# using Revise +# using REopt +# using JSON +# using DelimitedFiles +# using PlotlyJS +# using Dates +# using Test +# using JuMP +# using HiGHS +# using DotEnv +# DotEnv.load!() + + +########### CHP production_factor_series test ############# +@testset "CHP with production_factor_series" begin + # Create a simple scenario with electric-only CHP and a custom production factor + input_data = Dict( + "Site" => Dict( + "latitude" => 39.7, + "longitude" => -104.9 + ), + "ElectricLoad" => Dict( + "loads_kw" => repeat([100.0], 8760), + "year" => 2025 + ), + "ElectricTariff" => Dict( + "urdb_label" => "5ed6c1a15457a3367add15ae" + ), + "CHP" => Dict( + "is_electric_only" => true, + "fuel_cost_per_mmbtu" => 4.0, + "max_kw" => 100.0, + "min_kw" => 100.0, + "production_factor_series" => vcat(repeat([0.5], 4380), repeat([1.0], 4380)) # 50% for first half year, 100% for second half + ) + ) + + s = Scenario(input_data) + p = REoptInputs(s) + + # Verify the production factor series was properly assigned + @test !isnothing(s.chps[1].production_factor_series) + @test length(s.chps[1].production_factor_series) == 8760 + @test s.chps[1].production_factor_series[1] ≈ 0.5 atol=0.001 + @test s.chps[1].production_factor_series[8760] ≈ 1.0 atol=0.001 + + # Run the optimization + m = Model(optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false, "log_to_console" => false, "mip_rel_gap" => 0.01)) + results = run_reopt(m, p) + + # Check that production matches the pattern (lower in first half, higher in second half) + first_half_avg = sum(results["CHP"]["electric_production_series_kw"][1:4380]) / 4380 + second_half_avg = sum(results["CHP"]["electric_production_series_kw"][4381:8760]) / 4380 + + println("CHP production_factor_series test passed!") + println("CHP size: ", results["CHP"]["size_kw"], " kW") + println("First half avg production: ", round(first_half_avg, digits=2), " kW") + println("Second half avg production: ", round(second_half_avg, digits=2), " kW") +end \ No newline at end of file diff --git a/test/test_multiple_chps.jl b/test/test_multiple_chps.jl new file mode 100644 index 000000000..94fb398c8 --- /dev/null +++ b/test/test_multiple_chps.jl @@ -0,0 +1,50 @@ +# using Revise +# using REopt +# using JSON +# using DelimitedFiles +# using PlotlyJS +# using Dates +# using Test +# using JuMP +# using HiGHS +# using DotEnv +# DotEnv.load!() + +############### Multiple CHPs Test ################### + +# With BAU Scenario +input_data = JSON.parsefile("./scenarios/multiple_chps.json") +s = Scenario(input_data) +inputs = REoptInputs(s) + +m1 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.05, "output_flag" => false, "log_to_console" => false)) +m2 = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.05, "output_flag" => false, "log_to_console" => false)) +results = run_reopt([m1,m2], inputs) + +@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"]) > 500.0 +@test sum(CHP_micro_turbine["electric_production_series_kw"]) > 500.0 + +# Test that LCC matches sum of discounted cashflows for optimal case (within 0.1% tolerance) +# Note, cashflows are negative where lcc is positive, so negating cashflows in comparison +optimal_lcc = results["Financial"]["lcc"] +optimal_cashflow_sum = -1*sum(results["Financial"]["offtaker_discounted_annual_free_cashflows"]) +@test isapprox(optimal_lcc, optimal_cashflow_sum, rtol=0.001) + +# Test that LCC matches sum of discounted cashflows for BAU case (within 0.1% tolerance) +bau_lcc = results["Financial"]["lcc_bau"] +bau_cashflow_sum = -1*sum(results["Financial"]["offtaker_discounted_annual_free_cashflows_bau"]) +@test isapprox(bau_lcc, bau_cashflow_sum, rtol=0.001) + +# Add a test to compare npv with difference in optimal and bau cashflows +npv_calculated = bau_cashflow_sum - optimal_cashflow_sum +npv_reported = results["Financial"]["npv"] +@test isapprox(npv_calculated, npv_reported, rtol=0.001) diff --git a/test/test_ramp.jl b/test/test_ramp.jl new file mode 100644 index 000000000..f1f96b1b7 --- /dev/null +++ b/test/test_ramp.jl @@ -0,0 +1,187 @@ +# using Revise +# using REopt +# using JSON +# using DelimitedFiles +# using PlotlyJS +# using Dates +# using Test +# using JuMP +# using HiGHS +# using DotEnv +# DotEnv.load!() + +# Notes to test for Off-Grid DC powered by Nuclear + Storage +# +# DONE: CHP electric only with ramp rate limit, no outage but force battery size to avoid utility serving the high ramp times +# Boiler + SteamTurbine (+ HighTempThermalStorage?) +# Add CHP and/or Boiler+ST to off-grid with SR requirements for nuclear +# Generator and Battery currently have unlimited SR **supply** +# PV has SR requirement - use that as proxy for nuclear and/or geothermal? +# Focus on analysis needs for project, but eventually consider GenericGenerator (+Heat?) with needed attributes + + +# Create a 15-minute interval load profile from hourly with sine wave sub-hourly variations +function upsample_load_variation(hourly_load; intra_hour_variation=0.2) + # Create 15-minute load profile with sine wave sub-hourly variations + fifteen_min_loads = Float64[] + for (hour_idx, hourly_load) in enumerate(hourly_load) + # Create 4 sub-hourly values per hour with sine wave variation + # Base load with ±10% sine wave variation within each hour + for quarter_hour in 1:4 + # Intra-hour phase angle for this quarter hour (0, π/2, π, 3π/2) + intra_hour_phase = (quarter_hour - 1) * π / 2 + # Sine wave variation: ±20% of hourly load + variation = intra_hour_variation * sin(intra_hour_phase + 2π * hour_idx / 24) # Daily cycle component + sub_hourly_load = hourly_load * (1.0 + variation) + push!(fifteen_min_loads, max(0.0, sub_hourly_load)) # Ensure non-negative + end + end + + return fifteen_min_loads +end + + +############### Nuclear as CHP with ramp rate and battery to support ################### + +# scenario = 1 # "Specify CHP+battery sizes 15-min load with variation, no outages" +# scenario = 2 # "Least-cost sizes CHP+battery for year long outage hourly load from CRB" +# scenario = 3 # "Force CHP/Nuclear to peak load size, but still size battery for ramp rate limit" +scenario = 3 +# Load scenario with CHP (electric only) cost and performance params similar to nuclear +input_data = JSON.parsefile("./scenarios/nuclear_battery.json") + +# New ramp rate input/constraint (fraction of capacity per hour) +input_data["CHP"]["ramp_rate_fraction_per_hour"] = 0.1 + +if scenario == 1 + # Create 15-minute load profile from hourly with sine wave intra-hour variation + sim_input = Dict( + "load_type" => "electric", + "doe_reference_name" => "Hospital", + "annual_kwh" => 4.0e6, + "latitude" => input_data["Site"]["latitude"], + "longitude" => input_data["Site"]["longitude"], + "year" => 2023) + hourly_loads_kw = simulated_load(sim_input)["loads_kw"] + fifteen_min_loads = upsample_load_variation(hourly_loads_kw; intra_hour_variation=0.05) + input_data["Settings"] = Dict("time_steps_per_hour" => 4) + input_data["ElectricLoad"] = Dict("loads_kw" => fifteen_min_loads, "year" => 2023) + + # You can adjust the fixed size here if desired + fixed_chp_size_kw = 800.0 + input_data["CHP"]["min_kw"] = fixed_chp_size_kw + input_data["CHP"]["max_kw"] = fixed_chp_size_kw + + # See if REopt sizes battery to support nuclear ramp rate limitation + input_data["ElectricStorage"] = Dict( + "min_kw" => 800.0, + "min_kwh" => 800.0, + "max_kw" => 800.0, + "max_kwh" => 800.0 + ) +elseif scenario == 2 || scenario == 3 + # Year-long outage scenario from CHP test scenarios + input_data["ElectricLoad"] = Dict("doe_reference_name" => "Hospital", "annual_kwh" => 4.0e6) + input_data["ElectricLoad"]["critical_load_fraction"] = 1.0 + input_data["ElectricUtility"] = Dict("outage_start_time_step" => 1, + "outage_end_time_step" => 8760) + input_data["ElectricStorage"] = Dict() + if scenario == 3 + # Force CHP to peak load size, but still size battery for ramp rate limit + # Estimate peak load from annual energy assuming capacity factor + estimated_peak_load_kw = 740.0 # Assuming 50% capacity factor + input_data["CHP"]["min_kw"] = estimated_peak_load_kw + input_data["CHP"]["max_kw"] = estimated_peak_load_kw + end +end + +# This was somehow being set to 0 in dictkeys_tosymbols in utils.py which was causing an error, +# but function was updated to avoid this +# input_data["Settings"]["off_grid_flag"] = false + +# Create scenario and inputs +s = Scenario(input_data) +inputs = REoptInputs(s) +ts_per_hour = s.settings.time_steps_per_hour + +# Run optimization with single model (no BAU comparison needed for fixed sizing) +m = Model(optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.01, "output_flag" => false, "log_to_console" => false)) +results = run_reopt(m, inputs) + +# Print key results +println("\n===== Nuclear as CHP with ramp rate and battery to support =====") +println("\nCHP Results:") +println("CHP Size (kW): ", results["CHP"]["size_kw"]) +println("Annual CHP Production (kWh): ", round(results["CHP"]["annual_electric_production_kwh"], digits=0)) +println("CHP to Load (kWh): ", round(sum(results["CHP"]["electric_to_load_series_kw"])/ts_per_hour, digits=0)) +println("CHP to Battery (kWh): ", round(sum(results["CHP"]["electric_to_storage_series_kw"])/ts_per_hour, digits=0)) +println("CHP to Grid (kWh): ", round(sum(results["CHP"]["electric_to_grid_series_kw"])/ts_per_hour, digits=0)) +println("CHP Curtailed (kWh): ", round(sum(results["CHP"]["electric_curtailed_series_kw"])/ts_per_hour, digits=0)) +println("\nBattery Results:") +println("Battery Energy Size (kWh): ", results["ElectricStorage"]["size_kwh"]) +println("Battery Power Size (kW): ", results["ElectricStorage"]["size_kw"]) +println("Battery to Load (kWh): ", round(sum(results["ElectricStorage"]["storage_to_load_series_kw"])/ts_per_hour, digits=0)) +println("\nGrid Results:") +println("Utility to Load (kWh): ", round(sum(results["ElectricUtility"]["annual_energy_supplied_kwh"]), digits=0)) +# println("\nFinancial Results:") +# println("NPV (\$): ", round(results["Financial"]["npv"], digits=0)) +# println("Simple Payback (years): ", round(results["Financial"]["simple_payback_years"], digits=2)) +println("==========================================\n") + +# # Create stacked area chart showing how CHP and battery serve the load +# load_series = results["ElectricLoad"]["load_series_kw"] +# chp_to_load = results["CHP"]["electric_to_load_series_kw"] +# batt_to_load = results["ElectricStorage"]["storage_to_load_series_kw"] +# grid_to_load = results["ElectricUtility"]["electric_to_load_series_kw"] + +# # Create time steps array (in hours for x-axis) +# time_hours = collect(1:length(load_series)) ./ ts_per_hour + +# # Create stacked area chart +# plt = plot( +# [ +# scatter( +# x=time_hours, +# y=grid_to_load, +# mode="lines", +# name="Grid to Load", +# fill="tozeroy", +# line=attr(width=0.5, color="rgb(128,128,128)"), +# fillcolor="rgba(128,128,128,0.5)" +# ), +# scatter( +# x=time_hours, +# y=grid_to_load .+ chp_to_load, +# mode="lines", +# name="CHP to Load", +# fill="tonexty", +# line=attr(width=0.5, color="rgb(255,128,0)"), +# fillcolor="rgba(255,128,0,0.6)" +# ), +# scatter( +# x=time_hours, +# y=grid_to_load .+ chp_to_load .+ batt_to_load, +# mode="lines", +# name="Battery to Load", +# fill="tonexty", +# line=attr(width=0.5, color="rgb(0,128,255)"), +# fillcolor="rgba(0,128,255,0.6)" +# ), +# scatter( +# x=time_hours, +# y=load_series, +# mode="lines", +# name="Total Load", +# line=attr(width=2, color="rgb(0,0,0)", dash="dot") +# ) +# ], +# Layout( +# title="Electric Load Service Breakdown - CHP with Ramp Rate and Battery Support", +# xaxis=attr(title="Time (hours)", range=[0, 8760]), +# yaxis=attr(title="Power (kW)"), +# showlegend=true, +# hovermode="x unified", +# legend=attr(x=1.02, y=1, xanchor="left", yanchor="top") +# ) +# ) +# display(plt) \ No newline at end of file diff --git a/test/test_with_xpress.jl b/test/test_with_xpress.jl index b9847d7ff..b171b662f 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