From d29a2cfe9b985857e25c76c8a8f40be7be1c42aa Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Thu, 29 Jan 2026 15:24:10 -0700
Subject: [PATCH 01/38] starting a module for initialization

---
 pyomo/contrib/initialization/__init__.py | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 create mode 100644 pyomo/contrib/initialization/__init__.py

diff --git a/pyomo/contrib/initialization/__init__.py b/pyomo/contrib/initialization/__init__.py
new file mode 100644
index 00000000000..e69de29bb2d

From 8914a638203b048d28c06ef2d2dac7a6529c3396 Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Thu, 29 Jan 2026 16:13:14 -0700
Subject: [PATCH 02/38] starting a module for initialization

---
 pyomo/contrib/initialization/README.md        |  1 +
 .../contrib/initialization/bounds/__init__.py |  0
 .../initialization/bounds/bound_variables.py  | 10 +++
 pyomo/contrib/initialization/initialize.py    | 42 ++++++++++++
 pyomo/contrib/initialization/pwl_init.py      | 68 +++++++++++++++++++
 pyomo/contrib/initialization/utils.py         | 38 +++++++++++
 6 files changed, 159 insertions(+)
 create mode 100644 pyomo/contrib/initialization/README.md
 create mode 100644 pyomo/contrib/initialization/bounds/__init__.py
 create mode 100644 pyomo/contrib/initialization/bounds/bound_variables.py
 create mode 100644 pyomo/contrib/initialization/initialize.py
 create mode 100644 pyomo/contrib/initialization/pwl_init.py
 create mode 100644 pyomo/contrib/initialization/utils.py

diff --git a/pyomo/contrib/initialization/README.md b/pyomo/contrib/initialization/README.md
new file mode 100644
index 00000000000..24ff0a0de86
--- /dev/null
+++ b/pyomo/contrib/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/contrib/initialization/bounds/__init__.py b/pyomo/contrib/initialization/bounds/__init__.py
new file mode 100644
index 00000000000..e69de29bb2d
diff --git a/pyomo/contrib/initialization/bounds/bound_variables.py b/pyomo/contrib/initialization/bounds/bound_variables.py
new file mode 100644
index 00000000000..5eb59ec3c35
--- /dev/null
+++ b/pyomo/contrib/initialization/bounds/bound_variables.py
@@ -0,0 +1,10 @@
+from pyomo.core.base.block import BlockData
+
+
+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.
+    """
+    raise NotImplementedError('not done yet')
diff --git a/pyomo/contrib/initialization/initialize.py b/pyomo/contrib/initialization/initialize.py
new file mode 100644
index 00000000000..e20b18d51df
--- /dev/null
+++ b/pyomo/contrib/initialization/initialize.py
@@ -0,0 +1,42 @@
+from pyomo.core.base.block import BlockData
+from enum import Enum
+from pyomo.contrib.initialization.utils import get_vars, shallow_clone
+from pyomo.common.collections import ComponentMap
+from pyomo.contrib.initialization.pwl_init import _initialize_with_piecewise_linear_approximation
+
+
+class InitializationMethod(Enum):
+    pwl_approximation = "pwl_approximation"
+
+
+def initialize_nlp(
+    nlp: BlockData, 
+    method: InitializationMethod = InitializationMethod.pwl_approximation
+):
+    # 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
+    )
+
+    # create a shallow clone of the model so that the initialization method can 
+    # can work with the original variables but not make any other 
+    # modifications to the model
+    nlp = shallow_clone(nlp)
+
+    # run the initialization
+    if method == InitializationMethod.pwl_approximation:
+        _initialize_with_piecewise_linear_approximation(nlp)
+    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()
diff --git a/pyomo/contrib/initialization/pwl_init.py b/pyomo/contrib/initialization/pwl_init.py
new file mode 100644
index 00000000000..a5b0ee8aa84
--- /dev/null
+++ b/pyomo/contrib/initialization/pwl_init.py
@@ -0,0 +1,68 @@
+from pyomo.core.base.block import BlockData
+import pyomo.environ as pe
+from pyomo.contrib.initialization.bounds.bound_variables import bound_all_nonlinear_variables
+from pyomo.contrib.initialization.utils import fix_vars_with_equal_bounds
+
+
+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 _initialize_with_piecewise_linear_approximation(nlp: BlockData):
+    # 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(nlp, aggressive_substitution=True)
+
+    # now we need to try to get bounds on all of the nonlinear variables
+    bound_all_nonlinear_variables(nlp)
+
+    # Now, we need to fix variables with equal (or nearly equal) bounds.
+    # Otherwise, the PWL transformation complains
+    fix_vars_with_equal_bounds(nlp)
+
+    # now we modify the model by introducing slacks to make sure the PWL
+    # approximatin is feasible
+    _minimize_infeasibility(nlp)
+
+    raise NotImplementedError('not done yet')
diff --git a/pyomo/contrib/initialization/utils.py b/pyomo/contrib/initialization/utils.py
new file mode 100644
index 00000000000..792509583f0
--- /dev/null
+++ b/pyomo/contrib/initialization/utils.py
@@ -0,0 +1,38 @@
+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(c.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))

From 5ebfdd4dc3dda7407f89b2032cbfac939b865804 Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Tue, 3 Feb 2026 08:26:01 -0700
Subject: [PATCH 03/38] initialization

---
 .../initialization/bounds/bound_variables.py  | 16 +++++++++-
 pyomo/contrib/initialization/pwl_init.py      | 32 +++++++++++++++----
 2 files changed, 41 insertions(+), 7 deletions(-)

diff --git a/pyomo/contrib/initialization/bounds/bound_variables.py b/pyomo/contrib/initialization/bounds/bound_variables.py
index 5eb59ec3c35..0c243b8eb33 100644
--- a/pyomo/contrib/initialization/bounds/bound_variables.py
+++ b/pyomo/contrib/initialization/bounds/bound_variables.py
@@ -1,4 +1,10 @@
 from pyomo.core.base.block import BlockData
+from pyomo.contrib.fbbt.fbbt import fbbt
+from pyomo.contrib.initialization.utils import get_vars
+import logging
+
+
+logger = logging.getLogger(__name__)
 
 
 def bound_all_nonlinear_variables(m: BlockData, default_bound: float = 1.0e8):
@@ -7,4 +13,12 @@ def bound_all_nonlinear_variables(m: BlockData, default_bound: float = 1.0e8):
     constraints in the model, m. If variable bounds cannot be obtained, 
     we use default_bound.
     """
