[pull] master from industrial-optimization-group:master

desdeo/problem/cvxpy_evaluator.py

@@ -34,6 +34,7 @@ class CVXPYEvaluator:
     variables: dict[str, cp.Variable]
     parameters: dict[str, cp.Parameter]
     constraints: dict[str, cp.Constraint]
+    constraint_expressions: dict[str, cp.Expression]
     objective_functions: dict[str, cp.Expression]
     scalarizations: dict[str, cp.Expression]
     extra_functions: dict[str, cp.Expression]
@@ -55,6 +56,7 @@ def __init__(self, problem: Problem):
         self.variables = {}
         self.parameters = {}
         self.constraints = {}
+        self.constraint_expressions = {}
         self.objective_expr = None
         self.problem_model = None
 
@@ -98,65 +100,9 @@ def init_variables(self, problem: Problem) -> dict[str, cp.Variable]:
         Returns:
             dict[str, cp.Variable]: the variables for the problem.
         """
-        variables = {}
         for var in problem.variables:
-            if isinstance(var, Variable):
-                # handle regular variables
-                lowerbound = var.lowerbound
-                upperbound = var.upperbound
-
-                # Set bounds
-                bounds = None
-                if lowerbound is not None or upperbound is not None:
-                    lb = lowerbound if lowerbound is not None else -np.inf
-                    ub = upperbound if upperbound is not None else np.inf
-                    bounds = [lb, ub]
-
-                # figure out the variable type
-                match var.variable_type:
-                    case VariableTypeEnum.integer:
-                        # variable is integer
-                        cv_var = cp.Variable(integer=True, bounds=bounds, name=var.symbol)
-                    case VariableTypeEnum.real:
-                        # variable is real
-                        cv_var = cp.Variable(bounds=bounds, name=var.symbol)
-                    case VariableTypeEnum.binary:
-                        cv_var = cp.Variable(boolean=True, bounds=bounds, name=var.symbol)
-                    case _:
-                        msg = f"Could not figure out the type for variable {var}."
-                        raise CVXPYEvaluatorError(msg)
-
-                variables[var.symbol] = cv_var
-
-            elif isinstance(var, TensorVariable):
-                # handle tensor variables, i.e., vectors etc..
-                lowerbounds = var.get_lowerbound_values() if var.lowerbounds is not None else None
-                upperbounds = var.get_upperbound_values() if var.upperbounds is not None else None
-
-                # Set bounds
-                bounds = None
-                if lowerbounds is not None or upperbounds is not None:
-                    lb = np.array(lowerbounds) if lowerbounds is not None else -np.inf
-                    ub = np.array(upperbounds) if upperbounds is not None else np.inf
-                    bounds = [lb, ub]
-
-                # figure out the variable type
-                match var.variable_type:
-                    case VariableTypeEnum.integer:
-                        # variable is integer
-                        cv_var = cp.Variable(shape=tuple(var.shape), integer=True, bounds=bounds, name=var.symbol)
-                    case VariableTypeEnum.real:
-                        # variable is real
-                        cv_var = cp.Variable(shape=tuple(var.shape), bounds=bounds, name=var.symbol)
-                    case VariableTypeEnum.binary:
-                        cv_var = cp.Variable(shape=tuple(var.shape), boolean=True, bounds=bounds, name=var.symbol)
-                    case _:
-                        msg = f"Could not figure out the type for variable {var}."
-                        raise CVXPYEvaluatorError(msg)
-
-                variables[var.symbol] = cv_var
-
-        return variables
+            self.add_variable(var)
+        return self.variables
 
     def init_parameters(
         self, problem: Problem
@@ -221,22 +167,9 @@ def init_objectives(self, problem: Problem) -> dict[str, cp.Expression]:
         Returns:
             dict[str, cp.Expression]: dict containing the objective functions.
         """
-        objective_functions: dict[str, cp.Expression] = {}
         for obj in problem.objectives:
