Custom python solver for pyomo

Question

I am responsible for a series of exercises for nonlinear optimization. I thought it would be cool to start with some examples of optimization problems and solve them with pyomo + some black box solvers.

However, as the students learn more about optimization algorithms I wanted them to also implement some simple methods and test there implementation for the same examples. I hoped there would be an "easy" way to add a custom solver to pyomo however I cannot find any information about this.

Basically that would allow the students to check their implementation by just changing a single line in there code and compare to a well tested solver.

I would also try to implement a simple wrapper myself but I do not know anything about the pyomo internals.

Q: Can I add my own solvers written in python to pyomo? Solver could have an interface like the ones of scipy.optimize.

Ty for reading,

Franz

---

Related:

• Pyomo-Solver Communication
• Call scipy.optimize inside pyomo

Answer 1

Yes, we have a number of ways to do this in Pyomo and lots of examples in various modules in pyomo.contrib. I think the easiest way for students would be to leverage a Pyomo extension module called PyNumero which was designed for writing nonlinear optimization algorithms in Python and can be found here: https://github.com/Pyomo/pyomo/tree/main/pyomo/contrib/pynumero

Here are a couple implementation examples using PyNumero:

SQP --> https://github.com/Pyomo/pyomo/blob/main/pyomo/contrib/pynumero/examples/sqp.py

Interior Point --> https://github.com/Pyomo/pyomo/tree/main/pyomo/contrib/interior_point

Interfacing to Scipy square solvers --> https://github.com/Pyomo/pyomo/blob/main/pyomo/contrib/pynumero/algorithms/solvers/scipy_solvers.py

Answer 2

The short answer is, "Yes: Pyomo is designed to be extensible and user-defined solvers are definitely supported."

The details are a bit more involved. While Pyomo defines several "solver" class hierarchies, they are by no means required. The core Pyomo relies on duck-typing and not inheritance when registering and working with solvers. At the moment, the simplest solver interface needs only to implement the following methods:

class GenericSolverInterface(object):
    def solve(self,
              model: _BlockData,
              tee: bool = False,
              load_solutions: bool = True,
              logfile: Optional[str] = None,
              solnfile: Optional[str] = None,
              timelimit: Optional[float] = None,
              report_timing: bool = False,
              solver_io: Optional[str] = None,
              suffixes: Optional[Sequence] = None,
              options: Optional[Dict] = None,
              keepfiles: bool = False,
              symbolic_solver_labels: bool = False):
        pass

    def available(self, exception_flag=True) -> bool:
        pass

    def license_is_valid(self) -> bool:
        pass

    def version(self) -> Tuple:
        pass

    @property
    def options(self) -> ConfigDict:
        pass

    @options.setter
    def options(self, val):
        pass

Solvers are "registered" with the SolverFactory by either decorating the class:

from pyomo.opt import SolverFactory

@SolverFactory.register('demo', doc='DEMO Solver interface')
class DemoSolver(GenericSolverInterface):
    #...

or by explicitly calling the registration function

SolverFactory.register('demo', doc='DEMO Solver Interface')(DemoSolver)

The real trick is now how to implement the solve() method. Pyomo will hand solve() the model (or Block) that the user wants to solve. solve() then needs to convert the Pyomo model/Block into the format required by the solver. This step is (obviously) very solver-specific, although there are somewhat standardized approaches for some common solver interface standards (e.g., LP files, MPS files, NL files, BAR files, etc.). There are also other "in-development" interfaces that are more performant for specific use cases (direct in-memory solver interfaces, APPSI, PyNumero).

The simplest (and not very performant) approach could be to define a helper class that maps the Pyomo data structures to a simple set of lists:

from pyomo.core.base import Objective, Constraint
from pyomo.core.expr.visitor import identify_variables

class model_interface(object):
    def __init__(self, model):
        self.model = model
        self._obj = list(model.component_data_objects(Objective, active=True))
        self._con = list(model.component_data_objects(Constraint, active=True))
        self._var = {}
        for c in self._con:
            self._var.update({id(v): v for v in identify_variables(c.body)})
        self._var = list(self._var.values())

    @property
    def x(self):
        """Return the current list of variable values"""
        return [v.value for v in self._var]

    @x.setter
    def x(self, values):
        """Set the variables to new values"""
        for v, val in zip(self._var, values):
            v.set_value(val)

    @property
    def x_lb(self):
        """Return the list of variable lower bounds (may include None)"""
        return [v.lb for v in self._var]

    @property
    def x_ub(self):
        """Return the list of variable upper bounds (may include None)"""
        return [v.ub for v in self._var]

    @property
    def obj(self):
        """Return the list of objective values (computed using the current
        variable values)"""
        return [value(o) for o in self._obj]

    @property
    def con(self):
        """Return the list of constraint 'body' values (computed using the
        current variable values)
        """
        return [value(c) for c in self._con]

    @property
    def con_lb(self):
        """Return the list of constraint lower bounds (may include None)"""
        return [c.lb for c in self._con]

    @property
    def con_ub(self):
        """Return the list of constraint upper bounds (may include None)"""
        return [c.ub for c in self._con]