-    raise NotImplementedError('not done yet')
+    fbbt(m)
+    for v in get_vars(m):
+        if v.lb is None or v.lb < -default_bound:
+            logger.warning(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.warning(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/contrib/initialization/pwl_init.py b/pyomo/contrib/initialization/pwl_init.py
index a5b0ee8aa84..969adbecd99 100644
--- a/pyomo/contrib/initialization/pwl_init.py
+++ b/pyomo/contrib/initialization/pwl_init.py
@@ -1,7 +1,7 @@
 from pyomo.core.base.block import BlockData
 import pyomo.environ as pe
 from pyomo.contrib.initialization.bounds.bound_variables import bound_all_nonlinear_variables
-from pyomo.contrib.initialization.utils import fix_vars_with_equal_bounds
+from pyomo.contrib.initialization.utils import fix_vars_with_equal_bounds, shallow_clone
 
 
 def _minimize_infeasibility(m):
@@ -48,21 +48,41 @@ def _minimize_infeasibility(m):
     m.slack_obj = pe.Objective(expr=10*sum(m.slacks.values()) + obj_expr)
 
 
-def _initialize_with_piecewise_linear_approximation(nlp: BlockData):
+def _refine_pwl_approx(m):
+
+
+def _initialize_with_piecewise_linear_approximation(nlp: BlockData, default_bound=1.0e8):
+    pwl = shallow_clone(nlp)
+
     # 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(nlp, aggressive_substitution=True)
+    trans.apply_to(pwl, aggressive_substitution=True)
 
     # now we need to try to get bounds on all of the nonlinear variables
-    bound_all_nonlinear_variables(nlp)
+    bound_all_nonlinear_variables(pwl, default_bound=default_bound)
 
     # Now, we need to fix variables with equal (or nearly equal) bounds.
     # Otherwise, the PWL transformation complains
-    fix_vars_with_equal_bounds(nlp)
+    fix_vars_with_equal_bounds(pwl)
 
     # now we modify the model by introducing slacks to make sure the PWL
     # approximatin is feasible
-    _minimize_infeasibility(nlp)
+    # all of the slacks appear linearly, so we don't need to worry about 
+    # upper bounds for them
+    _minimize_infeasibility(pwl)
+
+    # build the PWL approximation
+    trans = pe.TransformationFactory('contrib.piecewise.nonlinear_to_pwl')
+    trans.apply_to(pwl, num_points=num_points, additively_decompose=False)
+
+    """
+    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
+    """
 
     raise NotImplementedError('not done yet')

From 488a539925fecda3e57f08e894650e7e6d80cf0c Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Tue, 3 Feb 2026 18:57:28 -0700
Subject: [PATCH 04/38] pwl initialization is almost working

---
 pyomo/contrib/initialization/initialize.py    |  17 +-
 pyomo/contrib/initialization/pwl_init.py      | 237 +++++++++++++++++-
 .../piecewise/transform/nonlinear_to_pwl.py   |  16 +-
 3 files changed, 259 insertions(+), 11 deletions(-)

diff --git a/pyomo/contrib/initialization/initialize.py b/pyomo/contrib/initialization/initialize.py
index e20b18d51df..7b7e8af7338 100644
--- a/pyomo/contrib/initialization/initialize.py
+++ b/pyomo/contrib/initialization/initialize.py
@@ -3,6 +3,7 @@
 from pyomo.contrib.initialization.utils import get_vars, shallow_clone
 from pyomo.common.collections import ComponentMap
 from pyomo.contrib.initialization.pwl_init import _initialize_with_piecewise_linear_approximation
+from pyomo.contrib.solver.common.base import SolverBase
 
 
 class InitializationMethod(Enum):
@@ -11,7 +12,12 @@ class InitializationMethod(Enum):
 
 def initialize_nlp(
     nlp: BlockData, 
-    method: InitializationMethod = InitializationMethod.pwl_approximation
+    mip_solver: SolverBase,
+    nlp_solver: SolverBase,
+    method: InitializationMethod = InitializationMethod.pwl_approximation,
+    default_bound=1.0e8,
+    max_pwl_refinement_iter=100,
+    num_pwl_cons_to_refine_per_iter=5,
 ):
     # get all variable bounds, domains, etc. to restore them later
     orig_vars = get_vars(nlp)
@@ -26,7 +32,14 @@ def initialize_nlp(
 
     # run the initialization
     if method == InitializationMethod.pwl_approximation:
-        _initialize_with_piecewise_linear_approximation(nlp)
+        _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,
+        )
     else:
         raise ValueError(f'unexpected initialization method: {method}')    
 
diff --git a/pyomo/contrib/initialization/pwl_init.py b/pyomo/contrib/initialization/pwl_init.py
index 969adbecd99..d6fdae92c25 100644
--- a/pyomo/contrib/initialization/pwl_init.py
+++ b/pyomo/contrib/initialization/pwl_init.py
@@ -1,7 +1,44 @@
 from pyomo.core.base.block import BlockData
 import pyomo.environ as pe
 from pyomo.contrib.initialization.bounds.bound_variables import bound_all_nonlinear_variables
-from pyomo.contrib.initialization.utils import fix_vars_with_equal_bounds, shallow_clone
+from pyomo.contrib.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
+
+
+logger = logging.getLogger(__name__)
 
 
 def _minimize_infeasibility(m):
@@ -48,33 +85,184 @@ def _minimize_infeasibility(m):
     m.slack_obj = pe.Objective(expr=10*sum(m.slacks.values()) + obj_expr)
 
 
-def _refine_pwl_approx(m):
+def _get_pwl_constraints(m: BlockData) -> MutableMapping[
+    PiecewiseLinearExpression, 
+    List[ConstraintData]
+]:
+    comp_types = set()
+    comp_types.add(PiecewiseLinearExpression)
+    pwl_expr_to_con_map = ComponentMap()
+    for con in m.component_data_objects(pe.Constraint, active=True, descend_into=True):
+        pwl_exprs = list(identify_components(con.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(con)
+    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)
+        self.named_expr_map[node] = new_expr
+        self.substitution[child] = new_expr.expr
+        return False, new_expr.expr
 
-def _initialize_with_piecewise_linear_approximation(nlp: BlockData, default_bound=1.0e8):
+
+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 = tuple(i.value for i in expr.args)
+        if any(i is None for i in var_vals):
+            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')
+    
+    if math.isclose(violations[0][0], 0):
+        raise RuntimeError('All of the original nonlinear functions are satisfied!')
+
+    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,
+):
+    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=True)
+    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
     # approximatin 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=num_points, additively_decompose=False)
+    trans.apply_to(pwl, num_points=2, additively_decompose=False)
+    logger.info('replaced nonlinear expressions with piecewise linear expressions')
 
     """
     Now we want to 
@@ -84,5 +272,44 @@ def _initialize_with_piecewise_linear_approximation(nlp: BlockData, default_boun
     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
+    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)
+        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.value = nv.value
+
+        # 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)
+        logger.info(f'solved NLP: {res.solution_status}, {res.termination_condition}')
+        if res.incumbent_objective is not None:
+            solved = True
+            res.solution_loader.load_vars()
+            # break
 
-    raise NotImplementedError('not done yet')
+    if not solved:
+        raise RuntimeError('no feasible solution found')
diff --git a/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py b/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py
index 65936eda109..e6239a1dcda 100644
--- a/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py
+++ b/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py
@@ -347,6 +347,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 "
@@ -742,10 +753,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)

From e99362341fc69dec82359d12071e549f00a54d13 Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Tue, 3 Feb 2026 20:58:24 -0700
Subject: [PATCH 05/38] debugging

---
 pyomo/contrib/initialization/initialize.py | 2 +-
 pyomo/contrib/initialization/pwl_init.py   | 6 +++++-
 2 files changed, 6 insertions(+), 2 deletions(-)

diff --git a/pyomo/contrib/initialization/initialize.py b/pyomo/contrib/initialization/initialize.py
index 7b7e8af7338..6bc0c9e9ff6 100644
--- a/pyomo/contrib/initialization/initialize.py
+++ b/pyomo/contrib/initialization/initialize.py
@@ -28,7 +28,7 @@ def initialize_nlp(
     # create a shallow clone of the model so that the initialization method can 
     # can work with the original variables but not make any other 
     # modifications to the model
-    nlp = shallow_clone(nlp)
+    # nlp = shallow_clone(nlp)
 
     # run the initialization
     if method == InitializationMethod.pwl_approximation:
diff --git a/pyomo/contrib/initialization/pwl_init.py b/pyomo/contrib/initialization/pwl_init.py
index d6fdae92c25..b70b82a2838 100644
--- a/pyomo/contrib/initialization/pwl_init.py
+++ b/pyomo/contrib/initialization/pwl_init.py
@@ -179,6 +179,8 @@ def beforeChild(self, node, child, child_idx):
         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.value = val
         self.named_expr_map[node] = new_expr
         self.substitution[child] = new_expr.expr
         return False, new_expr.expr
@@ -196,6 +198,8 @@ def _refine_pwl_approx(
     for expr in pwl_expr_to_con_map.keys():
         func = expr.pw_linear_function
         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):
             continue
         approx_value = func(*var_vals)
@@ -304,7 +308,7 @@ def _initialize_with_piecewise_linear_approximation(
         logger.info('refined PWL approximation')
 
         # try solving the NLP
-        res = nlp_solver.solve(nlp, load_solutions=False, raise_exception_on_nonoptimal_result=False)
+        res = nlp_solver.solve(nlp, tee=True, load_solutions=False, raise_exception_on_nonoptimal_result=False)
         logger.info(f'solved NLP: {res.solution_status}, {res.termination_condition}')
         if res.incumbent_objective is not None:
             solved = True

From aa450e1cedb4f69fbc0b5a9ddd8f8f84416e5138 Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Tue, 3 Feb 2026 21:07:56 -0700
Subject: [PATCH 06/38] initialization

---
 pyomo/contrib/initialization/pwl_init.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyomo/contrib/initialization/pwl_init.py b/pyomo/contrib/initialization/pwl_init.py
index b70b82a2838..65885f0ac5a 100644
--- a/pyomo/contrib/initialization/pwl_init.py
+++ b/pyomo/contrib/initialization/pwl_init.py
@@ -313,7 +313,7 @@ def _initialize_with_piecewise_linear_approximation(
         if res.incumbent_objective is not None:
             solved = True
             res.solution_loader.load_vars()
-            # break
+            break
 
     if not solved:
         raise RuntimeError('no feasible solution found')

From 043f8c8a6da16541116fb97e2ca21b586edece4b Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Mon, 23 Mar 2026 11:25:11 -0600
Subject: [PATCH 07/38] some updates to initialization module

---
 .../initialization/bounds/bound_variables.py  |   4 +-
 pyomo/contrib/initialization/initialize.py    |  12 +-
 .../contrib/initialization/lp_approx_init.py  | 182 ++++++++++++++++++
 pyomo/contrib/initialization/pwl_init.py      |  27 ++-
 4 files changed, 217 insertions(+), 8 deletions(-)
 create mode 100644 pyomo/contrib/initialization/lp_approx_init.py

diff --git a/pyomo/contrib/initialization/bounds/bound_variables.py b/pyomo/contrib/initialization/bounds/bound_variables.py
index 0c243b8eb33..0d0232f53b8 100644
--- a/pyomo/contrib/initialization/bounds/bound_variables.py
+++ b/pyomo/contrib/initialization/bounds/bound_variables.py
@@ -16,9 +16,9 @@ def bound_all_nonlinear_variables(m: BlockData, default_bound: float = 1.0e8):
     fbbt(m)
     for v in get_vars(m):
         if v.lb is None or v.lb < -default_bound:
-            logger.warning(f'Could not obtain a lower bound for {str(v)} better than {-default_bound}; setting the lower bound to {-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.warning(f'Could not obtain an upper bound for {str(v)} better than {default_bound}; setting the upper bound to {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/contrib/initialization/initialize.py b/pyomo/contrib/initialization/initialize.py
index 6bc0c9e9ff6..baf5accaec6 100644
--- a/pyomo/contrib/initialization/initialize.py
+++ b/pyomo/contrib/initialization/initialize.py
@@ -3,11 +3,13 @@
 from pyomo.contrib.initialization.utils import get_vars, shallow_clone
 from pyomo.common.collections import ComponentMap
 from pyomo.contrib.initialization.pwl_init import _initialize_with_piecewise_linear_approximation
+from pyomo.contrib.initialization.lp_approx_init import _initialize_with_LP_approximation
 from pyomo.contrib.solver.common.base import SolverBase
 
 
 class InitializationMethod(Enum):
     pwl_approximation = "pwl_approximation"
+    lp_approximation = "lp_approximation"
 
 
 def initialize_nlp(
@@ -32,7 +34,7 @@ def initialize_nlp(
 
     # run the initialization
     if method == InitializationMethod.pwl_approximation:
-        _initialize_with_piecewise_linear_approximation(
+        res = _initialize_with_piecewise_linear_approximation(
             nlp=nlp,
             mip_solver=mip_solver,
             nlp_solver=nlp_solver,
@@ -40,6 +42,12 @@ def initialize_nlp(
             max_iter=max_pwl_refinement_iter,
             num_cons_to_refine_per_iter=num_pwl_cons_to_refine_per_iter,
         )
+    elif method == InitializationMethod.lp_approximation:
+        res = _initialize_with_LP_approximation(
+            nlp=nlp,
+            lp_solver=mip_solver,
+            nlp_solver=nlp_solver,
+        )
     else:
         raise ValueError(f'unexpected initialization method: {method}')    
 
@@ -53,3 +61,5 @@ def initialize_nlp(
             assert v.fixed
         else:
             v.unfix()
+
+    return res
diff --git a/pyomo/contrib/initialization/lp_approx_init.py b/pyomo/contrib/initialization/lp_approx_init.py
new file mode 100644
index 00000000000..63d96011987
--- /dev/null
+++ b/pyomo/contrib/initialization/lp_approx_init.py
@@ -0,0 +1,182 @@
+from pyomo.core.base.block import BlockData
+import pyomo.environ as pe
+from pyomo.contrib.initialization.bounds.bound_variables import bound_all_nonlinear_variables
+from pyomo.contrib.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.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,
+):
+    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
+    trans = pe.TransformationFactory('contrib.piecewise.univariate_nonlinear_decomposition')
+    trans.apply_to(nlp, aggressive_substitution=False)
+    logger.info('applied the univariate_nonlinear_decomposition transformation')
+
+    # let's see if we can get bounds on the nonlinear variables
+    fbbt(nlp)
+    logger.info('ran FBBT')
+
+    # 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
+    # approximatin 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.incumbent_objective is not None:
+        nlp_res.solution_loader.load_vars()
+    else:
+        raise RuntimeError('no feasible solution found')
+
+    return nlp_res
diff --git a/pyomo/contrib/initialization/pwl_init.py b/pyomo/contrib/initialization/pwl_init.py
index 65885f0ac5a..30baa1d7cef 100644
--- a/pyomo/contrib/initialization/pwl_init.py
+++ b/pyomo/contrib/initialization/pwl_init.py
@@ -180,7 +180,7 @@ def beforeChild(self, node, child, child_idx):
         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.value = val
+            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
@@ -197,7 +197,20 @@ def _refine_pwl_approx(
     violations = []
     for expr in pwl_expr_to_con_map.keys():
         func = expr.pw_linear_function
-        var_vals = tuple(i.value for i in expr.args)
+        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):
@@ -278,6 +291,7 @@ def _initialize_with_piecewise_linear_approximation(
     """
     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}')
 
