test_solvers.py

# imports

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

import random

# *****************************************************************************
# *****************************************************************************


class TestSolvers:
    # *************************************************************************
    # *************************************************************************

    def test_solver_factory_arguments(self):
        # test a collection of problems using different solvers

        problem = self.problem_milp_feasible()

        # 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,
            # special option
            "tee": True,
        }

        solver_name = "scip"

        results, solver_interface = self.optimise(solver_name, solver_options, problem)

    # *************************************************************************
    # *************************************************************************

    def test_problems(self):
        # 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),
        ]

        # list of problem types

        list_problem_types = [
            SolverInterface.PROBLEM_QP,
            SolverInterface.PROBLEM_QP,
            SolverInterface.PROBLEM_LP,
            SolverInterface.PROBLEM_LP,
            SolverInterface.PROBLEM_LP,
            SolverInterface.PROBLEM_MILP,
            SolverInterface.PROBLEM_MILP,
            SolverInterface.PROBLEM_MILP,
            SolverInterface.PROBLEM_MILP,
            SolverInterface.PROBLEM_MILP,
            "unknown_problem_type",
        ]

        # expected

        list_problem_termination_conditions = [
            TerminationCondition.optimal,
            TerminationCondition.infeasible,
            TerminationCondition.unbounded,
            TerminationCondition.infeasible,
            TerminationCondition.optimal,
            TerminationCondition.unbounded,
            TerminationCondition.infeasible,
            TerminationCondition.optimal,
            None,  # if we don't know what to expect,
            None,  # if we don't know what to expect,
            None,  # if we don't know what to expect
        ]

        list_problem_optimisation_sucess = [
            True,
            True,
            False,
            False,
            True,
            False,
            False,
            True,
            True,
            True,
            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

        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

                # *************************************************************
                # *************************************************************

                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
                    )
                ):
                    pass

                else:
                    # check if the solver status matches the one one would expect
                    # if the termination condition predicted was obtained

                    assert (
                        TerminationCondition.to_solver_status(exp_term_cond)
                        == results.solver.status
                    )

                    # if valid, the solver status is correct despite other issues

                # *************************************************************
                # *************************************************************

                # make sure the optimisation went as expected

                exp_optim_result = list_problem_optimisation_sucess[problem_index]

                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)

                    pass

                else:
                    optim_result = solver_interface.was_optimisation_sucessful(
                        results, problem_type
                    )

                # *************************************************************
                # *************************************************************

                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

                else:
                    assert optim_result == exp_optim_result

                # *************************************************************
                # *************************************************************

                # test additional scenarios

                if optim_result == False:
                    continue

                # force unknown solver status error

                results.solver.status = "false_solver_status"

                try:
                    _ = solver_interface.was_optimisation_sucessful(
                        results, problem_type
                    )

                except solver_interface.UnknownSolverStatusError:
                    assert True

                # force unknown termination condition error

                results.solver.termination_condition = "false_termin_condition"

                try:
                    _ = solver_interface.was_optimisation_sucessful(
                        results, problem_type
                    )

                except solver_interface.UnknownTerminationConditionError:
                    assert True

                # force an InconsistentSolverStatusError

                results.solver.termination_condition = TerminationCondition.optimal

                results.solver.status = TerminationCondition.to_solver_status(
                    results.solver.termination_condition
                )

                results.solver.termination_condition = TerminationCondition.unknown

                try:
                    _ = solver_interface.was_optimisation_sucessful(
                        results, problem_type
                    )

                except solver_interface.InconsistentSolverStatusError:
                    assert True

                # force an InconsistentProblemTypeAndSolverError

                if problem_type == SolverInterface.PROBLEM_LP and solver_name == "glpk":
                    problem_type = SolverInterface.PROBLEM_QP

                    try:
                        _ = solver_interface.was_optimisation_sucessful(
                            results, problem_type
                        )

                    except solver_interface.InconsistentProblemTypeAndSolverError:
                        assert True

        # *********************************************************************
        # *********************************************************************

    # *************************************************************************
    # *************************************************************************

    # carry out optimisations

    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)

        # 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)

        # return
        return results, solver_interface

    # *************************************************************************
    # *************************************************************************

    def problem_qp_optimal(self):
        model = pyo.ConcreteModel("qp_optimal")

        model.x = pyo.Var(within=pyo.NonNegativeReals)
        model.y = pyo.Var(within=pyo.NonNegativeReals)

        def constraint_rule(model):
            return model.x + model.y >= 10

        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)
            )

        model.objective = pyo.Objective(rule=objective_rule, sense=pyo.minimize)

        return model

    # *************************************************************************
    # *************************************************************************

    def problem_qp_infeasible(self):
        model = pyo.ConcreteModel("qp_infeasible")

        # model.x = pyo.Var(within=pyo.NonNegativeReals, bounds=(0,5))
        # model.y = pyo.Var(within=pyo.NonNegativeReals, bounds=(0,4))

        model.x = pyo.Var(bounds=(0, 5))
        model.y = pyo.Var(bounds=(0, 4))

        def constraint_rule(model):
            return model.x + model.y >= 10

        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)
            )

        model.objective = pyo.Objective(rule=objective_rule, sense=pyo.minimize)

        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])

        model.Constraint1 = pyo.Constraint(expr=3 * model.x[1] + 4 * model.x[2] >= 1)

        return model

    # *************************************************************************
    # *************************************************************************

    def problem_lp_infeasible(self):
        model = pyo.ConcreteModel("lp_infeasible")

        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)

        return model

    # *************************************************************************
    # *************************************************************************

    def problem_lp_unbounded(self):
        model = pyo.ConcreteModel("lp_unbounded")

        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
        )

        model.Constraint1 = pyo.Constraint(expr=3 * model.x[1] + 4 * model.x[2] >= 1)

        return model

    # *************************************************************************
    # *************************************************************************

    def problem_milp_optimal(self):
        model = pyo.ConcreteModel("milp_optimal")

        model.x = pyo.Var([1, 2], domain=pyo.Binary)

        model.OBJ = pyo.Objective(expr=2.15 * model.x[1] + 3.8 * model.x[2])

        model.Constraint1 = pyo.Constraint(expr=3 * model.x[1] + 4 * model.x[2] >= 1)

        return model

    # *************************************************************************
    # *************************************************************************

    def problem_milp_infeasible(self):
        model = pyo.ConcreteModel("milp_infeasible")

        model.x = pyo.Var([1, 2], domain=pyo.Binary)

        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)

        return model

    # *************************************************************************
    # *************************************************************************

    def problem_milp_unbounded(self):
        model = pyo.ConcreteModel("milp_unbounded")

        model.x = pyo.Var([1, 2], domain=pyo.NonNegativeReals)

        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
        )

        model.Constraint1 = pyo.Constraint(expr=3 * model.x[1] + 4 * model.x[2] >= 1)

        return model

    # *************************************************************************
    # *************************************************************************

    def problem_milp_feasible(self, number_binary_variables=25, seed_number=None):
        if seed_number != None:
            random.seed(seed_number)

        model = pyo.ConcreteModel("milp_feasible")

        # a knapsack-type problem

        model.Y = pyo.RangeSet(number_binary_variables)

        model.y = pyo.Var(model.Y, domain=pyo.Binary)

        model.OBJ = pyo.Objective(
            expr=sum(model.y[j] * random.random() for j in model.Y), sense=pyo.maximize
        )

        model.Constraint1 = pyo.Constraint(
            expr=sum(model.y[j] * random.random() for j in model.Y)
            <= round(number_binary_variables / 5)
        )

        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

    # *************************************************************************
    # *************************************************************************

    def test_inexistent_solver(self):
        fake_solver = "fake_solver"
        good_solver = "glpk"
        # solver_options: dict = None

        # try using a fake solver and a problem incompatible with another solver

        # 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(),
        ]

        # problem types

        list_problem_types = [
            SolverInterface.PROBLEM_LP,
            SolverInterface.PROBLEM_LP,
            SolverInterface.PROBLEM_QP,
            "fake_problem_type",
        ]

        # list of solvers: one fake, one real

        list_solvers = [fake_solver, good_solver]

        # 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,
        }

        # *********************************************************************
        # *********************************************************************

        for solver_name in list_solvers:
            for index, problem in enumerate(list_problems):
                # optimise

                try:
                    # test problem-solver compatibility

                    problem_type = list_problem_types[index]

                    if (
                        SolverInterface.problem_and_solver_are_compatible(
                            solver_name, problem_type
                        )
                        == False
                    ):
                        continue

                except SolverInterface.UnknownSolverError:
                    assert True

                except SolverInterface.UnknownProblemTypeError:
                    assert True

                # test the solver interface

                try:
                    # configure common solver interface

                    _ = SolverInterface(solver_name=solver_name, **solver_options)

                except SolverInterface.UnknownSolverError:
                    assert True

    # **************************************************************************
    # **************************************************************************


# ******************************************************************************
# ******************************************************************************