diff --git a/src/constraints/chp_constraints.jl b/src/constraints/chp_constraints.jl index 0f25bf54f..67af76526 100644 --- a/src/constraints/chp_constraints.jl +++ b/src/constraints/chp_constraints.jl @@ -1,17 +1,11 @@ # REopt®, Copyright (c) Alliance for Sustainable Energy, LLC. See also https://github.com/NREL/REopt.jl/blob/master/LICENSE. -function add_chp_fuel_burn_constraints(m, p; _n="") - # Fuel burn slope and intercept - fuel_burn_slope, fuel_burn_intercept = fuel_slope_and_intercept(; - electric_efficiency_full_load = p.s.chp.electric_efficiency_full_load, - electric_efficiency_half_load = p.s.chp.electric_efficiency_half_load, - fuel_higher_heating_value_kwh_per_unit=1 - ) - +function add_chp_fuel_burn_constraints(m, p; _n="", binary_created=true, fuel_burn_slope, fuel_burn_intercept) # Fuel cost m[:TotalCHPFuelCosts] = @expression(m, sum(p.pwf_fuel[t] * m[:dvFuelUsage][t, ts] * p.fuel_cost_per_kwh[t][ts] for t in p.techs.chp, ts in p.time_steps) ) # Conditionally add dvFuelBurnYIntercept if coefficient p.FuelBurnYIntRate is greater than ~zero + # Note: if intercept exists, binary_created will always be true (checked in add_chp_constraints) if abs(fuel_burn_intercept) > 1.0E-7 dv = "dvFuelBurnYIntercept"*_n m[Symbol(dv)] = @variable(m, [p.techs.chp, p.time_steps], base_name=dv) @@ -40,15 +34,9 @@ function add_chp_fuel_burn_constraints(m, p; _n="") end end -function add_chp_thermal_production_constraints(m, p; _n="") - # Thermal production slope and intercept - thermal_prod_full_load = 1.0 / p.s.chp.electric_efficiency_full_load * p.s.chp.thermal_efficiency_full_load # [kWt/kWe] - thermal_prod_half_load = 0.5 / p.s.chp.electric_efficiency_half_load * p.s.chp.thermal_efficiency_half_load # [kWt/kWe] - thermal_prod_slope = (thermal_prod_full_load - thermal_prod_half_load) / (1.0 - 0.5) # [kWt/kWe] - thermal_prod_intercept = thermal_prod_full_load - thermal_prod_slope * 1.0 # [kWt/kWe_rated - - +function add_chp_thermal_production_constraints(m, p; _n="", binary_created=true, thermal_prod_slope, thermal_prod_intercept) # Conditionally add dvHeatingProductionYIntercept if coefficient p.s.chpThermalProdIntercept is greater than ~zero + # Note: if intercept exists, binary_created will always be true (checked in add_chp_constraints) if abs(thermal_prod_intercept) > 1.0E-7 dv = "dvHeatingProductionYIntercept"*_n m[Symbol(dv)] = @variable(m, [p.techs.chp, p.time_steps], base_name=dv) @@ -92,16 +80,20 @@ Used by add_chp_constraints to add supplementary firing constraints if p.s.chp.supplementary_firing_max_steam_ratio > 1.0 to add CHP supplementary firing operating constraints. Else, the supplementary firing dispatch and size decision variables are set to zero. """ -function add_chp_supplementary_firing_constraints(m, p; _n="") - thermal_prod_full_load = 1.0 / p.s.chp.electric_efficiency_full_load * p.s.chp.thermal_efficiency_full_load # [kWt/kWe] - thermal_prod_half_load = 0.5 / p.s.chp.electric_efficiency_half_load * p.s.chp.thermal_efficiency_half_load # [kWt/kWe] - thermal_prod_slope = (thermal_prod_full_load - thermal_prod_half_load) / (1.0 - 0.5) # [kWt/kWe] - +function add_chp_supplementary_firing_constraints(m, p; _n="", thermal_prod_slope, thermal_prod_intercept) # Constrain upper limit of dvSupplementaryThermalProduction, using auxiliary variable for (size * useSupplementaryFiring) - @constraint(m, CHPSupplementaryFireCon[t in p.techs.chp, ts in p.time_steps], - m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] <= - (p.s.chp.supplementary_firing_max_steam_ratio - 1.0) * p.production_factor[t,ts] * (thermal_prod_slope * m[Symbol("dvSupplementaryFiringSize"*_n)][t] + m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts]) - ) + # Note: dvHeatingProductionYIntercept only exists when thermal_prod_intercept > 0 + if abs(thermal_prod_intercept) > 1.0E-7 + @constraint(m, CHPSupplementaryFireCon[t in p.techs.chp, ts in p.time_steps], + m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] <= + (p.s.chp.supplementary_firing_max_steam_ratio - 1.0) * p.production_factor[t,ts] * (thermal_prod_slope * m[Symbol("dvSupplementaryFiringSize"*_n)][t] + m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts]) + ) + else + @constraint(m, CHPSupplementaryFireCon[t in p.techs.chp, ts in p.time_steps], + m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] <= + (p.s.chp.supplementary_firing_max_steam_ratio - 1.0) * p.production_factor[t,ts] * thermal_prod_slope * m[Symbol("dvSupplementaryFiringSize"*_n)][t] + ) + end if solver_is_compatible_with_indicator_constraints(p.s.settings.solver_name) # Constrain lower limit of 0 if CHP tech is off @constraint(m, NoCHPSupplementaryFireOffCon[t in p.techs.chp, ts in p.time_steps], @@ -174,13 +166,32 @@ end Used in src/reopt.jl to add_chp_constraints if !isempty(p.techs.chp) to add CHP operating constraints and cost expressions. """ -function add_chp_constraints(m, p; _n="") - # TODO if chp.min_turn_down_fraction is 0.0, and there is no fuel burn or thermal y-intercept, we don't need the binary below - @warn """Adding binary variable to model CHP. - Some solvers are very slow with integer variables""" - @variables m begin - binCHPIsOnInTS[p.techs.chp, p.time_steps], Bin # 1 If technology t is operating in time step; 0 otherwise - end +function add_chp_constraints(m, p; _n="") + # Calculate fuel burn and thermal production slopes and intercepts which are passed into some constraint functions + fuel_burn_slope, fuel_burn_intercept = fuel_slope_and_intercept(; + electric_efficiency_full_load = p.s.chp.electric_efficiency_full_load, + electric_efficiency_half_load = p.s.chp.electric_efficiency_half_load, + fuel_higher_heating_value_kwh_per_unit=1 + ) + + thermal_prod_full_load = 1.0 / p.s.chp.electric_efficiency_full_load * p.s.chp.thermal_efficiency_full_load + thermal_prod_half_load = 0.5 / p.s.chp.electric_efficiency_half_load * p.s.chp.thermal_efficiency_half_load + thermal_prod_slope = (thermal_prod_full_load - thermal_prod_half_load) / (1.0 - 0.5) + thermal_prod_intercept = thermal_prod_full_load - thermal_prod_slope * 1.0 + + # Check if binary variable is needed which creates binCHPIsOnInTS binary decision variable and activates other constraints + binary_needed = (abs(fuel_burn_intercept) > 1.0E-7) || + (abs(thermal_prod_intercept) > 1.0E-7) || + (p.s.chp.min_turn_down_fraction > 1.0E-7) || + (p.s.chp.om_cost_per_hr_per_kw_rated > 1.0E-7) + + if binary_needed + @warn """Adding binary variable binCHPIsOnInTS to model CHP. + Some solvers are very slow with integer variables""" + @variables m begin + binCHPIsOnInTS[p.techs.chp, p.time_steps], Bin # 1 If technology t is operating in time step; 0 otherwise + end + end m[:TotalHourlyCHPOMCosts] = 0 m[:TotalCHPFuelCosts] = 0 @@ -189,17 +200,20 @@ function add_chp_constraints(m, p; _n="") m[:dvRatedProduction][t, ts] for t in p.techs.chp, ts in p.time_steps) ) - if p.s.chp.om_cost_per_hr_per_kw_rated > 1.0E-7 - add_chp_hourly_om_charges(m, p) - end - - add_chp_fuel_burn_constraints(m, p; _n=_n) - add_chp_thermal_production_constraints(m, p; _n=_n) - add_binCHPIsOnInTS_constraints(m, p; _n=_n) + # These constraints are always needed + add_chp_fuel_burn_constraints(m, p; _n=_n, binary_created=binary_needed, fuel_burn_slope=fuel_burn_slope, fuel_burn_intercept=fuel_burn_intercept) + add_chp_thermal_production_constraints(m, p; _n=_n, binary_created=binary_needed, thermal_prod_slope=thermal_prod_slope, thermal_prod_intercept=thermal_prod_intercept) add_chp_rated_prod_constraint(m, p; _n=_n) + # These constraints and decision variables created within the functions are only needed if binary_needed or more specific parameters is non-zero + if binary_needed + add_binCHPIsOnInTS_constraints(m, p; _n=_n) + end + if p.s.chp.om_cost_per_hr_per_kw_rated > 1.0E-7 + add_chp_hourly_om_charges(m, p; _n=_n) + end if p.s.chp.supplementary_firing_max_steam_ratio > 1.0 - add_chp_supplementary_firing_constraints(m,p; _n=_n) + add_chp_supplementary_firing_constraints(m, p; _n=_n, thermal_prod_slope=thermal_prod_slope, thermal_prod_intercept=thermal_prod_intercept) else for t in p.techs.chp fix(m[Symbol("dvSupplementaryFiringSize"*_n)][t], 0.0, force=true) diff --git a/test/test_avoid_chp_binary.jl b/test/test_avoid_chp_binary.jl new file mode 100644 index 000000000..199683cd8 --- /dev/null +++ b/test/test_avoid_chp_binary.jl @@ -0,0 +1,205 @@ + +using Revise +using REopt +using JSON +using DelimitedFiles +using PlotlyJS +using Dates +using Test +using JuMP +using HiGHS +using DotEnv +DotEnv.load!() + +############### CHP Binary Creation Tests ################### +# Tests to verify that binCHPIsOnInTS is only created when needed + +@testset "CHP Binary Creation Tests" begin + + @testset "CHP with min_turn_down_fraction > 0 (binary SHOULD be created)" begin + # Test that binCHPIsOnInTS gets created when min_turn_down_fraction > 0 + data = JSON.parsefile("./scenarios/chp_unavailability_outage.json") + + # Set min_turn_down_fraction to a non-zero value + data["CHP"]["min_turn_down_fraction"] = 0.5 + + # Remove outage to simplify the test + delete!(data, "ElectricUtility") + data["ElectricUtility"] = Dict("net_metering_limit_kw" => 0.0, "co2_from_avert" => true) + + s = Scenario(data) + inputs = REoptInputs(s) + m = Model(optimizer_with_attributes(HiGHS.Optimizer, + "output_flag" => false, + "log_to_console" => false, + "mip_rel_gap" => 0.01)) + results = run_reopt(m, inputs) + + # Verify the binary variable was created + @test haskey(m.obj_dict, :binCHPIsOnInTS) + binary_var = m[:binCHPIsOnInTS] + + # Verify model solved successfully + @test results["CHP"]["size_kw"] > 0.0 + + println("✓ Test 1 passed: Binary created with min_turn_down_fraction = 0.5") + finalize(backend(m)) + empty!(m) + end + + @testset "CHP with min_turn_down_fraction = 0 (binary should NOT be created)" begin + # Test that binCHPIsOnInTS is NOT created when all conditions are zero/equal + data = JSON.parsefile("./scenarios/chp_unavailability_outage.json") + + # Set min_turn_down_fraction to zero + data["CHP"]["min_turn_down_fraction"] = 0.0 + + # Ensure no intercepts by making efficiencies equal (no part-load curve) + # When electric_efficiency_full_load == electric_efficiency_half_load, intercept = 0 + data["CHP"]["electric_efficiency_full_load"] = 0.35 + data["CHP"]["electric_efficiency_half_load"] = 0.35 + data["CHP"]["thermal_efficiency_full_load"] = 0.45 + data["CHP"]["thermal_efficiency_half_load"] = 0.45 + + # Ensure no hourly O&M cost + data["CHP"]["om_cost_per_hr_per_kw_rated"] = 0.0 + + # Remove outage to simplify + delete!(data, "ElectricUtility") + data["ElectricUtility"] = Dict("net_metering_limit_kw" => 0.0, "co2_from_avert" => true) + + s = Scenario(data) + inputs = REoptInputs(s) + m = Model(optimizer_with_attributes(HiGHS.Optimizer, + "output_flag" => false, + "log_to_console" => false, + "mip_rel_gap" => 0.01)) + results = run_reopt(m, inputs) + + # Verify the binary variable was NOT created + @test !haskey(m.obj_dict, :binCHPIsOnInTS) + + # Verify model still solved successfully + @test results["CHP"]["size_kw"] > 0.0 + + println("✓ Test 2 passed: Binary NOT created with min_turn_down_fraction = 0 and no intercepts") + finalize(backend(m)) + empty!(m) + end + + @testset "CHP with fuel_burn_intercept > 0 (binary SHOULD be created)" begin + # Test that binary is created when fuel burn has non-zero intercept + data = JSON.parsefile("./scenarios/chp_unavailability_outage.json") + + # Set min_turn_down_fraction to zero + data["CHP"]["min_turn_down_fraction"] = 0.0 + + # Create non-zero fuel burn intercept by having different efficiencies + data["CHP"]["electric_efficiency_full_load"] = 0.40 + data["CHP"]["electric_efficiency_half_load"] = 0.30 # Different from full load + data["CHP"]["thermal_efficiency_full_load"] = 0.45 + data["CHP"]["thermal_efficiency_half_load"] = 0.45 # Keep thermal same to isolate fuel burn + + # Ensure no hourly O&M cost + data["CHP"]["om_cost_per_hr_per_kw_rated"] = 0.0 + + delete!(data, "ElectricUtility") + data["ElectricUtility"] = Dict("net_metering_limit_kw" => 0.0, "co2_from_avert" => true) + + s = Scenario(data) + inputs = REoptInputs(s) + m = Model(optimizer_with_attributes(HiGHS.Optimizer, + "output_flag" => false, + "log_to_console" => false, + "mip_rel_gap" => 0.01)) + results = run_reopt(m, inputs) + + # Verify the binary variable was created + @test haskey(m.obj_dict, :binCHPIsOnInTS) + + # verify model solved successfully + @test results["Financial"]["lcc"] != 0.0 + + println("✓ Test 3 passed: Binary created with non-zero fuel_burn_intercept") + finalize(backend(m)) + empty!(m) + end + + @testset "CHP with thermal_prod_intercept > 0 (binary SHOULD be created)" begin + # Test that binary is created when thermal production has non-zero intercept + data = JSON.parsefile("./scenarios/chp_unavailability_outage.json") + + # Set min_turn_down_fraction to zero + data["CHP"]["min_turn_down_fraction"] = 0.0 + + # Create non-zero thermal prod intercept by having different thermal efficiencies + data["CHP"]["electric_efficiency_full_load"] = 0.35 + data["CHP"]["electric_efficiency_half_load"] = 0.35 # Keep electric same + data["CHP"]["thermal_efficiency_full_load"] = 0.40 + data["CHP"]["thermal_efficiency_half_load"] = 0.50 # Different from full load + + # Ensure no hourly O&M cost + data["CHP"]["om_cost_per_hr_per_kw_rated"] = 0.0 + + delete!(data, "ElectricUtility") + data["ElectricUtility"] = Dict("net_metering_limit_kw" => 0.0, "co2_from_avert" => true) + + s = Scenario(data) + inputs = REoptInputs(s) + m = Model(optimizer_with_attributes(HiGHS.Optimizer, + "output_flag" => false, + "log_to_console" => false, + "mip_rel_gap" => 0.01)) + results = run_reopt(m, inputs) + + # Verify the binary variable was created + @test haskey(m.obj_dict, :binCHPIsOnInTS) + + # Verify model solved successfully + @test results["Financial"]["lcc"] != 0.0 + + println("✓ Test 4 passed: Binary created with non-zero thermal_prod_intercept") + finalize(backend(m)) + empty!(m) + end + + @testset "CHP with om_cost_per_hr_per_kw_rated > 0 (binary SHOULD be created)" begin + # Test that binary is created when hourly O&M cost is non-zero + data = JSON.parsefile("./scenarios/chp_unavailability_outage.json") + + # Set min_turn_down_fraction to zero + data["CHP"]["min_turn_down_fraction"] = 0.0 + + # Make efficiencies equal (no intercepts) + data["CHP"]["electric_efficiency_full_load"] = 0.35 + data["CHP"]["electric_efficiency_half_load"] = 0.35 + data["CHP"]["thermal_efficiency_full_load"] = 0.45 + data["CHP"]["thermal_efficiency_half_load"] = 0.45 + + # Set non-zero hourly O&M cost + data["CHP"]["om_cost_per_hr_per_kw_rated"] = 0.01 + + delete!(data, "ElectricUtility") + data["ElectricUtility"] = Dict("net_metering_limit_kw" => 0.0, "co2_from_avert" => true) + + s = Scenario(data) + inputs = REoptInputs(s) + m = Model(optimizer_with_attributes(HiGHS.Optimizer, + "output_flag" => false, + "log_to_console" => false, + "mip_rel_gap" => 0.01)) + results = run_reopt(m, inputs) + + # Verify the binary variable was created + @test haskey(m.obj_dict, :binCHPIsOnInTS) + + # Verify model solved successfully + @test results["Financial"]["lcc"] != 0.0 + + println("✓ Test 5 passed: Binary created with om_cost_per_hr_per_kw_rated = 0.01") + finalize(backend(m)) + empty!(m) + end + + println("\n===== All CHP Binary Creation Tests Passed! =====\n") +end \ No newline at end of file