Compute conflicting constraints for infeasible models

CHANGELOG.md

@@ -24,6 +24,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Maintenance formulation for thermal-commit plants (#556).
 - Add new tests for GenX: three-zone, multi-stage, electrolyzer, VRE+storage, 
   piecewise_fuel+CO2, and TDR (#563 and #578).
+- Add functions to compute conflicting constraints when model is infeasible if supported by the solver (#624).
 
 
 ### Fixed

src/case_runners/case_runner.jl

@@ -67,19 +67,21 @@ function run_genx_case_simple!(case::AbstractString, mysetup::Dict, optimizer::A
     myinputs["solve_time"] = solve_time # Store the model solve time in myinputs
 
     # Run MGA if the MGA flag is set to 1 else only save the least cost solution
-    println("Writing Output")
-    outputs_path = get_default_output_folder(case)
-    elapsed_time = @elapsed outputs_path = write_outputs(EP, outputs_path, mysetup, myinputs)
-    println("Time elapsed for writing is")
-    println(elapsed_time)
-    if mysetup["ModelingToGenerateAlternatives"] == 1
-        println("Starting Model to Generate Alternatives (MGA) Iterations")
-        mga(EP, case, mysetup, myinputs, outputs_path)
-    end
+    if has_values(EP)
+        println("Writing Output")
+        outputs_path = get_default_output_folder(case)
+        elapsed_time = @elapsed outputs_path = write_outputs(EP, outputs_path, mysetup, myinputs)
+        println("Time elapsed for writing is")
+        println(elapsed_time)
+        if mysetup["ModelingToGenerateAlternatives"] == 1
+            println("Starting Model to Generate Alternatives (MGA) Iterations")
+            mga(EP, case, mysetup, myinputs, outputs_path)
+        end
 
-    if mysetup["MethodofMorris"] == 1
-        println("Starting Global sensitivity analysis with Method of Morris")
-        morris(EP, case, mysetup, myinputs, outputs_path, OPTIMIZER)
+        if mysetup["MethodofMorris"] == 1
+            println("Starting Global sensitivity analysis with Method of Morris")
+            morris(EP, case, mysetup, myinputs, outputs_path, OPTIMIZER)
+        end
     end
 end
 

src/configure_settings/configure_settings.jl

@@ -21,8 +21,9 @@ function default_settings()
         "MultiStage" => 0,
         "MethodofMorris" => 0,
         "IncludeLossesInESR" => 0,
+        "HydrogenHourlyMatching" => 0,
         "EnableJuMPStringNames" => false,
-        "HydrogenHourlyMatching" => 0
+        "ComputeConflicts" => 0
     )
 end
 
@@ -49,4 +50,8 @@ function validate_settings!(settings::Dict{Any,Any})
         Please see the Methods page in the documentation.""" maxlog=1
     end
 
+    if settings["EnableJuMPStringNames"]==0 && settings["ComputeConflicts"]==1
+        settings["EnableJuMPStringNames"]=1;
+    end
+
 end

src/model/solve_model.jl

@@ -62,22 +62,63 @@ function solve_model(EP::Model, setup::Dict)
 	## Solve Model
 	optimize!(EP)
 
-	if has_duals(EP) # fully linear model
-		println("LP solved for primal")
-	else
-		println("MILP solved for primal")
-	end
+	if has_values(EP)
 
-	if !has_duals(EP) && setup["WriteShadowPrices"] == 1
-		# function to fix integers and linearize problem
-		fix_integers(EP)
-		# re-solve statement for LP solution
-		println("Solving LP solution for duals")
-		optimize!(EP)
-	end
+		if has_duals(EP) # fully linear model
+			println("LP solved for primal")
+		else
+			println("MILP solved for primal")
+		end
+
+		if !has_duals(EP) && setup["WriteShadowPrices"] == 1
+			# function to fix integers and linearize problem
+			fix_integers(EP)
+			# re-solve statement for LP solution
+			println("Solving LP solution for duals")
+			optimize!(EP)
+		end
+
+		## Record solver time
+		solver_time = time() - solver_start_time
+	elseif setup["ComputeConflicts"]==0
+
+		@info "No model solution. You can try to set ComputeConflicts to 1 in the genx_settings.yml file to compute conflicting constraints."
 
-	## Record solver time
-	solver_time = time() - solver_start_time
+	elseif setup["ComputeConflicts"]==1
 
+		@info "No model solution. Trying to identify conflicting constriants..."
+
+		try
+			compute_conflict!(EP)
+		catch e
+			if isa(e, JuMP.ArgumentError)
+				@warn "$(solver_name(EP)) does not support computing conflicting constraints. This is available using either Gurobi or CPLEX."
+				solver_time = time() - solver_start_time
+				return EP, solver_time
+			else
+			 	rethrow(e)
+			end
+		end
+
+		list_of_conflicting_constraints = ConstraintRef[]
+		if get_attribute(EP, MOI.ConflictStatus()) == MOI.CONFLICT_FOUND
+			for (F, S) in list_of_constraint_types(EP)
+				for con in all_constraints(EP, F, S)
+					if get_attribute(con, MOI.ConstraintConflictStatus()) == MOI.IN_CONFLICT
+						push!(list_of_conflicting_constraints, con)
+					end
+				end
+			end
+			display(list_of_conflicting_constraints)
+			solver_time = time() - solver_start_time
+			return EP, solver_time, list_of_conflicting_constraints
+		else
+			@info "Conflicts computation failed."
+			solver_time = time() - solver_start_time
+			return EP, solver_time, list_of_conflicting_constraints
+		end
+
+	end
+
 	return EP, solver_time
-end # END solve_model()
+end # END solve_model()

test/ComputeConflicts/CO2_cap.csv

@@ -0,0 +1,2 @@
+,Network_zones,CO_2_Cap_Zone_1,CO_2_Max_tons_MWh_1,CO_2_Max_Mtons_1
+NE,z1,1,0,0

test/ComputeConflicts/Capacity_reserve_margin.csv

@@ -0,0 +1,2 @@
+,Network_zones,CapRes_1
+NE,z1,0.156

test/ComputeConflicts/Demand_data.csv

@@ -0,0 +1,25 @@
+Voll,Demand_Segment,Cost_of_Demand_Curtailment_per_MW,Max_Demand_Curtailment,Rep_Periods,Timesteps_per_Rep_Period,Sub_Weights,Time_Index,Demand_MW_z1
+50000,1,1,0,1,24,8760,1,11162
+,,,,,,,2,10556
+,,,,,,,3,10105
+,,,,,,,4,9878
+,,,,,,,5,9843
+,,,,,,,6,10017
+,,,,,,,7,10390
+,,,,,,,8,10727
+,,,,,,,9,11298
+,,,,,,,10,11859
+,,,,,,,11,12196
+,,,,,,,12,12321
+,,,,,,,13,12381
+,,,,,,,14,12270
+,,,,,,,15,12149
+,,,,,,,16,12219
+,,,,,,,17,13410
+,,,,,,,18,14539
+,,,,,,,19,14454
+,,,,,,,20,14012
+,,,,,,,21,13494
+,,,,,,,22,12772
+,,,,,,,23,11877
+,,,,,,,24,10874

test/ComputeConflicts/Fuels_data.csv

@@ -0,0 +1,26 @@
+Time_Index,NG
+0,0.05306
+1,5.28
+2,5.28
+3,5.28
+4,5.28
+5,5.28
+6,5.28
+7,5.28
+8,5.28
+9,5.28
+10,5.28
+11,5.28
+12,5.28
+13,5.28
+14,5.28
+15,5.28
+16,5.28
+17,5.28
+18,5.28
+19,5.28
+20,5.28
+21,5.28
+22,5.28
+23,5.28
+24,5.28

test/ComputeConflicts/Generators_data.csv

@@ -0,0 +1,5 @@
+Resource,Zone,THERM,MUST_RUN,STOR,FLEX,HYDRO,VRE,LDS,Num_VRE_Bins,New_Build,Can_Retire,Existing_Cap_MW,Existing_Cap_MWh,Existing_Charge_Cap_MW,Max_Cap_MW,Max_Cap_MWh,Max_Charge_Cap_MW,Min_Cap_MW,Min_Cap_MWh,Min_Charge_Cap_MW,Inv_Cost_per_MWyr,Inv_Cost_per_MWhyr,Inv_Cost_Charge_per_MWyr,Fixed_OM_Cost_per_MWyr,Fixed_OM_Cost_per_MWhyr,Fixed_OM_Cost_Charge_per_MWyr,Var_OM_Cost_per_MWh,Var_OM_Cost_per_MWh_In,Heat_Rate_MMBTU_per_MWh,Fuel,Cap_Size,Start_Cost_per_MW,Start_Fuel_MMBTU_per_MW,Up_Time,Down_Time,Ramp_Up_Percentage,Ramp_Dn_Percentage,Min_Power,Self_Disch,Eff_Up,Eff_Down,Min_Duration,Max_Duration,MaxCapTag_1,MaxCapTag_2,MaxCapTag_3,Resource_Type,CapRes_1,region,cluster
+natural_gas_combined_cycle,1,1,0,0,0,0,0,0,0,1,0,0,0,0,-1,-1,-1,0,0,0,65400,0,0,10287,0,0,3.55,0,7.43,NG,250,91,2,6,6,0.64,0.64,0.468,0,1,1,0,0,0,0,0,natural_gas_fired_combined_cycle,0.93,NE,1
+solar_pv,1,0,0,0,0,0,1,0,1,1,0,0,0,0,-1,-1,-1,0,0,0,85300,0,0,18760,0,0,0,0,9.13,None,0,0,0,0,0,1,1,0,0,1,1,0,0,1,0,0,solar_photovoltaic,0.8,NE,1
+onshore_wind,1,0,0,0,0,0,1,0,1,1,0,0,0,0,-1,-1,-1,0,0,0,97200,0,0,43205,0,0,0.1,0,9.12,None,0,0,0,0,0,1,1,0,0,1,1,0,0,0,1,0,onshore_wind_turbine,0.8,NE,1
+battery,1,0,0,1,0,0,0,0,0,1,0,0,0,0,-1,-1,-1,0,0,0,19584,22494,0,4895,5622,0,0.15,0.15,0,None,0,0,0,0,0,1,1,0,0,0.92,0.92,1,10,0,0,1,battery_mid,0.95,NE,0

test/ComputeConflicts/Generators_variability.csv

@@ -0,0 +1,25 @@
+Time_Index,solar_pv,onshore_wind
+1,0,0.889717042
+2,0,0.877715468
+3,0,0.903424203
+4,0,0.895153165
+5,0,0.757258117
+6,0,0.630928695
+7,0,0.557177782
+8,0,0.6072492
+9,0.1779,0.423417866
+10,0.429,0.007470775
+11,0.5748,0.002535942
+12,0.6484,0.002153709
+13,0.6208,0.00445132
+14,0.596,0.007711587
+15,0.5013,0.100848213
+16,0.3311,0.201802149
+17,0.0642,0.141933054
+18,0,0.567022562
+19,0,0.946024895
+20,0,0.923394203
+21,0,0.953386247
+22,0,0.929205418
+23,0,0.849528909
+24,0,0.665570974

test/ComputeConflicts/Maximum_capacity_requirement.csv

@@ -0,0 +1,4 @@
+MaxCapReqConstraint,ConstraintDescription,Max_MW
+1,PV,0
+2,Wind,0
+3,Batteries,0

test/ComputeConflicts/highs_settings.yml

@@ -0,0 +1,13 @@
+# HiGHS Solver Parameters
+# Common solver settings
+Feasib_Tol: 1.0e-05        # Primal feasibility tolerance # [type: double, advanced: false, range: [1e-10, inf], default: 1e-07]
+Optimal_Tol: 1.0e-05       # Dual feasibility tolerance # [type: double, advanced: false, range: [1e-10, inf], default: 1e-07]
+TimeLimit: 1.0e23             # Time limit # [type: double, advanced: false, range: [0, inf], default: inf]
+Pre_Solve: choose          # Presolve option: "off", "choose" or "on" # [type: string, advanced: false, default: "choose"]
+Method: choose             #HiGHS-specific solver settings # Solver option: "simplex", "choose" or "ipm" # [type: string, advanced: false, default: "choose"]
+
+#highs-specific solver settings
+
+# run the crossover routine for ipx
+# [type: string, advanced: "on", range: {"off", "on"}, default: "off"]
+run_crossover: "on"

test/runtests.jl

@@ -42,4 +42,8 @@ end
     @testset "Multi Stage" begin
         include("test_multistage.jl")
     end
+
+    @testset "Compute Conflicts" begin
+        include("test_compute_conflicts.jl")
+    end
 end

test/test_compute_conflicts.jl

@@ -0,0 +1,24 @@
+module TestConflicts
+
+using Test
+
+include(joinpath(@__DIR__, "utilities.jl"))
+test_path = joinpath(@__DIR__,"ComputeConflicts");
+
+# Define test inputs
+genx_setup =    Dict{Any,Any}(
+    "Trans_Loss_Segments" => 1,
+    "CO2Cap" => 1,
+    "StorageLosses" => 1,
+    "MaxCapReq" => 1,
+    "ComputeConflicts" => 1
+)
+
+genxoutput = redirect_stdout(devnull) do
+    run_genx_case_conflict_testing(test_path, genx_setup)
+end
+
+test_result = @test  length(genxoutput)==2
+write_testlog(test_path,"Testing that the infeasible model is correctly handled",test_result)
+
+end

test/utilities.jl

@@ -36,6 +36,29 @@ function run_genx_case_testing(
     return EP, inputs, OPTIMIZER
 end
 
+function run_genx_case_conflict_testing(
+    test_path::AbstractString,
+    test_setup::Dict,
+    optimizer::Any = HiGHS.Optimizer,
+)
+
+    # Merge the genx_setup with the default settings
+    settings = GenX.default_settings()
+    merge!(settings, test_setup)
+
+    @assert settings["MultiStage"] ∈ [0, 1]
+    # Create a ConsoleLogger that prints any log messages with level >= Error to stderr
+    error_logger = ConsoleLogger(stderr, Logging.Error)
+
+    output = with_logger(error_logger) do
+        OPTIMIZER = configure_solver(test_path, optimizer)
+        inputs = load_inputs(settings, test_path)
+        EP = generate_model(settings, inputs, OPTIMIZER)
+        solve_model(EP, settings)
+    end
+    return output
+end
+
 function run_genx_case_simple_testing(
     test_path::AbstractString,
     genx_setup::Dict,