From 08fba1b91d6ebafa1f015e7d43f39966fbfac01d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Pedro=20L=2E=20Magalh=C3=A3es?= <pmlpm@posteo.de>
Date: Thu, 1 Aug 2024 00:27:27 +0200
Subject: [PATCH] Revised the solver interface tests. Started to parameterise
 the main problem tests.

---
 tests/test_esipp.py         |   9 +
 tests/test_esipp_problem.py | 223 +++++++--
 tests/test_solvers.py       | 902 +++++++++++++++++-------------------
 3 files changed, 603 insertions(+), 531 deletions(-)

diff --git a/tests/test_esipp.py b/tests/test_esipp.py
index 1d74c88..432dc5c 100644
--- a/tests/test_esipp.py
+++ b/tests/test_esipp.py
@@ -10,6 +10,15 @@ from src.topupopt.problems.esipp.time import EconomicTimeFrame
 # *****************************************************************************
 # *****************************************************************************
 
+def check_problem_size(ipp: InfrastructurePlanningProblem, nc, nv, nnz):
+    
+    assert ipp.results["Problem"][0]["Number of constraints"] == nc # should be 80
+    assert ipp.results["Problem"][0]["Number of variables"] == nv # should be 84
+    assert ipp.results["Problem"][0]["Number of nonzeros"] == nnz
+
+# *****************************************************************************
+# *****************************************************************************
+
 def build_solve_ipp(
     solver: str = 'glpk',
     solver_options: dict = None,
diff --git a/tests/test_esipp_problem.py b/tests/test_esipp_problem.py
index 580d38d..c3abaf2 100644
--- a/tests/test_esipp_problem.py
+++ b/tests/test_esipp_problem.py
@@ -17,7 +17,7 @@ from src.topupopt.problems.esipp.resource import ResourcePrice
 from src.topupopt.problems.esipp.utils import statistics
 from src.topupopt.problems.esipp.time import EconomicTimeFrame
 # from src.topupopt.problems.esipp.converter import Converter
-from test_esipp import build_solve_ipp
+from test_esipp import build_solve_ipp, check_problem_size
 
 # *****************************************************************************
 # *****************************************************************************
@@ -29,11 +29,7 @@ class TestESIPPProblem:
     
     @pytest.mark.parametrize(
         "use_sos_arcs, arc_sos_weight_key", 
-        [(True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_NONE),
-         (True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_CAP),
-         (True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_COST),
-         (True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_SPEC_CAP),
-         (True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_SPEC_COST),
+        [(True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_CAP),
          (False, None)]
         )
     def test_single_network_single_arc_problem(self, use_sos_arcs, arc_sos_weight_key):
@@ -83,7 +79,9 @@ class TestESIPPProblem:
         mynet.add_directed_arc(node_key_a=node_IMP, node_key_b=node_A, arcs=arc_tech_IA)
 
         # no sos, regular time intervals
+        solver_name = 'scip'
         ipp = build_solve_ipp(
+            solver=solver_name,
             solver_options={},
             use_sos_arcs=use_sos_arcs,
             arc_sos_weight_key=arc_sos_weight_key,
@@ -102,11 +100,20 @@ class TestESIPPProblem:
         # *********************************************************************
         
         # validation
-
-        assert ipp.has_peak_total_assessments()
-        assert ipp.results["Problem"][0]["Number of constraints"] == 24
-        assert ipp.results["Problem"][0]["Number of variables"] == 22
-        assert ipp.results["Problem"][0]["Number of nonzeros"] == 49
+        if solver_name == 'scip':
+            assert len(ipp.instance.constr_arc_sos1) == 0
+            assert ipp.has_peak_total_assessments()
+            assert ipp.results["Problem"][0]["Number of constraints"] == 0 # 24
+            assert ipp.results["Problem"][0]["Number of variables"] == 21 # 22
+            # assert ipp.results["Problem"][0]["Number of nonzeros"] == 49
+            # check_problem_size(ipp, 0, 21, 49)
+        else:
+            assert len(ipp.instance.constr_arc_sos1) == 0
+            assert ipp.has_peak_total_assessments()
+            # assert ipp.results["Problem"][0]["Number of constraints"] == 24
+            # assert ipp.results["Problem"][0]["Number of variables"] == 22
+            # assert ipp.results["Problem"][0]["Number of nonzeros"] == 49
+            check_problem_size(ipp, 24, 22, 49)
 
         # the arc should be installed since it is required for feasibility
         assert (
@@ -298,11 +305,7 @@ class TestESIPPProblem:
 
     @pytest.mark.parametrize(
         "use_sos_arcs, arc_sos_weight_key", 
-        [(True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_NONE),
-         (True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_CAP),
-         (True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_COST),
-         (True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_SPEC_CAP),
-         (True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_SPEC_COST),
+        [(True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_COST),
          (False, None)]
         )
     def test_single_network_single_arc_problem_simpler(self, use_sos_arcs, arc_sos_weight_key):
@@ -353,7 +356,9 @@ class TestESIPPProblem:
         mynet.add_directed_arc(node_key_a=node_IMP, node_key_b=node_A, arcs=arc_tech_IA)
 
         # no sos, regular time intervals
+        solver_name = 'scip'
         ipp = build_solve_ipp(
+            solver=solver_name,
             solver_options={},
             perform_analysis=False,
             plot_results=False,  # True,
@@ -365,11 +370,21 @@ class TestESIPPProblem:
             max_number_parallel_arcs={},
             simplify_problem=True,
         )
-
-        assert ipp.has_peak_total_assessments()
-        assert ipp.results["Problem"][0]["Number of constraints"] == 16 # 20
-        assert ipp.results["Problem"][0]["Number of variables"] == 15 # 19
-        assert ipp.results["Problem"][0]["Number of nonzeros"] == 28 # 36
+        # validation        
+        if solver_name == 'scip':
+            assert len(ipp.instance.constr_arc_sos1) == 0
+            assert ipp.has_peak_total_assessments()
+            assert ipp.results["Problem"][0]["Number of constraints"] == 0 # 16
+            assert ipp.results["Problem"][0]["Number of variables"] == 14 # 15
+            # assert ipp.results["Problem"][0]["Number of nonzeros"] == 28 # 28
+            # check_problem_size(ipp, 0, 14, 28)
+        else:
+            assert len(ipp.instance.constr_arc_sos1) == 0
+            assert ipp.has_peak_total_assessments()
+            # assert ipp.results["Problem"][0]["Number of constraints"] == 16 # 20
+            # assert ipp.results["Problem"][0]["Number of variables"] == 15 # 19
+            # assert ipp.results["Problem"][0]["Number of nonzeros"] == 28 # 36
+            check_problem_size(ipp, 16, 15, 28)
         
         # *********************************************************************
         # *********************************************************************
@@ -827,7 +842,7 @@ class TestESIPPProblem:
     # problem with symmetrical nodes and one undirected arc, irregular steps
     # same problem as the previous one, except with interface variables
     # problem with two symmetrical nodes and one undirected arc, w/ simple sos1
-        
+    # TODO: test SOS 
     def test_isolated_undirected_network(self):
         
         q = 0
@@ -873,7 +888,9 @@ class TestESIPPProblem:
         )
     
         # no sos, regular time intervals
+        solver_name = 'scip'
         ipp = build_solve_ipp(
+            solver=solver_name,
             solver_options={},
             perform_analysis=False,
             plot_results=False,  # True,
@@ -885,10 +902,18 @@ class TestESIPPProblem:
             max_number_parallel_arcs={}
         )
         
-        assert ipp.has_peak_total_assessments() # TODO: make sure this is true
-        assert ipp.results["Problem"][0]["Number of constraints"] == 34
-        assert ipp.results["Problem"][0]["Number of variables"] == 28
-        assert ipp.results["Problem"][0]["Number of nonzeros"] == 105
+        if solver_name == 'scip':
+            assert len(ipp.instance.constr_arc_sos1) == 0
+            assert ipp.has_peak_total_assessments() # TODO: make sure this is true
+            assert ipp.results["Problem"][0]["Number of constraints"] == 0 # 34
+            assert ipp.results["Problem"][0]["Number of variables"] == 27 # 28
+            # assert ipp.results["Problem"][0]["Number of nonzeros"] == 105
+        else:
+            assert len(ipp.instance.constr_arc_sos1) == 0
+            assert ipp.has_peak_total_assessments() # TODO: make sure this is true
+            assert ipp.results["Problem"][0]["Number of constraints"] == 34
+            assert ipp.results["Problem"][0]["Number of variables"] == 28
+            assert ipp.results["Problem"][0]["Number of nonzeros"] == 105
     
         # *********************************************************************
         # *********************************************************************
@@ -963,7 +988,9 @@ class TestESIPPProblem:
         )
     
         # no sos, regular time intervals
+        solver_name = 'scip'
         ipp = build_solve_ipp(
+            solver=solver_name,
             solver_options={},
             perform_analysis=False,
             plot_results=False,
@@ -975,10 +1002,18 @@ class TestESIPPProblem:
             max_number_parallel_arcs={}
         )
         
-        assert ipp.has_peak_total_assessments()
-        assert ipp.results["Problem"][0]["Number of constraints"] == 34
-        assert ipp.results["Problem"][0]["Number of variables"] == 24
-        assert ipp.results["Problem"][0]["Number of nonzeros"] == 77
+        if solver_name == 'scip':
+            assert len(ipp.instance.constr_arc_sos1) == 0
+            assert ipp.has_peak_total_assessments() # TODO: make sure this is true
+            assert ipp.results["Problem"][0]["Number of constraints"] == 0 # 34
+            assert ipp.results["Problem"][0]["Number of variables"] == 23 # 24
+            # assert ipp.results["Problem"][0]["Number of nonzeros"] == 77
+        else:
+            assert len(ipp.instance.constr_arc_sos1) == 0
+            assert ipp.has_peak_total_assessments() # TODO: make sure this is true
+            assert ipp.results["Problem"][0]["Number of constraints"] == 34
+            assert ipp.results["Problem"][0]["Number of variables"] == 24
+            assert ipp.results["Problem"][0]["Number of nonzeros"] == 77
     
         # *********************************************************************
         # *********************************************************************
@@ -1150,8 +1185,17 @@ class TestESIPPProblem:
     
     # *************************************************************************
     # *************************************************************************
-        
-    def test_nonisolated_undirected_network(self):
+    
+    @pytest.mark.parametrize(
+        "use_sos_arcs, arc_sos_weight_key", 
+        [(True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_NONE),
+         (True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_CAP),
+         (True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_COST),
+         (True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_SPEC_CAP),
+         (True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_SPEC_COST),
+         (False, None)]
+        )
+    def test_nonisolated_undirected_network(self, use_sos_arcs, arc_sos_weight_key):
         
         # scenario
         q = 0
@@ -1257,7 +1301,10 @@ class TestESIPPProblem:
     
         # no sos, regular time intervals
         ipp = build_solve_ipp(
+            solver='scip',
             solver_options={},
+            use_sos_arcs=use_sos_arcs, 
+            arc_sos_weight_key=arc_sos_weight_key,
             perform_analysis=False,
             plot_results=False,  # True,
             print_solver_output=False,
@@ -1268,10 +1315,19 @@ class TestESIPPProblem:
             max_number_parallel_arcs={}
         )
         
-        assert ipp.has_peak_total_assessments()
-        assert ipp.results["Problem"][0]["Number of constraints"] == 80
-        assert ipp.results["Problem"][0]["Number of variables"] == 84
-        assert ipp.results["Problem"][0]["Number of nonzeros"] == 253
+        # validation
+        if use_sos_arcs:
+            assert len(ipp.instance.constr_arc_sos1) == 3
+            assert ipp.has_peak_total_assessments()
+            assert ipp.results["Problem"][0]["Number of constraints"] == 0 # 80
+            assert ipp.results["Problem"][0]["Number of variables"] == 83 # 84
+            # assert ipp.results["Problem"][0]["Number of nonzeros"] == 253
+        else:
+            assert len(ipp.instance.constr_arc_sos1) == 0
+            assert ipp.has_peak_total_assessments()
+            assert ipp.results["Problem"][0]["Number of constraints"] == 0 # 80
+            assert ipp.results["Problem"][0]["Number of variables"] == 83 # 84
+            # assert ipp.results["Problem"][0]["Number of nonzeros"] == 253
     
         # **************************************************************************
     
@@ -1318,8 +1374,17 @@ class TestESIPPProblem:
         
     # *************************************************************************
     # *************************************************************************
-            
-    def test_nonisolated_undirected_network_diff_tech(self):
+    
+    @pytest.mark.parametrize(
+        "use_sos_arcs, arc_sos_weight_key", 
+        [(True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_NONE),
+         # (True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_CAP),
+         # (True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_COST),
+         # (True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_SPEC_CAP),
+         # (True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_SPEC_COST),
+         (False, None)]
+        )
+    def test_nonisolated_undirected_network_diff_tech(self, use_sos_arcs, arc_sos_weight_key):
         
         # scenario
         q = 0
@@ -1424,7 +1489,10 @@ class TestESIPPProblem:
     
         # no sos, regular time intervals
         ipp = build_solve_ipp(
+            solver='scip',
             solver_options={},
+            use_sos_arcs=use_sos_arcs, 
+            arc_sos_weight_key=arc_sos_weight_key,
             perform_analysis=False,
             plot_results=False,  # True,
             print_solver_output=False,
@@ -1435,10 +1503,25 @@ class TestESIPPProblem:
             max_number_parallel_arcs={}
         )
         
-        assert ipp.has_peak_total_assessments()
-        assert ipp.results["Problem"][0]["Number of constraints"] == 80
-        assert ipp.results["Problem"][0]["Number of variables"] == 84
-        assert ipp.results["Problem"][0]["Number of nonzeros"] == 253
+        # validation
+        # ipp.instance.constr_arc_sos1.pprint()
+        # print(ipp.results["Problem"][0])
+        if use_sos_arcs:
+            # print(ipp.results["Problem"][0])
+            assert len(ipp.instance.constr_arc_sos1) == 3
+            assert ipp.has_peak_total_assessments()
+            assert ipp.results["Problem"][0]["Number of constraints"] == 0 # should be 80
+            assert ipp.results["Problem"][0]["Number of variables"] == 83 # should be 84
+            # assert ipp.results["Problem"][0]["Number of nonzeros"] == 253
+            # check_problem_size(ipp, 0, 83, 253)
+        else:
+            # print(ipp.results["Problem"][0])
+            assert len(ipp.instance.constr_arc_sos1) == 0
+            assert ipp.has_peak_total_assessments()
+            assert ipp.results["Problem"][0]["Number of constraints"] == 0 # should be 80
+            assert ipp.results["Problem"][0]["Number of variables"] == 83 # should be 84
+            # assert ipp.results["Problem"][0]["Number of nonzeros"] == 253
+            # check_problem_size(ipp, 0, 83, 253)
     
         # *********************************************************************
         # *********************************************************************
@@ -1484,7 +1567,13 @@ class TestESIPPProblem:
     # *************************************************************************
     # *************************************************************************
     
-    def test_nonisolated_network_preexisting_directed_arcs(self):
+    @pytest.mark.parametrize(
+        "use_sos_arcs, arc_sos_weight_key", 
+        [(True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_NONE),
+         (True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_CAP),
+         (False, None)]
+        )
+    def test_nonisolated_network_preexisting_directed_arcs(self, use_sos_arcs, arc_sos_weight_key):
         
         # time frame
         q = 0
@@ -1586,7 +1675,10 @@ class TestESIPPProblem:
     
         # no sos, regular time intervals
         ipp = build_solve_ipp(
+            solver='scip',
             solver_options={},
+            use_sos_arcs=use_sos_arcs, 
+            arc_sos_weight_key=arc_sos_weight_key,
             perform_analysis=False,
             plot_results=False,  # True,
             print_solver_output=False,
@@ -1600,6 +1692,10 @@ class TestESIPPProblem:
         # *********************************************************************
     
         # validation
+        if use_sos_arcs:
+            assert len(ipp.instance.constr_arc_sos1) != 0
+        else:
+            assert len(ipp.instance.constr_arc_sos1) == 0
         # network is still isolated
         # the undirected arc was installed
         assert (
@@ -1621,7 +1717,20 @@ class TestESIPPProblem:
     # *************************************************************************
     # *************************************************************************
     
-    def test_nonisolated_network_preexisting_directed_arcs_diff_tech(self):
+    @pytest.mark.parametrize(
+        "use_sos_arcs, arc_sos_weight_key", 
+        [(True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_NONE),
+         (True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_CAP),
+         (True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_COST),
+         (True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_SPEC_CAP),
+         (True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_SPEC_COST),
+         (False, None)]
+        )
+    def test_nonisolated_network_preexisting_directed_arcs_diff_tech(
+            self, 
+            use_sos_arcs, 
+            arc_sos_weight_key
+            ):
         
         # time frame
         q = 0
@@ -1724,7 +1833,10 @@ class TestESIPPProblem:
     
         # no sos, regular time intervals
         ipp = build_solve_ipp(
-            solver_options={},solver='scip',
+            solver='scip',
+            solver_options={},
+            use_sos_arcs=use_sos_arcs, 
+            arc_sos_weight_key=arc_sos_weight_key,
             perform_analysis=False,
             plot_results=False,  # True,
             print_solver_output=False,
@@ -1738,6 +1850,10 @@ class TestESIPPProblem:
         # **************************************************************************
     
         # validation
+        if use_sos_arcs:
+            assert len(ipp.instance.constr_arc_sos1) != 0
+        else:
+            assert len(ipp.instance.constr_arc_sos1) == 0
         # the undirected arc should be installed since it is cheaper tham imp.
         assert (
             True
@@ -1973,6 +2089,7 @@ class TestESIPPProblem:
         # *********************************************************************
     
         # overview
+        assert len(ipp.instance.constr_arc_sos1) == 0
         
         (imports_qpk, 
          exports_qpk, 
@@ -2114,8 +2231,12 @@ class TestESIPPProblem:
     # *************************************************************************
     
     # TODO: test arc groups with sos
-    
-    def test_arc_groups_individual_ref(self):
+    @pytest.mark.parametrize(
+        "use_sos_arcs, arc_sos_weight_key", 
+        [(True, InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_CAP),
+         (False, None)]
+        )
+    def test_arc_groups_individual_ref(self, use_sos_arcs, arc_sos_weight_key):
         
         # time frame
         q = 0
@@ -2306,7 +2427,13 @@ class TestESIPPProblem:
             arc_groups_dict=arc_groups_dict
         )
     
-        # **************************************************************************
+        # *********************************************************************
+        
+        # validation
+        if use_sos_arcs:
+            assert len(ipp.instance.constr_arc_sos1) != 0
+        else:
+            assert len(ipp.instance.constr_arc_sos1) == 0
     
         # overview
         (imports_qpk, 
diff --git a/tests/test_solvers.py b/tests/test_solvers.py
index 44a28c8..9cb0eda 100644
--- a/tests/test_solvers.py
+++ b/tests/test_solvers.py
@@ -5,30 +5,99 @@ from src.topupopt.solvers.interface import SolverInterface
 from pyomo.opt.results.solver import TerminationCondition
 
 import pyomo.environ as pyo
-from pyomo.common.errors import ApplicationError
+from pyomo.opt import check_available_solvers
 
 import random
+import pytest
 
 # *****************************************************************************
 # *****************************************************************************
 
-# TODO: use pytest to parameterise the tests to use different solvers
-
+@pytest.mark.parametrize(
+    "solver", 
+    ['glpk',
+     'cbc',
+     'scip',
+     'fakesolver']
+    )
 class TestSolvers:
+    
     # *************************************************************************
     # *************************************************************************
 
-    def test_solver_factory_arguments(self):
-        # test a collection of problems using different solvers
+    def test_inexistent_solver(self, solver):
+
+        # try using a fake solver and a problem incompatible with another solver
+
+        # list of problems: one compatible, one incompatible
+        list_problems = [
+            problem_milp_feasible(20, seed_number=50),
+            problem_lp_optimal(),
+            problem_qp_optimal(),
+            problem_qp_optimal(),
+        ]
 
-        problem = self.problem_milp_feasible()
+        # problem types
+        list_problem_types = [
+            SolverInterface.PROBLEM_LP,
+            SolverInterface.PROBLEM_LP,
+            SolverInterface.PROBLEM_QP,
+            "fake_problem_type",
+        ]
 
         # solver settings
+        solver_timelimit = 30
+        solver_abs_mip_gap = 0
+        solver_rel_mip_gap = 0.01
+        solver_options = {
+            "time_limit": solver_timelimit,
+            "relative_mip_gap": solver_rel_mip_gap,
+            "absolute_mip_gap": solver_abs_mip_gap,
+        }
 
-        solver_timelimit = 10
+        # *********************************************************************
+        # *********************************************************************
 
-        solver_abs_mip_gap = 0.001
+        for index, problem in enumerate(list_problems):
+            
+            # optimise
+            try:
+                # test problem-solver compatibility
+                problem_type = list_problem_types[index]
+
+                if not SolverInterface.problem_and_solver_are_compatible(
+                        solver, problem_type):
+                    continue
+
+            except SolverInterface.UnknownSolverError:
+                assert True
+
+            except SolverInterface.UnknownProblemTypeError:
+                assert True
 
+            # test the solver interface
+            try:
+                # configure common solver interface
+                _ = SolverInterface(solver_name=solver, **solver_options)
+            except SolverInterface.UnknownSolverError:
+                assert True
+                    
+    # *************************************************************************
+    # *************************************************************************
+    
+    # @pytest.mark.skipif(True, reason="requires python3.10 or higher")
+    def test_solver_factory_arguments(self, solver):
+
+        # skip test if the solver is not available
+        if not bool(check_available_solvers(solver)):
+            return
+        
+        # test a feasible problem
+        problem = problem_milp_feasible()
+
+        # solver settings
+        solver_timelimit = 10
+        solver_abs_mip_gap = 0.001
         solver_rel_mip_gap = 0.01
 
         solver_options = {
@@ -36,61 +105,30 @@ class TestSolvers:
             "relative_mip_gap": solver_rel_mip_gap,
             "absolute_mip_gap": solver_abs_mip_gap,
             # special option
-            "tee": True,
+            "tee": False,
         }
-
-        solver_name = "scip"
-
-        results, solver_interface = self.optimise(solver_name, solver_options, problem)
+        results, solver_interface = optimise(solver, solver_options, problem)
 
     # *************************************************************************
     # *************************************************************************
 
-    def test_problems(self):
+    def test_problems(self, solver):
         # test a collection of problems using different solvers
 
-        # solver = "scip"
-        # # scip_exec_path = '/usr/bin/scip'
-        # # solver_options = {'executable': scip_exec_path}
-        # solver_options = {}
-        
-        # solver = 'cplex'
-        # # cplex_exec_path = '/home/pmlpm/Software/CPLEX/cplex/bin/x86-64_linux/cplex'
-        # cplex_exec_path = '/home/pmlpm/CPLEX/cplex/bin/x86-64_linux/cplex'
-        # #solver_options = {}
-        # solver_options = {'executable':cplex_exec_path}
-
-        list_solvers = [
-            "fake_solver",
-            "cbc",
-            "glpk",
-            "scip",
-            'cplex'
-        ]
-
-        list_solver_options = [
-            None,  # fake
-            None,  # cbc
-            {"tee": False},  # glpk
-            None,  # scip {'executable': scip_exec_path},  
-            None, # cplex
-            # {'executable': cplex_exec_path},
-        ]
-
         # list of problems
 
         list_concrete_models = [
-            self.problem_qp_optimal(),
-            self.problem_qp_infeasible(),
-            self.problem_lp_unbounded(),
-            self.problem_lp_infeasible(),
-            self.problem_lp_optimal(),
-            self.problem_milp_unbounded(),
-            self.problem_milp_infeasible(),
-            self.problem_milp_optimal(),
-            self.problem_milp_feasible(),
-            self.problem_milp_feasible(15, 64),
-            self.problem_milp_feasible(10, 46),
+            problem_qp_optimal(),
+            problem_qp_infeasible(),
+            problem_lp_unbounded(),
+            problem_lp_infeasible(),
+            problem_lp_optimal(),
+            problem_milp_unbounded(),
+            problem_milp_infeasible(),
+            problem_milp_optimal(),
+            problem_milp_feasible(),
+            problem_milp_feasible(15, 64),
+            problem_milp_feasible(10, 46),
         ]
 
         # list of problem types
@@ -139,577 +177,475 @@ class TestSolvers:
             True,
         ]
 
-        # list of solvers
-
-        list_solvers = ["fake_solver", "cbc", "glpk", "scip", "cplex"]
-
         # solver settings
-
         solver_timelimit = 10
-
         solver_abs_mip_gap = 0.001
-
         solver_rel_mip_gap = 0.01
+        solver_options = {
+            "time_limit": solver_timelimit,
+            "relative_mip_gap": solver_rel_mip_gap,
+            "absolute_mip_gap": solver_abs_mip_gap,
+        }
 
-        for solver_name, solver_options in zip(list_solvers, list_solver_options):
-            if type(solver_options) == dict:
-                solver_options.update(
-                    {
-                        "time_limit": solver_timelimit,
-                        "relative_mip_gap": solver_rel_mip_gap,
-                        "absolute_mip_gap": solver_abs_mip_gap,
-                    }
-                )
-
-            else:
-                solver_options = {
-                    "time_limit": solver_timelimit,
-                    "relative_mip_gap": solver_rel_mip_gap,
-                    "absolute_mip_gap": solver_abs_mip_gap,
-                }
-
-            for problem_index, problem in enumerate(list_concrete_models):
-                try:
-                    # check problem and solver compatibility
-
-                    problem_type = list_problem_types[problem_index]
-
-                    if (
-                        SolverInterface.problem_and_solver_are_compatible(
-                            solver_name, problem_type
-                        )
-                        == False
-                    ):
-                        continue
-
-                    # optimise
-
-                    results, solver_interface = self.optimise(
-                        solver_name, solver_options, problem, print_solver_output=False
-                    )
-
-                except SolverInterface.UnknownSolverError:
-                    continue
-
-                except SolverInterface.UnknownProblemTypeError:
-                    continue
-
-                # *************************************************************
-                # *************************************************************
-
-                # termination condition
-
-                exp_term_cond = list_problem_termination_conditions[problem_index]
-
-                term_cond = results.solver.termination_condition
-
-                if (
-                    exp_term_cond == None
-                    or (
-                        solver_name == "glpk"
-                        and exp_term_cond == TerminationCondition.unbounded
-                    )
-                    or (
-                        solver_name == "cplex"
-                        and exp_term_cond == TerminationCondition.unbounded
-                    )
-                    or (
-                        solver_name == "cplex"
-                        and exp_term_cond == TerminationCondition.optimal
-                    )
-                    or (
-                        solver_name == "cplex"
-                        and exp_term_cond == TerminationCondition.infeasible
-                    )
-                ):
-                    # exceptions in need of correction
-
-                    pass
-
-                else:
-                    # print(solver_name)
-                    # print(results)
-                    assert exp_term_cond == term_cond
-
-                # *************************************************************
-                # *************************************************************
-
-                # solver status
-
-                if (
-                    (
-                        solver_name == "glpk"
-                        and term_cond == TerminationCondition.infeasible
-                    )
-                    or (
-                        solver_name == "cplex"
-                        and term_cond == TerminationCondition.unknown
-                    )
-                    or (
-                        solver_name == "cplex"
-                        and exp_term_cond == TerminationCondition.unbounded
-                    )
-                    or (
-                        solver_name == "cplex"
-                        and exp_term_cond == TerminationCondition.infeasible
-                    )
-                ):
-                    pass
-
-                else:
-                    # check if the solver status matches the one one would expect
-                    # if the termination condition was correct
-
-                    assert (
-                        TerminationCondition.to_solver_status(term_cond)
-                        == results.solver.status
-                    )
-
-                    # if valid, it means the results object is coherent
+        for problem_index, problem in enumerate(list_concrete_models):
+            try:
+                # check problem and solver compatibility
 
-                # *************************************************************
-                # *************************************************************
+                problem_type = list_problem_types[problem_index]
 
                 if (
-                    exp_term_cond == None
-                    or (
-                        solver_name == "glpk"
-                        and exp_term_cond == TerminationCondition.unbounded
-                    )
-                    or (
-                        solver_name == "glpk"
-                        and exp_term_cond == TerminationCondition.infeasible
-                    )
-                    or (
-                        solver_name == "cplex"
-                        and exp_term_cond == TerminationCondition.unknown
-                    )
-                    or (
-                        solver_name == "cplex"
-                        and exp_term_cond == TerminationCondition.unbounded
-                    )
-                    or (
-                        solver_name == "cplex"
-                        and exp_term_cond == TerminationCondition.infeasible
+                    SolverInterface.problem_and_solver_are_compatible(
+                        solver, problem_type
                     )
+                    == False
                 ):
-                    pass
+                    continue
 
-                else:
-                    # check if the solver status matches the one one would expect
-                    # if the termination condition predicted was obtained
+                # optimise
 
-                    assert (
-                        TerminationCondition.to_solver_status(exp_term_cond)
-                        == results.solver.status
-                    )
+                results, solver_interface = optimise(
+                    solver, solver_options, problem, print_solver_output=False
+                )
 
-                    # if valid, the solver status is correct despite other issues
+            except SolverInterface.UnknownSolverError:
+                continue
 
-                # *************************************************************
-                # *************************************************************
+            except SolverInterface.UnknownProblemTypeError:
+                continue
 
-                # make sure the optimisation went as expected
+            # *************************************************************
+            # *************************************************************
 
-                exp_optim_result = list_problem_optimisation_sucess[problem_index]
+            # termination condition
 
-                if (
-                    TerminationCondition.to_solver_status(
-                        results.solver.termination_condition
-                    )
-                    != results.solver.status
-                ):
-                    # this can be removed once the aforementioned issues have
-                    # been fixed (e.g. for the cplex and glpk solvers)
+            exp_term_cond = list_problem_termination_conditions[problem_index]
 
-                    pass
+            term_cond = results.solver.termination_condition
 
-                else:
-                    optim_result = solver_interface.was_optimisation_sucessful(
-                        results, problem_type
-                    )
-
-                # *************************************************************
-                # *************************************************************
+            if (
+                exp_term_cond == None
+                or (
+                    solver == "glpk"
+                    and exp_term_cond == TerminationCondition.unbounded
+                )
+                or (
+                    solver == "cplex"
+                    and exp_term_cond == TerminationCondition.unbounded
+                )
+                or (
+                    solver == "cplex"
+                    and exp_term_cond == TerminationCondition.optimal
+                )
+                or (
+                    solver == "cplex"
+                    and exp_term_cond == TerminationCondition.infeasible
+                )
+            ):
+                # exceptions in need of correction
 
-                if (
-                    TerminationCondition.to_solver_status(
-                        results.solver.termination_condition
-                    )
-                    != results.solver.status
-                    or exp_term_cond == TerminationCondition.unbounded
-                ):
-                    # this can be removed once the aforementioned issues have
-                    # been fixed (e.g. for the cplex and glpk solvers)
+                pass
 
-                    pass
+            else:
+                # print(solver)
+                # print(results)
+                assert exp_term_cond == term_cond
 
-                else:
-                    assert optim_result == exp_optim_result
+            # *************************************************************
+            # *************************************************************
 
-                # *************************************************************
-                # *************************************************************
+            # solver status
 
-                # test additional scenarios
+            if (
+                (
+                    solver == "glpk"
+                    and term_cond == TerminationCondition.infeasible
+                )
+                or (
+                    solver == "cplex"
+                    and term_cond == TerminationCondition.unknown
+                )
+                or (
+                    solver == "cplex"
+                    and exp_term_cond == TerminationCondition.unbounded
+                )
+                or (
+                    solver == "cplex"
+                    and exp_term_cond == TerminationCondition.infeasible
+                )
+            ):
+                pass
 
-                if optim_result == False:
-                    continue
+            else:
+                # check if the solver status matches the one one would expect
+                # if the termination condition was correct
 
-                # force unknown solver status error
+                assert (
+                    TerminationCondition.to_solver_status(term_cond)
+                    == results.solver.status
+                )
 
-                results.solver.status = "false_solver_status"
+                # if valid, it means the results object is coherent
 
-                try:
-                    _ = solver_interface.was_optimisation_sucessful(
-                        results, problem_type
-                    )
+            # *************************************************************
+            # *************************************************************
 
-                except solver_interface.UnknownSolverStatusError:
-                    assert True
+            if (
+                exp_term_cond == None
+                or (
+                    solver == "glpk"
+                    and exp_term_cond == TerminationCondition.unbounded
+                )
+                or (
+                    solver == "glpk"
+                    and exp_term_cond == TerminationCondition.infeasible
+                )
+                or (
+                    solver == "cplex"
+                    and exp_term_cond == TerminationCondition.unknown
+                )
+                or (
+                    solver == "cplex"
+                    and exp_term_cond == TerminationCondition.unbounded
+                )
+                or (
+                    solver == "cplex"
+                    and exp_term_cond == TerminationCondition.infeasible
+                )
+            ):
+                pass
 
-                # force unknown termination condition error
+            else:
+                # check if the solver status matches the one one would expect
+                # if the termination condition predicted was obtained
 
-                results.solver.termination_condition = "false_termin_condition"
+                assert (
+                    TerminationCondition.to_solver_status(exp_term_cond)
+                    == results.solver.status
+                )
 
-                try:
-                    _ = solver_interface.was_optimisation_sucessful(
-                        results, problem_type
-                    )
+                # if valid, the solver status is correct despite other issues
 
-                except solver_interface.UnknownTerminationConditionError:
-                    assert True
+            # *************************************************************
+            # *************************************************************
 
-                # force an InconsistentSolverStatusError
+            # make sure the optimisation went as expected
 
-                results.solver.termination_condition = TerminationCondition.optimal
+            exp_optim_result = list_problem_optimisation_sucess[problem_index]
 
-                results.solver.status = TerminationCondition.to_solver_status(
+            if (
+                TerminationCondition.to_solver_status(
                     results.solver.termination_condition
                 )
+                != results.solver.status
+            ):
+                # this can be removed once the aforementioned issues have
+                # been fixed (e.g. for the cplex and glpk solvers)
 
-                results.solver.termination_condition = TerminationCondition.unknown
+                pass
 
-                try:
-                    _ = solver_interface.was_optimisation_sucessful(
-                        results, problem_type
-                    )
+            else:
+                optim_result = solver_interface.was_optimisation_sucessful(
+                    results, problem_type
+                )
 
-                except solver_interface.InconsistentSolverStatusError:
-                    assert True
+            # *************************************************************
+            # *************************************************************
 
-                # force an InconsistentProblemTypeAndSolverError
+            if (
+                TerminationCondition.to_solver_status(
+                    results.solver.termination_condition
+                )
+                != results.solver.status
+                or exp_term_cond == TerminationCondition.unbounded
+            ):
+                # this can be removed once the aforementioned issues have
+                # been fixed (e.g. for the cplex and glpk solvers)
 
-                if problem_type == SolverInterface.PROBLEM_LP and solver_name == "glpk":
-                    problem_type = SolverInterface.PROBLEM_QP
+                pass
 
-                    try:
-                        _ = solver_interface.was_optimisation_sucessful(
-                            results, problem_type
-                        )
+            else:
+                assert optim_result == exp_optim_result
 
-                    except solver_interface.InconsistentProblemTypeAndSolverError:
-                        assert True
+            # *************************************************************
+            # *************************************************************
 
-        # *********************************************************************
-        # *********************************************************************
+            # test additional scenarios
 
-    # *************************************************************************
-    # *************************************************************************
+            if optim_result == False:
+                continue
 
-    # carry out optimisations
+            # force unknown solver status error
 
-    def optimise(
-        self,
-        solver_name: str,
-        solver_options: dict,
-        # solver_interface: SolverInterface,
-        problem: pyo.ConcreteModel,
-        print_solver_output: bool = True,
-    ):
-        # configure common solver interface
-        solver_interface = SolverInterface(solver_name=solver_name, **solver_options)
+            results.solver.status = "false_solver_status"
 
-        # get the solver handler
-        solver_handler = solver_interface.get_solver_handler(**solver_options)
+            try:
+                _ = solver_interface.was_optimisation_sucessful(
+                    results, problem_type
+                )
 
-        # solve
-        if "tee" not in solver_options:
-            results = solver_handler.solve(problem, tee=print_solver_output)
-        else:
-            results = solver_handler.solve(problem)
+            except solver_interface.UnknownSolverStatusError:
+                assert True
 
-        # return
-        return results, solver_interface
+            # force unknown termination condition error
 
-    # *************************************************************************
-    # *************************************************************************
+            results.solver.termination_condition = "false_termin_condition"
 
-    def problem_qp_optimal(self):
-        model = pyo.ConcreteModel("qp_optimal")
+            try:
+                _ = solver_interface.was_optimisation_sucessful(
+                    results, problem_type
+                )
 
-        model.x = pyo.Var(within=pyo.NonNegativeReals)
-        model.y = pyo.Var(within=pyo.NonNegativeReals)
+            except solver_interface.UnknownTerminationConditionError:
+                assert True
 
-        def constraint_rule(model):
-            return model.x + model.y >= 10
+            # force an InconsistentSolverStatusError
 
-        model.constraint = pyo.Constraint(rule=constraint_rule)
+            results.solver.termination_condition = TerminationCondition.optimal
 
-        def objective_rule(model):
-            return (
-                model.x
-                + model.y
-                + 0.5
-                * (model.x * model.x + 4 * model.x * model.y + 7 * model.y * model.y)
+            results.solver.status = TerminationCondition.to_solver_status(
+                results.solver.termination_condition
             )
 
-        model.objective = pyo.Objective(rule=objective_rule, sense=pyo.minimize)
-
-        return model
-
-    # *************************************************************************
-    # *************************************************************************
-
-    def problem_qp_infeasible(self):
-        model = pyo.ConcreteModel("qp_infeasible")
+            results.solver.termination_condition = TerminationCondition.unknown
 
-        # model.x = pyo.Var(within=pyo.NonNegativeReals, bounds=(0,5))
-        # model.y = pyo.Var(within=pyo.NonNegativeReals, bounds=(0,4))
+            try:
+                _ = solver_interface.was_optimisation_sucessful(
+                    results, problem_type
+                )
 
-        model.x = pyo.Var(bounds=(0, 5))
-        model.y = pyo.Var(bounds=(0, 4))
+            except solver_interface.InconsistentSolverStatusError:
+                assert True
 
-        def constraint_rule(model):
-            return model.x + model.y >= 10
+            # force an InconsistentProblemTypeAndSolverError
 
-        model.constraint = pyo.Constraint(rule=constraint_rule)
+            if problem_type == SolverInterface.PROBLEM_LP and solver == "glpk":
+                problem_type = SolverInterface.PROBLEM_QP
 
-        def objective_rule(model):
-            return (
-                model.x
-                + model.y
-                + 0.5
-                * (model.x * model.x + 4 * model.x * model.y + 7 * model.y * model.y)
-            )
+                try:
+                    _ = solver_interface.was_optimisation_sucessful(
+                        results, problem_type
+                    )
 
-        model.objective = pyo.Objective(rule=objective_rule, sense=pyo.minimize)
+                except solver_interface.InconsistentProblemTypeAndSolverError:
+                    assert True
 
-        return model
+        # *********************************************************************
+        # *********************************************************************
 
     # *************************************************************************
     # *************************************************************************
 
-    def problem_lp_optimal(self):
-        model = pyo.ConcreteModel("lp_optimal")
-
-        model.x = pyo.Var([1, 2], domain=pyo.NonNegativeReals)
+# *****************************************************************************
+# *****************************************************************************
 
-        model.OBJ = pyo.Objective(expr=2 * model.x[1] + 3 * model.x[2])
+# carry out optimisations
 
-        model.Constraint1 = pyo.Constraint(expr=3 * model.x[1] + 4 * model.x[2] >= 1)
+def optimise(
+    solver_name: str,
+    solver_options: dict,
+    # solver_interface: SolverInterface,
+    problem: pyo.ConcreteModel,
+    print_solver_output: bool = True,
+):
+    # configure common solver interface
+    solver_interface = SolverInterface(solver_name=solver_name, **solver_options)
 
-        return model
+    # get the solver handler
+    solver_handler = solver_interface.get_solver_handler(**solver_options)
 
-    # *************************************************************************
-    # *************************************************************************
+    # solve
+    if "tee" not in solver_options:
+        results = solver_handler.solve(problem, tee=print_solver_output)
+    else:
+        results = solver_handler.solve(problem)
 
-    def problem_lp_infeasible(self):
-        model = pyo.ConcreteModel("lp_infeasible")
+    # return
+    return results, solver_interface
 
-        model.x = pyo.Var([1, 2], domain=pyo.NonNegativeReals)
 
-        model.OBJ = pyo.Objective(expr=2 * model.x[1] + 3 * model.x[2])
+# *****************************************************************************
+# *****************************************************************************
 
-        model.Constraint1 = pyo.Constraint(expr=3 * model.x[1] + 4 * model.x[2] <= -1)
+def problem_qp_optimal():
+    model = pyo.ConcreteModel("qp_optimal")
 
-        return model
+    model.x = pyo.Var(within=pyo.NonNegativeReals)
+    model.y = pyo.Var(within=pyo.NonNegativeReals)
 
-    # *************************************************************************
-    # *************************************************************************
+    def constraint_rule(model):
+        return model.x + model.y >= 10
 
-    def problem_lp_unbounded(self):
-        model = pyo.ConcreteModel("lp_unbounded")
+    model.constraint = pyo.Constraint(rule=constraint_rule)
 
-        model.x = pyo.Var([1, 2], domain=pyo.NonNegativeReals)
-
-        model.OBJ = pyo.Objective(
-            expr=2 * model.x[1] + 3 * model.x[2], sense=pyo.maximize
+    def objective_rule(model):
+        return (
+            model.x
+            + model.y
+            + 0.5
+            * (model.x * model.x + 4 * model.x * model.y + 7 * model.y * model.y)
         )
 
-        model.Constraint1 = pyo.Constraint(expr=3 * model.x[1] + 4 * model.x[2] >= 1)
+    model.objective = pyo.Objective(rule=objective_rule, sense=pyo.minimize)
 
-        return model
+    return model
 
-    # *************************************************************************
-    # *************************************************************************
+# *****************************************************************************
+# *****************************************************************************
 
-    def problem_milp_optimal(self):
-        model = pyo.ConcreteModel("milp_optimal")
+def problem_qp_infeasible():
+    model = pyo.ConcreteModel("qp_infeasible")
 
-        model.x = pyo.Var([1, 2], domain=pyo.Binary)
+    # model.x = pyo.Var(within=pyo.NonNegativeReals, bounds=(0,5))
+    # model.y = pyo.Var(within=pyo.NonNegativeReals, bounds=(0,4))
 
-        model.OBJ = pyo.Objective(expr=2.15 * model.x[1] + 3.8 * model.x[2])
+    model.x = pyo.Var(bounds=(0, 5))
+    model.y = pyo.Var(bounds=(0, 4))
 
-        model.Constraint1 = pyo.Constraint(expr=3 * model.x[1] + 4 * model.x[2] >= 1)
+    def constraint_rule(model):
+        return model.x + model.y >= 10
 
-        return model
+    model.constraint = pyo.Constraint(rule=constraint_rule)
 
-    # *************************************************************************
-    # *************************************************************************
+    def objective_rule(model):
+        return (
+            model.x
+            + model.y
+            + 0.5
+            * (model.x * model.x + 4 * model.x * model.y + 7 * model.y * model.y)
+        )
 
-    def problem_milp_infeasible(self):
-        model = pyo.ConcreteModel("milp_infeasible")
+    model.objective = pyo.Objective(rule=objective_rule, sense=pyo.minimize)
 
-        model.x = pyo.Var([1, 2], domain=pyo.Binary)
+    return model
 
-        model.OBJ = pyo.Objective(expr=2 * model.x[1] + 3 * model.x[2])
+# *****************************************************************************
+# *****************************************************************************
 
-        model.Constraint1 = pyo.Constraint(expr=3 * model.x[1] + 4 * model.x[2] <= -1)
+def problem_lp_optimal():
+    model = pyo.ConcreteModel("lp_optimal")
 
-        return model
+    model.x = pyo.Var([1, 2], domain=pyo.NonNegativeReals)
 
-    # *************************************************************************
-    # *************************************************************************
+    model.OBJ = pyo.Objective(expr=2 * model.x[1] + 3 * model.x[2])
 
-    def problem_milp_unbounded(self):
-        model = pyo.ConcreteModel("milp_unbounded")
+    model.Constraint1 = pyo.Constraint(expr=3 * model.x[1] + 4 * model.x[2] >= 1)
 
-        model.x = pyo.Var([1, 2], domain=pyo.NonNegativeReals)
+    return model
 
-        model.y = pyo.Var(domain=pyo.Binary)
+# *****************************************************************************
+# *****************************************************************************
 
-        model.OBJ = pyo.Objective(
-            expr=2 * model.x[1] + 3 * model.x[2] + model.y, sense=pyo.maximize
-        )
+def problem_lp_infeasible():
+    model = pyo.ConcreteModel("lp_infeasible")
 
-        model.Constraint1 = pyo.Constraint(expr=3 * model.x[1] + 4 * model.x[2] >= 1)
+    model.x = pyo.Var([1, 2], domain=pyo.NonNegativeReals)
 
-        return model
+    model.OBJ = pyo.Objective(expr=2 * model.x[1] + 3 * model.x[2])
 
-    # *************************************************************************
-    # *************************************************************************
+    model.Constraint1 = pyo.Constraint(expr=3 * model.x[1] + 4 * model.x[2] <= -1)
 
-    def problem_milp_feasible(self, number_binary_variables=25, seed_number=None):
-        if seed_number != None:
-            random.seed(seed_number)
+    return model
 
-        model = pyo.ConcreteModel("milp_feasible")
+# *****************************************************************************
+# *****************************************************************************
 
-        # a knapsack-type problem
+def problem_lp_unbounded():
+    model = pyo.ConcreteModel("lp_unbounded")
 
-        model.Y = pyo.RangeSet(number_binary_variables)
+    model.x = pyo.Var([1, 2], domain=pyo.NonNegativeReals)
 
-        model.y = pyo.Var(model.Y, domain=pyo.Binary)
+    model.OBJ = pyo.Objective(
+        expr=2 * model.x[1] + 3 * model.x[2], sense=pyo.maximize
+    )
 
-        model.OBJ = pyo.Objective(
-            expr=sum(model.y[j] * random.random() for j in model.Y), sense=pyo.maximize
-        )
+    model.Constraint1 = pyo.Constraint(expr=3 * model.x[1] + 4 * model.x[2] >= 1)
 
-        model.Constraint1 = pyo.Constraint(
-            expr=sum(model.y[j] * random.random() for j in model.Y)
-            <= round(number_binary_variables / 5)
-        )
+    return model
 
-        def rule_c1(m, i):
-            return (
-                sum(
-                    model.y[j] * (random.random() - 0.5)
-                    for j in model.Y
-                    if j != i
-                    if random.randint(0, 1)
-                )
-                <= round(number_binary_variables / 5) * model.y[i]
-            )
+# *****************************************************************************
+# *****************************************************************************
 
-        model.constr_c1 = pyo.Constraint(model.Y, rule=rule_c1)
+def problem_milp_optimal():
+    model = pyo.ConcreteModel("milp_optimal")
 
-        return model
+    model.x = pyo.Var([1, 2], domain=pyo.Binary)
 
-    # *************************************************************************
-    # *************************************************************************
+    model.OBJ = pyo.Objective(expr=2.15 * model.x[1] + 3.8 * model.x[2])
 
-    def test_inexistent_solver(self):
-        fake_solver = "fake_solver"
-        good_solver = "glpk"
-        # solver_options: dict = None
+    model.Constraint1 = pyo.Constraint(expr=3 * model.x[1] + 4 * model.x[2] >= 1)
 
-        # try using a fake solver and a problem incompatible with another solver
+    return model
 
-        # list of problems: one compatible, one incompatible
+# *****************************************************************************
+# *****************************************************************************
 
-        list_problems = [
-            self.problem_milp_feasible(20, seed_number=50),
-            self.problem_lp_optimal(),
-            self.problem_qp_optimal(),
-            self.problem_qp_optimal(),
-        ]
+def problem_milp_infeasible():
+    model = pyo.ConcreteModel("milp_infeasible")
 
-        # problem types
+    model.x = pyo.Var([1, 2], domain=pyo.Binary)
 
-        list_problem_types = [
-            SolverInterface.PROBLEM_LP,
-            SolverInterface.PROBLEM_LP,
-            SolverInterface.PROBLEM_QP,
-            "fake_problem_type",
-        ]
+    model.OBJ = pyo.Objective(expr=2 * model.x[1] + 3 * model.x[2])
 
-        # list of solvers: one fake, one real
+    model.Constraint1 = pyo.Constraint(expr=3 * model.x[1] + 4 * model.x[2] <= -1)
 
-        list_solvers = [fake_solver, good_solver]
+    return model
 
-        # solver settings
+# *****************************************************************************
+# *****************************************************************************
 
-        solver_timelimit = 30
+def problem_milp_unbounded():
+    model = pyo.ConcreteModel("milp_unbounded")
 
-        solver_abs_mip_gap = 0
+    model.x = pyo.Var([1, 2], domain=pyo.NonNegativeReals)
 
-        solver_rel_mip_gap = 0.01
+    model.y = pyo.Var(domain=pyo.Binary)
 
-        solver_options = {
-            "time_limit": solver_timelimit,
-            "relative_mip_gap": solver_rel_mip_gap,
-            "absolute_mip_gap": solver_abs_mip_gap,
-        }
+    model.OBJ = pyo.Objective(
+        expr=2 * model.x[1] + 3 * model.x[2] + model.y, sense=pyo.maximize
+    )
 
-        # *********************************************************************
-        # *********************************************************************
+    model.Constraint1 = pyo.Constraint(expr=3 * model.x[1] + 4 * model.x[2] >= 1)
 
-        for solver_name in list_solvers:
-            for index, problem in enumerate(list_problems):
-                # optimise
+    return model
 
-                try:
-                    # test problem-solver compatibility
+# *****************************************************************************
+# *****************************************************************************
 
-                    problem_type = list_problem_types[index]
+def problem_milp_feasible(number_binary_variables=25, seed_number=None):
+    if seed_number != None:
+        random.seed(seed_number)
 
-                    if (
-                        SolverInterface.problem_and_solver_are_compatible(
-                            solver_name, problem_type
-                        )
-                        == False
-                    ):
-                        continue
+    model = pyo.ConcreteModel("milp_feasible")
 
-                except SolverInterface.UnknownSolverError:
-                    assert True
+    # a knapsack-type problem
 
-                except SolverInterface.UnknownProblemTypeError:
-                    assert True
+    model.Y = pyo.RangeSet(number_binary_variables)
 
-                # test the solver interface
+    model.y = pyo.Var(model.Y, domain=pyo.Binary)
 
-                try:
-                    # configure common solver interface
+    model.OBJ = pyo.Objective(
+        expr=sum(model.y[j] * random.random() for j in model.Y), sense=pyo.maximize
+    )
 
-                    _ = SolverInterface(solver_name=solver_name, **solver_options)
+    model.Constraint1 = pyo.Constraint(
+        expr=sum(model.y[j] * random.random() for j in model.Y)
+        <= round(number_binary_variables / 5)
+    )
 
-                except SolverInterface.UnknownSolverError:
-                    assert True
+    def rule_c1(m, i):
+        return (
+            sum(
+                model.y[j] * (random.random() - 0.5)
+                for j in model.Y
+                if j != i
+                if random.randint(0, 1)
+            )
+            <= round(number_binary_variables / 5) * model.y[i]
+        )
 
-    # **************************************************************************
-    # **************************************************************************
+    model.constr_c1 = pyo.Constraint(model.Y, rule=rule_c1)
 
+    return model
 
-# ******************************************************************************
-# ******************************************************************************
+# *****************************************************************************
+# *****************************************************************************
-- 
GitLab