Initialization Module

doc/OnlineDocs/conf.py

@@ -428,6 +428,7 @@ def check_output(self, want, got, optionflags):
 baron_available = bool(_opt.check_available_solvers('baron'))
 glpk_available = bool(_opt.check_available_solvers('glpk'))
 gurobipy_available = bool(_opt.check_available_solvers('gurobi_direct'))
+pyscipopt_available = bool(_opt.check_available_solvers('scip_direct'))
 
 baron = _opt.SolverFactory('baron')
 

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:

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 <pyomo.devel.initialization.initialize.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 <pyomo.devel.initialization.initialize.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 <pyomo.devel.initialization.initialize.initialize_nlp>`:
+
+.. testcode::
+   :skipif: not pyscipopt_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 <pyomo.contrib.solver.solvers.scip.scip_direct.ScipDirect>` and
+  :class:`persistent <pyomo.contrib.solver.solvers.scip.scip_direct.ScipPersistent>`)
+* :class:`Gurobi MINLP <pyomo.contrib.solver.solvers.gurobi.gurobi_direct_minlp.GurobiDirectMINLP>`
+
+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.

doc/OnlineDocs/explanation/index.rst

@@ -43,6 +43,7 @@ Explanations
       `Design of Experiments`
       `MPC`
       `AOS`
+      `NLP Initialization`
    `Modeling Utilities`
       `Latex Printer`
       `FME`

pyomo/contrib/piecewise/tests/test_nonlinear_to_pwl.py

@@ -116,8 +116,8 @@ def test_log_constraint_uniform_grid(self):
         self.assertEqual(len(pwlf), 1)
         pwlf = pwlf[0]
 
-        points = [(1.0009,), (5.5,), (9.9991,)]
-        x1, x2, x3 = 1.0009, 5.5, 9.9991
+        points = [(1,), (5.5,), (10,)]
+        x1, x2, x3 = 1, 5.5, 10
         self.check_pw_linear_log_x(m, pwlf, x1, x2, x3)
 
     @unittest.skipUnless(numpy_available, "Numpy is not available")
@@ -142,8 +142,8 @@ def test_clone_transformed_model(self):
         self.assertEqual(len(pwlf), 1)
         pwlf = pwlf[0]
 
-        points = [(1.0009,), (5.5,), (9.9991,)]
-        x1, x2, x3 = 1.0009, 5.5, 9.9991
+        points = [(1,), (5.5,), (10,)]
+        x1, x2, x3 = 1, 5.5, 10
 
         self.check_pw_linear_log_x(twin, pwlf, x1, x2, x3)
 
@@ -197,8 +197,8 @@ def test_do_not_transform_quadratic_constraint(self):
         self.assertEqual(len(pwlf), 1)
         pwlf = pwlf[0]
 
-        points = [(1.0009,), (5.5,), (9.9991,)]
-        x1, x2, x3 = 1.0009, 5.5, 9.9991
+        points = [(1,), (5.5,), (10,)]
+        x1, x2, x3 = 1, 5.5, 10
         self.check_pw_linear_log_x(m, pwlf, x1, x2, x3)
 
         # quad is not
@@ -228,8 +228,8 @@ def test_constraint_target(self):
         self.assertEqual(len(pwlf), 1)
         pwlf = pwlf[0]
 
-        points = [(1.0009,), (5.5,), (9.9991,)]
-        x1, x2, x3 = 1.0009, 5.5, 9.9991
+        points = [(1,), (5.5,), (10,)]
+        x1, x2, x3 = 1, 5.5, 10
         self.check_pw_linear_log_x(m, pwlf, x1, x2, x3)
 
         # quad is not
@@ -501,10 +501,10 @@ def test_paraboloid_objective_uniform_grid(self):
         self.assertEqual(len(pwlf), 1)
         pwlf = pwlf[0]
 
-        x1 = 0.00030000000000000003
-        x2 = 2.9997
-        y1 = 1.0006
-        y2 = 6.9994
+        x1 = 0
+        x2 = 3
+        y1 = 1
+        y2 = 7
 
         self.check_pw_linear_paraboloid(m, pwlf, x1, x2, y1, y2)
 
@@ -531,10 +531,10 @@ def test_multivariate_clone(self):
         self.assertEqual(len(pwlf), 1)
         pwlf = pwlf[0]
 
-        x1 = 0.00030000000000000003
-        x2 = 2.9997
-        y1 = 1.0006
-        y2 = 6.9994
+        x1 = 0
+        x2 = 3
+        y1 = 1
+        y2 = 7
 
         self.check_pw_linear_paraboloid(twin, pwlf, x1, x2, y1, y2)
 
@@ -562,10 +562,10 @@ def test_objective_target(self):
         self.assertEqual(len(pwlf), 1)
         pwlf = pwlf[0]
 
-        x1 = 0.00030000000000000003
-        x2 = 2.9997
-        y1 = 1.0006
-        y2 = 6.9994
+        x1 = 0
+        x2 = 3
+        y1 = 1
+        y2 = 7
 
         self.check_pw_linear_paraboloid(m, pwlf, x1, x2, y1, y2)
 

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)

pyomo/devel/initialization/README.md

@@ -0,0 +1 @@
+The purpose of this module is to provide methods for initializing nonlinear programming models.

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

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.
+# ____________________________________________________________________________________

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)

pyomo/devel/initialization/examples/__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.
+# ____________________________________________________________________________________