Source code for gemseo.algos.opt.lib_nlopt

# Copyright 2021 IRT Saint Exupéry,
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
# License version 3 as published by the Free Software Foundation.
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# Lesser General Public License for more details.
# You should have received a copy of the GNU Lesser General Public License
# along with this program; if not, write to the Free Software Foundation,
# Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
# Contributors:
#    INITIAL AUTHORS - initial API and implementation and/or initial
#                           documentation
#        :author: Damien Guenot
#         Francois Gallard : refactoring for v1, May 2016
"""NLopt library wrapper.

    If the objective, or a constraint, of the :class:`.OptimizationProblem`
    returns a value of type ``int``
    then ``nlopt.opt.optimize`` will terminate with
    ``ValueError: nlopt invalid argument``.

    This behavior has been identified as
    `a bug internal to NLopt 2.7.1 <>`_
    and has been fixed in the development version of NLopt.

    Until a new version of NLopt including the bugfix is released,
    the user of |g| shall provide objective and constraint functions
    that return values of type ``float`` and ``NDArray[float]``.

from __future__ import annotations

import logging
from dataclasses import dataclass
from typing import TYPE_CHECKING
from typing import Any
from typing import Callable
from typing import ClassVar
from typing import Union

import nlopt
from nlopt import RoundoffLimited
from numpy import atleast_1d
from numpy import atleast_2d
from numpy import ndarray

from gemseo.algos.opt.optimization_library import OptimizationAlgorithmDescription
from gemseo.algos.opt.optimization_library import OptimizationLibrary
from gemseo.algos.stop_criteria import TerminationCriterion
from gemseo.core.mdofunctions.mdo_function import MDOFunction

    from gemseo.algos.opt_problem import OptimizationProblem
    from gemseo.algos.opt_result import OptimizationResult

LOGGER = logging.getLogger(__name__)

NLoptOptionsType = Union[bool, int, float]

[docs] @dataclass class NLoptAlgorithmDescription(OptimizationAlgorithmDescription): """The description of an optimization algorithm from the NLopt library.""" library_name: str = "NLopt"
[docs] class Nlopt(OptimizationLibrary): """NLopt optimization library interface. See OptimizationLibrary. """ LIB_COMPUTE_GRAD = False INNER_MAXEVAL = "inner_maxeval" STOPVAL = "stopval" CTOL_ABS = "ctol_abs" INIT_STEP = "init_step" SUCCESS = "NLOPT_SUCCESS: Generic success return value" STOPVAL_REACHED = ( "NLOPT_STOPVAL_REACHED: Optimization stopped " "because stopval (above) was reached" ) FTOL_REACHED = ( "NLOPT_FTOL_REACHED: Optimization stopped " "because ftol_rel or ftol_abs (above) was reached" ) XTOL_REACHED = ( "NLOPT_XTOL_REACHED Optimization stopped " "because xtol_rel or xtol_abs (above) was reached" ) MAXEVAL_REACHED = ( "NLOPT_MAXEVAL_REACHED: Optimization stopped " "because maxeval (above) was reached" ) MAXTIME_REACHED = ( "NLOPT_MAXTIME_REACHED: Optimization stopped " "because maxtime (above) was reached" ) FAILURE = "NLOPT_FAILURE: Generic failure code" INVALID_ARGS = ( "NLOPT_INVALID_ARGS: Invalid arguments (e.g. lower " "bounds are bigger than upper bounds, an unknown" " algorithm was specified, etcetera)." ) OUT_OF_MEMORY = "OUT_OF_MEMORY: Ran out of memory" ROUNDOFF_LIMITED = ( "NLOPT_ROUNDOFF_LIMITED: Halted because " "roundoff errors limited progress. (In this " "case, the optimization still typically " "returns a useful result.)" ) FORCED_STOP = ( "NLOPT_FORCED_STOP: Halted because of a forced " "termination: the user called nlopt_force_stop" "(opt) on the optimization's nlopt_opt" " object opt from the user's objective " "function or constraints." ) NLOPT_MESSAGES: ClassVar[dict[int, str]] = { 1: SUCCESS, 2: STOPVAL_REACHED, 3: FTOL_REACHED, 4: XTOL_REACHED, 5: MAXEVAL_REACHED, 6: MAXTIME_REACHED, -1: FAILURE, -2: INVALID_ARGS, -3: OUT_OF_MEMORY, -4: ROUNDOFF_LIMITED, -5: FORCED_STOP, } LIBRARY_NAME = "NLopt" def __init__(self) -> None: # noqa:D107 super().__init__() nlopt_doc = "" self.descriptions = { "NLOPT_MMA": NLoptAlgorithmDescription( algorithm_name="MMA", description=( "Method of Moving Asymptotes (MMA)" "implemented in the NLOPT library" ), handle_inequality_constraints=True, internal_algorithm_name=nlopt.LD_MMA, require_gradient=True, website=f"{nlopt_doc}#mma-method-of-moving-asymptotes-and-ccsa", ), "NLOPT_COBYLA": NLoptAlgorithmDescription( algorithm_name="COBYLA", description=( "Constrained Optimization BY Linear " "Approximations (COBYLA) implemented " "in the NLOPT library" ), handle_equality_constraints=True, handle_inequality_constraints=True, internal_algorithm_name=nlopt.LN_COBYLA, website=( f"{nlopt_doc}#cobyla-constrained-optimization-by-linear-" "approximations" ), ), "NLOPT_SLSQP": NLoptAlgorithmDescription( algorithm_name="SLSQP", description=( "Sequential Least-Squares Quadratic " "Programming (SLSQP) implemented in " "the NLOPT library" ), handle_equality_constraints=True, handle_inequality_constraints=True, internal_algorithm_name=nlopt.LD_SLSQP, require_gradient=True, website=f"{nlopt_doc}#slsqp", ), "NLOPT_BOBYQA": NLoptAlgorithmDescription( algorithm_name="BOBYQA", description=( "Bound Optimization BY Quadratic " "Approximation (BOBYQA) implemented " "in the NLOPT library" ), internal_algorithm_name=nlopt.LN_BOBYQA, website=f"{nlopt_doc}#bobyqa", ), "NLOPT_BFGS": NLoptAlgorithmDescription( algorithm_name="BFGS", description=( "Broyden-Fletcher-Goldfarb-Shanno method " "(BFGS) implemented in the NLOPT library" ), internal_algorithm_name=nlopt.LD_LBFGS, require_gradient=True, website=f"{nlopt_doc}#low-storage-bfgs", ), # Does not work on Rastrigin => banned # 'NLOPT_ESCH': { Does not work on Rastrigin # self.INTERNAL_NAME: nlopt.GN_ESCH, # self.REQUIRE_GRAD: False, # self.HANDLE_EQ_CONS: False, # self.HANDLE_INEQ_CONS: False}, "NLOPT_NEWUOA": NLoptAlgorithmDescription( algorithm_name="NEWUOA", description=( "NEWUOA + bound constraints implemented in the NLOPT library" ), internal_algorithm_name=nlopt.LN_NEWUOA_BOUND, website=f"{nlopt_doc}#newuoa-bound-constraints", ), # Does not work on Rastrigin => banned # 'NLOPT_ISRES': { # self.INTERNAL_NAME: nlopt.GN_ISRES, # self.REQUIRE_GRAD: False, # self.HANDLE_EQ_CONS: True, # self.HANDLE_INEQ_CONS: True} } def _get_options( self, ftol_abs: float = 1e-14, # pylint: disable=W0221 xtol_abs: float = 1e-14, max_time: float = 0.0, max_iter: int = 999, ftol_rel: float = 1e-8, xtol_rel: float = 1e-8, ctol_abs: float = 1e-6, stopval: float | None = None, normalize_design_space: bool = True, eq_tolerance: float = 1e-2, ineq_tolerance: float = 1e-4, init_step: float = 0.25, kkt_tol_abs: float | None = None, kkt_tol_rel: float | None = None, stop_crit_n_x: int | None = None, **kwargs: Any, ) -> dict[str, NLoptOptionsType]: r"""Retrieve the options of the Nlopt library. Args: ftol_abs: The absolute tolerance on the objective function. xtol_abs: The absolute tolerance on the design parameters. max_time: The maximum runtime in seconds. The value 0 means no runtime limit. max_iter: The maximum number of iterations. ftol_rel: The relative tolerance on the objective function. xtol_rel: The relative tolerance on the design parameters. ctol_abs: The absolute tolerance on the constraints. stopval: The objective value at which the optimization will stop. Stop minimizing when an objective value :math:`\leq` stopval is found, or stop maximizing when a value :math:`\geq` stopval is found. If ``None``, this termination condition will not be active. kkt_tol_abs: The absolute tolerance on the KKT residual norm. If ``None`` this criterion is not activated. kkt_tol_rel: The relative tolerance on the KKT residual norm. If ``None`` this criterion is not activated. normalize_design_space: If ``True``, normalize the design variables between 0 and 1. eq_tolerance: The tolerance on the equality constraints. ineq_tolerance: The tolerance on the inequality constraints. init_step: The initial step size :math:`r` for derivative-free algorithms. Increasing the initial step size will make the initial DOE of size :math:`d+1` take wider steps in the design variables. In details, given a :math:`d`-length design vector initialized to :math:`x_0`, the first value of the design vector will be the initial one :math:`x^{(1)}=x_0`, the second one will be :math:`x^{(2)}=x^{(1)}+(r(\max(x_1)-\min(x_1)),0,\ldots,0)`, ..., the :math:`d+1`-th one will be :math:`x^{(d+1)}=x^{d}+(0,\ldots,0,r(\max(x_d)-\min(x_d)))`. Note that in a normalized design space, :math:`\min(x_i)=0` and :math:`\max(x_i)=1`. stop_crit_n_x: The minimum number of design vectors to take into account in the stopping criteria. **kwargs: The additional algorithm-specific options. Returns: The NLopt library options with their values. """ nds = normalize_design_space return self._process_options( ftol_rel=ftol_rel, ftol_abs=ftol_abs, xtol_rel=xtol_rel, xtol_abs=xtol_abs, stop_crit_n_x=stop_crit_n_x, max_time=max_time, max_iter=max_iter, ctol_abs=ctol_abs, normalize_design_space=nds, stopval=stopval, eq_tolerance=eq_tolerance, ineq_tolerance=ineq_tolerance, init_step=init_step, kkt_tol_abs=kkt_tol_abs, kkt_tol_rel=kkt_tol_rel, **kwargs, ) def __opt_objective_grad_nlopt( self, xn_vect: ndarray, grad: ndarray, ) -> float: """Evaluate the objective and gradient functions for NLopt. Args: xn_vect: The normalized design variables vector. grad: The gradient of the objective function. Returns: The evaluation of the objective function for the given `xn_vect`. """ obj_func = self.problem.objective if grad.size > 0: grad[:] = obj_func.jac(xn_vect) return float(obj_func.func(xn_vect).real) def __make_constraint( self, func: Callable[[ndarray], ndarray], jac: Callable[[ndarray], ndarray], index_cstr: int, ) -> Callable[[ndarray, ndarray], ndarray]: """Build NLopt-like constraints. No vector functions are allowed. The database will avoid multiple evaluations. Args: func: The function pointer. jac: The Jacobian pointer. index_cstr: The index of the constraint. Returns: The constraint function. """ def cstr_fun_grad( xn_vect: ndarray, grad: ndarray, ) -> ndarray: """Define the function to be given as a pointer to the optimizer. Used to compute constraints and constraints gradients if required. Args: xn_vect: The normalized design vector. grad: The gradient of the objective function. Returns: The result of evaluating the function for a given constraint. """ if self.descriptions[self.algo_name].require_gradient and grad.size > 0: cstr_jac = jac(xn_vect) grad[:] = atleast_2d(cstr_jac)[index_cstr,] return atleast_1d(func(xn_vect).real)[index_cstr] return cstr_fun_grad def __add_constraints( self, nlopt_problem: nlopt.opt, ctol: float = 0.0, ) -> None: """Add all the constraints to the optimization problem. Args: nlopt_problem: The optimization problem. ctol: The absolute tolerance on the constraints. """ for constraint in self.problem.constraints: f_type = constraint.f_type func = constraint.func jac = constraint.jac dim = constraint.dim for idim in range(dim): nl_fun = self.__make_constraint(func, jac, idim) if f_type == MDOFunction.ConstraintType.INEQ: nlopt_problem.add_inequality_constraint(nl_fun, ctol) elif f_type == MDOFunction.ConstraintType.EQ: nlopt_problem.add_equality_constraint(nl_fun, ctol) def __set_prob_options( self, nlopt_problem: nlopt.opt, **opt_options: Any, ) -> nlopt.opt: """Set the options for the NLopt algorithm. Args: nlopt_problem: The optimization problem from NLopt. **opt_options: The NLopt optimization options. Returns: The updated NLopt problem. """ # ALready 0 by default # nlopt_problem.set_xtol_abs(0.0) # nlopt_problem.set_xtol_rel(0.0) # nlopt_problem.set_ftol_rel(0.0) # nlopt_problem.set_ftol_abs(0.0) nlopt_problem.set_maxeval(int(1.5 * opt_options[self.MAX_ITER])) # anti-cycling nlopt_problem.set_maxtime(opt_options[self.MAX_TIME]) nlopt_problem.set_initial_step(opt_options[self.INIT_STEP]) if self.STOPVAL in opt_options: stopval = opt_options[self.STOPVAL] if stopval is not None: nlopt_problem.set_stopval(stopval) if self.INNER_MAXEVAL in opt_options: nlopt_problem.set_param(self.INNER_MAXEVAL, opt_options[self.INNER_MAXEVAL]) return nlopt_problem def _pre_run( self, problem: OptimizationProblem, algo_name: str, **options: NLoptOptionsType, ) -> None: """Set :attr:`.STOP_CRIT_NX` depending on the algorithm. The COBYLA and BOBYQA algorithms create sets of interpolation points of sizes ``N+1`` and ``2*N+1`` respectively at initialization, where ``N`` is the dimension of the design space. In some cases, a termination criterion can be matched during this phase, leading to a premature termination. In order to circumvent this, :attr:`.STOP_CRIT_NX` is set accordingly, depending on the algorithm used. It ensures that the termination criterion will not be triggered during this preliminary Design of Experiment phase of the algorithm. Args: problem: The optimization problem. algo_name: The name of the algorithm. **options: The options of the algorithm, see the associated JSON file. """ n_stop_crit_x = options[self.STOP_CRIT_NX] if algo_name == "NLOPT_COBYLA" and not n_stop_crit_x: design_space_dimension = self.problem.design_space.dimension options[self.STOP_CRIT_NX] = design_space_dimension + 1 elif algo_name == "NLOPT_BOBYQA" and not n_stop_crit_x: design_space_dimension = self.problem.design_space.dimension options[self.STOP_CRIT_NX] = 2 * design_space_dimension + 1 else: options[self.STOP_CRIT_NX] = n_stop_crit_x or 3 super()._pre_run(problem, algo_name, **options) def _run(self, **options: NLoptOptionsType) -> OptimizationResult: """Run the algorithm. Args: **options: The options for the algorithm, see associated JSON file. Returns: The optimization result. Raises: TerminationCriterion: If the driver stops for some reason. """ normalize_ds = options.pop(self.NORMALIZE_DESIGN_SPACE_OPTION, True) # Get the bounds anx x0 x_0, l_b, u_b = self.get_x0_and_bounds_vects(normalize_ds) nlopt_problem = nlopt.opt(self.internal_algo_name, x_0.shape[0]) # Set the normalized bounds: nlopt_problem.set_lower_bounds(l_b.real) nlopt_problem.set_upper_bounds(u_b.real) nlopt_problem.set_min_objective(self.__opt_objective_grad_nlopt) if self.CTOL_ABS in options: ctol = options[self.CTOL_ABS] self.__add_constraints(nlopt_problem, ctol) nlopt_problem = self.__set_prob_options(nlopt_problem, **options) try: nlopt_problem.optimize(x_0.real) except (RoundoffLimited, RuntimeError) as err: LOGGER.exception( "NLopt run failed: %s, %s", str(err.args[0]), str(err.__class__.__name__), ) raise TerminationCriterion from None message = self.NLOPT_MESSAGES[nlopt_problem.last_optimize_result()] status = nlopt_problem.last_optimize_result() return self.get_optimum_from_database(message, status)