Source code for pyomo.contrib.gdpopt.gloa

#  ___________________________________________________________________________
#
#  Pyomo: Python Optimization Modeling Objects
#  Copyright (c) 2008-2024
#  National Technology and Engineering Solutions of Sandia, LLC
#  Under the terms of Contract DE-NA0003525 with National Technology and
#  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
#  rights in this software.
#  This software is distributed under the 3-clause BSD License.
#  ___________________________________________________________________________

from pyomo.common.config import document_kwargs_from_configdict
from pyomo.common.errors import DeveloperError
from pyomo.common.modeling import unique_component_name
from pyomo.contrib.gdp_bounds.info import disjunctive_bounds
from pyomo.contrib.gdpopt.algorithm_base_class import _GDPoptAlgorithm
from pyomo.contrib.gdpopt.config_options import (
    _add_oa_configs,
    _add_mip_solver_configs,
    _add_nlp_solver_configs,
    _add_tolerance_configs,
)
from pyomo.contrib.gdpopt.create_oa_subproblems import (
    _get_discrete_problem_and_subproblem,
    add_constraints_by_disjunct,
    add_global_constraint_list,
)
from pyomo.contrib.gdpopt.cut_generation import add_no_good_cut
from pyomo.contrib.gdpopt.oa_algorithm_utils import _OAAlgorithmMixIn
from pyomo.contrib.gdpopt.solve_discrete_problem import solve_MILP_discrete_problem
from pyomo.contrib.gdpopt.util import (
    _add_bigm_constraint_to_transformed_model,
    time_code,
)
from pyomo.contrib.mcpp.pyomo_mcpp import McCormick as mc, MCPP_Error

from pyomo.core import Constraint, Block, NonNegativeIntegers, Objective, value
from pyomo.core.expr.numvalue import is_potentially_variable
from pyomo.core.expr.visitor import identify_variables
from pyomo.opt.base import SolverFactory