@@ -292,12 +306,12 @@ def _initialize_with_piecewise_linear_approximation(
         logger.info('applied the disaggregated logarithmic transformation')
 
         # solve the MILP
-        res = mip_solver.solve(_pwl, load_solutions=True)
+        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.value = nv.value
+            ov.set_value(nv.value, skip_validation=True)
 
         # refine the PWL approximation
         _refine_pwl_approx(
@@ -308,7 +322,8 @@ def _initialize_with_piecewise_linear_approximation(
         logger.info('refined PWL approximation')
 
         # try solving the NLP
-        res = nlp_solver.solve(nlp, tee=True, load_solutions=False, raise_exception_on_nonoptimal_result=False)
+        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.incumbent_objective is not None:
             solved = True
@@ -317,3 +332,5 @@ def _initialize_with_piecewise_linear_approximation(
 
     if not solved:
         raise RuntimeError('no feasible solution found')
+
+    return last_nlp_res

From c8a07305a16cba2abeec0333d89f43b9c9269ade Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Mon, 23 Mar 2026 15:45:32 -0600
Subject: [PATCH 08/38] initialize with global solvers

---
 pyomo/contrib/initialization/global_init.py | 28 +++++++++++++++++++++
 pyomo/contrib/initialization/initialize.py  | 15 +++++++++--
 2 files changed, 41 insertions(+), 2 deletions(-)
 create mode 100644 pyomo/contrib/initialization/global_init.py

diff --git a/pyomo/contrib/initialization/global_init.py b/pyomo/contrib/initialization/global_init.py
new file mode 100644
index 00000000000..e982570a2c9
--- /dev/null
+++ b/pyomo/contrib/initialization/global_init.py
@@ -0,0 +1,28 @@
+from pyomo.core.base.block import BlockData
+from pyomo.contrib.solver.common.base import SolverBase
+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,
+):
+    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=False, raise_exception_on_nonoptimal_result=False, solver_options=opts)
+    logger.info(f'solved NLP: {res.solution_status}, {res.termination_condition}')
+    if res.incumbent_objective is not None:
+        res.solution_loader.load_vars()
+    else:
+        raise RuntimeError('no feasible solution found')
+    
+    return res
diff --git a/pyomo/contrib/initialization/initialize.py b/pyomo/contrib/initialization/initialize.py
index baf5accaec6..b7e04d784d1 100644
--- a/pyomo/contrib/initialization/initialize.py
+++ b/pyomo/contrib/initialization/initialize.py
@@ -1,3 +1,4 @@
+from typing import Optional
 from pyomo.core.base.block import BlockData
 from enum import Enum
 from pyomo.contrib.initialization.utils import get_vars, shallow_clone
@@ -5,17 +6,20 @@
 from pyomo.contrib.initialization.pwl_init import _initialize_with_piecewise_linear_approximation
 from pyomo.contrib.initialization.lp_approx_init import _initialize_with_LP_approximation
 from pyomo.contrib.solver.common.base import SolverBase
+from pyomo.contrib.initialization.global_init import _initialize_with_global_solver
 
 
 class InitializationMethod(Enum):
     pwl_approximation = "pwl_approximation"
     lp_approximation = "lp_approximation"
+    global_opt = "global_opt"
 
 
 def initialize_nlp(
     nlp: BlockData, 
-    mip_solver: SolverBase,
-    nlp_solver: SolverBase,
+    mip_solver: Optional[SolverBase] = None,
+    nlp_solver: Optional[SolverBase] = None,
+    global_solver: Optional[SolverBase] = None,
     method: InitializationMethod = InitializationMethod.pwl_approximation,
     default_bound=1.0e8,
     max_pwl_refinement_iter=100,
@@ -34,6 +38,8 @@ def initialize_nlp(
 
     # run the initialization
     if method == InitializationMethod.pwl_approximation:
+        assert mip_solver is not None, f"mip_solver must be specified for {method}"
+        assert nlp_solver is not None, f"nlp_solver must be specified for {method}"
         res = _initialize_with_piecewise_linear_approximation(
             nlp=nlp,
             mip_solver=mip_solver,
@@ -43,11 +49,16 @@ def initialize_nlp(
             num_cons_to_refine_per_iter=num_pwl_cons_to_refine_per_iter,
         )
     elif method == InitializationMethod.lp_approximation:
+        assert mip_solver is not None, f"mip_solver must be specified for {method}"
+        assert nlp_solver is not None, f"nlp_solver must be specified for {method}"
         res = _initialize_with_LP_approximation(
             nlp=nlp,
             lp_solver=mip_solver,
             nlp_solver=nlp_solver,
         )
+    elif method == InitializationMethod.global_opt:
+        assert global_solver is not None, f"global_solver must be specified for {method}"
+        res = _initialize_with_global_solver(nlp=nlp, global_solver=global_solver)
     else:
         raise ValueError(f'unexpected initialization method: {method}')    
 

From e019bffc8f40fc2583623d3ff3cbd6a02b0b9233 Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Tue, 24 Mar 2026 07:15:46 -0600
Subject: [PATCH 09/38] initialization work

---
 pyomo/contrib/initialization/global_init.py |  7 ++--
 pyomo/contrib/initialization/initialize.py  | 40 ++++++++++++++-------
 2 files changed, 33 insertions(+), 14 deletions(-)

diff --git a/pyomo/contrib/initialization/global_init.py b/pyomo/contrib/initialization/global_init.py
index e982570a2c9..4f858a4b7c3 100644
--- a/pyomo/contrib/initialization/global_init.py
+++ b/pyomo/contrib/initialization/global_init.py
@@ -11,6 +11,7 @@
 def _initialize_with_global_solver(
     nlp: BlockData,
     global_solver: SolverBase,
+    nlp_solver: SolverBase,
 ):
     if isinstance(global_solver, (ScipDirect, ScipPersistent)):
         opts = {'limits/solutions': 1}
@@ -18,8 +19,10 @@ def _initialize_with_global_solver(
         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=False, raise_exception_on_nonoptimal_result=False, solver_options=opts)
-    logger.info(f'solved NLP: {res.solution_status}, {res.termination_condition}')
+    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.incumbent_objective is not None:
         res.solution_loader.load_vars()
     else:
diff --git a/pyomo/contrib/initialization/initialize.py b/pyomo/contrib/initialization/initialize.py
index b7e04d784d1..f8918f95528 100644
--- a/pyomo/contrib/initialization/initialize.py
+++ b/pyomo/contrib/initialization/initialize.py
@@ -7,6 +7,11 @@
 from pyomo.contrib.initialization.lp_approx_init import _initialize_with_LP_approximation
 from pyomo.contrib.solver.common.base import SolverBase
 from pyomo.contrib.initialization.global_init import _initialize_with_global_solver
+from pyomo.contrib.solver.common.factory import SolverFactory
+import logging
+
+
+logger = logging.getLogger(__name__)
 
 
 class InitializationMethod(Enum):
@@ -15,12 +20,21 @@ class InitializationMethod(Enum):
     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, 
     mip_solver: Optional[SolverBase] = None,
     nlp_solver: Optional[SolverBase] = None,
     global_solver: Optional[SolverBase] = None,
-    method: InitializationMethod = InitializationMethod.pwl_approximation,
+    method: InitializationMethod = InitializationMethod.global_opt,
     default_bound=1.0e8,
     max_pwl_refinement_iter=100,
     num_pwl_cons_to_refine_per_iter=5,
@@ -31,15 +45,12 @@ def initialize_nlp(
         (v, (v.lower, v.upper, v.domain, v.fixed, v.value)) for v in orig_vars
     )
 
-    # create a shallow clone of the model so that the initialization method can 
-    # can work with the original variables but not make any other 
-    # modifications to the model
-    # nlp = shallow_clone(nlp)
-
     # run the initialization
     if method == InitializationMethod.pwl_approximation:
-        assert mip_solver is not None, f"mip_solver must be specified for {method}"
-        assert nlp_solver is not None, f"nlp_solver must be specified for {method}"
+        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,
@@ -49,16 +60,21 @@ def initialize_nlp(
             num_cons_to_refine_per_iter=num_pwl_cons_to_refine_per_iter,
         )
     elif method == InitializationMethod.lp_approximation:
-        assert mip_solver is not None, f"mip_solver must be specified for {method}"
-        assert nlp_solver is not None, f"nlp_solver must be specified for {method}"
+        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:
-        assert global_solver is not None, f"global_solver must be specified for {method}"
-        res = _initialize_with_global_solver(nlp=nlp, global_solver=global_solver)
+        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}')    
 

From 9568bb0d6f9812bf3bbf0baeb35daa42bab25993 Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Tue, 24 Mar 2026 07:37:49 -0600
Subject: [PATCH 10/38] lp init should bound all variables

---
 pyomo/contrib/initialization/lp_approx_init.py | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/pyomo/contrib/initialization/lp_approx_init.py b/pyomo/contrib/initialization/lp_approx_init.py
index 63d96011987..0c72e82fbd7 100644
--- a/pyomo/contrib/initialization/lp_approx_init.py
+++ b/pyomo/contrib/initialization/lp_approx_init.py
@@ -137,6 +137,7 @@ 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')
@@ -149,9 +150,9 @@ def _initialize_with_LP_approximation(
     trans.apply_to(nlp, aggressive_substitution=False)
     logger.info('applied the univariate_nonlinear_decomposition transformation')
 
-    # let's see if we can get bounds on the nonlinear variables
-    fbbt(nlp)
-    logger.info('ran FBBT')
+    # 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

From 387a646a36aa607b3d4c388c6c19de41b1654251 Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Tue, 24 Mar 2026 07:58:48 -0600
Subject: [PATCH 11/38] move initialization module from contrib to devel

---
 pyomo/{contrib => devel}/initialization/README.md                 | 0
 pyomo/{contrib => devel}/initialization/__init__.py               | 0
 pyomo/{contrib => devel}/initialization/bounds/__init__.py        | 0
 pyomo/{contrib => devel}/initialization/bounds/bound_variables.py | 0
 pyomo/{contrib => devel}/initialization/global_init.py            | 0
 pyomo/{contrib => devel}/initialization/initialize.py             | 0
 pyomo/{contrib => devel}/initialization/lp_approx_init.py         | 0
 pyomo/{contrib => devel}/initialization/pwl_init.py               | 0
 pyomo/{contrib => devel}/initialization/utils.py                  | 0
 9 files changed, 0 insertions(+), 0 deletions(-)
 rename pyomo/{contrib => devel}/initialization/README.md (100%)
 rename pyomo/{contrib => devel}/initialization/__init__.py (100%)
 rename pyomo/{contrib => devel}/initialization/bounds/__init__.py (100%)
 rename pyomo/{contrib => devel}/initialization/bounds/bound_variables.py (100%)
 rename pyomo/{contrib => devel}/initialization/global_init.py (100%)
 rename pyomo/{contrib => devel}/initialization/initialize.py (100%)
 rename pyomo/{contrib => devel}/initialization/lp_approx_init.py (100%)
 rename pyomo/{contrib => devel}/initialization/pwl_init.py (100%)
 rename pyomo/{contrib => devel}/initialization/utils.py (100%)

diff --git a/pyomo/contrib/initialization/README.md b/pyomo/devel/initialization/README.md
similarity index 100%
rename from pyomo/contrib/initialization/README.md
rename to pyomo/devel/initialization/README.md
diff --git a/pyomo/contrib/initialization/__init__.py b/pyomo/devel/initialization/__init__.py
similarity index 100%
rename from pyomo/contrib/initialization/__init__.py
rename to pyomo/devel/initialization/__init__.py
diff --git a/pyomo/contrib/initialization/bounds/__init__.py b/pyomo/devel/initialization/bounds/__init__.py
similarity index 100%
rename from pyomo/contrib/initialization/bounds/__init__.py
rename to pyomo/devel/initialization/bounds/__init__.py
diff --git a/pyomo/contrib/initialization/bounds/bound_variables.py b/pyomo/devel/initialization/bounds/bound_variables.py
similarity index 100%
rename from pyomo/contrib/initialization/bounds/bound_variables.py
rename to pyomo/devel/initialization/bounds/bound_variables.py
diff --git a/pyomo/contrib/initialization/global_init.py b/pyomo/devel/initialization/global_init.py
similarity index 100%
rename from pyomo/contrib/initialization/global_init.py
rename to pyomo/devel/initialization/global_init.py
diff --git a/pyomo/contrib/initialization/initialize.py b/pyomo/devel/initialization/initialize.py
similarity index 100%
rename from pyomo/contrib/initialization/initialize.py
rename to pyomo/devel/initialization/initialize.py
diff --git a/pyomo/contrib/initialization/lp_approx_init.py b/pyomo/devel/initialization/lp_approx_init.py
similarity index 100%
rename from pyomo/contrib/initialization/lp_approx_init.py
rename to pyomo/devel/initialization/lp_approx_init.py
diff --git a/pyomo/contrib/initialization/pwl_init.py b/pyomo/devel/initialization/pwl_init.py
similarity index 100%
rename from pyomo/contrib/initialization/pwl_init.py
rename to pyomo/devel/initialization/pwl_init.py
diff --git a/pyomo/contrib/initialization/utils.py b/pyomo/devel/initialization/utils.py
similarity index 100%
rename from pyomo/contrib/initialization/utils.py
rename to pyomo/devel/initialization/utils.py

From 24dbd64bc36b26664ab3e4641a1bfc46ac288adf Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Tue, 24 Mar 2026 08:41:07 -0600
Subject: [PATCH 12/38] update imports

---
 pyomo/devel/initialization/__init__.py               |  1 +
 pyomo/devel/initialization/bounds/bound_variables.py |  2 +-
 pyomo/devel/initialization/global_init.py            |  3 ++-
 pyomo/devel/initialization/initialize.py             |  8 ++++----
 pyomo/devel/initialization/lp_approx_init.py         | 10 ++++++----
 pyomo/devel/initialization/pwl_init.py               |  7 ++++---
 6 files changed, 18 insertions(+), 13 deletions(-)

diff --git a/pyomo/devel/initialization/__init__.py b/pyomo/devel/initialization/__init__.py
index e69de29bb2d..099f15af824 100644
--- a/pyomo/devel/initialization/__init__.py
+++ b/pyomo/devel/initialization/__init__.py
@@ -0,0 +1 @@
+from pyomo.devel.initialization.initialize import initialize_nlp, InitializationMethod
\ No newline at end of file
diff --git a/pyomo/devel/initialization/bounds/bound_variables.py b/pyomo/devel/initialization/bounds/bound_variables.py
index 0d0232f53b8..8741b6ca979 100644
--- a/pyomo/devel/initialization/bounds/bound_variables.py
+++ b/pyomo/devel/initialization/bounds/bound_variables.py
@@ -1,6 +1,6 @@
 from pyomo.core.base.block import BlockData
 from pyomo.contrib.fbbt.fbbt import fbbt
-from pyomo.contrib.initialization.utils import get_vars
+from pyomo.devel.initialization.utils import get_vars
 import logging
 
 
diff --git a/pyomo/devel/initialization/global_init.py b/pyomo/devel/initialization/global_init.py
index 4f858a4b7c3..142d58ec802 100644
--- a/pyomo/devel/initialization/global_init.py
+++ b/pyomo/devel/initialization/global_init.py
@@ -1,5 +1,6 @@
 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
@@ -23,7 +24,7 @@ def _initialize_with_global_solver(
     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.incumbent_objective is not None:
+    if res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
         res.solution_loader.load_vars()
     else:
         raise RuntimeError('no feasible solution found')
diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py
index f8918f95528..f62bc9610b1 100644
--- a/pyomo/devel/initialization/initialize.py
+++ b/pyomo/devel/initialization/initialize.py
@@ -1,12 +1,12 @@
 from typing import Optional
 from pyomo.core.base.block import BlockData
 from enum import Enum
-from pyomo.contrib.initialization.utils import get_vars, shallow_clone
+from pyomo.devel.initialization.utils import get_vars, shallow_clone
 from pyomo.common.collections import ComponentMap
-from pyomo.contrib.initialization.pwl_init import _initialize_with_piecewise_linear_approximation
-from pyomo.contrib.initialization.lp_approx_init import _initialize_with_LP_approximation
+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.contrib.initialization.global_init import _initialize_with_global_solver
+from pyomo.devel.initialization.global_init import _initialize_with_global_solver
 from pyomo.contrib.solver.common.factory import SolverFactory
 import logging
 
diff --git a/pyomo/devel/initialization/lp_approx_init.py b/pyomo/devel/initialization/lp_approx_init.py
index 0c72e82fbd7..18752271c47 100644
--- a/pyomo/devel/initialization/lp_approx_init.py
+++ b/pyomo/devel/initialization/lp_approx_init.py
@@ -1,10 +1,11 @@
 from pyomo.core.base.block import BlockData
 import pyomo.environ as pe
-from pyomo.contrib.initialization.bounds.bound_variables import bound_all_nonlinear_variables
-from pyomo.contrib.initialization.utils import fix_vars_with_equal_bounds, shallow_clone, get_vars
+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
@@ -36,7 +37,7 @@
 from pyomo.contrib.solver.common.base import SolverBase
 import logging
 from pyomo.common.modeling import unique_component_name
-from pyomo.contrib.initialization.pwl_init import _minimize_infeasibility
+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
@@ -175,7 +176,8 @@ def _initialize_with_LP_approximation(
     # 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.incumbent_objective is not None:
+
+    if nlp_res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
         nlp_res.solution_loader.load_vars()
     else:
         raise RuntimeError('no feasible solution found')
diff --git a/pyomo/devel/initialization/pwl_init.py b/pyomo/devel/initialization/pwl_init.py
index 30baa1d7cef..d11d518608a 100644
--- a/pyomo/devel/initialization/pwl_init.py
+++ b/pyomo/devel/initialization/pwl_init.py
@@ -1,7 +1,7 @@
 from pyomo.core.base.block import BlockData
 import pyomo.environ as pe
-from pyomo.contrib.initialization.bounds.bound_variables import bound_all_nonlinear_variables
-from pyomo.contrib.initialization.utils import fix_vars_with_equal_bounds, shallow_clone, get_vars
+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
@@ -36,6 +36,7 @@
 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__)
@@ -325,7 +326,7 @@ def _initialize_with_piecewise_linear_approximation(
         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.incumbent_objective is not None:
+        if res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
             solved = True
             res.solution_loader.load_vars()
             break

From 261650220c608f6581ed500eff0a7d9f8f1bcbcd Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Tue, 24 Mar 2026 08:41:58 -0600
Subject: [PATCH 13/38] initialization example

---
 .../solver/solvers/scip/scip_direct.py        |  2 +-
 .../devel/initialization/examples/__init__.py |  0
 .../examples/init_polynomial_ex.py            | 33 +++++++++++++++++++
 3 files changed, 34 insertions(+), 1 deletion(-)
 create mode 100644 pyomo/devel/initialization/examples/__init__.py
 create mode 100644 pyomo/devel/initialization/examples/init_polynomial_ex.py

diff --git a/pyomo/contrib/solver/solvers/scip/scip_direct.py b/pyomo/contrib/solver/solvers/scip/scip_direct.py
index 022878f918d..1371ef61530 100644
--- a/pyomo/contrib/solver/solvers/scip/scip_direct.py
+++ b/pyomo/contrib/solver/solvers/scip/scip_direct.py
@@ -349,7 +349,7 @@ def load_vars(
         for v, val in self.get_vars(
             vars_to_load=vars_to_load, solution_id=solution_id
         ).items():
-            v.value = val
+            v.set_value(val, skip_validation=True)
 
     def get_vars(
         self, vars_to_load: Optional[Sequence[VarData]] = None, solution_id=None
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..fe99923a23c
--- /dev/null
+++ b/pyomo/devel/initialization/examples/init_polynomial_ex.py
@@ -0,0 +1,33 @@
+import pyomo.environ as pyo
+import pyomo.devel.initialization as ini
+from pyomo.contrib.solver.common.factory import SolverFactory
+import logging
+
+
+def build_model():
+    m = pyo.ConcreteModel()
+    m.x = pyo.Var()
+    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))

From 54d354e60383c66f4d385649160c5903afaa1087 Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Tue, 24 Mar 2026 12:30:46 -0600
Subject: [PATCH 14/38] initialization: starting tests

---
 .../examples/init_polynomial_ex.py            |  3 +-
 pyomo/devel/initialization/tests/__init__.py  |  0
 .../tests/test_initialization.py              | 29 +++++++++++++++++++
 3 files changed, 30 insertions(+), 2 deletions(-)
 create mode 100644 pyomo/devel/initialization/tests/__init__.py
 create mode 100644 pyomo/devel/initialization/tests/test_initialization.py

diff --git a/pyomo/devel/initialization/examples/init_polynomial_ex.py b/pyomo/devel/initialization/examples/init_polynomial_ex.py
index fe99923a23c..dc8dc88b31a 100644
--- a/pyomo/devel/initialization/examples/init_polynomial_ex.py
+++ b/pyomo/devel/initialization/examples/init_polynomial_ex.py
@@ -1,12 +1,11 @@
 import pyomo.environ as pyo
 import pyomo.devel.initialization as ini
 from pyomo.contrib.solver.common.factory import SolverFactory
-import logging
 
 
 def build_model():
     m = pyo.ConcreteModel()
-    m.x = pyo.Var()
+    m.x = pyo.Var(bounds=(-20, 20))
     m.c = pyo.Constraint(expr=(m.x+7)*(m.x+5)*(m.x-4) + 200 == 0)
     return m
 
diff --git a/pyomo/devel/initialization/tests/__init__.py b/pyomo/devel/initialization/tests/__init__.py
new file mode 100644
index 00000000000..e69de29bb2d
diff --git a/pyomo/devel/initialization/tests/test_initialization.py b/pyomo/devel/initialization/tests/test_initialization.py
new file mode 100644
index 00000000000..c1c6dd9bbd2
--- /dev/null
+++ b/pyomo/devel/initialization/tests/test_initialization.py
@@ -0,0 +1,29 @@
+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
+
+
+scip = SolverFactory('scip_direct')
+ipopt = SolverFactory('ipopt')
+
+
+@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)

From 056cc993a4664f9db22ef644c250708451ad3dfd Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Tue, 24 Mar 2026 14:22:09 -0600
Subject: [PATCH 15/38] initialization: add copyright statements

---
 pyomo/devel/initialization/__init__.py                 |  9 +++++++++
 pyomo/devel/initialization/bounds/__init__.py          |  8 ++++++++
 pyomo/devel/initialization/bounds/bound_variables.py   |  9 +++++++++
 .../initialization/examples/init_polynomial_ex.py      | 10 ++++++++++
 pyomo/devel/initialization/global_init.py              |  9 +++++++++
 pyomo/devel/initialization/initialize.py               |  9 +++++++++
 pyomo/devel/initialization/lp_approx_init.py           |  9 +++++++++
 pyomo/devel/initialization/pwl_init.py                 |  9 +++++++++
 pyomo/devel/initialization/tests/__init__.py           |  8 ++++++++
 .../devel/initialization/tests/test_initialization.py  |  9 +++++++++
 pyomo/devel/initialization/utils.py                    |  9 +++++++++
 11 files changed, 98 insertions(+)

diff --git a/pyomo/devel/initialization/__init__.py b/pyomo/devel/initialization/__init__.py
index 099f15af824..42e52b2d8ff 100644
--- a/pyomo/devel/initialization/__init__.py
+++ b/pyomo/devel/initialization/__init__.py
@@ -1 +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
\ No newline at end of file
diff --git a/pyomo/devel/initialization/bounds/__init__.py b/pyomo/devel/initialization/bounds/__init__.py
index e69de29bb2d..231b44987f6 100644
--- a/pyomo/devel/initialization/bounds/__init__.py
+++ 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
index 8741b6ca979..4402ea76510 100644
--- a/pyomo/devel/initialization/bounds/bound_variables.py
+++ b/pyomo/devel/initialization/bounds/bound_variables.py
@@ -1,3 +1,12 @@
+# ____________________________________________________________________________________
+#
+# 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
diff --git a/pyomo/devel/initialization/examples/init_polynomial_ex.py b/pyomo/devel/initialization/examples/init_polynomial_ex.py
index dc8dc88b31a..88b03aed0e0 100644
--- a/pyomo/devel/initialization/examples/init_polynomial_ex.py
+++ b/pyomo/devel/initialization/examples/init_polynomial_ex.py
@@ -1,3 +1,13 @@
+# ____________________________________________________________________________________
+#
+# 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
diff --git a/pyomo/devel/initialization/global_init.py b/pyomo/devel/initialization/global_init.py
index 142d58ec802..11533d9f516 100644
--- a/pyomo/devel/initialization/global_init.py
+++ b/pyomo/devel/initialization/global_init.py
@@ -1,3 +1,12 @@
+# ____________________________________________________________________________________
+#
+# 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
diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py
index f62bc9610b1..90eab0f3dad 100644
--- a/pyomo/devel/initialization/initialize.py
+++ b/pyomo/devel/initialization/initialize.py
@@ -2,6 +2,15 @@
 from pyomo.core.base.block import BlockData
 from enum import Enum
 from pyomo.devel.initialization.utils import get_vars, shallow_clone
+# ____________________________________________________________________________________
+#
+# 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.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
diff --git a/pyomo/devel/initialization/lp_approx_init.py b/pyomo/devel/initialization/lp_approx_init.py
index 18752271c47..4401a69e55d 100644
--- a/pyomo/devel/initialization/lp_approx_init.py
+++ b/pyomo/devel/initialization/lp_approx_init.py
@@ -1,3 +1,12 @@
+# ____________________________________________________________________________________
+#
+# 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
diff --git a/pyomo/devel/initialization/pwl_init.py b/pyomo/devel/initialization/pwl_init.py
index d11d518608a..bec34ddb612 100644
--- a/pyomo/devel/initialization/pwl_init.py
+++ b/pyomo/devel/initialization/pwl_init.py
@@ -1,3 +1,12 @@
+# ____________________________________________________________________________________
+#
+# 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
diff --git a/pyomo/devel/initialization/tests/__init__.py b/pyomo/devel/initialization/tests/__init__.py
index e69de29bb2d..231b44987f6 100644
--- a/pyomo/devel/initialization/tests/__init__.py
+++ 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
index c1c6dd9bbd2..9e7aa50753e 100644
--- a/pyomo/devel/initialization/tests/test_initialization.py
+++ b/pyomo/devel/initialization/tests/test_initialization.py
@@ -1,3 +1,12 @@
+# ____________________________________________________________________________________
+#
+# 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 pyo
 import pyomo.devel.initialization as ini
 from pyomo.devel.initialization.examples.init_polynomial_ex import main
diff --git a/pyomo/devel/initialization/utils.py b/pyomo/devel/initialization/utils.py
index 792509583f0..3c69ba1a55f 100644
--- a/pyomo/devel/initialization/utils.py
+++ b/pyomo/devel/initialization/utils.py
@@ -1,3 +1,12 @@
+# ____________________________________________________________________________________
+#
+# 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 d70e677a751af53a9a19025025aa46e44f2c0ab7 Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Tue, 24 Mar 2026 14:43:53 -0600
Subject: [PATCH 16/38] docstring

---
 pyomo/devel/initialization/initialize.py | 54 ++++++++++++++++++++++--
 1 file changed, 50 insertions(+), 4 deletions(-)

diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py
index 90eab0f3dad..c0cbcf35e68 100644
--- a/pyomo/devel/initialization/initialize.py
+++ b/pyomo/devel/initialization/initialize.py
@@ -17,6 +17,7 @@
 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
 
 
@@ -44,10 +45,55 @@ def initialize_nlp(
     nlp_solver: Optional[SolverBase] = None,
     global_solver: Optional[SolverBase] = None,
     method: InitializationMethod = InitializationMethod.global_opt,
-    default_bound=1.0e8,
-    max_pwl_refinement_iter=100,
-    num_pwl_cons_to_refine_per_iter=5,
-):
+    default_bound: float = 1.0e8,
+    max_pwl_refinement_iter: int = 100,
+    num_pwl_cons_to_refine_per_iter: int = 5,
+) -> 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.
+
+    Returns
+    -------
+    res: Results
+        The results object obtained the last time the nlp_solver was used to 
+        try and solve the model.
+    """
     # get all variable bounds, domains, etc. to restore them later
     orig_vars = get_vars(nlp)
     orig_var_data = ComponentMap(

From 823641d0dcf3207cb0a9c254d3b63f8acb07b336 Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Tue, 24 Mar 2026 20:22:49 -0600
Subject: [PATCH 17/38] initialization: docs

---
 doc/OnlineDocs/explanation/analysis/index.rst          |  1 +
 .../analysis/nlp_initialization/nlp_initialization.rst | 10 ++++++++++
 2 files changed, 11 insertions(+)
 create mode 100644 doc/OnlineDocs/explanation/analysis/nlp_initialization/nlp_initialization.rst

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..0e7c07664ac
--- /dev/null
+++ b/doc/OnlineDocs/explanation/analysis/nlp_initialization/nlp_initialization.rst
@@ -0,0 +1,10 @@
+.. _pyomo.devel.initialization
+
+NLP Initialization
+==================
+The initialization module within ``pyomo.devel.initialization`` is intended to provide methods to help initialize nonconvex nonlinear programs. Example usage is shown below.
+
+.. literalinclude:: /../../pyomo/devel/initialization/examples/ini_polynomial_ex.py
+    :start-after: # === Required imports ===
+
+The `:meth:~pyomo.devel.initialization.initialize.initialize_nlp` function 

From 0b02132ec10aec2739e11f6209092a8c116a4f6b Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Tue, 7 Apr 2026 07:26:55 -0600
Subject: [PATCH 18/38] working on docs for initialization

---
 .../nlp_initialization/nlp_initialization.rst | 131 +++++++++++++++++-
 doc/OnlineDocs/explanation/index.rst          |   1 +
 2 files changed, 127 insertions(+), 5 deletions(-)

diff --git a/doc/OnlineDocs/explanation/analysis/nlp_initialization/nlp_initialization.rst b/doc/OnlineDocs/explanation/analysis/nlp_initialization/nlp_initialization.rst
index 0e7c07664ac..c18f0f4ac0a 100644
--- a/doc/OnlineDocs/explanation/analysis/nlp_initialization/nlp_initialization.rst
+++ b/doc/OnlineDocs/explanation/analysis/nlp_initialization/nlp_initialization.rst
@@ -1,10 +1,131 @@
-.. _pyomo.devel.initialization
+.. _analysis_nlp_initialization:
 
 NLP Initialization
-==================
-The initialization module within ``pyomo.devel.initialization`` is intended to provide methods to help initialize nonconvex nonlinear programs. Example usage is shown below.
+******************
 
-.. literalinclude:: /../../pyomo/devel/initialization/examples/ini_polynomial_ex.py
+.. 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 `:meth:~pyomo.devel.initialization.initialize.initialize_nlp` function 
+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 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 <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 solutin. 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 implemenation 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`

From 235587037dd6328d551a2afdbfbf51302cc3ced8 Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Thu, 9 Apr 2026 09:47:34 -0600
Subject: [PATCH 19/38] bug

---
 pyomo/devel/initialization/utils.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyomo/devel/initialization/utils.py b/pyomo/devel/initialization/utils.py
index 3c69ba1a55f..b73da8625b6 100644
--- a/pyomo/devel/initialization/utils.py
+++ b/pyomo/devel/initialization/utils.py
@@ -19,7 +19,7 @@ def get_vars(m: BlockData):
     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(c.expr, include_fixed=False))
+        vset.update(identify_variables(o.expr, include_fixed=False))
     return vset
 
 

From e40f487286c1e70fb36d0d7e8b6e1de80b411148 Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Thu, 9 Apr 2026 10:26:38 -0600
Subject: [PATCH 20/38] try nlp before initialization

---
 pyomo/devel/initialization/initialize.py | 12 +++++++++++-
 1 file changed, 11 insertions(+), 1 deletion(-)

diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py
index c0cbcf35e68..2793cda79ca 100644
--- a/pyomo/devel/initialization/initialize.py
+++ b/pyomo/devel/initialization/initialize.py
@@ -19,6 +19,7 @@
 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__)
@@ -41,8 +42,8 @@ def _get_solver(sname, reason):
 
 def initialize_nlp(
     nlp: BlockData, 
+    nlp_solver: SolverBase,
     mip_solver: Optional[SolverBase] = None,
-    nlp_solver: Optional[SolverBase] = None,
     global_solver: Optional[SolverBase] = None,
     method: InitializationMethod = InitializationMethod.global_opt,
     default_bound: float = 1.0e8,
@@ -94,6 +95,15 @@ def initialize_nlp(
         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(

From 42d7fff2866e8e0effc8d6184f5922a5f76cea75 Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Thu, 9 Apr 2026 16:11:09 -0600
Subject: [PATCH 21/38] fix some docstrings

---
 .../analysis/nlp_initialization/nlp_initialization.rst |  2 +-
 pyomo/devel/initialization/initialize.py               | 10 +++++-----
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/doc/OnlineDocs/explanation/analysis/nlp_initialization/nlp_initialization.rst b/doc/OnlineDocs/explanation/analysis/nlp_initialization/nlp_initialization.rst
index c18f0f4ac0a..770f14369bf 100644
--- a/doc/OnlineDocs/explanation/analysis/nlp_initialization/nlp_initialization.rst
+++ b/doc/OnlineDocs/explanation/analysis/nlp_initialization/nlp_initialization.rst
@@ -84,7 +84,7 @@ Method ``pwl_approximation``
 ----------------------------
 
 This method builds a piecewise linear (PWL) approximation of the model, solves
-it, and initializes the NLP solver at the solutin. If the NLP solver does not 
+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. 
diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py
index 2793cda79ca..e240e46a12c 100644
--- a/pyomo/devel/initialization/initialize.py
+++ b/pyomo/devel/initialization/initialize.py
@@ -1,7 +1,3 @@
-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
 # ____________________________________________________________________________________
 #
 # Pyomo: Python Optimization Modeling Objects
@@ -11,6 +7,10 @@
 # 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
@@ -91,7 +91,7 @@ def initialize_nlp(
 
     Returns
     -------
-    res: Results
+    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.
     """

From fe98b7c72e86f32fcc2d4fc28c6eb55f08b7669e Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Sat, 11 Apr 2026 16:29:21 -0600
Subject: [PATCH 22/38] initialization: working on tests

---
 .../initialization/examples/init_polynomial_ex.py      |  2 +-
 pyomo/devel/initialization/global_init.py              |  2 +-
 pyomo/devel/initialization/lp_approx_init.py           | 10 ++++++----
 3 files changed, 8 insertions(+), 6 deletions(-)

diff --git a/pyomo/devel/initialization/examples/init_polynomial_ex.py b/pyomo/devel/initialization/examples/init_polynomial_ex.py
index 88b03aed0e0..111ad418363 100644
--- a/pyomo/devel/initialization/examples/init_polynomial_ex.py
+++ b/pyomo/devel/initialization/examples/init_polynomial_ex.py
@@ -15,7 +15,7 @@
 
 def build_model():
     m = pyo.ConcreteModel()
-    m.x = pyo.Var(bounds=(-20, 20))
+    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
 
diff --git a/pyomo/devel/initialization/global_init.py b/pyomo/devel/initialization/global_init.py
index 11533d9f516..b924fb61273 100644
--- a/pyomo/devel/initialization/global_init.py
+++ b/pyomo/devel/initialization/global_init.py
@@ -36,6 +36,6 @@ def _initialize_with_global_solver(
     if res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
         res.solution_loader.load_vars()
     else:
-        raise RuntimeError('no feasible solution found')
+        logger.warning('initialization was not successful via global optimization')
     
     return res
diff --git a/pyomo/devel/initialization/lp_approx_init.py b/pyomo/devel/initialization/lp_approx_init.py
index 4401a69e55d..91537436e10 100644
--- a/pyomo/devel/initialization/lp_approx_init.py
+++ b/pyomo/devel/initialization/lp_approx_init.py
@@ -156,9 +156,11 @@ def _initialize_with_LP_approximation(
 
     # 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(nlp, aggressive_substitution=False)
-    logger.info('applied the univariate_nonlinear_decomposition transformation')
+    # 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)
@@ -189,6 +191,6 @@ def _initialize_with_LP_approximation(
     if nlp_res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
         nlp_res.solution_loader.load_vars()
     else:
-        raise RuntimeError('no feasible solution found')
+        logger.warning('initialization was not successful via LP approximation')
 
     return nlp_res

From 5d9660b08f09d59c9f41a80951f4b5bee5929651 Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Sat, 11 Apr 2026 16:31:19 -0600
Subject: [PATCH 23/38] initialization: fixing bugs

---
 pyomo/devel/initialization/pwl_init.py | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/pyomo/devel/initialization/pwl_init.py b/pyomo/devel/initialization/pwl_init.py
index bec34ddb612..12f493f293b 100644
--- a/pyomo/devel/initialization/pwl_init.py
+++ b/pyomo/devel/initialization/pwl_init.py
@@ -63,7 +63,7 @@ def _minimize_infeasibility(m):
         if obj.sense == pe.minimize:
             obj_expr += 0.1*obj.expr
         else:
-            obj_expr -= 0.1*obj_expr
+            obj_expr -= 0.1*obj.expr
         obj.deactivate()
         found_obj = True
 
@@ -102,15 +102,17 @@ def _get_pwl_constraints(m: BlockData) -> MutableMapping[
     comp_types = set()
     comp_types.add(PiecewiseLinearExpression)
     pwl_expr_to_con_map = ComponentMap()
-    for con in m.component_data_objects(pe.Constraint, active=True, descend_into=True):
-        pwl_exprs = list(identify_components(con.expr, comp_types))
+    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(con)
+        pwl_expr_to_con_map[e].append(comp)
     return pwl_expr_to_con_map
 
 
@@ -224,6 +226,7 @@ def _refine_pwl_approx(
         # 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)
@@ -341,6 +344,6 @@ def _initialize_with_piecewise_linear_approximation(
             break
 
     if not solved:
-        raise RuntimeError('no feasible solution found')
+        logger.warning('initialization was not successful via PWL approximation')
 
     return last_nlp_res

From c08a92c4cfe40b9723b92120936f2849e6b24383 Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Sat, 11 Apr 2026 17:01:30 -0600
Subject: [PATCH 24/38] fix nudge

---
 pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py | 6 +++---
 pyomo/devel/initialization/pwl_init.py                | 7 +++++--
 2 files changed, 8 insertions(+), 5 deletions(-)

diff --git a/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py b/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py
index db3f4b21621..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)
 
 
diff --git a/pyomo/devel/initialization/pwl_init.py b/pyomo/devel/initialization/pwl_init.py
index 12f493f293b..786bcc2c699 100644
--- a/pyomo/devel/initialization/pwl_init.py
+++ b/pyomo/devel/initialization/pwl_init.py
@@ -237,8 +237,11 @@ def _refine_pwl_approx(
     if len(violations) == 0:
         raise RuntimeError('Did not find any piecewise linear functions with variable values')
     
-    if math.isclose(violations[0][0], 0):
-        raise RuntimeError('All of the original nonlinear functions are satisfied!')
+    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}')

From ec9f9b9496009fa013da3767da2d612655b0883d Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Sun, 12 Apr 2026 19:41:04 -0600
Subject: [PATCH 25/38] initialization: working on tests

---
 .../tests/test_initialization.py              | 154 +++++++++++++++++-
 1 file changed, 153 insertions(+), 1 deletion(-)

diff --git a/pyomo/devel/initialization/tests/test_initialization.py b/pyomo/devel/initialization/tests/test_initialization.py
index 9e7aa50753e..177f32cec76 100644
--- a/pyomo/devel/initialization/tests/test_initialization.py
+++ b/pyomo/devel/initialization/tests/test_initialization.py
@@ -7,18 +7,54 @@
 # 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
+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):
+            # print(v, v.value, val)
+            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):
@@ -36,3 +72,119 @@ 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.x1 = pyo.Var(bounds=(0.5, 1.5))
+        m.x2 = pyo.Var(bounds=(2.5, 3.5))
+        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-0, 1e-0),
+                2: ([1, 3], 1e-1, 1e-1),
+                3: ([1, 3], 1e-2, 1e-2),
+                4: ([1, 3], 1e-3, 1e-3),
+                5: ([1, 3], 1e-5, 1e-5),
+                6: ([1, 3], 1e-6, 1e-6),
+                7: ([1, 3], 1e-7, 1e-7),
+                8: ([1, 3], 1e-7, 1e-7),
+                9: ([1, 3], 1e-7, 1e-7),
+                10: ([1, 3], 1e-7, 1e-7),
+                11: ([1, 3], 1e-7, 1e-7),
+            },
+        )
+        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=10,
+        )
+
+
+if __name__ == '__main__':
+    import logging
+    logging.basicConfig(level=logging.INFO)
+    t = TestInit()
+    t.test_pwl_init()

From e52e1236921c09f0fadd84f78a2447740678e3c5 Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Sun, 12 Apr 2026 19:42:00 -0600
Subject: [PATCH 26/38] debugging pwl initialization

---
 pyomo/devel/initialization/pwl_init.py | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/pyomo/devel/initialization/pwl_init.py b/pyomo/devel/initialization/pwl_init.py
index 786bcc2c699..7e32b8b31e6 100644
--- a/pyomo/devel/initialization/pwl_init.py
+++ b/pyomo/devel/initialization/pwl_init.py
@@ -292,6 +292,10 @@ def _initialize_with_piecewise_linear_approximation(
     _minimize_infeasibility(pwl)
     logger.info('reformulated model to minimize infeasibility')
 
+    m_before_pwl = shallow_clone(pwl)
+    vset = list(ComponentSet(get_vars(m_before_pwl)))
+    vset.sort(key=lambda i: i.local_name)
+
     # build the PWL approximation
     trans = pe.TransformationFactory('contrib.piecewise.nonlinear_to_pwl')
     trans.apply_to(pwl, num_points=2, additively_decompose=False)
@@ -329,6 +333,16 @@ def _initialize_with_piecewise_linear_approximation(
         for ov, nv in zip(orig_vars, new_vars):
             ov.set_value(nv.value, skip_validation=True)
 
+        print('\n\n**********************')
+        print(m_before_pwl.obj.expr)
+        print('\nConstraints:')
+        for c in m_before_pwl.component_data_objects(pe.Constraint, active=True, descend_into=True):
+            print(c)
+            print('    ', c.expr)
+        for v in vset:
+            print(v, v.value)
+        m_before_pwl.display()
+
         # refine the PWL approximation
         _refine_pwl_approx(
             pwl,

From 5b4fe02706a685c537290a77e818b6ef0e3aaec4 Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Tue, 14 Apr 2026 08:19:16 -0600
Subject: [PATCH 27/38] comment out print statements

---
 pyomo/devel/initialization/pwl_init.py | 24 ++++++++++++------------
 1 file changed, 12 insertions(+), 12 deletions(-)

diff --git a/pyomo/devel/initialization/pwl_init.py b/pyomo/devel/initialization/pwl_init.py
index 7e32b8b31e6..3d2a37251ed 100644
--- a/pyomo/devel/initialization/pwl_init.py
+++ b/pyomo/devel/initialization/pwl_init.py
@@ -292,9 +292,9 @@ def _initialize_with_piecewise_linear_approximation(
     _minimize_infeasibility(pwl)
     logger.info('reformulated model to minimize infeasibility')
 
-    m_before_pwl = shallow_clone(pwl)
-    vset = list(ComponentSet(get_vars(m_before_pwl)))
-    vset.sort(key=lambda i: i.local_name)
+    # m_before_pwl = shallow_clone(pwl)
+    # vset = list(ComponentSet(get_vars(m_before_pwl)))
+    # vset.sort(key=lambda i: i.local_name)
 
     # build the PWL approximation
     trans = pe.TransformationFactory('contrib.piecewise.nonlinear_to_pwl')
@@ -333,15 +333,15 @@ def _initialize_with_piecewise_linear_approximation(
         for ov, nv in zip(orig_vars, new_vars):
             ov.set_value(nv.value, skip_validation=True)
 
-        print('\n\n**********************')
-        print(m_before_pwl.obj.expr)
-        print('\nConstraints:')
-        for c in m_before_pwl.component_data_objects(pe.Constraint, active=True, descend_into=True):
-            print(c)
-            print('    ', c.expr)
-        for v in vset:
-            print(v, v.value)
-        m_before_pwl.display()
+        # print('\n\n**********************')
+        # print(m_before_pwl.obj.expr)
+        # print('\nConstraints:')
+        # for c in m_before_pwl.component_data_objects(pe.Constraint, active=True, descend_into=True):
+        #     print(c)
+        #     print('    ', c.expr)
+        # for v in vset:
+        #     print(v, v.value)
+        # m_before_pwl.display()
 
         # refine the PWL approximation
         _refine_pwl_approx(

From 135b1527d70762aee4ccd0f6c2422c0eea0520fd Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Wed, 15 Apr 2026 05:54:18 -0600
Subject: [PATCH 28/38] initialization: testing

---
 pyomo/devel/initialization/initialize.py      |  2 +
 pyomo/devel/initialization/pwl_init.py        |  3 +-
 .../tests/test_initialization.py              | 54 ++++++++++++-------
 3 files changed, 38 insertions(+), 21 deletions(-)

diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py
index e240e46a12c..64439212b1e 100644
--- a/pyomo/devel/initialization/initialize.py
+++ b/pyomo/devel/initialization/initialize.py
@@ -49,6 +49,7 @@ def initialize_nlp(
     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``.
@@ -123,6 +124,7 @@ def initialize_nlp(
             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:
diff --git a/pyomo/devel/initialization/pwl_init.py b/pyomo/devel/initialization/pwl_init.py
index 3d2a37251ed..d3b6c3c25c6 100644
--- a/pyomo/devel/initialization/pwl_init.py
+++ b/pyomo/devel/initialization/pwl_init.py
@@ -265,6 +265,7 @@ def _initialize_with_piecewise_linear_approximation(
     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)
@@ -273,7 +274,7 @@ def _initialize_with_piecewise_linear_approximation(
     # 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=True)
+    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
diff --git a/pyomo/devel/initialization/tests/test_initialization.py b/pyomo/devel/initialization/tests/test_initialization.py
index 177f32cec76..0ba33881570 100644
--- a/pyomo/devel/initialization/tests/test_initialization.py
+++ b/pyomo/devel/initialization/tests/test_initialization.py
@@ -40,7 +40,6 @@ 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):
-            # print(v, v.value, val)
             assert v.value == pytest.approx(val, rel=rel_tol, abs=abs_tol)
     
     def solve(self, model, **kwds) -> Results:
@@ -149,28 +148,42 @@ def test_pwl_init(self):
         approximation improves as refinement is done.
         """
         m = pyo.ConcreteModel()
-        m.x1 = pyo.Var(bounds=(0.5, 1.5))
-        m.x2 = pyo.Var(bounds=(2.5, 3.5))
-        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)
+        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.x1, m.x2], 
+            varlist=[m.x], 
             sol_map={
-                0: ([None, None], 0, 0),
-                1: ([1, 3], 1e-0, 1e-0),
-                2: ([1, 3], 1e-1, 1e-1),
-                3: ([1, 3], 1e-2, 1e-2),
-                4: ([1, 3], 1e-3, 1e-3),
-                5: ([1, 3], 1e-5, 1e-5),
-                6: ([1, 3], 1e-6, 1e-6),
-                7: ([1, 3], 1e-7, 1e-7),
-                8: ([1, 3], 1e-7, 1e-7),
-                9: ([1, 3], 1e-7, 1e-7),
-                10: ([1, 3], 1e-7, 1e-7),
-                11: ([1, 3], 1e-7, 1e-7),
+                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')
@@ -179,7 +192,8 @@ def test_pwl_init(self):
             nlp_solver=nlp_solver,
             mip_solver=mip_solver,
             method=ini.InitializationMethod.pwl_approximation,
-            max_pwl_refinement_iter=10,
+            max_pwl_refinement_iter=27,
+            aggressive_substitution=False,
         )
 
 

From fc60afeb75d7f1e61a61b6b64c740a60ff6c0d95 Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Wed, 15 Apr 2026 05:57:13 -0600
Subject: [PATCH 29/38] initialization: update docstring; remove debugging code

---
 pyomo/devel/initialization/initialize.py |  4 ++++
 pyomo/devel/initialization/pwl_init.py   | 14 --------------
 2 files changed, 4 insertions(+), 14 deletions(-)

diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py
index 64439212b1e..1369db1a1f1 100644
--- a/pyomo/devel/initialization/initialize.py
+++ b/pyomo/devel/initialization/initialize.py
@@ -89,6 +89,10 @@ def initialize_nlp(
         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
     -------
diff --git a/pyomo/devel/initialization/pwl_init.py b/pyomo/devel/initialization/pwl_init.py
index d3b6c3c25c6..551e33a2f9a 100644
--- a/pyomo/devel/initialization/pwl_init.py
+++ b/pyomo/devel/initialization/pwl_init.py
@@ -293,10 +293,6 @@ def _initialize_with_piecewise_linear_approximation(
     _minimize_infeasibility(pwl)
     logger.info('reformulated model to minimize infeasibility')
 
-    # m_before_pwl = shallow_clone(pwl)
-    # vset = list(ComponentSet(get_vars(m_before_pwl)))
-    # vset.sort(key=lambda i: i.local_name)
-
     # build the PWL approximation
     trans = pe.TransformationFactory('contrib.piecewise.nonlinear_to_pwl')
     trans.apply_to(pwl, num_points=2, additively_decompose=False)
@@ -334,16 +330,6 @@ def _initialize_with_piecewise_linear_approximation(
         for ov, nv in zip(orig_vars, new_vars):
             ov.set_value(nv.value, skip_validation=True)
 
-        # print('\n\n**********************')
-        # print(m_before_pwl.obj.expr)
-        # print('\nConstraints:')
-        # for c in m_before_pwl.component_data_objects(pe.Constraint, active=True, descend_into=True):
-        #     print(c)
-        #     print('    ', c.expr)
-        # for v in vset:
-        #     print(v, v.value)
-        # m_before_pwl.display()
-
         # refine the PWL approximation
         _refine_pwl_approx(
             pwl,

From 48bae5223495427db0adfd840da74cefbca703a4 Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Wed, 15 Apr 2026 05:58:27 -0600
Subject: [PATCH 30/38] run black

---
 pyomo/devel/initialization/__init__.py        |  2 +-
 .../initialization/bounds/bound_variables.py  | 13 ++-
 .../examples/init_polynomial_ex.py            |  2 +-
 pyomo/devel/initialization/global_init.py     | 30 ++++--
 pyomo/devel/initialization/initialize.py      | 39 ++++----
 pyomo/devel/initialization/lp_approx_init.py  | 52 +++++++----
 pyomo/devel/initialization/pwl_init.py        | 93 ++++++++++++-------
 .../tests/test_initialization.py              | 64 ++++++-------
 pyomo/devel/initialization/utils.py           | 10 +-
 9 files changed, 187 insertions(+), 118 deletions(-)

diff --git a/pyomo/devel/initialization/__init__.py b/pyomo/devel/initialization/__init__.py
index 42e52b2d8ff..0c9fd0b1950 100644
--- a/pyomo/devel/initialization/__init__.py
+++ b/pyomo/devel/initialization/__init__.py
@@ -7,4 +7,4 @@
 # software.  This software is distributed under the 3-clause BSD License.
 # ____________________________________________________________________________________
 
-from pyomo.devel.initialization.initialize import initialize_nlp, InitializationMethod
\ No newline at end of file
+from pyomo.devel.initialization.initialize import initialize_nlp, InitializationMethod
diff --git a/pyomo/devel/initialization/bounds/bound_variables.py b/pyomo/devel/initialization/bounds/bound_variables.py
index 4402ea76510..ca6eabbc7f8 100644
--- a/pyomo/devel/initialization/bounds/bound_variables.py
+++ b/pyomo/devel/initialization/bounds/bound_variables.py
@@ -12,22 +12,25 @@
 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, 
+    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}')
+            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}')
+            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_polynomial_ex.py b/pyomo/devel/initialization/examples/init_polynomial_ex.py
index 111ad418363..d327c6c9e3c 100644
--- a/pyomo/devel/initialization/examples/init_polynomial_ex.py
+++ b/pyomo/devel/initialization/examples/init_polynomial_ex.py
@@ -16,7 +16,7 @@
 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)
+    m.c = pyo.Constraint(expr=(m.x + 7) * (m.x + 5) * (m.x - 4) + 200 == 0)
     return m
 
 
diff --git a/pyomo/devel/initialization/global_init.py b/pyomo/devel/initialization/global_init.py
index b924fb61273..c573e4e8162 100644
--- a/pyomo/devel/initialization/global_init.py
+++ b/pyomo/devel/initialization/global_init.py
@@ -14,28 +14,38 @@
 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,
+    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}')
+        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
index 1369db1a1f1..759a4887f88 100644
--- a/pyomo/devel/initialization/initialize.py
+++ b/pyomo/devel/initialization/initialize.py
@@ -12,7 +12,9 @@
 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.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
@@ -21,7 +23,6 @@
 import logging
 from pyomo.contrib.solver.common.results import SolutionStatus
 
-
 logger = logging.getLogger(__name__)
 
 
@@ -36,12 +37,14 @@ def _get_solver(sname, reason):
     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')
+        raise RuntimeError(
+            f'No solver was specified for {reason} and the default ({sname}) is not available'
+        )
     return opt
 
 
 def initialize_nlp(
-    nlp: BlockData, 
+    nlp: BlockData,
     nlp_solver: SolverBase,
     mip_solver: Optional[SolverBase] = None,
     global_solver: Optional[SolverBase] = None,
@@ -53,7 +56,7 @@ def initialize_nlp(
 ) -> 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 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
@@ -76,32 +79,34 @@ def initialize_nlp(
         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 
+        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 
+        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 
+        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 
+        Only used when method = InitializationMethod.pwl_approximation. This is
         passed along to the contrib.piecewise.univariate_nonlinear_decomposition
-        transformation. 
+        transformation.
 
     Returns
     -------
     res: pyomo.contrib.solver.common.results.Results
-        The results object obtained the last time the nlp_solver was used to 
+        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)
+    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:
@@ -136,18 +141,18 @@ def initialize_nlp(
         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,
+            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)
+        res = _initialize_with_global_solver(
+            nlp=nlp, global_solver=global_solver, nlp_solver=nlp_solver
+        )
     else:
-        raise ValueError(f'unexpected initialization method: {method}')    
+        raise ValueError(f'unexpected initialization method: {method}')
 
     # restore variable bounds, domain, etc.
     for v, (lb, ub, domain, fixed, value) in orig_var_data.items():
diff --git a/pyomo/devel/initialization/lp_approx_init.py b/pyomo/devel/initialization/lp_approx_init.py
index 91537436e10..583e467c912 100644
--- a/pyomo/devel/initialization/lp_approx_init.py
+++ b/pyomo/devel/initialization/lp_approx_init.py
@@ -9,10 +9,18 @@
 
 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.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_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
@@ -37,7 +45,11 @@
     NPV_SumExpression,
     NPV_UnaryFunctionExpression,
 )
-from pyomo.core.expr.relational_expr import EqualityExpression, InequalityExpression, RangedExpression
+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
@@ -53,7 +65,6 @@
 import numpy as np
 from scipy.stats import qmc
 
-
 logger = logging.getLogger(__name__)
 
 
@@ -105,10 +116,14 @@ def _build_lp_approx(nlp: BlockData) -> BlockData:
     lp.cons = pe.ConstraintList()
     visitor = LinearRepnVisitor(subexpression_cache={})
 
-    objs = list(nlp.component_data_objects(pe.Objective, active=True, descend_into=True))
+    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')
+            raise NotImplementedError(
+                'lp approximation does not support multiple objectives'
+            )
         obj = objs[0]
         repn = visitor.walk_expression(obj)
         assert repn.multiplier == 1
@@ -123,7 +138,9 @@ def _build_lp_approx(nlp: BlockData) -> BlockData:
             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):
+    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
@@ -144,17 +161,14 @@ def _build_lp_approx(nlp: BlockData) -> BlockData:
 
 
 def _initialize_with_LP_approximation(
-    nlp: BlockData,
-    lp_solver: SolverBase,
-    nlp_solver: SolverBase,
-    default_bound=1.0e8, 
+    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 
+    # 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
@@ -181,12 +195,18 @@ def _initialize_with_LP_approximation(
     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)
+    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}')
+    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()
diff --git a/pyomo/devel/initialization/pwl_init.py b/pyomo/devel/initialization/pwl_init.py
index 551e33a2f9a..f5513e34bd3 100644
--- a/pyomo/devel/initialization/pwl_init.py
+++ b/pyomo/devel/initialization/pwl_init.py
@@ -9,10 +9,18 @@
 
 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.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_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
@@ -36,7 +44,11 @@
     NPV_SumExpression,
     NPV_UnaryFunctionExpression,
 )
-from pyomo.core.expr.relational_expr import EqualityExpression, InequalityExpression, RangedExpression
+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
@@ -47,7 +59,6 @@
 from pyomo.common.modeling import unique_component_name
 from pyomo.contrib.solver.common.results import SolutionStatus
 
-
 logger = logging.getLogger(__name__)
 
 
@@ -61,9 +72,9 @@ def _minimize_infeasibility(m):
     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
+            obj_expr += 0.1 * obj.expr
         else:
-            obj_expr -= 0.1*obj.expr
+            obj_expr -= 0.1 * obj.expr
         obj.deactivate()
         found_obj = True
 
@@ -92,18 +103,21 @@ def _minimize_infeasibility(m):
             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)
+    m.slack_obj = pe.Objective(expr=10 * sum(m.slacks.values()) + obj_expr)
 
 
-def _get_pwl_constraints(m: BlockData) -> MutableMapping[
-    PiecewiseLinearExpression, 
-    List[ConstraintData]
-]:
+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))
+    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:
@@ -128,11 +142,11 @@ def _handle_node(node, data):
 for t in [float, int, VarData, ScalarVar, ParamData, ScalarParam, NumericConstant]:
     _handlers[t] = _handle_leaf
 for t in [
-    ProductExpression, 
-    SumExpression, 
-    DivisionExpression, 
-    PowExpression, 
-    MonomialTermExpression, 
+    ProductExpression,
+    SumExpression,
+    DivisionExpression,
+    PowExpression,
+    MonomialTermExpression,
     LinearExpression,
     ExpressionData,
     ScalarExpression,
@@ -170,7 +184,7 @@ def exitNode(self, node, data):
             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]
@@ -188,7 +202,9 @@ def beforeChild(self, node, child, child_idx):
         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')
+        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):
@@ -201,8 +217,7 @@ def beforeChild(self, node, child, child_idx):
 def _refine_pwl_approx(
     m,
     pwl_expr_to_con_map: MutableMapping[
-        PiecewiseLinearExpression, 
-        Sequence[ConstraintData],
+        PiecewiseLinearExpression, Sequence[ConstraintData]
     ],
     num_to_refine: int = 5,
 ):
@@ -235,12 +250,14 @@ def _refine_pwl_approx(
     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')
-    
+        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]:
@@ -259,11 +276,11 @@ def _refine_pwl_approx(
 
 
 def _initialize_with_piecewise_linear_approximation(
-    nlp: BlockData, 
+    nlp: BlockData,
     mip_solver: SolverBase,
     nlp_solver: SolverBase,
-    default_bound=1.0e8, 
-    max_iter=100, 
+    default_bound=1.0e8,
+    max_iter=100,
     num_cons_to_refine_per_iter=5,
     aggressive_substitution=True,
 ):
@@ -271,9 +288,11 @@ def _initialize_with_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 
+    # 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 = pe.TransformationFactory(
+        'contrib.piecewise.univariate_nonlinear_decomposition'
+    )
     trans.apply_to(pwl, aggressive_substitution=aggressive_substitution)
     logger.info('applied the univariate_nonlinear_decomposition transformation')
 
@@ -288,7 +307,7 @@ def _initialize_with_piecewise_linear_approximation(
 
     # now we modify the model by introducing slacks to make sure the PWL
     # approximatin is feasible
-    # all of the slacks appear linearly, so we don't need to worry about 
+    # 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')
@@ -323,10 +342,12 @@ def _initialize_with_piecewise_linear_approximation(
         logger.info('applied the disaggregated logarithmic transformation')
 
         # solve the MILP
-        res = mip_solver.solve(_pwl, load_solutions=True, raise_exception_on_nonoptimal_result=False)
+        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
+        # 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)
 
@@ -339,7 +360,9 @@ def _initialize_with_piecewise_linear_approximation(
         logger.info('refined PWL approximation')
 
         # try solving the NLP
-        res = nlp_solver.solve(nlp, load_solutions=False, raise_exception_on_nonoptimal_result=False)
+        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}:
diff --git a/pyomo/devel/initialization/tests/test_initialization.py b/pyomo/devel/initialization/tests/test_initialization.py
index 0ba33881570..bf78d3ee9f2 100644
--- a/pyomo/devel/initialization/tests/test_initialization.py
+++ b/pyomo/devel/initialization/tests/test_initialization.py
@@ -14,11 +14,14 @@
 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.results import (
+    SolutionStatus,
+    Results,
+    TerminationCondition,
+)
 from pyomo.contrib.solver.common.base import Availability, SolverBase
 import pytest
 
-
 scip = SolverFactory('scip_direct')
 ipopt = SolverFactory('ipopt')
 
@@ -32,16 +35,16 @@ def __init__(self, varlist, sol_map, **kwds) -> None:
 
     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()
@@ -95,17 +98,16 @@ def test_lp_init(self):
         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.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)
+        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),
-            },
+            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(
@@ -122,17 +124,16 @@ def test_global_init(self):
         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.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)
+        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),
-            },
+            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(
@@ -144,7 +145,7 @@ def test_global_init(self):
 
     def test_pwl_init(self):
         """
-        Here, we really just want to make sure that the 
+        Here, we really just want to make sure that the
         approximation improves as refinement is done.
         """
         m = pyo.ConcreteModel()
@@ -154,18 +155,18 @@ def test_pwl_init(self):
 
         # all the actual testing happens in the MockNLPSolver
         nlp_solver = MockNLPSolver(
-            varlist=[m.x], 
+            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),
+                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),
@@ -199,6 +200,7 @@ def test_pwl_init(self):
 
 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
index b73da8625b6..b8ce5095081 100644
--- a/pyomo/devel/initialization/utils.py
+++ b/pyomo/devel/initialization/utils.py
@@ -30,7 +30,9 @@ def shallow_clone(m1):
     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))
+    objlist = list(
+        m1.component_data_objects(pe.Objective, active=True, descend_into=True)
+    )
     assert len(objlist) <= 1
     if objlist:
         obj = objlist[0]
@@ -43,5 +45,9 @@ 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):
+        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))

From 0d1175f0afa360d7d17d20800a5ac6a83d3b009f Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Wed, 15 Apr 2026 07:09:50 -0600
Subject: [PATCH 31/38] fix typos

---
 .../analysis/nlp_initialization/nlp_initialization.rst          | 2 +-
 pyomo/devel/initialization/lp_approx_init.py                    | 2 +-
 pyomo/devel/initialization/pwl_init.py                          | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/doc/OnlineDocs/explanation/analysis/nlp_initialization/nlp_initialization.rst b/doc/OnlineDocs/explanation/analysis/nlp_initialization/nlp_initialization.rst
index 770f14369bf..1ce1ffaa0fc 100644
--- a/doc/OnlineDocs/explanation/analysis/nlp_initialization/nlp_initialization.rst
+++ b/doc/OnlineDocs/explanation/analysis/nlp_initialization/nlp_initialization.rst
@@ -102,7 +102,7 @@ Advantages
 Disadvantages
 ^^^^^^^^^^^^^
 
-* Current implemenation can be slow
+* Current implementation can be slow
 * Requires all nonlinear variables to be bounded
 
 
diff --git a/pyomo/devel/initialization/lp_approx_init.py b/pyomo/devel/initialization/lp_approx_init.py
index 583e467c912..44bdcb52a5c 100644
--- a/pyomo/devel/initialization/lp_approx_init.py
+++ b/pyomo/devel/initialization/lp_approx_init.py
@@ -186,7 +186,7 @@ def _initialize_with_LP_approximation(
     logger.info('fixed variables with equal bounds')
 
     # now we modify the model by introducing slacks to make sure the LP
-    # approximatin is feasible
+    # approximation is feasible
     _minimize_infeasibility(nlp)
     logger.info('reformulated model to minimize infeasibility')
 
diff --git a/pyomo/devel/initialization/pwl_init.py b/pyomo/devel/initialization/pwl_init.py
index f5513e34bd3..08ca99b9a75 100644
--- a/pyomo/devel/initialization/pwl_init.py
+++ b/pyomo/devel/initialization/pwl_init.py
@@ -306,7 +306,7 @@ def _initialize_with_piecewise_linear_approximation(
     logger.info('fixed variables with equal bounds')
 
     # now we modify the model by introducing slacks to make sure the PWL
-    # approximatin is feasible
+    # approximation is feasible
     # all of the slacks appear linearly, so we don't need to worry about
     # upper bounds for them
     _minimize_infeasibility(pwl)

From f5fe2373bd0ef981995b01916cb0be0b78138f88 Mon Sep 17 00:00:00 2001
From: Bethany Nicholson <blnicho@users.noreply.github.com>
Date: Tue, 19 May 2026 16:01:05 -0600
Subject: [PATCH 32/38] Fix ScipPersistent import

---
 pyomo/devel/initialization/global_init.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/pyomo/devel/initialization/global_init.py b/pyomo/devel/initialization/global_init.py
index c573e4e8162..cc558486ee2 100644
--- a/pyomo/devel/initialization/global_init.py
+++ b/pyomo/devel/initialization/global_init.py
@@ -10,7 +10,8 @@
 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.scip.scip_direct import ScipDirect
+from pyomo.contrib.solver.solvers.scip.scip_persistent import ScipPersistent
 from pyomo.contrib.solver.solvers.gurobi.gurobi_direct_minlp import GurobiDirectMINLP
 import logging
 

From 061e46f96a2d43073c98b78724fb913cf8214f1b Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Wed, 20 May 2026 05:21:53 -0600
Subject: [PATCH 33/38] update tests

---
 .../piecewise/tests/test_nonlinear_to_pwl.py  | 40 +++++++++----------
 1 file changed, 20 insertions(+), 20 deletions(-)

diff --git a/pyomo/contrib/piecewise/tests/test_nonlinear_to_pwl.py b/pyomo/contrib/piecewise/tests/test_nonlinear_to_pwl.py
index 61400eceadc..bf8b4ac2d61 100644
--- a/pyomo/contrib/piecewise/tests/test_nonlinear_to_pwl.py
+++ b/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)
 

From cba1223adc7cdd5f2a7e74ef2bcaefde9f1620e2 Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Wed, 20 May 2026 06:00:19 -0600
Subject: [PATCH 34/38] fix doctests

---
 doc/OnlineDocs/conf.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/doc/OnlineDocs/conf.py b/doc/OnlineDocs/conf.py
index 340af695bb0..c6f3662e519 100644
--- a/doc/OnlineDocs/conf.py
+++ b/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'))
+scip_available = bool(_opt.check_available_solvers('scip_direct'))
 
 baron = _opt.SolverFactory('baron')
 

From c85a86422ce8130f73bcb6f352f5c4c06c52213a Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Wed, 20 May 2026 06:09:38 -0600
Subject: [PATCH 35/38] use attempt_import

---
 pyomo/devel/initialization/lp_approx_init.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/pyomo/devel/initialization/lp_approx_init.py b/pyomo/devel/initialization/lp_approx_init.py
index 44bdcb52a5c..0ebacdc2539 100644
--- a/pyomo/devel/initialization/lp_approx_init.py
+++ b/pyomo/devel/initialization/lp_approx_init.py
@@ -62,8 +62,10 @@
 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
+from pyomo.common.dependencies import numpy as np
+from pyomo.common.dependencies import attempt_import
+
+qmc, _ = attempt_import('scipy.stats.qmc')
 
 logger = logging.getLogger(__name__)
 

From a33c665694b95db454a1041d4e7a09241dd10766 Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Wed, 20 May 2026 06:21:26 -0600
Subject: [PATCH 36/38] update tests

---
 pyomo/devel/initialization/tests/test_initialization.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/pyomo/devel/initialization/tests/test_initialization.py b/pyomo/devel/initialization/tests/test_initialization.py
index bf78d3ee9f2..08e5c99b6ff 100644
--- a/pyomo/devel/initialization/tests/test_initialization.py
+++ b/pyomo/devel/initialization/tests/test_initialization.py
@@ -24,6 +24,7 @@
 
 scip = SolverFactory('scip_direct')
 ipopt = SolverFactory('ipopt')
+highs = SolverFactory('highs')
 
 
 class MockNLPSolver(SolverBase):
@@ -77,6 +78,7 @@ def test_poly_lp(self):
 
 
 class TestInit(unittest.TestCase):
+    @unittest.skipUnless(highs.available(), 'highs is not available')
     def test_lp_init(self):
         """
         For this test, we will create a small linear program,
@@ -117,6 +119,7 @@ def test_lp_init(self):
             method=ini.InitializationMethod.lp_approximation,
         )
 
+    @unittest.skipUnless(scip.available(), 'scip is not available')
     def test_global_init(self):
         """
         Same as test_lp_init
@@ -143,6 +146,7 @@ def test_global_init(self):
             method=ini.InitializationMethod.global_opt,
         )
 
+    @unittest.skipUnless(highs.available(), 'highs is not available')
     def test_pwl_init(self):
         """
         Here, we really just want to make sure that the

From a800114e5d5e7f1f5457c1caa446b37fd94a8a75 Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Wed, 20 May 2026 07:27:29 -0600
Subject: [PATCH 37/38] guard tests

---
 pyomo/devel/initialization/lp_approx_init.py            | 2 +-
 pyomo/devel/initialization/tests/test_initialization.py | 2 ++
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/pyomo/devel/initialization/lp_approx_init.py b/pyomo/devel/initialization/lp_approx_init.py
index 0ebacdc2539..f61e34ad4b8 100644
--- a/pyomo/devel/initialization/lp_approx_init.py
+++ b/pyomo/devel/initialization/lp_approx_init.py
@@ -65,7 +65,7 @@
 from pyomo.common.dependencies import numpy as np
 from pyomo.common.dependencies import attempt_import
 
-qmc, _ = attempt_import('scipy.stats.qmc')
+qmc, qmc_avail = attempt_import('scipy.stats.qmc')
 
 logger = logging.getLogger(__name__)
 
diff --git a/pyomo/devel/initialization/tests/test_initialization.py b/pyomo/devel/initialization/tests/test_initialization.py
index 08e5c99b6ff..f7f3213a240 100644
--- a/pyomo/devel/initialization/tests/test_initialization.py
+++ b/pyomo/devel/initialization/tests/test_initialization.py
@@ -11,6 +11,7 @@
 
 import pyomo.environ as pyo
 import pyomo.devel.initialization as ini
+from pyomo.devel.initialization.lp_approx_init import qmc_avail
 from pyomo.devel.initialization.examples.init_polynomial_ex import main
 from pyomo.common import unittest
 from pyomo.contrib.solver.common.factory import SolverFactory
@@ -79,6 +80,7 @@ def test_poly_lp(self):
 
 class TestInit(unittest.TestCase):
     @unittest.skipUnless(highs.available(), 'highs is not available')
+    @unittest.skipUnless(qmc_avail, 'scipy.stats.qmc is not available')
     def test_lp_init(self):
         """
         For this test, we will create a small linear program,

From 2d4594a29e5656d20ba4ea7ffbacbffac333f962 Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Wed, 20 May 2026 07:56:05 -0600
Subject: [PATCH 38/38] try highs instead of scip

---
 pyomo/devel/initialization/examples/init_polynomial_ex.py | 2 +-
 pyomo/devel/initialization/tests/test_initialization.py   | 4 +++-
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/pyomo/devel/initialization/examples/init_polynomial_ex.py b/pyomo/devel/initialization/examples/init_polynomial_ex.py
index d327c6c9e3c..6550675f98f 100644
--- a/pyomo/devel/initialization/examples/init_polynomial_ex.py
+++ b/pyomo/devel/initialization/examples/init_polynomial_ex.py
@@ -24,7 +24,7 @@ def main(method: ini.InitializationMethod):
     m = build_model()
     nlp_solver = SolverFactory('ipopt')
     global_solver = SolverFactory('scip_direct')
-    mip_solver = SolverFactory('scip_direct')
+    mip_solver = SolverFactory('highs')
     results = ini.initialize_nlp(
         nlp=m,
         nlp_solver=nlp_solver,
diff --git a/pyomo/devel/initialization/tests/test_initialization.py b/pyomo/devel/initialization/tests/test_initialization.py
index f7f3213a240..11c8289948d 100644
--- a/pyomo/devel/initialization/tests/test_initialization.py
+++ b/pyomo/devel/initialization/tests/test_initialization.py
@@ -59,19 +59,21 @@ def solve(self, model, **kwds) -> Results:
         return res
 
 
-@unittest.skipUnless(scip.available(), 'scip is not available')
 @unittest.skipUnless(ipopt.available(), 'ipopt is not available')
 class TestExamples(unittest.TestCase):
+    @unittest.skipUnless(scip.available(), 'scip is not available')
     def test_poly_global(self):
         stat, x = main(method=ini.InitializationMethod.global_opt)
         self.assertEqual(stat, SolutionStatus.optimal)
         self.assertAlmostEqual(x, -9.920159607881597)
 
+    @unittest.skipUnless(highs.available(), 'highs is not available')
     def test_poly_pwl(self):
         stat, x = main(method=ini.InitializationMethod.pwl_approximation)
         self.assertEqual(stat, SolutionStatus.optimal)
         self.assertAlmostEqual(x, -9.920159607881597)
 
+    @unittest.skipUnless(highs.available(), 'highs is not available')
     def test_poly_lp(self):
         stat, x = main(method=ini.InitializationMethod.lp_approximation)
         self.assertEqual(stat, SolutionStatus.optimal)