diff --git a/doc/OnlineDocs/explanation/analysis/index.rst b/doc/OnlineDocs/explanation/analysis/index.rst index 0a8e3c3b416..5537636b1bb 100644 --- a/doc/OnlineDocs/explanation/analysis/index.rst +++ b/doc/OnlineDocs/explanation/analysis/index.rst @@ -12,6 +12,7 @@ Analysis in Pyomo mpc/index parmest/index sensitivity_toolbox + nlp_initialization/nlp_initialization .. Reorganization notes: diff --git a/doc/OnlineDocs/explanation/analysis/nlp_initialization/nlp_initialization.rst b/doc/OnlineDocs/explanation/analysis/nlp_initialization/nlp_initialization.rst new file mode 100644 index 00000000000..1ce1ffaa0fc --- /dev/null +++ b/doc/OnlineDocs/explanation/analysis/nlp_initialization/nlp_initialization.rst @@ -0,0 +1,131 @@ +.. _analysis_nlp_initialization: + +NLP Initialization +****************** + +.. warning:: + + This package lives in :mod:`pyomo.devel`. APIs, options, and behavior may + change without notice. + +The initialization module within ``pyomo.devel.initialization`` is intended to +provide methods to help initialize nonconvex nonlinear programs (NLPs). The +goal is to increase the chance of finding a local minimizer (i.e., decrease the +chance of getting stuck a point that locally minimizes infeasibility). If +you are already able to solve your problem with a local NLP solver, these +tools will not help you. Example usage is shown below. + +.. literalinclude:: /../../pyomo/devel/initialization/examples/init_polynomial_ex.py + :start-after: # === Required imports === + +The :func:`initialize_nlp ` +function uses the specified method to try to find a good starting point for the +NLP solver and then attempts to solve the problem with the given NLP solver. + +.. note:: + + Currently, this module only works with solvers from :mod:`pyomo.contrib.solver`. + + +Initialization Methods +====================== + +The initialization method is selected using the +:class:`InitializationMethod ` enum. + +.. note:: + + Not all of the methods described below require all nonlinear variables to be + bounded. However, all of the methods will perform better if all nonlinear + variables are bounded (the tighter the bounds, the better). + + +Method ``global_opt`` +--------------------- + +This method uses an MINLP solver to try to find a feasible solution. We +adjust the solver parameters so that the solver will stop as soon as any +feasible solution is found. We then initialize the NLP solver at that +feasible solution. Many MINLP solvers will default to a very large +time limit, so it can be useful to specify a time limit before +calling :func:`initialize_nlp `: + +.. testcode:: + :skipif: not scip_available + + import pyomo.environ as pyo + from pyomo.contrib.solver.common.factory import SolverFactory + + global_solver = SolverFactory('scip_direct') + global_solver.config.time_limit = 600 # 10 minutes + # now call initialize_nlp + +This method currently works with the following solver interfaces for MINLP solvers: + +* SCIP (:class:`direct ` and + :class:`persistent `) +* :class:`Gurobi MINLP ` + +Advantages +^^^^^^^^^^ + +* Currently, this is the method that is most likely to succeed in finding a + feasible solution. +* Does not strictly require variable bounds + +Disadvantages +^^^^^^^^^^^^^ + +* This method will only work if the model is completely algebraic. It will not + work with external functions. + + +Method ``pwl_approximation`` +---------------------------- + +This method builds a piecewise linear (PWL) approximation of the model, solves +it, and initializes the NLP solver at the solution. If the NLP solver does not +converge, then the PWL approximation will be refined by adding additional +"segments". This is repeated until either a feasible solution is found or +the iteration limit is reached. + +This method does not currently work as well as ``global_opt``, but it does +have a great deal of potential. We expect future versions of this method +to perform significantly better. + +Advantages +^^^^^^^^^^ + +* Does not require an MINLP solver +* Future versions will work with external functions + +Disadvantages +^^^^^^^^^^^^^ + +* Current implementation can be slow +* Requires all nonlinear variables to be bounded + + +Method ``lp_approximation`` +--------------------------- + +This method is similar to the PWL approximation method, but it builds +an LP approximation instead and does not do any refinement. Another +distinction is that the LP approximation uses a linear least-squares +fit, so the approximation may not equal the original function at the +variable bounds. This also means that variable bounds are not strictly +necessary, though they do help improve the approximation. + +Advantages +^^^^^^^^^^ + +* Fast +* Future versions will work with external functions +* Does not strictly require variable bounds +* Does not require an MINLP or even an MILP solver + +Disadvantages +^^^^^^^^^^^^^ + +* This method only attempts to initialize the problem once. If it does + not succeed, it is done. diff --git a/doc/OnlineDocs/explanation/index.rst b/doc/OnlineDocs/explanation/index.rst index 0121debcd32..27a3227de9d 100644 --- a/doc/OnlineDocs/explanation/index.rst +++ b/doc/OnlineDocs/explanation/index.rst @@ -43,6 +43,7 @@ Explanations `Design of Experiments` `MPC` `AOS` + `NLP Initialization` `Modeling Utilities` `Latex Printer` `FME` diff --git a/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py b/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py index 5aa04e3f0b8..6912811eb00 100644 --- a/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py +++ b/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py @@ -118,13 +118,13 @@ def _get_random_point_grid(bounds, n, func, config, seed=42): return list(itertools.product(*linspaces)) -def _get_uniform_point_grid(bounds, n, func, config): +def _get_uniform_point_grid(bounds, n, func, config, nudge_factor=0): # Generate non-randomized grid of points linspaces = [] for (lb, ub), is_integer in bounds: if not is_integer: # Issues happen when exactly using the boundary - nudge = (ub - lb) * 1e-4 + nudge = (ub - lb) * nudge_factor linspaces.append(np.linspace(lb + nudge, ub - nudge, n)) else: size = min(n, ub - lb + 1) @@ -139,7 +139,7 @@ def _get_points_lmt_random_sample(bounds, n, func, config, seed=42): def _get_points_lmt_uniform_sample(bounds, n, func, config, seed=42): - points = _get_uniform_point_grid(bounds, n, func, config) + points = _get_uniform_point_grid(bounds, n, func, config, nudge_factor=1e-4) return _get_points_lmt(points, bounds, func, config, seed) @@ -345,6 +345,17 @@ def _find_leaves(splits, leaves, input_node): return leaves_list +class _Evaluator: + def __init__(self, expr, expr_vars): + self.expr = expr + self.expr_vars = expr_vars + + def __call__(self, *args): + for i, v in enumerate(self.expr_vars): + v.value = args[i] + return value(self.expr) + + @TransformationFactory.register( 'contrib.piecewise.nonlinear_to_pwl', doc="Convert nonlinear constraints and objectives to piecewise-linear " @@ -755,10 +766,7 @@ def _approximate_expression( continue # else we approximate subexpr - def eval_expr(*args): - for i, v in enumerate(expr_vars): - v.value = args[i] - return value(subexpr) + eval_expr = _Evaluator(subexpr, expr_vars) pwlf = _get_pwl_function_approximation( eval_expr, config, self._get_bounds_list(expr_vars, obj) diff --git a/pyomo/devel/initialization/README.md b/pyomo/devel/initialization/README.md new file mode 100644 index 00000000000..24ff0a0de86 --- /dev/null +++ b/pyomo/devel/initialization/README.md @@ -0,0 +1 @@ +The purpose of this module is to provide methods for initializing nonlinear programming models. \ No newline at end of file diff --git a/pyomo/devel/initialization/__init__.py b/pyomo/devel/initialization/__init__.py new file mode 100644 index 00000000000..0c9fd0b1950 --- /dev/null +++ b/pyomo/devel/initialization/__init__.py @@ -0,0 +1,10 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 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.devel.initialization.initialize import initialize_nlp, InitializationMethod diff --git a/pyomo/devel/initialization/bounds/__init__.py b/pyomo/devel/initialization/bounds/__init__.py new file mode 100644 index 00000000000..231b44987f6 --- /dev/null +++ b/pyomo/devel/initialization/bounds/__init__.py @@ -0,0 +1,8 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 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. +# ____________________________________________________________________________________ diff --git a/pyomo/devel/initialization/bounds/bound_variables.py b/pyomo/devel/initialization/bounds/bound_variables.py new file mode 100644 index 00000000000..ca6eabbc7f8 --- /dev/null +++ b/pyomo/devel/initialization/bounds/bound_variables.py @@ -0,0 +1,36 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 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.core.base.block import BlockData +from pyomo.contrib.fbbt.fbbt import fbbt +from pyomo.devel.initialization.utils import get_vars +import logging + +logger = logging.getLogger(__name__) + + +def bound_all_nonlinear_variables(m: BlockData, default_bound: float = 1.0e8): + """ + Attempt to obtain valid bounds on all nonlinear variables based on the + constraints in the model, m. If variable bounds cannot be obtained, + we use default_bound. + """ + fbbt(m) + for v in get_vars(m): + if v.lb is None or v.lb < -default_bound: + logger.debug( + f'Could not obtain a lower bound for {str(v)} better than {-default_bound}; setting the lower bound to {-default_bound}' + ) + v.setlb(-default_bound) + if v.ub is None or v.ub > default_bound: + logger.debug( + f'Could not obtain an upper bound for {str(v)} better than {default_bound}; setting the upper bound to {default_bound}' + ) + v.setub(default_bound) + fbbt(m) diff --git a/pyomo/devel/initialization/examples/__init__.py b/pyomo/devel/initialization/examples/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/pyomo/devel/initialization/examples/init_polynomial_ex.py b/pyomo/devel/initialization/examples/init_polynomial_ex.py new file mode 100644 index 00000000000..d327c6c9e3c --- /dev/null +++ b/pyomo/devel/initialization/examples/init_polynomial_ex.py @@ -0,0 +1,42 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 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. +# ____________________________________________________________________________________ + +# === Required imports === +import pyomo.environ as pyo +import pyomo.devel.initialization as ini +from pyomo.contrib.solver.common.factory import SolverFactory + + +def build_model(): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(-20, 20), initialize=-3.6) + m.c = pyo.Constraint(expr=(m.x + 7) * (m.x + 5) * (m.x - 4) + 200 == 0) + return m + + +def main(method: ini.InitializationMethod): + m = build_model() + nlp_solver = SolverFactory('ipopt') + global_solver = SolverFactory('scip_direct') + mip_solver = SolverFactory('scip_direct') + results = ini.initialize_nlp( + nlp=m, + nlp_solver=nlp_solver, + mip_solver=mip_solver, + global_solver=global_solver, + method=method, + ) + + return results.solution_status, m.x.value + + +if __name__ == '__main__': + stat, x = main(ini.InitializationMethod.global_opt) + print(stat) + print(round(x, 4)) diff --git a/pyomo/devel/initialization/global_init.py b/pyomo/devel/initialization/global_init.py new file mode 100644 index 00000000000..c573e4e8162 --- /dev/null +++ b/pyomo/devel/initialization/global_init.py @@ -0,0 +1,51 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 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.core.base.block import BlockData +from pyomo.contrib.solver.common.base import SolverBase +from pyomo.contrib.solver.common.results import SolutionStatus +from pyomo.contrib.solver.solvers.scip.scip_direct import ScipDirect, ScipPersistent +from pyomo.contrib.solver.solvers.gurobi.gurobi_direct_minlp import GurobiDirectMINLP +import logging + +logger = logging.getLogger(__name__) + + +def _initialize_with_global_solver( + nlp: BlockData, global_solver: SolverBase, nlp_solver: SolverBase +): + if isinstance(global_solver, (ScipDirect, ScipPersistent)): + opts = {'limits/solutions': 1} + elif isinstance(global_solver, (GurobiDirectMINLP,)): + opts = {'SolutionLimit': 1} + else: + raise NotImplementedError( + 'Currently, the initialization module only works with new solver interface, so the global solvers are limited to ScipDirect, ScipPersistent, and GurobiDirectMINLP.' + ) + res = global_solver.solve( + nlp, + load_solutions=True, + raise_exception_on_nonoptimal_result=False, + solver_options=opts, + ) + logger.info( + f'solved NLP with {global_solver.name}: {res.solution_status}, {res.termination_condition}' + ) + res = nlp_solver.solve( + nlp, load_solutions=False, raise_exception_on_nonoptimal_result=False + ) + logger.info( + f'solved NLP with {nlp_solver.name}: {res.solution_status}, {res.termination_condition}' + ) + if res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}: + res.solution_loader.load_vars() + else: + logger.warning('initialization was not successful via global optimization') + + return res diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py new file mode 100644 index 00000000000..759a4887f88 --- /dev/null +++ b/pyomo/devel/initialization/initialize.py @@ -0,0 +1,168 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 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 typing import Optional +from pyomo.core.base.block import BlockData +from enum import Enum +from pyomo.devel.initialization.utils import get_vars, shallow_clone +from pyomo.common.collections import ComponentMap +from pyomo.devel.initialization.pwl_init import ( + _initialize_with_piecewise_linear_approximation, +) +from pyomo.devel.initialization.lp_approx_init import _initialize_with_LP_approximation +from pyomo.contrib.solver.common.base import SolverBase +from pyomo.devel.initialization.global_init import _initialize_with_global_solver +from pyomo.contrib.solver.common.factory import SolverFactory +from pyomo.contrib.solver.common.results import Results +import logging +from pyomo.contrib.solver.common.results import SolutionStatus + +logger = logging.getLogger(__name__) + + +class InitializationMethod(Enum): + pwl_approximation = "pwl_approximation" + lp_approximation = "lp_approximation" + global_opt = "global_opt" + + +def _get_solver(sname, reason): + opt = SolverFactory(sname) + if opt.available(): + logger.info(f'Using {sname} for {reason} because a solver was not specified') + else: + raise RuntimeError( + f'No solver was specified for {reason} and the default ({sname}) is not available' + ) + return opt + + +def initialize_nlp( + nlp: BlockData, + nlp_solver: SolverBase, + mip_solver: Optional[SolverBase] = None, + global_solver: Optional[SolverBase] = None, + method: InitializationMethod = InitializationMethod.global_opt, + default_bound: float = 1.0e8, + max_pwl_refinement_iter: int = 100, + num_pwl_cons_to_refine_per_iter: int = 5, + aggressive_substitution: bool = True, +) -> Results: + """ + Attempt to initialize and subsequently solve the model given by ``nlp``. + The basic idea is to apply some method to find good initial values for + the variables and then try to solve the problem with ``nlp_solver``. + + Parameters + ---------- + nlp: BlockData + The pyomo model to be initialized. + mip_solver: Optional[SolverBase] + A solver interface appropriate for LPs and MILPs. + Needed for the following methods: + - pwl_approximation + - lp_approximation + Default: gurobi_persistent + nlp_solver: Optional[SolverBase] + A solver interface appropriate for NLPs. + Default: ipopt + global_solver: Optional[SolverBase] + A solver interface appropriate for global solution of NLPs + Default: gurobi_direct_minlp + method: InitializationMethod + The method used to initialize the model. + default_bound: float + Some initialize methods require all nonlinear variables to be bounded. + For these methods, all unbounded variables will be given lower and + upper bounds equal to default_bound. + Needed for the following methods: + - pwl_approximation + - lp_approximation + max_pwl_refinement_iter: int + Only used when method = InitializationMethod.pwl_approximation. This is + the maximum number of iterations used to refine the piecewise linear + approximation. + num_pwl_cons_to_refine_per_iter: int + Only used when method = InitializationMethod.pwl_approximation. This is + the maximum number of constraints to be refined with additional + segments in the piecewise linear approximation each iteration. + aggressive_substitution: bool + Only used when method = InitializationMethod.pwl_approximation. This is + passed along to the contrib.piecewise.univariate_nonlinear_decomposition + transformation. + + Returns + ------- + res: pyomo.contrib.solver.common.results.Results + The results object obtained the last time the nlp_solver was used to + try and solve the model. + """ + # in all cases, try to solve the nlp before doing extra work + res = nlp_solver.solve( + nlp, load_solutions=False, raise_exception_on_nonoptimal_result=False + ) + logger.info(f'solved NLP: {res.solution_status}, {res.termination_condition}') + + if res.solution_status == SolutionStatus.optimal: + res.solution_loader.load_vars() + logger.info('NLP solved without any initialization') + return res + + # get all variable bounds, domains, etc. to restore them later + orig_vars = get_vars(nlp) + orig_var_data = ComponentMap( + (v, (v.lower, v.upper, v.domain, v.fixed, v.value)) for v in orig_vars + ) + + # run the initialization + if method == InitializationMethod.pwl_approximation: + if mip_solver is None: + mip_solver = _get_solver('gurobi_persistent', 'MILP solver') + if nlp_solver is None: + nlp_solver = _get_solver('ipopt', 'local NLP solver') + res = _initialize_with_piecewise_linear_approximation( + nlp=nlp, + mip_solver=mip_solver, + nlp_solver=nlp_solver, + default_bound=default_bound, + max_iter=max_pwl_refinement_iter, + num_cons_to_refine_per_iter=num_pwl_cons_to_refine_per_iter, + aggressive_substitution=aggressive_substitution, + ) + elif method == InitializationMethod.lp_approximation: + if mip_solver is None: + mip_solver = _get_solver('gurobi_persistent', 'MILP solver') + if nlp_solver is None: + nlp_solver = _get_solver('ipopt', 'local NLP solver') + res = _initialize_with_LP_approximation( + nlp=nlp, lp_solver=mip_solver, nlp_solver=nlp_solver + ) + elif method == InitializationMethod.global_opt: + if global_solver is None: + global_solver = _get_solver('gurobi_direct_minlp', 'global NLP solver') + if nlp_solver is None: + nlp_solver = _get_solver('ipopt', 'local NLP solver') + res = _initialize_with_global_solver( + nlp=nlp, global_solver=global_solver, nlp_solver=nlp_solver + ) + else: + raise ValueError(f'unexpected initialization method: {method}') + + # restore variable bounds, domain, etc. + for v, (lb, ub, domain, fixed, value) in orig_var_data.items(): + v.setlb(lb) + v.setub(ub) + v.domain = domain + if fixed: + assert v.value == value + assert v.fixed + else: + v.unfix() + + return res diff --git a/pyomo/devel/initialization/lp_approx_init.py b/pyomo/devel/initialization/lp_approx_init.py new file mode 100644 index 00000000000..44bdcb52a5c --- /dev/null +++ b/pyomo/devel/initialization/lp_approx_init.py @@ -0,0 +1,216 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 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.core.base.block import BlockData +import pyomo.environ as pe +from pyomo.devel.initialization.bounds.bound_variables import ( + bound_all_nonlinear_variables, +) +from pyomo.devel.initialization.utils import ( + fix_vars_with_equal_bounds, + shallow_clone, + get_vars, +) +from pyomo.core.expr.visitor import identify_components +from pyomo.contrib.piecewise.piecewise_linear_expression import ( + PiecewiseLinearExpression, +) +from pyomo.contrib.piecewise.piecewise_linear_function import PiecewiseLinearFunction +from pyomo.contrib.solver.common.results import SolutionStatus +from pyomo.common.collections import ComponentMap, ComponentSet +from typing import MutableMapping, Sequence, List +from pyomo.core.base.constraint import ConstraintData +from pyomo.core.expr.visitor import StreamBasedExpressionVisitor +from pyomo.common.numeric_types import native_numeric_types +from pyomo.core.expr.numvalue import NumericConstant +from pyomo.core.expr.numeric_expr import ( + NegationExpression, + PowExpression, + ProductExpression, + MonomialTermExpression, + DivisionExpression, + SumExpression, + LinearExpression, + UnaryFunctionExpression, + NPV_NegationExpression, + NPV_PowExpression, + NPV_ProductExpression, + NPV_DivisionExpression, + NPV_SumExpression, + NPV_UnaryFunctionExpression, +) +from pyomo.core.expr.relational_expr import ( + EqualityExpression, + InequalityExpression, + RangedExpression, +) +from pyomo.repn.util import ExitNodeDispatcher +from pyomo.core.base.var import ScalarVar, VarData +from pyomo.core.base.param import ScalarParam, ParamData +from pyomo.core.base.expression import ScalarExpression, ExpressionData +import math +from pyomo.contrib.solver.common.base import SolverBase +import logging +from pyomo.common.modeling import unique_component_name +from pyomo.devel.initialization.pwl_init import _minimize_infeasibility +from pyomo.contrib.fbbt.fbbt import fbbt +from pyomo.repn.linear import LinearRepnVisitor, LinearRepn +from pyomo.core.expr.visitor import identify_variables +import numpy as np +from scipy.stats import qmc + +logger = logging.getLogger(__name__) + + +def _replace_expression_with_linear_approx(expr, num_samples=100): + vset = ComponentSet(identify_variables(expr, include_fixed=False)) + vlist = list(vset) + n_vars = len(vlist) + bnds_list = [] + for v in vlist: + if v.lb is None: + lb = -1e6 + else: + lb = v.lb + if v.ub is None: + ub = 1e6 + else: + ub = v.ub + bnds_list.append((lb, ub)) + sampler = qmc.LatinHypercube(d=n_vars) + sample = sampler.random(n=num_samples) + l_bounds = [i[0] for i in bnds_list] + u_bounds = [i[1] for i in bnds_list] + sample = qmc.scale(sample, l_bounds, u_bounds) + + # we have our samples + # now we want to build the matrix and the right hand side + # we have a linear coefficient for each variable plus a constant + n_coefs = n_vars + 1 + A = np.zeros((num_samples, n_coefs), dtype=float) + b = np.zeros(num_samples, dtype=float) + A[:, :n_vars] = sample + A[:, n_vars] = 1 + for sample_ndx in range(num_samples): + for v, val in zip(vlist, sample[sample_ndx, :]): + v.value = float(val) + b[sample_ndx] = pe.value(expr) + coefs = np.linalg.solve(A.transpose().dot(A), A.transpose().dot(b)) + coefs = [float(i) for i in coefs] + + new_expr = 0 + for c, v in zip(coefs[:n_vars], vlist): + new_expr += c * v + new_expr += coefs[-1] + return new_expr + + +def _build_lp_approx(nlp: BlockData) -> BlockData: + lp = pe.Block(concrete=True) + lp.cons = pe.ConstraintList() + visitor = LinearRepnVisitor(subexpression_cache={}) + + objs = list( + nlp.component_data_objects(pe.Objective, active=True, descend_into=True) + ) + if objs: + if len(objs) > 1: + raise NotImplementedError( + 'lp approximation does not support multiple objectives' + ) + obj = objs[0] + repn = visitor.walk_expression(obj) + assert repn.multiplier == 1 + linear_part = LinearRepn() + linear_part.multiplier = 1 + linear_part.constant = repn.constant + linear_part.linear = repn.linear + linear_part.nonlinear = None + new_obj_expr = linear_part.to_expression(visitor=visitor) + if repn.nonlinear is not None: + replacement = _replace_expression_with_linear_approx(repn.nonlinear) + new_obj_expr += replacement + lp.obj = pe.Objective(expr=new_obj_expr, sense=obj.sense) + + for con in nlp.component_data_objects( + pe.Constraint, active=True, descend_into=True + ): + lb, body, ub = con.to_bounded_expression() + repn = visitor.walk_expression(body) + assert repn.multiplier == 1 + linear_part = LinearRepn() + linear_part.multiplier = 1 + linear_part.constant = repn.constant + linear_part.linear = repn.linear + linear_part.nonlinear = None + new_body = linear_part.to_expression(visitor=visitor) + if repn.nonlinear is not None: + replacement = _replace_expression_with_linear_approx(repn.nonlinear) + new_body += replacement + if lb == ub: + lp.cons.add(new_body == lb) + else: + lp.cons.add((lb, new_body, ub)) + return lp + + +def _initialize_with_LP_approximation( + nlp: BlockData, lp_solver: SolverBase, nlp_solver: SolverBase, default_bound=1.0e8 +): + orig_nlp = nlp + logger.info('Starting initialization using a linear programming approximation') + nlp = shallow_clone(nlp) + logger.info('created a shallow clone of the model') + + # first introduce auxiliary variables so that we don't try to + # approximate any functions of more than two variables + # actually, this is not necessary for this method + # we will just comment this out for now + # trans = pe.TransformationFactory('contrib.piecewise.univariate_nonlinear_decomposition') + # trans.apply_to(nlp, aggressive_substitution=False) + # logger.info('applied the univariate_nonlinear_decomposition transformation') + + # bounds on the nonlinear variables + bound_all_nonlinear_variables(nlp, default_bound=default_bound) + logger.info('bounded nonlinear variables') + + # Now, we need to fix variables with equal (or nearly equal) bounds. + # Otherwise, the PWL transformation complains + fix_vars_with_equal_bounds(nlp) + logger.info('fixed variables with equal bounds') + + # now we modify the model by introducing slacks to make sure the LP + # approximation is feasible + _minimize_infeasibility(nlp) + logger.info('reformulated model to minimize infeasibility') + + # build the LP approximation + lp = _build_lp_approx(nlp) + logger.info('replaced nonlinear expressions with linear approximations') + + # solve the LP + lp_res = lp_solver.solve( + lp, load_solutions=True, raise_exception_on_nonoptimal_result=False + ) + logger.info(f'solved LP: {lp_res.solution_status}, {lp_res.termination_condition}') + + # try solving the NLP + nlp_res = nlp_solver.solve( + orig_nlp, load_solutions=False, raise_exception_on_nonoptimal_result=False + ) + logger.info( + f'solved NLP: {nlp_res.solution_status}, {nlp_res.termination_condition}' + ) + + if nlp_res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}: + nlp_res.solution_loader.load_vars() + else: + logger.warning('initialization was not successful via LP approximation') + + return nlp_res diff --git a/pyomo/devel/initialization/pwl_init.py b/pyomo/devel/initialization/pwl_init.py new file mode 100644 index 00000000000..08ca99b9a75 --- /dev/null +++ b/pyomo/devel/initialization/pwl_init.py @@ -0,0 +1,376 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 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.core.base.block import BlockData +import pyomo.environ as pe +from pyomo.devel.initialization.bounds.bound_variables import ( + bound_all_nonlinear_variables, +) +from pyomo.devel.initialization.utils import ( + fix_vars_with_equal_bounds, + shallow_clone, + get_vars, +) +from pyomo.core.expr.visitor import identify_components +from pyomo.contrib.piecewise.piecewise_linear_expression import ( + PiecewiseLinearExpression, +) +from pyomo.contrib.piecewise.piecewise_linear_function import PiecewiseLinearFunction +from pyomo.common.collections import ComponentMap, ComponentSet +from typing import MutableMapping, Sequence, List +from pyomo.core.base.constraint import ConstraintData +from pyomo.core.expr.visitor import StreamBasedExpressionVisitor +from pyomo.common.numeric_types import native_numeric_types +from pyomo.core.expr.numvalue import NumericConstant +from pyomo.core.expr.numeric_expr import ( + NegationExpression, + PowExpression, + ProductExpression, + MonomialTermExpression, + DivisionExpression, + SumExpression, + LinearExpression, + UnaryFunctionExpression, + NPV_NegationExpression, + NPV_PowExpression, + NPV_ProductExpression, + NPV_DivisionExpression, + NPV_SumExpression, + NPV_UnaryFunctionExpression, +) +from pyomo.core.expr.relational_expr import ( + EqualityExpression, + InequalityExpression, + RangedExpression, +) +from pyomo.repn.util import ExitNodeDispatcher +from pyomo.core.base.var import ScalarVar, VarData +from pyomo.core.base.param import ScalarParam, ParamData +from pyomo.core.base.expression import ScalarExpression, ExpressionData +import math +from pyomo.contrib.solver.common.base import SolverBase +import logging +from pyomo.common.modeling import unique_component_name +from pyomo.contrib.solver.common.results import SolutionStatus + +logger = logging.getLogger(__name__) + + +def _minimize_infeasibility(m): + m.slacks = pe.VarList() + m.extra_cons = pe.ConstraintList() + + obj_expr = 0 + + found_obj = False + for obj in m.component_data_objects(pe.Objective, active=True, descend_into=True): + assert not found_obj + if obj.sense == pe.minimize: + obj_expr += 0.1 * obj.expr + else: + obj_expr -= 0.1 * obj.expr + obj.deactivate() + found_obj = True + + for con in m.component_data_objects(pe.Constraint, active=True, descend_into=True): + lb, body, ub = con.to_bounded_expression(evaluate_bounds=True) + if lb == ub: + ps = m.slacks.add() + ns = m.slacks.add() + ps.setlb(0) + ns.setlb(0) + con.set_value(body - lb - ps + ns == 0) + elif lb is None: + ps = m.slacks.add() + ps.setlb(0) + con.set_value(body - ub - ps <= 0) + elif ub is None: + ns = m.slacks.add() + ns.setlb(0) + con.set_value(body - lb + ns >= 0) + else: + con.deactivate() + ps = m.slacks.add() + ns = m.slacks.add() + ps.setlb(0) + ns.setlb(0) + m.extra_cons.add(body - ub - ps <= 0) + m.extra_cons.add(body - lb + ns >= 0) + + m.slack_obj = pe.Objective(expr=10 * sum(m.slacks.values()) + obj_expr) + + +def _get_pwl_constraints( + m: BlockData, +) -> MutableMapping[PiecewiseLinearExpression, List[ConstraintData]]: + comp_types = set() + comp_types.add(PiecewiseLinearExpression) + pwl_expr_to_con_map = ComponentMap() + con_list = list( + m.component_data_objects(pe.Constraint, active=True, descend_into=True) + ) + obj_list = list( + m.component_data_objects(pe.Objective, active=True, descend_into=True) + ) + for comp in con_list + obj_list: + pwl_exprs = list(identify_components(comp.expr, comp_types)) + if not pwl_exprs: + continue + assert len(pwl_exprs) == 1 + e = pwl_exprs[0] + if e not in pwl_expr_to_con_map: + pwl_expr_to_con_map[e] = [] + pwl_expr_to_con_map[e].append(comp) + return pwl_expr_to_con_map + + +def _handle_leaf(node, data): + return node + + +def _handle_node(node, data): + return node.create_node_with_local_data(data) + + +_handlers = ExitNodeDispatcher() +for t in [float, int, VarData, ScalarVar, ParamData, ScalarParam, NumericConstant]: + _handlers[t] = _handle_leaf +for t in [ + ProductExpression, + SumExpression, + DivisionExpression, + PowExpression, + MonomialTermExpression, + LinearExpression, + ExpressionData, + ScalarExpression, + NegationExpression, + UnaryFunctionExpression, + NPV_NegationExpression, + NPV_PowExpression, + NPV_ProductExpression, + NPV_DivisionExpression, + NPV_SumExpression, + NPV_UnaryFunctionExpression, + EqualityExpression, + InequalityExpression, + RangedExpression, +]: + _handlers[t] = _handle_node + + +class _PWLRefinementVisitor(StreamBasedExpressionVisitor): + def __init__(self, m, pwl_exprs, **kwds): + self.m = m + self.pwl_exprs = ComponentSet(pwl_exprs) + self.substitution = ComponentMap() + self.named_expr_map = ComponentMap() + super().__init__(**kwds) + + def exitNode(self, node, data): + if node in self.named_expr_map: + return self.named_expr_map[node] + nt = type(node) + if nt in _handlers: + return _handlers[nt](node, data) + elif nt in native_numeric_types: + _handlers[nt] = _handle_leaf + return _handle_leaf(node, data) + else: + raise NotImplementedError(f'unrecognized expression type: {nt}') + + def beforeChild(self, node, child, child_idx): + if child in self.substitution: + return False, self.substitution[child] + + if child not in self.pwl_exprs: + return True, None + + old_func = child.pw_linear_function + _func = old_func._func + points = list(old_func._points) + variables = child.args + var_values = tuple(i.value for i in variables) + points.append(var_values) + points.sort() + if len(points[0]) == 1: + points = [i[0] for i in points] + new_func = PiecewiseLinearFunction(points=points, function=_func) + fname = unique_component_name( + self.m.auxiliary._pyomo_contrib_nonlinear_to_pwl, 'f' + ) + setattr(self.m.auxiliary._pyomo_contrib_nonlinear_to_pwl, fname, new_func) + new_expr = new_func(*variables) + for v, val in zip(variables, var_values): + v.set_value(val, skip_validation=True) + self.named_expr_map[node] = new_expr + self.substitution[child] = new_expr.expr + return False, new_expr.expr + + +def _refine_pwl_approx( + m, + pwl_expr_to_con_map: MutableMapping[ + PiecewiseLinearExpression, Sequence[ConstraintData] + ], + num_to_refine: int = 5, +): + violations = [] + for expr in pwl_expr_to_con_map.keys(): + func = expr.pw_linear_function + var_vals = [] + for v in expr.args: + if math.isclose(v.lb, v.ub, rel_tol=1e-6, abs_tol=1e-6): + val = 0.5 * (v.lb + v.ub) + elif v.value is None: + val = None + else: + val = v.value + if val <= v.lb + 1e-6 + 1e-6 * abs(v.lb): + val += 1e-5 + if val >= v.ub - 1e-6 - 1e-6 * abs(v.ub): + val -= 1e-5 + var_vals.append(val) + # var_vals = tuple(i.value for i in expr.args) + # for v, val in zip(expr.args, var_vals): + # print(f'{str(v):<20}{val:<20.5f}{v.lb:<20.5f}{v.ub:<20.5f}{id(v):<20}') + if any(i is None for i in var_vals): + logger.info(f'missing variable values for {expr}') + continue + approx_value = func(*var_vals) + true_value = func._func(*var_vals) + err = abs(true_value - approx_value) + violations.append((err, expr)) + violations.sort(key=lambda i: i[0], reverse=True) + + if len(violations) == 0: + raise RuntimeError( + 'Did not find any piecewise linear functions with variable values' + ) + + tol = 1e-5 + if math.isclose(violations[0][0], 0, abs_tol=tol): + logger.info('All of the original nonlinear functions are satisfied!') + + violations = [i for i in violations if i[0] > tol] + + for err, expr in violations[:num_to_refine]: + logger.info(f'refining {expr.pw_linear_function._func.expr} with error {err}') + + funcs_to_refine = ComponentSet(i[1] for i in violations[:num_to_refine]) + visitor = _PWLRefinementVisitor(m, funcs_to_refine) + + for expr in funcs_to_refine: + for con in pwl_expr_to_con_map[expr]: + con.set_value(visitor.walk_expression(con.expr)) + + for e1, e2 in visitor.substitution.items(): + cons = pwl_expr_to_con_map.pop(e1) + pwl_expr_to_con_map[e2] = cons + + +def _initialize_with_piecewise_linear_approximation( + nlp: BlockData, + mip_solver: SolverBase, + nlp_solver: SolverBase, + default_bound=1.0e8, + max_iter=100, + num_cons_to_refine_per_iter=5, + aggressive_substitution=True, +): + logger.info('Starting initialization using a piecewise linear approximation') + pwl = shallow_clone(nlp) + logger.info('created a shallow clone of the model') + + # first introduce auxiliary variables so that we don't try to + # approximate any functions of more than two variables + trans = pe.TransformationFactory( + 'contrib.piecewise.univariate_nonlinear_decomposition' + ) + trans.apply_to(pwl, aggressive_substitution=aggressive_substitution) + logger.info('applied the univariate_nonlinear_decomposition transformation') + + # now we need to try to get bounds on all of the nonlinear variables + bound_all_nonlinear_variables(pwl, default_bound=default_bound) + logger.info('bounded nonlinear variables') + + # Now, we need to fix variables with equal (or nearly equal) bounds. + # Otherwise, the PWL transformation complains + fix_vars_with_equal_bounds(pwl) + logger.info('fixed variables with equal bounds') + + # now we modify the model by introducing slacks to make sure the PWL + # approximation is feasible + # all of the slacks appear linearly, so we don't need to worry about + # upper bounds for them + _minimize_infeasibility(pwl) + logger.info('reformulated model to minimize infeasibility') + + # build the PWL approximation + trans = pe.TransformationFactory('contrib.piecewise.nonlinear_to_pwl') + trans.apply_to(pwl, num_points=2, additively_decompose=False) + logger.info('replaced nonlinear expressions with piecewise linear expressions') + + """ + Now we want to + 1. solve the PWL approximation + 2. Initialize the NLP to the solution + 3. Try solving the NLP + 4. If the NLP converges => done + 5. If the NLP does not converge, refine the PWL approximation and repeat + """ + pwl_expr_to_con_map = _get_pwl_constraints(pwl) + solved = False + last_nlp_res = None + for _iter in range(max_iter): + logger.info(f'PWL initialization: iter {_iter}') + + # PWL transformation (and map the variables) + orig_vars = list(get_vars(pwl)) + pwl.orig_vars = orig_vars + trans = pe.TransformationFactory('contrib.piecewise.disaggregated_logarithmic') + _pwl = trans.create_using(pwl) + new_vars = _pwl.orig_vars + del pwl.orig_vars + del _pwl.orig_vars + logger.info('applied the disaggregated logarithmic transformation') + + # solve the MILP + res = mip_solver.solve( + _pwl, load_solutions=True, raise_exception_on_nonoptimal_result=False + ) + logger.info(f'solved MILP: {res.solution_status}, {res.termination_condition}') + + # load the variable values back into orig_vars + for ov, nv in zip(orig_vars, new_vars): + ov.set_value(nv.value, skip_validation=True) + + # refine the PWL approximation + _refine_pwl_approx( + pwl, + pwl_expr_to_con_map=pwl_expr_to_con_map, + num_to_refine=num_cons_to_refine_per_iter, + ) + logger.info('refined PWL approximation') + + # try solving the NLP + res = nlp_solver.solve( + nlp, load_solutions=False, raise_exception_on_nonoptimal_result=False + ) + last_nlp_res = res + logger.info(f'solved NLP: {res.solution_status}, {res.termination_condition}') + if res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}: + solved = True + res.solution_loader.load_vars() + break + + if not solved: + logger.warning('initialization was not successful via PWL approximation') + + return last_nlp_res diff --git a/pyomo/devel/initialization/tests/__init__.py b/pyomo/devel/initialization/tests/__init__.py new file mode 100644 index 00000000000..231b44987f6 --- /dev/null +++ b/pyomo/devel/initialization/tests/__init__.py @@ -0,0 +1,8 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 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. +# ____________________________________________________________________________________ diff --git a/pyomo/devel/initialization/tests/test_initialization.py b/pyomo/devel/initialization/tests/test_initialization.py new file mode 100644 index 00000000000..bf78d3ee9f2 --- /dev/null +++ b/pyomo/devel/initialization/tests/test_initialization.py @@ -0,0 +1,206 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 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 typing import Tuple + +import pyomo.environ as pyo +import pyomo.devel.initialization as ini +from pyomo.devel.initialization.examples.init_polynomial_ex import main +from pyomo.common import unittest +from pyomo.contrib.solver.common.factory import SolverFactory +from pyomo.contrib.solver.common.results import ( + SolutionStatus, + Results, + TerminationCondition, +) +from pyomo.contrib.solver.common.base import Availability, SolverBase +import pytest + +scip = SolverFactory('scip_direct') +ipopt = SolverFactory('ipopt') + + +class MockNLPSolver(SolverBase): + def __init__(self, varlist, sol_map, **kwds) -> None: + super().__init__(**kwds) + self.varlist = varlist + self.sol_map = sol_map + self.iter = 0 + + def available(self) -> Availability: + return Availability.FullLicense + + def version(self) -> Tuple: + return (1, 0, 0) + + def check_solution(self): + expected, rel_tol, abs_tol = self.sol_map[self.iter] + self.iter += 1 + for v, val in zip(self.varlist, expected): + assert v.value == pytest.approx(val, rel=rel_tol, abs=abs_tol) + + def solve(self, model, **kwds) -> Results: + self.check_solution() + res = Results() + res.termination_condition = TerminationCondition.error + res.solution_status = SolutionStatus.noSolution + res.incumbent_objective = None + res.objective_bound = None + res.solver_name = 'MockNLPSolver' + res.solver_version = self.version() + return res + + +@unittest.skipUnless(scip.available(), 'scip is not available') +@unittest.skipUnless(ipopt.available(), 'ipopt is not available') +class TestExamples(unittest.TestCase): + def test_poly_global(self): + stat, x = main(method=ini.InitializationMethod.global_opt) + self.assertEqual(stat, SolutionStatus.optimal) + self.assertAlmostEqual(x, -9.920159607881597) + + def test_poly_pwl(self): + stat, x = main(method=ini.InitializationMethod.pwl_approximation) + self.assertEqual(stat, SolutionStatus.optimal) + self.assertAlmostEqual(x, -9.920159607881597) + + def test_poly_lp(self): + stat, x = main(method=ini.InitializationMethod.lp_approximation) + self.assertEqual(stat, SolutionStatus.optimal) + self.assertAlmostEqual(x, -9.920159607881597) + + +class TestInit(unittest.TestCase): + def test_lp_init(self): + """ + For this test, we will create a small linear program, + but we will make it look nonlinear. Then, the linear + approximation should be exact. The LP is + + max 3*x1 + 2*x2 + s.t. + x1 + x2 <= 4 + 2*x1 + x2 <= 5 + x1 >= 0 + x2 >= 0 + + The solution is + + x1 = 1 + x2 = 3 + """ + m = pyo.ConcreteModel() + m.x1 = pyo.Var(bounds=(0, 100)) + m.x2 = pyo.Var(bounds=(0, 100)) + m.obj = pyo.Objective( + expr=(3 * m.x1 * m.x1 + 2 * m.x2 * m.x1) / m.x1, sense=pyo.maximize + ) + m.c1 = pyo.Constraint(expr=pyo.exp(pyo.log(m.x1 + m.x2)) <= 4) + m.c2 = pyo.Constraint(expr=((2 * m.x1 + m.x2) ** 2) ** 0.5 <= 5) + + # all the actual testing happens in the MockNLPSolver + nlp_solver = MockNLPSolver( + varlist=[m.x1, m.x2], + sol_map={0: ([None, None], 0, 0), 1: ([1, 3], 1e-6, 1e-6)}, + ) + mip_solver = SolverFactory('highs') + results = ini.initialize_nlp( + nlp=m, + nlp_solver=nlp_solver, + mip_solver=mip_solver, + method=ini.InitializationMethod.lp_approximation, + ) + + def test_global_init(self): + """ + Same as test_lp_init + """ + m = pyo.ConcreteModel() + m.x1 = pyo.Var(bounds=(0, 100)) + m.x2 = pyo.Var(bounds=(0, 100)) + m.obj = pyo.Objective( + expr=(3 * m.x1 * m.x1 + 2 * m.x2 * m.x1) / m.x1, sense=pyo.maximize + ) + m.c1 = pyo.Constraint(expr=pyo.exp(pyo.log(m.x1 + m.x2)) <= 4) + m.c2 = pyo.Constraint(expr=((2 * m.x1 + m.x2) ** 2) ** 0.5 <= 5) + + # all the actual testing happens in the MockNLPSolver + nlp_solver = MockNLPSolver( + varlist=[m.x1, m.x2], + sol_map={0: ([None, None], 0, 0), 1: ([1, 3], 1e-6, 1e-6)}, + ) + global_solver = SolverFactory('scip_direct') + results = ini.initialize_nlp( + nlp=m, + nlp_solver=nlp_solver, + global_solver=global_solver, + method=ini.InitializationMethod.global_opt, + ) + + def test_pwl_init(self): + """ + Here, we really just want to make sure that the + approximation improves as refinement is done. + """ + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(-15, 5)) + m.c = pyo.Constraint(expr=(m.x + 7) * (m.x + 5) * (m.x - 4) + 200 == 0) + m.obj = pyo.Objective(expr=m.x) + + # all the actual testing happens in the MockNLPSolver + nlp_solver = MockNLPSolver( + varlist=[m.x], + sol_map={ + 0: ([None], 0, 0), + 1: ([1.0975609756097562], 1e-6, 1e-6), + 2: ([0.4346767574185112], 1e-6, 1e-6), + 3: ([-0.19286405313201946], 1e-6, 1e-6), + 4: ([-0.8653073960726083], 1e-6, 1e-6), + 5: ([-1.6404750700409576], 1e-6, 1e-6), + 6: ([-2.5676344169949443], 1e-6, 1e-6), + 7: ([-3.6759614495828297], 1e-6, 1e-6), + 8: ([-4.942429761325623], 1e-6, 1e-6), + 9: ([-6.259703235160286], 1e-6, 1e-6), + 10: ([-7.457220752001633], 1e-6, 1e-6), + 11: ([-8.393746738936832], 1e-6, 1e-6), + 12: ([-9.032852172775847], 1e-6, 1e-6), + 13: ([-9.426202540402329], 1e-6, 1e-6), + 14: ([-9.652335186743512], 1e-6, 1e-6), + 15: ([-9.777115808257673], 1e-6, 1e-6), + 16: ([-9.844390507666596], 1e-6, 1e-6), + 17: ([-9.880203709976758], 1e-6, 1e-6), + 18: ([-9.899139197799068], 1e-6, 1e-6), + 19: ([-9.90911480313665], 1e-6, 1e-6), + 20: ([-9.914360132504347], 1e-6, 1e-6), + 21: ([-9.91711543776506], 1e-6, 1e-6), + 22: ([-9.918562000338856], 1e-6, 1e-6), + 23: ([-9.919321249329018], 1e-6, 1e-6), + 24: ([-9.919719693944147], 1e-6, 1e-6), + 25: ([-9.91992877683681], 1e-6, 1e-6), + 26: ([-9.920038488200985], 1e-6, 1e-6), + 27: ([-9.920096055464825], 1e-6, 1e-6), + }, + ) + mip_solver = SolverFactory('highs') + results = ini.initialize_nlp( + nlp=m, + nlp_solver=nlp_solver, + mip_solver=mip_solver, + method=ini.InitializationMethod.pwl_approximation, + max_pwl_refinement_iter=27, + aggressive_substitution=False, + ) + + +if __name__ == '__main__': + import logging + + logging.basicConfig(level=logging.INFO) + t = TestInit() + t.test_pwl_init() diff --git a/pyomo/devel/initialization/utils.py b/pyomo/devel/initialization/utils.py new file mode 100644 index 00000000000..b8ce5095081 --- /dev/null +++ b/pyomo/devel/initialization/utils.py @@ -0,0 +1,53 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 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. +# ____________________________________________________________________________________ + +import pyomo.environ as pe +from pyomo.common.collections import ComponentSet +from pyomo.core.base.block import BlockData +from pyomo.core.expr.visitor import identify_variables +import math + + +def get_vars(m: BlockData): + vset = ComponentSet() + for c in m.component_data_objects(pe.Constraint, active=True, descend_into=True): + vset.update(identify_variables(c.body, include_fixed=False)) + for o in m.component_data_objects(pe.Objective, active=True, descend_into=True): + vset.update(identify_variables(o.expr, include_fixed=False)) + return vset + + +def shallow_clone(m1): + m2 = pe.ConcreteModel() + m2.cons = pe.ConstraintList() + + for con in m1.component_data_objects(pe.Constraint, active=True, descend_into=True): + m2.cons.add(con.expr) + + objlist = list( + m1.component_data_objects(pe.Objective, active=True, descend_into=True) + ) + assert len(objlist) <= 1 + if objlist: + obj = objlist[0] + m2.obj = pe.Objective(expr=obj.expr, sense=obj.sense) + + return m2 + + +def fix_vars_with_equal_bounds(m): + for v in get_vars(m): + if v.fixed: + continue + if ( + v.lb is not None + and v.ub is not None + and math.isclose(v.lb, v.ub, abs_tol=1e-4, rel_tol=1e-4) + ): + v.fix(0.5 * (v.lb + v.ub))