[docs]@SolverFactory.register( 'gdpopt.gloa', doc="The GLOA (global logic-based outer approximation) Generalized " "Disjunctive Programming (GDP) solver", ) class GDP_GLOA_Solver(_GDPoptAlgorithm, _OAAlgorithmMixIn): """The GDPopt (Generalized Disjunctive Programming optimizer) global logic-based outer approximation (GLOA) solver. Accepts models that can include nonlinear, continuous variables and constraints, as well as logical conditions. """ CONFIG = _GDPoptAlgorithm.CONFIG() _add_oa_configs(CONFIG) _add_mip_solver_configs(CONFIG) _add_nlp_solver_configs(CONFIG, default_solver='couenne') _add_tolerance_configs(CONFIG) algorithm = 'GLOA' # Override solve() to customize the docstring for this solver
[docs] @document_kwargs_from_configdict(CONFIG, doc=_GDPoptAlgorithm.solve.__doc__) def solve(self, model, **kwds): return super().solve(model, **kwds)
def _log_citation(self, config): config.logger.info( "\n" + """- GLOA algorithm: Lee, S; Grossmann, IE. A Global Optimization Algorithm for Nonconvex Generalized Disjunctive Programming and Applications to Process Systems. Comp. and Chem. Eng. 2001, 25, 1675-1697. DOI: 10.1016/S0098-1354(01)00732-3. """.strip() ) def _solve_gdp(self, original_model, config): logger = config.logger # we need to gather a map of Disjuncts to their active Constraints # before we call any GDP transformations, as we will need this # information for cut generation later add_constraints_by_disjunct(self.original_util_block) # We also save these in advance because we know only linear logical # constraints will be added by the transformation to a MIP, so these are # all we'll ever need. add_global_constraint_list(self.original_util_block) (discrete_problem_util_block, subproblem_util_block) = ( _get_discrete_problem_and_subproblem(self, config) ) discrete = discrete_problem_util_block.parent_block() subproblem = subproblem_util_block.parent_block() discrete_obj = next( discrete.component_data_objects(Objective, active=True, descend_into=True) ) self._log_header(logger) # main loop while not config.iterlim or self.iteration < config.iterlim: self.iteration += 1 # solve linear discrete problem with time_code(self.timing, 'mip'): mip_feasible = solve_MILP_discrete_problem( discrete_problem_util_block, self, config ) self._update_bounds_after_discrete_problem_solve( mip_feasible, discrete_obj, logger ) # Check termination conditions if self.any_termination_criterion_met(config): break with time_code(self.timing, 'nlp'): self._fix_discrete_soln_solve_subproblem_and_add_cuts( discrete_problem_util_block, subproblem_util_block, config ) # Add integer cut with time_code(self.timing, "integer cut generation"): add_no_good_cut(discrete_problem_util_block, config) # Check termination conditions if self.any_termination_criterion_met(config): break def _add_cuts_to_discrete_problem( self, subproblem_util_block, discrete_problem_util_block, objective_sense, config, timing, ): """Add affine cuts""" m = discrete_problem_util_block.parent_block() if hasattr(discrete_problem_util_block, "aff_utils_blocks"): aff_utils_blocks = discrete_problem_util_block.aff_utils_blocks else: aff_utils_blocks = discrete_problem_util_block.aff_utils_blocks = dict() config.logger.debug("Adding affine cuts.") counter = 0 for discrete_var, subprob_var in zip( discrete_problem_util_block.algebraic_variable_list, subproblem_util_block.algebraic_variable_list, ): val = subprob_var.value if val is not None and not discrete_var.fixed: discrete_var.set_value(val, skip_validation=True) for constr in self._get_active_untransformed_constraints( discrete_problem_util_block, config ): disjunctive_var_bounds = disjunctive_bounds(constr.parent_block()) if constr.body.polynomial_degree() in (1, 0): continue vars_in_constr = list(identify_variables(constr.body)) if any(var.value is None for var in vars_in_constr): continue # a variable has no values # mcpp stuff try: mc_eqn = mc(constr.body, disjunctive_var_bounds) except MCPP_Error as e: config.logger.debug( "Skipping constraint %s due to MCPP " "error %s" % (constr.name, str(e)) ) continue # skip to the next constraint ccSlope = mc_eqn.subcc() cvSlope = mc_eqn.subcv() ccStart = mc_eqn.concave() cvStart = mc_eqn.convex() ub_int = ( min(value(constr.upper), mc_eqn.upper()) if constr.has_ub() else mc_eqn.upper() ) lb_int = ( max(value(constr.lower), mc_eqn.lower()) if constr.has_lb() else mc_eqn.lower() ) parent_block = constr.parent_block() # Create a block on which to put outer approximation cuts. aff_utils = aff_utils_blocks.get(parent_block) if aff_utils is None: aff_utils = Block(doc="Block holding affine constraints") nm = unique_component_name(parent_block, "GDPopt_aff") parent_block.add_component(nm, aff_utils) aff_utils_blocks[parent_block] = aff_utils aff_utils.GDPopt_aff_cons = Constraint(NonNegativeIntegers) aff_cuts = aff_utils.GDPopt_aff_cons cut_body = sum( ccSlope[var] * (var - var.value) for var in vars_in_constr if not var.fixed ) if not is_potentially_variable(cut_body): if ( cut_body + ccStart >= lb_int - config.constraint_tolerance and cut_body + cvStart <= ub_int + config.constraint_tolerance ): # We won't add them, but nothing is wrong--they hold config.logger.debug("Affine cut is trivially True.") else: # something went wrong. raise DeveloperError("One of the affine cuts is trivially False.") else: concave_cut = cut_body + ccStart >= lb_int convex_cut = cut_body + cvStart <= ub_int idx = len(aff_cuts) aff_cuts[idx] = concave_cut aff_cuts[idx + 1] = convex_cut _add_bigm_constraint_to_transformed_model(m, aff_cuts[idx], aff_cuts) _add_bigm_constraint_to_transformed_model( m, aff_cuts[idx + 1], aff_cuts ) counter += 2 config.logger.debug("Added %s affine cuts" % counter)