-            expr = self.parse(obj.func, callback=self.get_expression_by_name)
-            if isinstance(expr, (int, float)):
-                warnings.warn(
-                    "One or more of the problem objectives seems to be a constant.",
-                    CVXPYEvaluatorWarning,
-                    stacklevel=2,
-                )
-
-            objective_functions[obj.symbol] = expr
-
-            # the obj.symbol_min objectives are used when optimizing and building scalarizations etc...
-            objective_functions[f"{obj.symbol}_min"] = -expr if obj.maximize else expr
-
-        return objective_functions
+            self.add_objective(obj)
+        return self.objective_functions
 
     def init_constraints(self, problem: Problem) -> dict[str, cp.Constraint]:
         """Add constraint expressions to a CVXPY problem.
@@ -250,23 +183,9 @@ def init_constraints(self, problem: Problem) -> dict[str, cp.Constraint]:
         Returns:
             dict[str, cp.Constraint]: dict of constraints keyed by symbol.
         """
-        constraints = {}
         for cons in problem.constraints:
-            expr = self.parse(cons.func, callback=self.get_expression_by_name)
-
-            match con_type := cons.cons_type:
-                case ConstraintTypeEnum.LTE:
-                    # constraints in DESDEO are defined such that they must be less than zero
-                    constraint = expr <= 0
-                case ConstraintTypeEnum.EQ:
-                    constraint = expr == 0
-                case _:
-                    msg = f"Constraint type of {con_type} not supported. Must be one of {ConstraintTypeEnum}."
-                    raise CVXPYEvaluatorError(msg)
-
-            constraints[cons.symbol] = constraint
-
-        return constraints
+            self.add_constraint(cons)
+        return self.constraints
 
     def init_scalarizations(self, problem: Problem) -> dict[str, cp.Expression]:
         """Add scalarization expressions to a CVXPY evaluator.
@@ -300,15 +219,15 @@ def add_constraint(self, constraint: Constraint) -> cp.Constraint:
             cp.Constraint: The CVXPY constraint that was added.
         """
         expr = self.parse(constraint.func, self.get_expression_by_name)
+        self.constraint_expressions[constraint.symbol] = expr
 
-        match con_type := constraint.cons_type:
+        match constraint.cons_type:
             case ConstraintTypeEnum.LTE:
-                # constraints in DESDEO are defined such that they must be less than zero
                 cvxpy_constraint = expr <= 0
             case ConstraintTypeEnum.EQ:
                 cvxpy_constraint = expr == 0
             case _:
-                msg = f"Constraint type of {con_type} not supported. Must be one of {ConstraintTypeEnum}."
+                msg = f"Constraint type of {constraint.cons_type} not supported. Must be one of {ConstraintTypeEnum}."
                 raise CVXPYEvaluatorError(msg)
 
         self.constraints[constraint.symbol] = cvxpy_constraint
@@ -369,7 +288,7 @@ def add_variable(self, var: Variable | TensorVariable) -> cp.Variable:
             if lowerbound is not None or upperbound is not None:
                 lb = lowerbound if lowerbound is not None else -cp.inf
                 ub = upperbound if upperbound is not None else cp.inf
-                bounds = (lb, ub)
+                bounds = [lb, ub]
 
             # figure out the variable type
             match var.variable_type:
@@ -395,7 +314,7 @@ def add_variable(self, var: Variable | TensorVariable) -> cp.Variable:
             if lowerbounds is not None or upperbounds is not None:
                 lb = np.array(lowerbounds) if lowerbounds is not None else -np.inf
                 ub = np.array(upperbounds) if upperbounds is not None else np.inf
-                bounds = (lb, ub)
+                bounds = [lb, ub]
 
             # figure out the variable type
             match var.variable_type:
@@ -457,7 +376,7 @@ def get_values(self) -> dict[str, float | int | bool | list[float] | list[int] |
         Raises:
             CVXPYEvaluatorError: if the problem has not been solved yet.
         """
-        if self.problem_model is None or self.problem_model.status not in ["optimal", "optimal_inaccurate"]:
+        if self.problem_model is None or self.problem_model.status not in [cp.OPTIMAL, cp.OPTIMAL_INACCURATE]:
             msg = "Problem has not been solved yet or did not achieve optimal status."
             raise CVXPYEvaluatorError(msg)
 
@@ -487,8 +406,7 @@ def get_values(self) -> dict[str, float | int | bool | list[float] | list[int] |
 
         if self.problem.constraints is not None:
             for con in self.problem.constraints:
-                constraint_expr = self.parse(con.func, callback=self.get_expression_by_name)
-                result_dict[con.symbol] = constraint_expr.value
+                result_dict[con.symbol] = self.constraint_expressions[con.symbol].value
 
         return result_dict
 
@@ -497,8 +415,6 @@ def set_optimization_target(self, target: str):
 
         Args:
             target (str): a str representing a symbol. Needs to match an objective function or scalarization
-                function already found in the evaluator.
-            maximize (bool): If true, the target function is maximized instead of minimized.
 
         Raises:
             CVXPYEvaluatorError: the given target was not found in the evaluator.
@@ -509,7 +425,7 @@ def set_optimization_target(self, target: str):
         elif target in self.scalarizations:
             maximize = False
         else:
-            msg = f"The gurobipy model has no objective or scalarization named {target}."
+            msg = f"The CVXPY model has no objective or scalarization named {target}."
             raise CVXPYEvaluatorError(msg)
 
         obj_expr = self.get_expression_by_name(target)
@@ -532,7 +448,8 @@ def solve(self, **kwargs):
         if self.problem_model.is_dcp():
             self.problem_model.solve(**kwargs)
         elif self.problem_model.is_dgp():
-            self.problem_model.solve(gp=True, **kwargs)
+            kwargs["gp"] = True
+            self.problem_model.solve(**kwargs)
         else:
             warnings.warn(
                 "The problem does not appear to be DCP or DGP. CVXPY may not be able to solve it.",

desdeo/problem/testproblems/simple_problem.py

@@ -250,6 +250,8 @@ def simple_constrained_quadratic_tensor_test_problem(dqp=False) -> Problem:
         constants=[mmult, bvector],
         constraints=[cons],
         objectives=[obj],
+        is_twice_differentiable=True,
+        is_convex=True,
     )
 
 

desdeo/tools/cvxpy_solver_interfaces.py

@@ -1,5 +1,6 @@
 """Defines solver interfaces for cvxpy."""
 
+import cvxpy as cp
 from pydantic import BaseModel, Field
 
 from desdeo.problem import (
@@ -79,18 +80,18 @@ def parse_cvxpy_optimizer_results(problem: Problem, evaluator: CVXPYEvaluator) -
     )
     lagrange_multipliers = None
 
-    success = evaluator.problem_model.status in {"optimal", "optimal_inaccurate"}
-    if evaluator.problem_model.status == "optimal":
+    success = evaluator.problem_model.status in {cp.OPTIMAL, cp.OPTIMAL_INACCURATE}
+    if evaluator.problem_model.status == cp.OPTIMAL:
         status = "Optimal solution found."
-    elif evaluator.problem_model.status == "optimal_inaccurate":
+    elif evaluator.problem_model.status == cp.OPTIMAL_INACCURATE:
         status = "Optimal solution found (inaccurate)."
-    elif evaluator.problem_model.status == "infeasible":
+    elif evaluator.problem_model.status == cp.INFEASIBLE:
         status = "Problem is infeasible."
-    elif evaluator.problem_model.status == "unbounded":
+    elif evaluator.problem_model.status == cp.UNBOUNDED:
         status = "Problem is unbounded."
-    elif evaluator.problem_model.status == "infeasible_inaccurate":
+    elif evaluator.problem_model.status == cp.INFEASIBLE_INACCURATE:
         status = "Problem is infeasible (inaccurate)."
-    elif evaluator.problem_model.status == "unbounded_inaccurate":
+    elif evaluator.problem_model.status == cp.UNBOUNDED_INACCURATE:
         status = "Problem is unbounded (inaccurate)."
     else:
         status = f"Optimization ended with status: {evaluator.problem_model.status}"
@@ -144,3 +145,18 @@ def solve(self, target: str) -> SolverResults:
         self.evaluator.set_optimization_target(target)
         self.evaluator.solve(**self.solve_options)
         return parse_cvxpy_optimizer_results(self.problem, self.evaluator)
+
+
+def check_cvxpy_suitability(problem: Problem) -> bool:
+    """Checks whether a problem is suitable for being solved with CVXPY."""
+    try:
+        evaluator = CVXPYEvaluator(problem)
+        for obj in problem.objectives:
+            evaluator.set_optimization_target(obj.symbol)
+            if not (evaluator.problem_model.is_dcp() or evaluator.problem_model.is_dgp()):
+                return False
+                break
+        else:
+            return True
+    except Exception:
+        return False

desdeo/tools/gurobipy_solver_interfaces.py

@@ -1,5 +1,8 @@
 """Defines solver interfaces for gurobipy."""
 
+import io
+import sys
+
 import gurobipy as gp
 
 from desdeo.problem import (
@@ -276,3 +279,28 @@ def solve(self, target: str) -> SolverResults:
         self.evaluator.set_optimization_target(target)
         self.evaluator.model.optimize()
         return parse_gurobipy_optimizer_results(self.problem, self.evaluator)
+
+
+def check_gurobi_license():
+    """Check if Gurobi is using a full license (not trial).
+
+    Returns:
+        True if using full academic/commercial license
+        False if using trial license or no license found
+    """
+    captured_output = io.StringIO()
+    original_stdout = sys.stdout
+
+    try:
+        sys.stdout = captured_output
+        with gp.Env(empty=True) as env:
+            env.setParam("OutputFlag", 1)
+            env.start()
+        sys.stdout = original_stdout
+
+        output = captured_output.getvalue()
+        return "Restricted license - for non-production use only" not in output  # noqa: TRY300
+
+    except Exception:
+        sys.stdout = original_stdout
+        return False

desdeo/tools/utils.py

@@ -14,8 +14,9 @@
     numpy_array_to_objective_dict,
     variable_dimension_enumerate,
 )
+from desdeo.tools.cvxpy_solver_interfaces import CVXPYSolver, CVXPYSolverOptions, check_cvxpy_suitability
 from desdeo.tools.generics import BaseSolver
-from desdeo.tools.gurobipy_solver_interfaces import GurobipySolver, PersistentGurobipySolver
+from desdeo.tools.gurobipy_solver_interfaces import GurobipySolver, PersistentGurobipySolver, check_gurobi_license
 from desdeo.tools.ng_solver_interfaces import NevergradGenericOptions, NevergradGenericSolver
 from desdeo.tools.proximal_solver import ProximalSolver
 from desdeo.tools.pyomo_solver_interfaces import (
@@ -69,6 +70,10 @@
         "constructor": PersistentGurobipySolver,
         "options": None,
     },
+    "cvxpy": {
+        "constructor": CVXPYSolver,
+        "options": CVXPYSolverOptions,
+    },
 }
 
 
@@ -122,8 +127,12 @@ def find_compatible_solvers(problem: Problem) -> list[BaseSolver]:
     ):
         solvers.append(available_solvers["pyomo_ipopt"]["constructor"])  # ipopt has to be installed
 
+    # check if the problem is convex or log-log convex
+    if check_cvxpy_suitability(problem):
+        solvers.append(available_solvers["cvxpy"]["constructor"])
+
     # check if the problem is linear
-    if problem.is_linear:
+    if problem.is_linear and check_gurobi_license():
         solvers.append(available_solvers["gurobipy"]["constructor"])
     if problem.is_linear and shutil.which("gurobi"):
         solvers.append(available_solvers["pyomo_gurobi"]["constructor"])  # gurobi has to be installed
@@ -197,9 +206,13 @@ def guess_best_solver(problem: Problem) -> BaseSolver:
         return available_solvers["proximal"]["constructor"]
 
     # check if the problem is linear
-    if problem.is_linear:
+    if problem.is_linear and check_gurobi_license():
         return available_solvers["gurobipy"]["constructor"]
 
+    # check if the problem is convex or log-log convex
+    if check_cvxpy_suitability(problem):
+        return available_solvers["cvxpy"]["constructor"]
+
     # check if the problem is differentiable and if it is mixed integer
     if (
         problem.is_twice_differentiable