From 7eb0cbe205b3774c70e509c7cf3787f8198cca5b Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Thu, 2 Jul 2026 13:07:42 -0600
Subject: [PATCH 01/26] Add time limit warning to global solver

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

diff --git a/pyomo/devel/initialization/global_init.py b/pyomo/devel/initialization/global_init.py
index 83b8afba98e..8a90b373f94 100644
--- a/pyomo/devel/initialization/global_init.py
+++ b/pyomo/devel/initialization/global_init.py
@@ -31,6 +31,12 @@ def _initialize_with_global_solver(
             'interfaces, so the global solvers are limited to ScipDirect, '
             'ScipPersistent, and GurobiDirectMINLP.'
         )
+    # Check if time limit is provided for global solver
+    if global_solver.config.time_limit is None:
+        logger.warning('No time limit set for global optimizer. ' 
+                       'For a large model, this may take a long time. ' 
+                       'Consider setting a time limit using global_solver.config.time_limit.')
+
     res = global_solver.solve(
         nlp,
         load_solutions=True,

From 0f97cb0274e1e723c630225d68bcad24bb872da7 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Wed, 8 Jul 2026 08:41:30 -0600
Subject: [PATCH 02/26] Ran black

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

diff --git a/pyomo/devel/initialization/global_init.py b/pyomo/devel/initialization/global_init.py
index 8a90b373f94..10b33a9a0cb 100644
--- a/pyomo/devel/initialization/global_init.py
+++ b/pyomo/devel/initialization/global_init.py
@@ -33,9 +33,11 @@ def _initialize_with_global_solver(
         )
     # Check if time limit is provided for global solver
     if global_solver.config.time_limit is None:
-        logger.warning('No time limit set for global optimizer. ' 
-                       'For a large model, this may take a long time. ' 
-                       'Consider setting a time limit using global_solver.config.time_limit.')
+        logger.warning(
+            'No time limit set for global optimizer. '
+            'For a large model, this may take a long time. '
+            'Consider setting a time limit using global_solver.config.time_limit.'
+        )
 
     res = global_solver.solve(
         nlp,

From 352c828c54b09972577b42bdc5bf93707cc33280 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Wed, 8 Jul 2026 09:29:06 -0600
Subject: [PATCH 03/26] Adjusted order to solve unbounded nlp at end of method

---
 pyomo/devel/initialization/global_init.py    | 10 -----
 pyomo/devel/initialization/initialize.py     | 41 +++++++++++++++++++-
 pyomo/devel/initialization/lp_approx_init.py | 15 +------
 3 files changed, 40 insertions(+), 26 deletions(-)

diff --git a/pyomo/devel/initialization/global_init.py b/pyomo/devel/initialization/global_init.py
index 10b33a9a0cb..3f3ecb9fee0 100644
--- a/pyomo/devel/initialization/global_init.py
+++ b/pyomo/devel/initialization/global_init.py
@@ -48,15 +48,5 @@ 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.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 7b1117cf5b9..fd481644b30 100644
--- a/pyomo/devel/initialization/initialize.py
+++ b/pyomo/devel/initialization/initialize.py
@@ -155,8 +155,22 @@ def initialize_with_piecewise_linear_approximation(
         )
     finally:
         _cleanup(orig_var_data)
+        # Try final nlp solve
 
-    return res
+        # solve the original problem from the initialized solution
+        nlp_res = nlp_solver.solve(
+            nlp, load_solutions=False, raise_exception_on_nonoptimal_result=False
+        )
+        logger.info(
+            f'solved NLP: {nlp_res.solution_status}, {nlp_res.termination_condition}'
+        )
+
+        if nlp_res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
+            nlp_res.solution_loader.load_vars()
+        else:
+            logger.warning('initialization was not successful via LP approximation')
+
+    return nlp_res
 
 
 def initialize_with_LP_approximation(
@@ -239,7 +253,20 @@ def initialize_with_LP_approximation(
     finally:
         _cleanup(orig_var_data)
 
-    return res
+        # solve the original problem from the initialized solution
+        nlp_res = nlp_solver.solve(
+            nlp, load_solutions=False, raise_exception_on_nonoptimal_result=False
+        )
+        logger.info(
+            f'solved NLP: {nlp_res.solution_status}, {nlp_res.termination_condition}'
+        )
+
+        if nlp_res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
+            nlp_res.solution_loader.load_vars()
+        else:
+            logger.warning('initialization was not successful via LP approximation')
+
+    return nlp_res
 
 
 def initialize_with_global_opt(
@@ -293,4 +320,14 @@ def initialize_with_global_opt(
     finally:
         _cleanup(orig_var_data)
 
+        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/lp_approx_init.py b/pyomo/devel/initialization/lp_approx_init.py
index c8f7c2ff4fc..d20c92c0203 100644
--- a/pyomo/devel/initialization/lp_approx_init.py
+++ b/pyomo/devel/initialization/lp_approx_init.py
@@ -218,17 +218,4 @@ def _initialize_with_LP_approximation(
     )
     logger.info(f'solved LP: {lp_res.solution_status}, {lp_res.termination_condition}')
 
-    # try solving the NLP
-    nlp_res = nlp_solver.solve(
-        orig_nlp, load_solutions=False, raise_exception_on_nonoptimal_result=False
-    )
-    logger.info(
-        f'solved NLP: {nlp_res.solution_status}, {nlp_res.termination_condition}'
-    )
-
-    if nlp_res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
-        nlp_res.solution_loader.load_vars()
-    else:
-        logger.warning('initialization was not successful via LP approximation')
-
-    return nlp_res
+    return lp_res

From 8c0a358441100f400f07e935fd0d2cfbd73969ce Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Wed, 8 Jul 2026 10:30:53 -0600
Subject: [PATCH 04/26] Made naming consistent, commented out pwl while
 debugging tests

---
 pyomo/devel/initialization/initialize.py | 40 +++++++++++++-----------
 1 file changed, 22 insertions(+), 18 deletions(-)

diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py
index fd481644b30..1371ab077af 100644
--- a/pyomo/devel/initialization/initialize.py
+++ b/pyomo/devel/initialization/initialize.py
@@ -155,22 +155,25 @@ def initialize_with_piecewise_linear_approximation(
         )
     finally:
         _cleanup(orig_var_data)
-        # Try final nlp solve
 
-        # solve the original problem from the initialized solution
-        nlp_res = nlp_solver.solve(
-            nlp, load_solutions=False, raise_exception_on_nonoptimal_result=False
-        )
-        logger.info(
-            f'solved NLP: {nlp_res.solution_status}, {nlp_res.termination_condition}'
-        )
+        # Commented out while I fix testing
+        # # Try final nlp solve
 
-        if nlp_res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
-            nlp_res.solution_loader.load_vars()
-        else:
-            logger.warning('initialization was not successful via LP approximation')
+        # # solve the original problem from the initialized solution
+        # nlp_res = nlp_solver.solve(
+        #     nlp, load_solutions=False, raise_exception_on_nonoptimal_result=False
+        # )
+        # logger.info(
+        #     f'solved NLP: {nlp_res.solution_status}, {nlp_res.termination_condition}'
+        # )
 
-    return nlp_res
+        # if nlp_res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
+        #     nlp_res.solution_loader.load_vars()
+        # else:
+        #     logger.warning('initialization was not successful via LP approximation')
+
+    # return nlp_res
+    return res
 
 
 def initialize_with_LP_approximation(
@@ -320,14 +323,15 @@ def initialize_with_global_opt(
     finally:
         _cleanup(orig_var_data)
 
-        res = nlp_solver.solve(
+        nlp_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}'
+            f'solved NLP with {nlp_solver.name}: {nlp_res.solution_status}, {nlp_res.termination_condition}'
         )
-        if res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
-            res.solution_loader.load_vars()
+        if nlp_res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
+            nlp_res.solution_loader.load_vars()
         else:
             logger.warning('initialization was not successful via global optimization')
-    return res
+
+    return nlp_res

From b6947a0c66446d83687d53649a97ea68b2f89e09 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Wed, 8 Jul 2026 10:52:51 -0600
Subject: [PATCH 05/26] Fixed test, ran black

---
 pyomo/devel/initialization/initialize.py      | 29 +++++++++----------
 .../tests/test_initialization.py              |  5 ++--
 2 files changed, 17 insertions(+), 17 deletions(-)

diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py
index 1371ab077af..b227d7819d4 100644
--- a/pyomo/devel/initialization/initialize.py
+++ b/pyomo/devel/initialization/initialize.py
@@ -156,21 +156,20 @@ def initialize_with_piecewise_linear_approximation(
     finally:
         _cleanup(orig_var_data)
 
-        # Commented out while I fix testing
-        # # Try final nlp solve
-
-        # # solve the original problem from the initialized solution
-        # nlp_res = nlp_solver.solve(
-        #     nlp, load_solutions=False, raise_exception_on_nonoptimal_result=False
-        # )
-        # logger.info(
-        #     f'solved NLP: {nlp_res.solution_status}, {nlp_res.termination_condition}'
-        # )
-
-        # if nlp_res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
-        #     nlp_res.solution_loader.load_vars()
-        # else:
-        #     logger.warning('initialization was not successful via LP approximation')
+        # Try final nlp solve
+
+        # solve the original problem from the initialized solution
+        nlp_res = nlp_solver.solve(
+            nlp, load_solutions=False, raise_exception_on_nonoptimal_result=False
+        )
+        logger.info(
+            f'solved NLP: {nlp_res.solution_status}, {nlp_res.termination_condition}'
+        )
+
+        if nlp_res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
+            nlp_res.solution_loader.load_vars()
+        else:
+            logger.warning('initialization was not successful via LP approximation')
 
     # return nlp_res
     return res
diff --git a/pyomo/devel/initialization/tests/test_initialization.py b/pyomo/devel/initialization/tests/test_initialization.py
index 67a260d9904..b9851888f39 100644
--- a/pyomo/devel/initialization/tests/test_initialization.py
+++ b/pyomo/devel/initialization/tests/test_initialization.py
@@ -196,6 +196,8 @@ def test_pwl_init(self):
                 25: ([-9.91992877683681], 1e-6, 1e-6),
                 26: ([-9.920038488200985], 1e-6, 1e-6),
                 27: ([-9.920096055464825], 1e-6, 1e-6),
+                # For new second solve, repeat last value.
+                28: ([-9.920096055464825], 1e-6, 1e-6),
             },
         )
         mip_solver = SolverFactory('highs')
@@ -230,5 +232,4 @@ def test_pwl_ineq(self):
     import logging
 
     logging.basicConfig(level=logging.INFO)
-    t = TestInit()
-    t.test_pwl_init()
+    t = TestInit()
\ No newline at end of file

From 3cb15856b34415431aa58b59326ca88ec0447b16 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Wed, 8 Jul 2026 10:55:10 -0600
Subject: [PATCH 06/26] Forgot comment, ran black

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

diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py
index b227d7819d4..815cd94669e 100644
--- a/pyomo/devel/initialization/initialize.py
+++ b/pyomo/devel/initialization/initialize.py
@@ -157,7 +157,6 @@ def initialize_with_piecewise_linear_approximation(
         _cleanup(orig_var_data)
 
         # Try final nlp solve
-
         # solve the original problem from the initialized solution
         nlp_res = nlp_solver.solve(
             nlp, load_solutions=False, raise_exception_on_nonoptimal_result=False
@@ -171,8 +170,7 @@ def initialize_with_piecewise_linear_approximation(
         else:
             logger.warning('initialization was not successful via LP approximation')
 
-    # return nlp_res
-    return res
+    return nlp_res
 
 
 def initialize_with_LP_approximation(
diff --git a/pyomo/devel/initialization/tests/test_initialization.py b/pyomo/devel/initialization/tests/test_initialization.py
index b9851888f39..042bb890828 100644
--- a/pyomo/devel/initialization/tests/test_initialization.py
+++ b/pyomo/devel/initialization/tests/test_initialization.py
@@ -232,4 +232,4 @@ def test_pwl_ineq(self):
     import logging
 
     logging.basicConfig(level=logging.INFO)
-    t = TestInit()
\ No newline at end of file
+    t = TestInit()

From cc40909598114cf598898506bf8614436adc88ef Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Wed, 8 Jul 2026 13:53:46 -0600
Subject: [PATCH 07/26] Initial setup progress for multistart initialization

---
 pyomo/contrib/multistart/multi.py             | 23 ++++--
 pyomo/devel/initialization/initialize.py      | 75 +++++++++++++++++++
 pyomo/devel/initialization/multistart_init.py | 40 ++++++++++
 3 files changed, 130 insertions(+), 8 deletions(-)
 create mode 100644 pyomo/devel/initialization/multistart_init.py

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index 32fb8ce3cff..46dfe12c749 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -120,6 +120,13 @@ class MultiStart:
             description="Tolerance on HCS objective value equality. Defaults to Python float equality precision.",
         ),
     )
+    CONFIG.declare(
+        "seed",
+        ConfigValue(
+            default=None,
+            description="Seed for reproducibility in random sampling methods"
+        )
+    )
 
     def available(self, exception_flag=True):
         """Check if solver is available.
@@ -149,14 +156,14 @@ def solve(self, model, **kwds):
             raise RuntimeError(
                 "Multistart solver is unable to handle model with multiple active objectives."
             )
-        if obj is None:
-            raise RuntimeError(
-                "Multistart solver is unable to handle model with no active objective."
-            )
-        if obj.polynomial_degree() == 0:
-            raise RuntimeError(
-                "Multistart solver received model with constant objective"
-            )
+        # if obj is None:
+        #     raise RuntimeError(
+        #         "Multistart solver is unable to handle model with no active objective."
+        #     )
+        # if obj.polynomial_degree() == 0:
+        #     raise RuntimeError(
+        #         "Multistart solver received model with constant objective"
+        #     )
 
         # store objective values and objective/result information for best
         # solution obtained
diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py
index 815cd94669e..9d2f633c021 100644
--- a/pyomo/devel/initialization/initialize.py
+++ b/pyomo/devel/initialization/initialize.py
@@ -18,10 +18,14 @@
 from pyomo.devel.initialization.lp_approx_init import _initialize_with_LP_approximation
 from pyomo.contrib.solver.common.base import SolverBase
 from pyomo.devel.initialization.global_init import _initialize_with_global_solver
+from pyomo.devel.initialization.multistart_init import (
+    _initialize_with_multistart_solver,
+)
 from pyomo.contrib.solver.common.factory import SolverFactory
 from pyomo.contrib.solver.common.results import Results
 import logging
 from pyomo.contrib.solver.common.results import SolutionStatus
+import pyomo.environ as pyo
 
 logger = logging.getLogger(__name__)
 
@@ -332,3 +336,74 @@ def initialize_with_global_opt(
             logger.warning('initialization was not successful via global optimization')
 
     return nlp_res
+
+
+def initialize_with_multistart_opt(
+    nlp: BlockData,
+    nlp_solver: SolverBase | None = None,
+    multistart_solver = None,
+    skip_initial_nlp_solve: bool = False,
+    default_bound: float = 1e8,
+    seed=0,
+
+) -> 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.
+    nlp_solver: Optional[SolverBase]
+        A solver interface appropriate for NLPs.
+        Default: ipopt
+    multistart_solver: Optional[SolverFactory]
+        A configured multistart solver object for performing multistart optimization
+    skip_initial_nlp_solve: bool
+        If True, the initial attempt at solving the NLP without initialization
+        will be skipped.
+    seed: 0
+        Set reproducibility seed to make result deterministic.
+
+    Returns
+    -------
+    res: pyomo.contrib.solver.common.results.Results
+        The results object obtained the last time the nlp_solver was used to
+        try and solve the model.
+    """
+
+    if nlp_solver is None:
+        nlp_solver = _get_solver('ipopt', 'local NLP solver')
+
+    if multistart_solver is None:
+        multistart_solver = pyo.SolverFactory("multistart")
+
+    if not skip_initial_nlp_solve:
+        res = _try_nlp_solve(nlp, nlp_solver)
+        if res.solution_status == SolutionStatus.optimal:
+            return res
+
+    orig_var_data = _setup(nlp)
+
+    try:
+        res = _initialize_with_multistart_solver(
+            nlp=nlp, multistart_solver=multistart_solver,
+            default_bound=default_bound, seed=seed
+        )
+    finally:
+        _cleanup(orig_var_data)
+
+        nlp_res = nlp_solver.solve(
+            nlp, load_solutions=False, raise_exception_on_nonoptimal_result=False
+        )
+        logger.info(
+            f'solved NLP with {nlp_solver.name}: {nlp_res.solution_status}, {nlp_res.termination_condition}'
+        )
+        if nlp_res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
+            nlp_res.solution_loader.load_vars()
+        else:
+            logger.warning('initialization was not successful via global optimization')
+
+    return nlp_res
diff --git a/pyomo/devel/initialization/multistart_init.py b/pyomo/devel/initialization/multistart_init.py
new file mode 100644
index 00000000000..d3e169133be
--- /dev/null
+++ b/pyomo/devel/initialization/multistart_init.py
@@ -0,0 +1,40 @@
+# ____________________________________________________________________________________
+#
+# 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.multistart.multi import Multistart
+import pyomo.environ as pyo
+from pyomo.contrib.solver.common.results import SolutionStatus
+from pyomo.devel.initialization.bounds.bound_variables import (
+    bound_all_nonlinear_variables,
+)
+from pyomo.devel.initialization.utils import shallow_clone
+import logging
+
+logger = logging.getLogger(__name__)
+
+
+def _initialize_with_multistart_solver(
+    nlp: BlockData, 
+    multistart_solver,
+    default_bound=1e6,
+    seed = None,
+    ):
+
+    # Make a shallow clone
+    nlp = shallow_clone(nlp)
+    # bounds on the nonlinear variables
+    bound_all_nonlinear_variables(nlp, default_bound=default_bound)
+    res = multistart_solver.solve(nlp)
+    logger.info(
+        'solved multistart run'
+        )
+
+    return res

From 4e46f906edecc892e63a31a124a444f4eae820b5 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Wed, 8 Jul 2026 15:03:56 -0600
Subject: [PATCH 08/26] Work on lhs, in progress

---
 pyomo/contrib/multistart/reinit.py       | 32 ++++++++++++++++++++++++
 pyomo/devel/initialization/initialize.py |  2 +-
 2 files changed, 33 insertions(+), 1 deletion(-)

diff --git a/pyomo/contrib/multistart/reinit.py b/pyomo/contrib/multistart/reinit.py
index a3b52a8d611..5a7ecec8e31 100644
--- a/pyomo/contrib/multistart/reinit.py
+++ b/pyomo/contrib/multistart/reinit.py
@@ -11,6 +11,10 @@
 
 import logging
 import random
+from pyomo.common.dependencies.scipy import stats
+from pyomo.core.expr.visitor import (
+    identify_variables,
+)
 
 from pyomo.core import Var
 
@@ -20,6 +24,26 @@
 def rand(val, lb, ub):
     return random.uniform(lb, ub)  # uniform distribution between lb and ub
 
+def latin_hypercube(val, lb, ub, sampler):
+    sample = sampler.random(n=1)
+    sample = stats.qmc.scale(sample, lb, ub)
+    return sample
+
+def _generate_lhs_sample(vlist, config):
+    n_vars = len(vlist)
+    bnds_list = []
+    for v in vlist:
+        # the bounds should not be None because we
+        # set the bounds to default_bound in
+        # bound_all_nonlinear_variables
+        lb = v.lb
+        ub = v.ub
+        bnds_list.append((lb, ub))
+    sampler = stats.qmc.LatinHypercube(d=n_vars, seed=config.seed)
+    sample = sampler.random(n=config.seed)
+    l_bounds = [i[0] for i in bnds_list]
+    u_bounds = [i[1] for i in bnds_list]
+    sample = stats.qmc.scale(sample, l_bounds, u_bounds)
 
 def midpoint_guess_and_bound(val, lb, ub):
     """Midpoint between current value and farthest bound."""
@@ -54,6 +78,7 @@ def linspace(lower, upper, n):
     "rand_guess_and_bound": rand_guess_and_bound,
     "rand_distributed": rand_distributed,
     "midpoint": simple_midpoint,
+    "latin_hypercube": latin_hypercube
 }
 
 
@@ -63,6 +88,13 @@ def reinitialize_variables(model, config):
     Excludes fixed, noncontinuous, and unbounded variables.
 
     """
+    if config.strategy is "latin_hypercube":
+        vlist = list(identify_variables(model, include_fixed=False))
+        _
+        for v in vlist:
+
+        
+    else:
     for var in model.component_data_objects(ctype=Var, descend_into=True):
         if var.is_fixed() or not var.is_continuous():
             continue
diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py
index 9d2f633c021..c93545b02a9 100644
--- a/pyomo/devel/initialization/initialize.py
+++ b/pyomo/devel/initialization/initialize.py
@@ -404,6 +404,6 @@ def initialize_with_multistart_opt(
         if nlp_res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
             nlp_res.solution_loader.load_vars()
         else:
-            logger.warning('initialization was not successful via global optimization')
+            logger.warning('initialization was not successful via multistart optimization')
 
     return nlp_res

From fb66e6c34e92ff8b5237ba1e737fb6aeaec67f56 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Wed, 8 Jul 2026 16:18:19 -0600
Subject: [PATCH 09/26] Semi working version of multistart!

---
 pyomo/contrib/multistart/multi.py             | 15 ++++++++++-
 pyomo/contrib/multistart/reinit.py            | 27 ++++++++++---------
 pyomo/devel/initialization/initialize.py      |  1 +
 pyomo/devel/initialization/multistart_init.py |  5 ++--
 4 files changed, 32 insertions(+), 16 deletions(-)

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index 46dfe12c749..682347bb764 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -17,6 +17,7 @@
     document_kwargs_from_configdict,
 )
 from pyomo.common.modeling import unique_component_name
+from pyomo.common.dependencies import numpy as np
 from pyomo.contrib.multistart.high_conf_stop import should_stop
 from pyomo.contrib.multistart.reinit import reinitialize_variables, strategies
 from pyomo.core import Objective, Var, minimize, value
@@ -127,6 +128,14 @@ class MultiStart:
             description="Seed for reproducibility in random sampling methods"
         )
     )
+    CONFIG.declare(
+        "rng",
+        ConfigValue(
+            default=None,
+            description="Random number generator for reproducibility in random sampling methods." \
+                "Preferred over seed."
+        )
+    )
 
     def available(self, exception_flag=True):
         """Check if solver is available.
@@ -144,6 +153,9 @@ def solve(self, model, **kwds):
         # initialize keyword args
         config = self.CONFIG(kwds.pop('options', {}))
         config.set_value(kwds)
+        
+        if config.rng is None:
+            config.rng = np.random.default_rng(config.seed)
 
         # initialize the solver
         solver = SolverFactory(config.solver)
@@ -211,11 +223,12 @@ def solve(self, model, **kwds):
                 ):
                     HCS_completed = True
                     break
+                print(f"num_iter: {num_iter}\n")
                 num_iter += 1
                 # at first iteration, solve the originally passed model
                 m = model.clone() if num_iter > 1 else model
                 reinitialize_variables(m, config)
-                result = solver.solve(m, **config.solver_args)
+                result = solver.solve(m, **config.solver_args, ) # tee = True)
                 if (
                     result.solver.status is SolverStatus.ok
                     and result.solver.termination_condition is tc.optimal
diff --git a/pyomo/contrib/multistart/reinit.py b/pyomo/contrib/multistart/reinit.py
index 5a7ecec8e31..6e66010995c 100644
--- a/pyomo/contrib/multistart/reinit.py
+++ b/pyomo/contrib/multistart/reinit.py
@@ -11,6 +11,7 @@
 
 import logging
 import random
+from pyomo.common.dependencies import numpy as np
 from pyomo.common.dependencies.scipy import stats
 from pyomo.core.expr.visitor import (
     identify_variables,
@@ -21,8 +22,10 @@
 logger = logging.getLogger('pyomo.contrib.multistart')
 
 
-def rand(val, lb, ub):
-    return random.uniform(lb, ub)  # uniform distribution between lb and ub
+def rand(val, lb, ub, rng):
+    sample = rng.uniform(lb, ub) # uniform distribution between lb and ub
+    print(f"sample={sample})\n")
+    return sample
 
 def latin_hypercube(val, lb, ub, sampler):
     sample = sampler.random(n=1)
@@ -51,16 +54,16 @@ def midpoint_guess_and_bound(val, lb, ub):
     return (far_bound + val) / 2
 
 
-def rand_guess_and_bound(val, lb, ub):
+def rand_guess_and_bound(val, lb, ub, rng):
     """Random choice between current value and farthest bound."""
     far_bound = ub if ((ub - val) >= (val - lb)) else lb  # farther bound
-    return random.uniform(val, far_bound)
+    return rng.uniform(val, far_bound)
 
 
-def rand_distributed(val, lb, ub, divisions=9):
+def rand_distributed(val, lb, ub, rng, divisions=9):
     """Random choice among evenly distributed set of values between bounds."""
     set_distributed_vals = linspace(lb, ub, divisions)
-    return random.choice(set_distributed_vals)
+    return rng.choice(set_distributed_vals)
 
 
 def simple_midpoint(val, lb, ub):
@@ -88,13 +91,9 @@ def reinitialize_variables(model, config):
     Excludes fixed, noncontinuous, and unbounded variables.
 
     """
-    if config.strategy is "latin_hypercube":
-        vlist = list(identify_variables(model, include_fixed=False))
-        _
-        for v in vlist:
-
+    # if config.strategy == "latin_hypercube":
+    #     vlist = list(identify_variables(model, include_fixed=False))
         
-    else:
     for var in model.component_data_objects(ctype=Var, descend_into=True):
         if var.is_fixed() or not var.is_continuous():
             continue
@@ -108,7 +107,9 @@ def reinitialize_variables(model, config):
                 )
             continue
         val = var.value if var.value is not None else (var.lb + var.ub) / 2
+        print(f"val = {val}\n")
         # apply reinitialization strategy to variable
         var.set_value(
-            strategies[config.strategy](val, var.lb, var.ub), skip_validation=True
+            strategies[config.strategy](val, var.lb, var.ub, config.rng), skip_validation=True
+
         )
diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py
index c93545b02a9..5fa895af2c6 100644
--- a/pyomo/devel/initialization/initialize.py
+++ b/pyomo/devel/initialization/initialize.py
@@ -379,6 +379,7 @@ def initialize_with_multistart_opt(
 
     if multistart_solver is None:
         multistart_solver = pyo.SolverFactory("multistart")
+        multistart_solver.CONFIG.seed=seed
 
     if not skip_initial_nlp_solve:
         res = _try_nlp_solve(nlp, nlp_solver)
diff --git a/pyomo/devel/initialization/multistart_init.py b/pyomo/devel/initialization/multistart_init.py
index d3e169133be..bfa28542180 100644
--- a/pyomo/devel/initialization/multistart_init.py
+++ b/pyomo/devel/initialization/multistart_init.py
@@ -9,7 +9,7 @@
 
 from pyomo.core.base.block import BlockData
 from pyomo.contrib.solver.common.base import SolverBase
-# from pyomo.contrib.multistart.multi import Multistart
+from pyomo.common.dependencies import numpy as np
 import pyomo.environ as pyo
 from pyomo.contrib.solver.common.results import SolutionStatus
 from pyomo.devel.initialization.bounds.bound_variables import (
@@ -27,11 +27,12 @@ def _initialize_with_multistart_solver(
     default_bound=1e6,
     seed = None,
     ):
-
+    
     # Make a shallow clone
     nlp = shallow_clone(nlp)
     # bounds on the nonlinear variables
     bound_all_nonlinear_variables(nlp, default_bound=default_bound)
+
     res = multistart_solver.solve(nlp)
     logger.info(
         'solved multistart run'

From 86cbd736a53f6493d4326ceb079503d4c9fcd43b Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Thu, 9 Jul 2026 07:16:29 -0600
Subject: [PATCH 10/26] Changes to resolve testing errors, in progress

---
 pyomo/contrib/multistart/multi.py             | 39 +++++++++++++++----
 pyomo/devel/initialization/initialize.py      |  1 +
 pyomo/devel/initialization/multistart_init.py |  2 +-
 3 files changed, 33 insertions(+), 9 deletions(-)

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index 682347bb764..4c473a0473c 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -125,15 +125,23 @@ class MultiStart:
         "seed",
         ConfigValue(
             default=None,
-            description="Seed for reproducibility in random sampling methods"
+            description="Seed for reproducibility in random sampling methods."
         )
     )
     CONFIG.declare(
         "rng",
         ConfigValue(
             default=None,
-            description="Random number generator for reproducibility in random sampling methods." \
-                "Preferred over seed."
+            description="Random number generator for reproducibility in random sampling methods. \
+                Preferred over seed."
+        )
+    )
+    CONFIG.declare(
+        "new_solvers_bool",
+        ConfigValue(
+            default=False,
+            description="Boolean option for whether to use the new solver interface, default to no \
+                until solver testing complete (?)"
         )
     )
 
@@ -158,6 +166,12 @@ def solve(self, model, **kwds):
             config.rng = np.random.default_rng(config.seed)
 
         # initialize the solver
+        if config.new_solvers_bool == True:
+            from pyomo.contrib.solver.common.factory import SolverFactory
+            from pyomo.contrib.solver.common.results import Results
+            from pyomo.contrib.solver.common.results import SolutionStatus
+
+
         solver = SolverFactory(config.solver)
 
         # Model sense
@@ -195,10 +209,15 @@ def solve(self, model, **kwds):
             )
 
             best_result = result = solver.solve(model, **config.solver_args)
+            if best_result.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
+                best_result.solution_loader.load_vars()
+                logger.info(f'solved NLP: {best_result.solution_status}, {best_result.termination_condition}')
+            
             if (
-                result.solver.status is SolverStatus.ok
-                and result.solver.termination_condition is tc.optimal
+                result.solution_status is SolverStatus.ok
+                and result.termination_condition is tc.optimal
             ):
+                            
                 obj_val = value(obj.expr)
                 best_objective = obj_val
                 objectives.append(obj_val)
@@ -228,10 +247,14 @@ def solve(self, model, **kwds):
                 # at first iteration, solve the originally passed model
                 m = model.clone() if num_iter > 1 else model
                 reinitialize_variables(m, config)
-                result = solver.solve(m, **config.solver_args, ) # tee = True)
+                result = solver.solve(m, **config.solver_args) #, tee=True)
+
+                if result.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
+                    result.solution_loader.load_vars()
+                    logger.info(f'solved NLP: {result.solution_status}, {result.termination_condition}')   
                 if (
-                    result.solver.status is SolverStatus.ok
-                    and result.solver.termination_condition is tc.optimal
+                    result.solution_status is SolverStatus.ok
+                    and result.termination_condition is tc.optimal
                 ):
                     model_objectives = m.component_data_objects(Objective, active=True)
                     mobj = next(model_objectives)
diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py
index 5fa895af2c6..953dc2a607a 100644
--- a/pyomo/devel/initialization/initialize.py
+++ b/pyomo/devel/initialization/initialize.py
@@ -380,6 +380,7 @@ def initialize_with_multistart_opt(
     if multistart_solver is None:
         multistart_solver = pyo.SolverFactory("multistart")
         multistart_solver.CONFIG.seed=seed
+        multistart_solver.CONFIG.new_solvers_bool=True
 
     if not skip_initial_nlp_solve:
         res = _try_nlp_solve(nlp, nlp_solver)
diff --git a/pyomo/devel/initialization/multistart_init.py b/pyomo/devel/initialization/multistart_init.py
index bfa28542180..70307f35fea 100644
--- a/pyomo/devel/initialization/multistart_init.py
+++ b/pyomo/devel/initialization/multistart_init.py
@@ -35,7 +35,7 @@ def _initialize_with_multistart_solver(
 
     res = multistart_solver.solve(nlp)
     logger.info(
-        'solved multistart run'
+        'Finished multistart optimization iterations.'
         )
 
     return res

From 2f17c59550a1af92283587ea41e3eed4df80f8e5 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Thu, 9 Jul 2026 10:15:55 -0600
Subject: [PATCH 11/26] Added more loggers to debug stuff

---
 pyomo/contrib/multistart/high_conf_stop.py |  5 +++++
 pyomo/contrib/multistart/multi.py          | 18 +++++++++++-------
 2 files changed, 16 insertions(+), 7 deletions(-)

diff --git a/pyomo/contrib/multistart/high_conf_stop.py b/pyomo/contrib/multistart/high_conf_stop.py
index b9efdfd65f9..ae9558ef05f 100644
--- a/pyomo/contrib/multistart/high_conf_stop.py
+++ b/pyomo/contrib/multistart/high_conf_stop.py
@@ -17,6 +17,7 @@
 
 from collections import Counter
 from math import log, sqrt
+import logger
 
 
 def num_one_occurrences(observed_obj_vals, tolerance):
@@ -59,4 +60,8 @@ def should_stop(solutions, stopping_mass, stopping_delta, tolerance):
     d = stopping_delta
     c = stopping_mass
     confidence = f / n + (2 * sqrt(2) + sqrt(3)) * sqrt(log(3 / d) / n)
+    # Add temporary logger
+    logger.info(f"Number of solutions [n]:{n}; Optima viewed once [f]:{f}; \
+    Confidence:{confidence}"
+    )
     return confidence < c
diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index 4c473a0473c..86d535ee7c2 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -209,15 +209,17 @@ def solve(self, model, **kwds):
             )
 
             best_result = result = solver.solve(model, **config.solver_args)
-            if best_result.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
-                best_result.solution_loader.load_vars()
-                logger.info(f'solved NLP: {best_result.solution_status}, {best_result.termination_condition}')
+            # only use one condition to record results, this might be error causing
+            # if best_result.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
+            #     best_result.solution_loader.load_vars()
+            #     logger.info(f'solved NLP: {best_result.solution_status}, {best_result.termination_condition}')
             
             if (
                 result.solution_status is SolverStatus.ok
                 and result.termination_condition is tc.optimal
             ):
-                            
+                best_result.solution_loader.load_vars()
+                logger.info(f'solved NLP: {best_result.solution_status}, {best_result.termination_condition}')
                 obj_val = value(obj.expr)
                 best_objective = obj_val
                 objectives.append(obj_val)
@@ -249,13 +251,15 @@ def solve(self, model, **kwds):
                 reinitialize_variables(m, config)
                 result = solver.solve(m, **config.solver_args) #, tee=True)
 
-                if result.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
-                    result.solution_loader.load_vars()
-                    logger.info(f'solved NLP: {result.solution_status}, {result.termination_condition}')   
+                # if result.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
+                #     result.solution_loader.load_vars()
+                #     logger.info(f'solved NLP: {result.solution_status}, {result.termination_condition}')   
                 if (
                     result.solution_status is SolverStatus.ok
                     and result.termination_condition is tc.optimal
                 ):
+                    result.solution_loader.load_vars()
+                    logger.info(f'solved NLP: {result.solution_status}, {result.termination_condition}') 
                     model_objectives = m.component_data_objects(Objective, active=True)
                     mobj = next(model_objectives)
                     obj_val = value(mobj.expr)

From 65f73986123d587adfe14144f655c43328f33a9e Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Thu, 9 Jul 2026 10:22:52 -0600
Subject: [PATCH 12/26] Fix typo

---
 pyomo/contrib/multistart/high_conf_stop.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pyomo/contrib/multistart/high_conf_stop.py b/pyomo/contrib/multistart/high_conf_stop.py
index ae9558ef05f..2373cf06e0f 100644
--- a/pyomo/contrib/multistart/high_conf_stop.py
+++ b/pyomo/contrib/multistart/high_conf_stop.py
@@ -17,7 +17,7 @@
 
 from collections import Counter
 from math import log, sqrt
-import logger
+import logging
 
 
 def num_one_occurrences(observed_obj_vals, tolerance):
@@ -61,7 +61,7 @@ def should_stop(solutions, stopping_mass, stopping_delta, tolerance):
     c = stopping_mass
     confidence = f / n + (2 * sqrt(2) + sqrt(3)) * sqrt(log(3 / d) / n)
     # Add temporary logger
-    logger.info(f"Number of solutions [n]:{n}; Optima viewed once [f]:{f}; \
+    logging.info(f"Number of solutions [n]:{n}; Optima viewed once [f]:{f}; \
     Confidence:{confidence}"
     )
     return confidence < c

From 55a1e883150ed09b670a94dd4ca9b240555c132c Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Thu, 9 Jul 2026 10:35:34 -0600
Subject: [PATCH 13/26] Solution loading too strict

---
 pyomo/contrib/multistart/multi.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index 86d535ee7c2..ea9d5f1378f 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -213,13 +213,14 @@ def solve(self, model, **kwds):
             # if best_result.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
             #     best_result.solution_loader.load_vars()
             #     logger.info(f'solved NLP: {best_result.solution_status}, {best_result.termination_condition}')
+            logger.info(f'solved NLP: {best_result.solution_status}, {best_result.termination_condition}')
             
             if (
                 result.solution_status is SolverStatus.ok
                 and result.termination_condition is tc.optimal
             ):
-                best_result.solution_loader.load_vars()
-                logger.info(f'solved NLP: {best_result.solution_status}, {best_result.termination_condition}')
+                # best_result.solution_loader.load_vars()
+
                 obj_val = value(obj.expr)
                 best_objective = obj_val
                 objectives.append(obj_val)

From c86cbc08aca76ca388be6bd140a3cbd2757a9401 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Thu, 9 Jul 2026 10:42:25 -0600
Subject: [PATCH 14/26] Caused other issues, putting this back

---
 pyomo/contrib/multistart/multi.py | 20 +++++++-------------
 1 file changed, 7 insertions(+), 13 deletions(-)

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index ea9d5f1378f..55d60bd46ac 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -209,18 +209,14 @@ def solve(self, model, **kwds):
             )
 
             best_result = result = solver.solve(model, **config.solver_args)
-            # only use one condition to record results, this might be error causing
-            # if best_result.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
-            #     best_result.solution_loader.load_vars()
-            #     logger.info(f'solved NLP: {best_result.solution_status}, {best_result.termination_condition}')
-            logger.info(f'solved NLP: {best_result.solution_status}, {best_result.termination_condition}')
+            if best_result.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
+                best_result.solution_loader.load_vars()
+                logger.info(f'solved NLP: {best_result.solution_status}, {best_result.termination_condition}')
             
             if (
                 result.solution_status is SolverStatus.ok
                 and result.termination_condition is tc.optimal
-            ):
-                # best_result.solution_loader.load_vars()
-
+            ):      
                 obj_val = value(obj.expr)
                 best_objective = obj_val
                 objectives.append(obj_val)
@@ -252,15 +248,13 @@ def solve(self, model, **kwds):
                 reinitialize_variables(m, config)
                 result = solver.solve(m, **config.solver_args) #, tee=True)
 
-                # if result.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
-                #     result.solution_loader.load_vars()
-                #     logger.info(f'solved NLP: {result.solution_status}, {result.termination_condition}')   
+                if result.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
+                    result.solution_loader.load_vars()
+                    logger.info(f'solved NLP: {result.solution_status}, {result.termination_condition}')   
                 if (
                     result.solution_status is SolverStatus.ok
                     and result.termination_condition is tc.optimal
                 ):
-                    result.solution_loader.load_vars()
-                    logger.info(f'solved NLP: {result.solution_status}, {result.termination_condition}') 
                     model_objectives = m.component_data_objects(Objective, active=True)
                     mobj = next(model_objectives)
                     obj_val = value(mobj.expr)

From 6aa9ce0e9a80d0eaf40a1f31dc7102bc5092a898 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Thu, 9 Jul 2026 13:17:05 -0600
Subject: [PATCH 15/26] define helper fcn to retry solve, removed from finally
 block

---
 pyomo/devel/initialization/initialize.py | 55 +++++++++---------------
 1 file changed, 20 insertions(+), 35 deletions(-)

diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py
index 815cd94669e..6d449fa8e54 100644
--- a/pyomo/devel/initialization/initialize.py
+++ b/pyomo/devel/initialization/initialize.py
@@ -76,6 +76,21 @@ def _try_nlp_solve(nlp: BlockData, nlp_solver: SolverBase):
         logger.info('NLP solved without any initialization')
     return res
 
+def _retry_nlp_solve(nlp: BlockData, nlp_solver: SolverBase):
+    # retry to solve the original nlp after using an initialization method
+    nlp_res = nlp_solver.solve(
+        nlp, load_solutions=False, raise_exception_on_nonoptimal_result=False
+    )
+    logger.info(
+        f'solved NLP with {nlp_solver.name}: {nlp_res.solution_status}, {nlp_res.termination_condition}'
+    )
+    if nlp_res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
+        nlp_res.solution_loader.load_vars()
+    else:
+        logger.warning('initialization did not find feasible solution')
+
+    return nlp_res
+
 
 def initialize_with_piecewise_linear_approximation(
     nlp: BlockData,
@@ -156,19 +171,7 @@ def initialize_with_piecewise_linear_approximation(
     finally:
         _cleanup(orig_var_data)
 
-        # Try final nlp solve
-        # solve the original problem from the initialized solution
-        nlp_res = nlp_solver.solve(
-            nlp, load_solutions=False, raise_exception_on_nonoptimal_result=False
-        )
-        logger.info(
-            f'solved NLP: {nlp_res.solution_status}, {nlp_res.termination_condition}'
-        )
-
-        if nlp_res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
-            nlp_res.solution_loader.load_vars()
-        else:
-            logger.warning('initialization was not successful via LP approximation')
+    nlp_res = _retry_nlp_solve(nlp, nlp_solver)
 
     return nlp_res
 
@@ -253,18 +256,7 @@ def initialize_with_LP_approximation(
     finally:
         _cleanup(orig_var_data)
 
-        # solve the original problem from the initialized solution
-        nlp_res = nlp_solver.solve(
-            nlp, load_solutions=False, raise_exception_on_nonoptimal_result=False
-        )
-        logger.info(
-            f'solved NLP: {nlp_res.solution_status}, {nlp_res.termination_condition}'
-        )
-
-        if nlp_res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
-            nlp_res.solution_loader.load_vars()
-        else:
-            logger.warning('initialization was not successful via LP approximation')
+    nlp_res = _retry_nlp_solve(nlp, nlp_solver)
 
     return nlp_res
 
@@ -320,15 +312,8 @@ def initialize_with_global_opt(
     finally:
         _cleanup(orig_var_data)
 
-        nlp_res = nlp_solver.solve(
-            nlp, load_solutions=False, raise_exception_on_nonoptimal_result=False
-        )
-        logger.info(
-            f'solved NLP with {nlp_solver.name}: {nlp_res.solution_status}, {nlp_res.termination_condition}'
-        )
-        if nlp_res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
-            nlp_res.solution_loader.load_vars()
-        else:
-            logger.warning('initialization was not successful via global optimization')
+    nlp_res = _retry_nlp_solve(nlp, nlp_solver)
 
     return nlp_res
+
+

From 3753f4b81923c46b3537d7897d0eee8030e26f07 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Thu, 9 Jul 2026 13:18:44 -0600
Subject: [PATCH 16/26] Ran black

---
 pyomo/devel/initialization/initialize.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py
index 6d449fa8e54..2fd4f5c145c 100644
--- a/pyomo/devel/initialization/initialize.py
+++ b/pyomo/devel/initialization/initialize.py
@@ -76,6 +76,7 @@ def _try_nlp_solve(nlp: BlockData, nlp_solver: SolverBase):
         logger.info('NLP solved without any initialization')
     return res
 
+
 def _retry_nlp_solve(nlp: BlockData, nlp_solver: SolverBase):
     # retry to solve the original nlp after using an initialization method
     nlp_res = nlp_solver.solve(
@@ -315,5 +316,3 @@ def initialize_with_global_opt(
     nlp_res = _retry_nlp_solve(nlp, nlp_solver)
 
     return nlp_res
-
-

From 9dd6c48616a06b182de26e9e32c8e674828987b9 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Thu, 9 Jul 2026 13:28:23 -0600
Subject: [PATCH 17/26] Moved string, Ran black

---
 pyomo/devel/initialization/initialize.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py
index 2fd4f5c145c..58f8f9c0a97 100644
--- a/pyomo/devel/initialization/initialize.py
+++ b/pyomo/devel/initialization/initialize.py
@@ -82,9 +82,8 @@ def _retry_nlp_solve(nlp: BlockData, nlp_solver: SolverBase):
     nlp_res = nlp_solver.solve(
         nlp, load_solutions=False, raise_exception_on_nonoptimal_result=False
     )
-    logger.info(
-        f'solved NLP with {nlp_solver.name}: {nlp_res.solution_status}, {nlp_res.termination_condition}'
-    )
+    logger.info(f'resolved NLP with {nlp_solver.name}: {nlp_res.solution_status}, \
+              {nlp_res.termination_condition}')
     if nlp_res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
         nlp_res.solution_loader.load_vars()
     else:

From 2870aff0f1a3640d525eec283981554e2f3fecd7 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Thu, 9 Jul 2026 13:44:13 -0600
Subject: [PATCH 18/26] Made multistart also use nlp solve

---
 pyomo/devel/initialization/initialize.py | 11 +----------
 1 file changed, 1 insertion(+), 10 deletions(-)

diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py
index 6531348f8e7..c793e14c1c1 100644
--- a/pyomo/devel/initialization/initialize.py
+++ b/pyomo/devel/initialization/initialize.py
@@ -380,15 +380,6 @@ def initialize_with_multistart_opt(
     finally:
         _cleanup(orig_var_data)
 
-        nlp_res = nlp_solver.solve(
-            nlp, load_solutions=False, raise_exception_on_nonoptimal_result=False
-        )
-        logger.info(
-            f'solved NLP with {nlp_solver.name}: {nlp_res.solution_status}, {nlp_res.termination_condition}'
-        )
-        if nlp_res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
-            nlp_res.solution_loader.load_vars()
-        else:
-            logger.warning('initialization was not successful via multistart optimization')
+    nlp_res = _retry_nlp_solve(nlp, nlp_solver)
 
     return nlp_res

From e84f9fb0b58d86941ea55bd8c3f07bb7a7e88030 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Mon, 13 Jul 2026 07:38:38 -0600
Subject: [PATCH 19/26] Adding new features for initialization, in progress

---
 pyomo/contrib/multistart/multi.py             | 17 ++++++++++++++++-
 pyomo/contrib/multistart/reinit.py            | 12 ++++++++----
 pyomo/devel/initialization/multistart_init.py |  2 +-
 3 files changed, 25 insertions(+), 6 deletions(-)

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index 55d60bd46ac..806b8304355 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -121,6 +121,21 @@ class MultiStart:
             description="Tolerance on HCS objective value equality. Defaults to Python float equality precision.",
         ),
     )
+    CONFIG.declare(
+        "break_when_optimal",
+        ConfigValue(
+            default=False,
+            description="Condition to break if a feasible or optimal solution is found. Defaults to False."
+        )
+    )
+    CONFIG.declare(
+        "sampling_method",
+        ConfigValue(
+            default="random_uniform",
+            description="Method for sampling random starting points for reinitialization step. Supported options are \
+            'random_uniform', 'latin_hypercube', and 'sobol_sampling'"
+        )
+    )
     CONFIG.declare(
         "seed",
         ConfigValue(
@@ -241,7 +256,7 @@ def solve(self, model, **kwds):
                 ):
                     HCS_completed = True
                     break
-                print(f"num_iter: {num_iter}\n")
+                logger.info(f"num_iter: {num_iter}\n")
                 num_iter += 1
                 # at first iteration, solve the originally passed model
                 m = model.clone() if num_iter > 1 else model
diff --git a/pyomo/contrib/multistart/reinit.py b/pyomo/contrib/multistart/reinit.py
index 6e66010995c..84091ab8887 100644
--- a/pyomo/contrib/multistart/reinit.py
+++ b/pyomo/contrib/multistart/reinit.py
@@ -22,9 +22,13 @@
 logger = logging.getLogger('pyomo.contrib.multistart')
 
 
-def rand(val, lb, ub, rng):
-    sample = rng.uniform(lb, ub) # uniform distribution between lb and ub
-    print(f"sample={sample})\n")
+def rand(val, lb, ub, rng, sampling="random_uniform"):
+    # sample = rng.uniform(lb, ub) # uniform distribution between lb and ub
+    # print(f"sample={sample})\n")
+
+    # Changing to other style
+    # Basic layout
+    # sample = _generate_sample()
     return sample
 
 def latin_hypercube(val, lb, ub, sampler):
@@ -32,7 +36,7 @@ def latin_hypercube(val, lb, ub, sampler):
     sample = stats.qmc.scale(sample, lb, ub)
     return sample
 
-def _generate_lhs_sample(vlist, config):
+def _generate_sample(vlist, config):
     n_vars = len(vlist)
     bnds_list = []
     for v in vlist:
diff --git a/pyomo/devel/initialization/multistart_init.py b/pyomo/devel/initialization/multistart_init.py
index 70307f35fea..394decb6cc1 100644
--- a/pyomo/devel/initialization/multistart_init.py
+++ b/pyomo/devel/initialization/multistart_init.py
@@ -24,7 +24,7 @@
 def _initialize_with_multistart_solver(
     nlp: BlockData, 
     multistart_solver,
-    default_bound=1e6,
+    default_bound=1.0e8,
     seed = None,
     ):
     

From 64ea87acfb3c80ab7f6e34690ab5882675bc5ab6 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Mon, 13 Jul 2026 10:00:37 -0600
Subject: [PATCH 20/26] First draft of additional sampling support and
 break_on_solution

---
 pyomo/contrib/multistart/multi.py  | 87 +++++++++++++++++++++++++-----
 pyomo/contrib/multistart/reinit.py | 72 ++++++++++++-------------
 2 files changed, 108 insertions(+), 51 deletions(-)

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index 806b8304355..ca4bf0dc88a 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -23,6 +23,8 @@
 from pyomo.core import Objective, Var, minimize, value
 from pyomo.opt import SolverFactory, SolverStatus
 from pyomo.opt import TerminationCondition as tc
+from pyomo.common.dependencies.scipy import stats
+from pyomo.common.dependencies import numpy as np
 
 logger = logging.getLogger('pyomo.contrib.multistart')
 
@@ -122,7 +124,7 @@ class MultiStart:
         ),
     )
     CONFIG.declare(
-        "break_when_optimal",
+        "break_on_solution",
         ConfigValue(
             default=False,
             description="Condition to break if a feasible or optimal solution is found. Defaults to False."
@@ -176,9 +178,6 @@ def solve(self, model, **kwds):
         # initialize keyword args
         config = self.CONFIG(kwds.pop('options', {}))
         config.set_value(kwds)
-        
-        if config.rng is None:
-            config.rng = np.random.default_rng(config.seed)
 
         # initialize the solver
         if config.new_solvers_bool == True:
@@ -186,6 +185,9 @@ def solve(self, model, **kwds):
             from pyomo.contrib.solver.common.results import Results
             from pyomo.contrib.solver.common.results import SolutionStatus
 
+        # Create centralized sampler once
+        sampler = SamplingManager(method=config.sampling_method, 
+                                  rng=config.rng, seed=config.seed)
 
         solver = SolverFactory(config.solver)
 
@@ -224,17 +226,19 @@ def solve(self, model, **kwds):
             )
 
             best_result = result = solver.solve(model, **config.solver_args)
+            # Check the solution status before loading variables into the model.
             if best_result.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
                 best_result.solution_loader.load_vars()
                 logger.info(f'solved NLP: {best_result.solution_status}, {best_result.termination_condition}')
+                # If we are looking for the first feasible solution, then return immediately
+                if config.break_on_solution:
+                    return best_result
             
-            if (
-                result.solution_status is SolverStatus.ok
-                and result.termination_condition is tc.optimal
-            ):      
+            if result.termination_condition is tc.optimal:      
                 obj_val = value(obj.expr)
                 best_objective = obj_val
                 objectives.append(obj_val)
+
             num_iter = 0
             max_iter = config.iterations
             # if HCS rule is specified, reinitialize completely randomly until
@@ -260,16 +264,19 @@ def solve(self, model, **kwds):
                 num_iter += 1
                 # at first iteration, solve the originally passed model
                 m = model.clone() if num_iter > 1 else model
-                reinitialize_variables(m, config)
+                reinitialize_variables(m, config, sampler)
                 result = solver.solve(m, **config.solver_args) #, tee=True)
 
                 if result.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
                     result.solution_loader.load_vars()
-                    logger.info(f'solved NLP: {result.solution_status}, {result.termination_condition}')   
-                if (
-                    result.solution_status is SolverStatus.ok
-                    and result.termination_condition is tc.optimal
-                ):
+                    logger.info(f'solved NLP: {result.solution_status}, {result.termination_condition}')
+                    if config.break_on_solution:
+                        best_model = m
+                        best_result = result
+                        break
+
+
+                if result.termination_condition is tc.optimal:
                     model_objectives = m.component_data_objects(Objective, active=True)
                     mobj = next(model_objectives)
                     obj_val = value(mobj.expr)
@@ -315,3 +322,55 @@ def __enter__(self):
 
     def __exit__(self, t, v, traceback):
         pass
+
+# Sampling class to organize and configure random samplers
+
+class SamplingManager:
+    def __init__(self, method="uniform", rng=None, seed=None):
+        aliases = {
+            "random_uniform": "uniform",
+            "uniform": "uniform",
+            "latin_hypercube": "lhs",
+            "lhs": "lhs",
+            "sobol_sampling": "sobol",
+            "sobol": "sobol",
+        }
+        self.method = aliases[method.lower()]
+        
+        self.seed = seed
+
+        # Define or create a random number generator
+        # All 
+
+        if rng is not None:
+            self.rng = rng
+        else:
+            self.rng = np.random.default_rng(seed)
+
+        self.qmc_sampler = None
+
+    def _ensure_qmc(self, dim):
+        if self.qmc_sampler is not None:
+            return
+
+        if self.method == "lhs":
+            self.qmc_sampler = stats.qmc.LatinHypercube(d=dim, seed=self.seed)
+        elif self.method == "sobol":
+            self.qmc_sampler = stats.qmc.Sobol(d=dim, scramble=True, seed=self.seed)
+        else:
+            raise ValueError(f"QMC sampler not valid for method '{self.method}'")
+
+    def sample_vector(self, lower, upper):
+        """Vector sample for uniform/lhs/sobol over all vars at once."""
+        lower = np.asarray(lower, dtype=float)
+        upper = np.asarray(upper, dtype=float)
+
+        if self.method == "uniform":
+            return self.rng.uniform(lower, upper)
+
+        if self.method in ("lhs", "sobol"):
+            self._ensure_qmc(dim=len(lower))
+            x = self.qmc_sampler.random(n=1)   # shape (1, d)
+            return stats.qmc.scale(x, lower, upper)[0]
+
+        raise ValueError(f"Unknown sampling method '{self.method}'")
\ No newline at end of file
diff --git a/pyomo/contrib/multistart/reinit.py b/pyomo/contrib/multistart/reinit.py
index 84091ab8887..557e49b650b 100644
--- a/pyomo/contrib/multistart/reinit.py
+++ b/pyomo/contrib/multistart/reinit.py
@@ -22,37 +22,11 @@
 logger = logging.getLogger('pyomo.contrib.multistart')
 
 
-def rand(val, lb, ub, rng, sampling="random_uniform"):
-    # sample = rng.uniform(lb, ub) # uniform distribution between lb and ub
-    # print(f"sample={sample})\n")
-
-    # Changing to other style
-    # Basic layout
-    # sample = _generate_sample()
-    return sample
-
-def latin_hypercube(val, lb, ub, sampler):
-    sample = sampler.random(n=1)
-    sample = stats.qmc.scale(sample, lb, ub)
+def rand(val, lb, ub, rng):
+    sample = rng.uniform(lb, ub) # uniform distribution between lb and ub
     return sample
 
-def _generate_sample(vlist, config):
-    n_vars = len(vlist)
-    bnds_list = []
-    for v in vlist:
-        # the bounds should not be None because we
-        # set the bounds to default_bound in
-        # bound_all_nonlinear_variables
-        lb = v.lb
-        ub = v.ub
-        bnds_list.append((lb, ub))
-    sampler = stats.qmc.LatinHypercube(d=n_vars, seed=config.seed)
-    sample = sampler.random(n=config.seed)
-    l_bounds = [i[0] for i in bnds_list]
-    u_bounds = [i[1] for i in bnds_list]
-    sample = stats.qmc.scale(sample, l_bounds, u_bounds)
-
-def midpoint_guess_and_bound(val, lb, ub):
+def midpoint_guess_and_bound(val, lb, ub, rng=None):
     """Midpoint between current value and farthest bound."""
     far_bound = ub if ((ub - val) >= (val - lb)) else lb  # farther bound
     return (far_bound + val) / 2
@@ -70,7 +44,7 @@ def rand_distributed(val, lb, ub, rng, divisions=9):
     return rng.choice(set_distributed_vals)
 
 
-def simple_midpoint(val, lb, ub):
+def simple_midpoint(val, lb, ub, rng=None):
     return (lb + ub) * 0.5
 
 
@@ -85,20 +59,20 @@ def linspace(lower, upper, n):
     "rand_guess_and_bound": rand_guess_and_bound,
     "rand_distributed": rand_distributed,
     "midpoint": simple_midpoint,
-    "latin_hypercube": latin_hypercube
+
 }
 
 
-def reinitialize_variables(model, config):
+def reinitialize_variables(model, config, sampler):
     """Reinitializes all variable values in the model.
 
     Excludes fixed, noncontinuous, and unbounded variables.
 
     """
-    # if config.strategy == "latin_hypercube":
-    #     vlist = list(identify_variables(model, include_fixed=False))
-        
-    for var in model.component_data_objects(ctype=Var, descend_into=True):
+
+    eligible_vars = []
+
+    for var in model.component_data_objects(ctype=Var, descend_into=True):     
         if var.is_fixed() or not var.is_continuous():
             continue
         if var.lb is None or var.ub is None:
@@ -110,10 +84,34 @@ def reinitialize_variables(model, config):
                     'suppress_unbounded_warning flag.' % (var.name, var.lb, var.ub)
                 )
             continue
+
+        eligible_vars.append(var)
+
+    # Sample for new methods as a vector
+    if sampler.method in {"uniform", "lhs", "sobol"}:
+        if len(eligible_vars) == 0:
+                raise ValueError("No eligible variables to reinitialize." \
+                "Please add bounds.")
+        
+        # Collect lower and upper bounds for sampler
+        lowers = [v.lb for v in eligible_vars]
+        uppers = [v.ub for v in eligible_vars]
+        
+        # Generate vector of samples using sampler
+        samples = sampler.sample_vector(lowers, uppers)
+
+        # assign samples to variables
+        for var, sample in zip(eligible_vars, samples):
+                var.set_value(sample, skip_validation=True)
+
+        return
+    
+    # Otherwise use strategies to maintain original functionality
+    for var in eligible_vars:
         val = var.value if var.value is not None else (var.lb + var.ub) / 2
         print(f"val = {val}\n")
         # apply reinitialization strategy to variable
         var.set_value(
-            strategies[config.strategy](val, var.lb, var.ub, config.rng), skip_validation=True
+            strategies[config.strategy](val, var.lb, var.ub, sampler.rng), skip_validation=True
 
         )

From 4424e25ddd1d1d372ec6952e3b57ef6833ccd843 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Mon, 13 Jul 2026 10:01:24 -0600
Subject: [PATCH 21/26] Ran black

---
 pyomo/contrib/multistart/high_conf_stop.py |  3 +-
 pyomo/contrib/multistart/multi.py          | 60 +++++++++++++---------
 pyomo/contrib/multistart/reinit.py         | 27 +++++-----
 3 files changed, 50 insertions(+), 40 deletions(-)

diff --git a/pyomo/contrib/multistart/high_conf_stop.py b/pyomo/contrib/multistart/high_conf_stop.py
index 2373cf06e0f..91fd0244775 100644
--- a/pyomo/contrib/multistart/high_conf_stop.py
+++ b/pyomo/contrib/multistart/high_conf_stop.py
@@ -62,6 +62,5 @@ def should_stop(solutions, stopping_mass, stopping_delta, tolerance):
     confidence = f / n + (2 * sqrt(2) + sqrt(3)) * sqrt(log(3 / d) / n)
     # Add temporary logger
     logging.info(f"Number of solutions [n]:{n}; Optima viewed once [f]:{f}; \
-    Confidence:{confidence}"
-    )
+    Confidence:{confidence}")
     return confidence < c
diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index ca4bf0dc88a..513c9f77ab5 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -127,39 +127,39 @@ class MultiStart:
         "break_on_solution",
         ConfigValue(
             default=False,
-            description="Condition to break if a feasible or optimal solution is found. Defaults to False."
-        )
+            description="Condition to break if a feasible or optimal solution is found. Defaults to False.",
+        ),
     )
     CONFIG.declare(
         "sampling_method",
         ConfigValue(
             default="random_uniform",
             description="Method for sampling random starting points for reinitialization step. Supported options are \
-            'random_uniform', 'latin_hypercube', and 'sobol_sampling'"
-        )
+            'random_uniform', 'latin_hypercube', and 'sobol_sampling'",
+        ),
     )
     CONFIG.declare(
         "seed",
         ConfigValue(
             default=None,
-            description="Seed for reproducibility in random sampling methods."
-        )
+            description="Seed for reproducibility in random sampling methods.",
+        ),
     )
     CONFIG.declare(
         "rng",
         ConfigValue(
             default=None,
             description="Random number generator for reproducibility in random sampling methods. \
-                Preferred over seed."
-        )
+                Preferred over seed.",
+        ),
     )
     CONFIG.declare(
         "new_solvers_bool",
         ConfigValue(
             default=False,
             description="Boolean option for whether to use the new solver interface, default to no \
-                until solver testing complete (?)"
-        )
+                until solver testing complete (?)",
+        ),
     )
 
     def available(self, exception_flag=True):
@@ -186,8 +186,9 @@ def solve(self, model, **kwds):
             from pyomo.contrib.solver.common.results import SolutionStatus
 
         # Create centralized sampler once
-        sampler = SamplingManager(method=config.sampling_method, 
-                                  rng=config.rng, seed=config.seed)
+        sampler = SamplingManager(
+            method=config.sampling_method, rng=config.rng, seed=config.seed
+        )
 
         solver = SolverFactory(config.solver)
 
@@ -227,14 +228,19 @@ def solve(self, model, **kwds):
 
             best_result = result = solver.solve(model, **config.solver_args)
             # Check the solution status before loading variables into the model.
-            if best_result.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
+            if best_result.solution_status in {
+                SolutionStatus.feasible,
+                SolutionStatus.optimal,
+            }:
                 best_result.solution_loader.load_vars()
-                logger.info(f'solved NLP: {best_result.solution_status}, {best_result.termination_condition}')
+                logger.info(
+                    f'solved NLP: {best_result.solution_status}, {best_result.termination_condition}'
+                )
                 # If we are looking for the first feasible solution, then return immediately
                 if config.break_on_solution:
                     return best_result
-            
-            if result.termination_condition is tc.optimal:      
+
+            if result.termination_condition is tc.optimal:
                 obj_val = value(obj.expr)
                 best_objective = obj_val
                 objectives.append(obj_val)
@@ -265,17 +271,21 @@ def solve(self, model, **kwds):
                 # at first iteration, solve the originally passed model
                 m = model.clone() if num_iter > 1 else model
                 reinitialize_variables(m, config, sampler)
-                result = solver.solve(m, **config.solver_args) #, tee=True)
+                result = solver.solve(m, **config.solver_args)  # , tee=True)
 
-                if result.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
+                if result.solution_status in {
+                    SolutionStatus.feasible,
+                    SolutionStatus.optimal,
+                }:
                     result.solution_loader.load_vars()
-                    logger.info(f'solved NLP: {result.solution_status}, {result.termination_condition}')
+                    logger.info(
+                        f'solved NLP: {result.solution_status}, {result.termination_condition}'
+                    )
                     if config.break_on_solution:
                         best_model = m
                         best_result = result
                         break
 
-
                 if result.termination_condition is tc.optimal:
                     model_objectives = m.component_data_objects(Objective, active=True)
                     mobj = next(model_objectives)
@@ -323,8 +333,10 @@ def __enter__(self):
     def __exit__(self, t, v, traceback):
         pass
 
+
 # Sampling class to organize and configure random samplers
 
+
 class SamplingManager:
     def __init__(self, method="uniform", rng=None, seed=None):
         aliases = {
@@ -336,11 +348,11 @@ def __init__(self, method="uniform", rng=None, seed=None):
             "sobol": "sobol",
         }
         self.method = aliases[method.lower()]
-        
+
         self.seed = seed
 
         # Define or create a random number generator
-        # All 
+        # All
 
         if rng is not None:
             self.rng = rng
@@ -370,7 +382,7 @@ def sample_vector(self, lower, upper):
 
         if self.method in ("lhs", "sobol"):
             self._ensure_qmc(dim=len(lower))
-            x = self.qmc_sampler.random(n=1)   # shape (1, d)
+            x = self.qmc_sampler.random(n=1)  # shape (1, d)
             return stats.qmc.scale(x, lower, upper)[0]
 
-        raise ValueError(f"Unknown sampling method '{self.method}'")
\ No newline at end of file
+        raise ValueError(f"Unknown sampling method '{self.method}'")
diff --git a/pyomo/contrib/multistart/reinit.py b/pyomo/contrib/multistart/reinit.py
index 557e49b650b..9eaee477cc5 100644
--- a/pyomo/contrib/multistart/reinit.py
+++ b/pyomo/contrib/multistart/reinit.py
@@ -13,9 +13,7 @@
 import random
 from pyomo.common.dependencies import numpy as np
 from pyomo.common.dependencies.scipy import stats
-from pyomo.core.expr.visitor import (
-    identify_variables,
-)
+from pyomo.core.expr.visitor import identify_variables
 
 from pyomo.core import Var
 
@@ -23,9 +21,10 @@
 
 
 def rand(val, lb, ub, rng):
-    sample = rng.uniform(lb, ub) # uniform distribution between lb and ub
+    sample = rng.uniform(lb, ub)  # uniform distribution between lb and ub
     return sample
 
+
 def midpoint_guess_and_bound(val, lb, ub, rng=None):
     """Midpoint between current value and farthest bound."""
     far_bound = ub if ((ub - val) >= (val - lb)) else lb  # farther bound
@@ -59,7 +58,6 @@ def linspace(lower, upper, n):
     "rand_guess_and_bound": rand_guess_and_bound,
     "rand_distributed": rand_distributed,
     "midpoint": simple_midpoint,
-
 }
 
 
@@ -72,7 +70,7 @@ def reinitialize_variables(model, config, sampler):
 
     eligible_vars = []
 
-    for var in model.component_data_objects(ctype=Var, descend_into=True):     
+    for var in model.component_data_objects(ctype=Var, descend_into=True):
         if var.is_fixed() or not var.is_continuous():
             continue
         if var.lb is None or var.ub is None:
@@ -90,28 +88,29 @@ def reinitialize_variables(model, config, sampler):
     # Sample for new methods as a vector
     if sampler.method in {"uniform", "lhs", "sobol"}:
         if len(eligible_vars) == 0:
-                raise ValueError("No eligible variables to reinitialize." \
-                "Please add bounds.")
-        
+            raise ValueError(
+                "No eligible variables to reinitialize." "Please add bounds."
+            )
+
         # Collect lower and upper bounds for sampler
         lowers = [v.lb for v in eligible_vars]
         uppers = [v.ub for v in eligible_vars]
-        
+
         # Generate vector of samples using sampler
         samples = sampler.sample_vector(lowers, uppers)
 
         # assign samples to variables
         for var, sample in zip(eligible_vars, samples):
-                var.set_value(sample, skip_validation=True)
+            var.set_value(sample, skip_validation=True)
 
         return
-    
+
     # Otherwise use strategies to maintain original functionality
     for var in eligible_vars:
         val = var.value if var.value is not None else (var.lb + var.ub) / 2
         print(f"val = {val}\n")
         # apply reinitialization strategy to variable
         var.set_value(
-            strategies[config.strategy](val, var.lb, var.ub, sampler.rng), skip_validation=True
-
+            strategies[config.strategy](val, var.lb, var.ub, sampler.rng),
+            skip_validation=True,
         )

From ead5c802843c7a510c042767788f5777b914babf Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Mon, 13 Jul 2026 11:08:41 -0600
Subject: [PATCH 22/26] Trying to support old and new solver interfaces, in
 progress.

---
 pyomo/contrib/multistart/multi.py | 139 +++++++++++++++++++++---------
 1 file changed, 97 insertions(+), 42 deletions(-)

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index 513c9f77ab5..f5df1dfec0e 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -173,17 +173,25 @@ def available(self, exception_flag=True):
 
     def license_is_valid(self):
         return True
+    
+    def _get_solver_api(self, use_new):
+        if use_new:
+            from pyomo.contrib.solver.common.factory import SolverFactory
+            from pyomo.contrib.solver.common.results import SolutionStatus
+            return SolverFactory, SolutionStatus, None
+        else:
+            from pyomo.opt import SolverFactory
+            from pyomo.opt import TerminationCondition as tc
+            return SolverFactory, None, tc
+
 
     def solve(self, model, **kwds):
         # initialize keyword args
         config = self.CONFIG(kwds.pop('options', {}))
         config.set_value(kwds)
 
-        # initialize the solver
-        if config.new_solvers_bool == True:
-            from pyomo.contrib.solver.common.factory import SolverFactory
-            from pyomo.contrib.solver.common.results import Results
-            from pyomo.contrib.solver.common.results import SolutionStatus
+        # initialize the solver and get accurate api
+        SolverFactory, SolutionStatus, tc = self._get_solver_api(config.new_solvers_bool)
 
         # Create centralized sampler once
         sampler = SamplingManager(
@@ -228,22 +236,32 @@ def solve(self, model, **kwds):
 
             best_result = result = solver.solve(model, **config.solver_args)
             # Check the solution status before loading variables into the model.
-            if best_result.solution_status in {
-                SolutionStatus.feasible,
-                SolutionStatus.optimal,
-            }:
-                best_result.solution_loader.load_vars()
-                logger.info(
-                    f'solved NLP: {best_result.solution_status}, {best_result.termination_condition}'
-                )
-                # If we are looking for the first feasible solution, then return immediately
-                if config.break_on_solution:
-                    return best_result
+            if config.new_solvers_bool:
+                if result.solution_status in {
+                    SolutionStatus.feasible,
+                    SolutionStatus.optimal,
+                }:
+                    result.solution_loader.load_vars()
+                    logger.info(
+                        f'solved NLP: {result.solution_status}, {result.termination_condition}'
+                    )
+                    
+                if best_result.solution_status is SolutionStatus.optimal:
+                    obj_val = value(obj.expr)
+                    best_objective = obj_val
+                    objectives.append(obj_val)
 
-            if result.termination_condition is tc.optimal:
-                obj_val = value(obj.expr)
-                best_objective = obj_val
-                objectives.append(obj_val)
+            else:
+                if result.termination_condition in {
+                    tc.feasible,
+                    tc.optimal,
+                }:
+                    result.solution_loader.load_vars()
+
+                if result.termination_condition is tc.optimal:
+                    obj_val = value(obj.expr)
+                    best_objective = obj_val
+                    objectives.append(obj_val)
 
             num_iter = 0
             max_iter = config.iterations
@@ -273,29 +291,66 @@ def solve(self, model, **kwds):
                 reinitialize_variables(m, config, sampler)
                 result = solver.solve(m, **config.solver_args)  # , tee=True)
 
-                if result.solution_status in {
-                    SolutionStatus.feasible,
-                    SolutionStatus.optimal,
-                }:
-                    result.solution_loader.load_vars()
-                    logger.info(
-                        f'solved NLP: {result.solution_status}, {result.termination_condition}'
-                    )
-                    if config.break_on_solution:
-                        best_model = m
-                        best_result = result
-                        break
-
-                if result.termination_condition is tc.optimal:
-                    model_objectives = m.component_data_objects(Objective, active=True)
-                    mobj = next(model_objectives)
-                    obj_val = value(mobj.expr)
-                    objectives.append(obj_val)
-                    if obj_val * obj_sign < obj_sign * best_objective:
-                        # objective has improved
+                # Check the solution status before loading variables into the model.
+                if config.new_solvers_bool:
+                    if result.solution_status in {
+                        SolutionStatus.feasible,
+                        SolutionStatus.optimal,
+                    }:
+                        result.solution_loader.load_vars()
+                        logger.info(
+                            f'solved NLP: {result.solution_status}, {result.termination_condition}'
+                        )
+                        # If we are looking for the first feasible solution, then return immediately
+                        if config.break_on_solution:
+                            return best_result
+                        
+                    if best_result.solution_status is SolutionStatus.optimal:
+                        obj_val = value(obj.expr)
                         best_objective = obj_val
-                        best_model = m
-                        best_result = result
+                        objectives.append(obj_val)
+
+                else:
+                    if result.termination_condition in {
+                        tc.feasible,
+                        tc.optimal,
+                    }:
+                        result.solution_loader.load_vars()
+
+                    if result.termination_condition is tc.optimal:
+                        model_objectives = m.component_data_objects(Objective, active=True)
+                        mobj = next(model_objectives)
+                        obj_val = value(mobj.expr)
+                        objectives.append(obj_val)
+                        if obj_val * obj_sign < obj_sign * best_objective:
+                            # objective has improved
+                            best_objective = obj_val
+                            best_model = m
+                            best_result = result
+
+                # if result.solution_status in {
+                #     SolutionStatus.feasible,
+                #     SolutionStatus.optimal,
+                # }:
+                #     result.solution_loader.load_vars()
+                #     logger.info(
+                #         f'solved NLP: {result.solution_status}, {result.termination_condition}'
+                #     )
+                #     if config.break_on_solution:
+                #         best_model = m
+                #         best_result = result
+                #         break
+
+                # if result.termination_condition is tc.optimal:
+                #     model_objectives = m.component_data_objects(Objective, active=True)
+                #     mobj = next(model_objectives)
+                #     obj_val = value(mobj.expr)
+                #     objectives.append(obj_val)
+                #     if obj_val * obj_sign < obj_sign * best_objective:
+                #         # objective has improved
+                #         best_objective = obj_val
+                #         best_model = m
+                #         best_result = result
                 if num_iter == 1:
                     # if it's the first iteration, set the best_model and
                     # best_result regardless of solution status in case the

From 1561d175cf61a5466290871d48f9f34f25cb2168 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Mon, 13 Jul 2026 11:20:11 -0600
Subject: [PATCH 23/26] Ran black

---
 pyomo/contrib/multistart/multi.py | 27 +++++++++++++--------------
 1 file changed, 13 insertions(+), 14 deletions(-)

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index f5df1dfec0e..47274e4a767 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -173,17 +173,18 @@ def available(self, exception_flag=True):
 
     def license_is_valid(self):
         return True
-    
+
     def _get_solver_api(self, use_new):
         if use_new:
             from pyomo.contrib.solver.common.factory import SolverFactory
             from pyomo.contrib.solver.common.results import SolutionStatus
+
             return SolverFactory, SolutionStatus, None
         else:
             from pyomo.opt import SolverFactory
             from pyomo.opt import TerminationCondition as tc
-            return SolverFactory, None, tc
 
+            return SolverFactory, None, tc
 
     def solve(self, model, **kwds):
         # initialize keyword args
@@ -191,7 +192,9 @@ def solve(self, model, **kwds):
         config.set_value(kwds)
 
         # initialize the solver and get accurate api
-        SolverFactory, SolutionStatus, tc = self._get_solver_api(config.new_solvers_bool)
+        SolverFactory, SolutionStatus, tc = self._get_solver_api(
+            config.new_solvers_bool
+        )
 
         # Create centralized sampler once
         sampler = SamplingManager(
@@ -245,17 +248,14 @@ def solve(self, model, **kwds):
                     logger.info(
                         f'solved NLP: {result.solution_status}, {result.termination_condition}'
                     )
-                    
+
                 if best_result.solution_status is SolutionStatus.optimal:
                     obj_val = value(obj.expr)
                     best_objective = obj_val
                     objectives.append(obj_val)
 
             else:
-                if result.termination_condition in {
-                    tc.feasible,
-                    tc.optimal,
-                }:
+                if result.termination_condition in {tc.feasible, tc.optimal}:
                     result.solution_loader.load_vars()
 
                 if result.termination_condition is tc.optimal:
@@ -304,21 +304,20 @@ def solve(self, model, **kwds):
                         # If we are looking for the first feasible solution, then return immediately
                         if config.break_on_solution:
                             return best_result
-                        
+
                     if best_result.solution_status is SolutionStatus.optimal:
                         obj_val = value(obj.expr)
                         best_objective = obj_val
                         objectives.append(obj_val)
 
                 else:
-                    if result.termination_condition in {
-                        tc.feasible,
-                        tc.optimal,
-                    }:
+                    if result.termination_condition in {tc.feasible, tc.optimal}:
                         result.solution_loader.load_vars()
 
                     if result.termination_condition is tc.optimal:
-                        model_objectives = m.component_data_objects(Objective, active=True)
+                        model_objectives = m.component_data_objects(
+                            Objective, active=True
+                        )
                         mobj = next(model_objectives)
                         obj_val = value(mobj.expr)
                         objectives.append(obj_val)

From 207bd130cbcaecb54cca64ff2d805155663708a0 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Mon, 13 Jul 2026 11:24:57 -0600
Subject: [PATCH 24/26] Ran black on devel/initialization

---
 pyomo/devel/initialization/initialize.py      | 13 +++++++------
 pyomo/devel/initialization/multistart_init.py | 13 ++++---------
 2 files changed, 11 insertions(+), 15 deletions(-)

diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py
index c793e14c1c1..e3ace56edc6 100644
--- a/pyomo/devel/initialization/initialize.py
+++ b/pyomo/devel/initialization/initialize.py
@@ -324,11 +324,10 @@ def initialize_with_global_opt(
 def initialize_with_multistart_opt(
     nlp: BlockData,
     nlp_solver: SolverBase | None = None,
-    multistart_solver = None,
+    multistart_solver=None,
     skip_initial_nlp_solve: bool = False,
     default_bound: float = 1e8,
     seed=0,
-
 ) -> Results:
     """
     Attempt to initialize and subsequently solve the model given by ``nlp``.
@@ -362,8 +361,8 @@ def initialize_with_multistart_opt(
 
     if multistart_solver is None:
         multistart_solver = pyo.SolverFactory("multistart")
-        multistart_solver.CONFIG.seed=seed
-        multistart_solver.CONFIG.new_solvers_bool=True
+        multistart_solver.CONFIG.seed = seed
+        multistart_solver.CONFIG.new_solvers_bool = True
 
     if not skip_initial_nlp_solve:
         res = _try_nlp_solve(nlp, nlp_solver)
@@ -374,8 +373,10 @@ def initialize_with_multistart_opt(
 
     try:
         res = _initialize_with_multistart_solver(
-            nlp=nlp, multistart_solver=multistart_solver,
-            default_bound=default_bound, seed=seed
+            nlp=nlp,
+            multistart_solver=multistart_solver,
+            default_bound=default_bound,
+            seed=seed,
         )
     finally:
         _cleanup(orig_var_data)
diff --git a/pyomo/devel/initialization/multistart_init.py b/pyomo/devel/initialization/multistart_init.py
index 394decb6cc1..46ea2ff7e80 100644
--- a/pyomo/devel/initialization/multistart_init.py
+++ b/pyomo/devel/initialization/multistart_init.py
@@ -22,20 +22,15 @@
 
 
 def _initialize_with_multistart_solver(
-    nlp: BlockData, 
-    multistart_solver,
-    default_bound=1.0e8,
-    seed = None,
-    ):
-    
+    nlp: BlockData, multistart_solver, default_bound=1.0e8, seed=None
+):
+
     # Make a shallow clone
     nlp = shallow_clone(nlp)
     # bounds on the nonlinear variables
     bound_all_nonlinear_variables(nlp, default_bound=default_bound)
 
     res = multistart_solver.solve(nlp)
-    logger.info(
-        'Finished multistart optimization iterations.'
-        )
+    logger.info('Finished multistart optimization iterations.')
 
     return res

From 9e734b1c96425879dde4b6729f7a9aeb60196dfc Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Mon, 20 Jul 2026 14:33:55 -0600
Subject: [PATCH 25/26] Switching over to new solver factory interface, ran
 black

---
 pyomo/contrib/multistart/multi.py            | 143 +++++--------------
 pyomo/contrib/multistart/tests/test_multi.py |  25 ++--
 2 files changed, 50 insertions(+), 118 deletions(-)

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index 47274e4a767..55ac5906b18 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -22,7 +22,8 @@
 from pyomo.contrib.multistart.reinit import reinitialize_variables, strategies
 from pyomo.core import Objective, Var, minimize, value
 from pyomo.opt import SolverFactory, SolverStatus
-from pyomo.opt import TerminationCondition as tc
+from pyomo.contrib.solver.common.factory import SolverFactory as NewSolverFactory
+from pyomo.contrib.solver.common.results import SolutionStatus
 from pyomo.common.dependencies.scipy import stats
 from pyomo.common.dependencies import numpy as np
 
@@ -153,14 +154,6 @@ class MultiStart:
                 Preferred over seed.",
         ),
     )
-    CONFIG.declare(
-        "new_solvers_bool",
-        ConfigValue(
-            default=False,
-            description="Boolean option for whether to use the new solver interface, default to no \
-                until solver testing complete (?)",
-        ),
-    )
 
     def available(self, exception_flag=True):
         """Check if solver is available.
@@ -174,34 +167,17 @@ def available(self, exception_flag=True):
     def license_is_valid(self):
         return True
 
-    def _get_solver_api(self, use_new):
-        if use_new:
-            from pyomo.contrib.solver.common.factory import SolverFactory
-            from pyomo.contrib.solver.common.results import SolutionStatus
-
-            return SolverFactory, SolutionStatus, None
-        else:
-            from pyomo.opt import SolverFactory
-            from pyomo.opt import TerminationCondition as tc
-
-            return SolverFactory, None, tc
-
     def solve(self, model, **kwds):
         # initialize keyword args
         config = self.CONFIG(kwds.pop('options', {}))
         config.set_value(kwds)
 
-        # initialize the solver and get accurate api
-        SolverFactory, SolutionStatus, tc = self._get_solver_api(
-            config.new_solvers_bool
-        )
-
         # Create centralized sampler once
         sampler = SamplingManager(
             method=config.sampling_method, rng=config.rng, seed=config.seed
         )
 
-        solver = SolverFactory(config.solver)
+        solver = NewSolverFactory(config.solver)
 
         # Model sense
         objectives = model.component_data_objects(Objective, active=True)
@@ -239,29 +215,19 @@ def solve(self, model, **kwds):
 
             best_result = result = solver.solve(model, **config.solver_args)
             # Check the solution status before loading variables into the model.
-            if config.new_solvers_bool:
-                if result.solution_status in {
-                    SolutionStatus.feasible,
-                    SolutionStatus.optimal,
-                }:
-                    result.solution_loader.load_vars()
-                    logger.info(
-                        f'solved NLP: {result.solution_status}, {result.termination_condition}'
-                    )
-
-                if best_result.solution_status is SolutionStatus.optimal:
-                    obj_val = value(obj.expr)
-                    best_objective = obj_val
-                    objectives.append(obj_val)
-
-            else:
-                if result.termination_condition in {tc.feasible, tc.optimal}:
-                    result.solution_loader.load_vars()
+            if result.solution_status in {
+                SolutionStatus.feasible,
+                SolutionStatus.optimal,
+            }:
+                result.solution_loader.load_vars()
+                logger.info(
+                    f'solved NLP: {result.solution_status}, {result.termination_condition}'
+                )
 
-                if result.termination_condition is tc.optimal:
-                    obj_val = value(obj.expr)
-                    best_objective = obj_val
-                    objectives.append(obj_val)
+            if best_result.solution_status is SolutionStatus.optimal:
+                obj_val = value(obj.expr)
+                best_objective = obj_val
+                objectives.append(obj_val)
 
             num_iter = 0
             max_iter = config.iterations
@@ -292,64 +258,29 @@ def solve(self, model, **kwds):
                 result = solver.solve(m, **config.solver_args)  # , tee=True)
 
                 # Check the solution status before loading variables into the model.
-                if config.new_solvers_bool:
-                    if result.solution_status in {
-                        SolutionStatus.feasible,
-                        SolutionStatus.optimal,
-                    }:
-                        result.solution_loader.load_vars()
-                        logger.info(
-                            f'solved NLP: {result.solution_status}, {result.termination_condition}'
-                        )
-                        # If we are looking for the first feasible solution, then return immediately
-                        if config.break_on_solution:
-                            return best_result
-
-                    if best_result.solution_status is SolutionStatus.optimal:
-                        obj_val = value(obj.expr)
+                if result.solution_status in {
+                    SolutionStatus.feasible,
+                    SolutionStatus.optimal,
+                }:
+                    result.solution_loader.load_vars()
+                    logger.info(
+                        f'solved NLP: {result.solution_status}, {result.termination_condition}'
+                    )
+                    # If we are looking for the first feasible solution, then return immediately
+                    if config.break_on_solution:
+                        return best_result
+
+                if best_result.solution_status is SolutionStatus.optimal:
+                    model_objectives = m.component_data_objects(Objective, active=True)
+                    mobj = next(model_objectives)
+                    obj_val = value(mobj.expr)
+                    objectives.append(obj_val)
+                    if obj_val * obj_sign < obj_sign * best_objective:
+                        # objective has improved
                         best_objective = obj_val
-                        objectives.append(obj_val)
-
-                else:
-                    if result.termination_condition in {tc.feasible, tc.optimal}:
-                        result.solution_loader.load_vars()
-
-                    if result.termination_condition is tc.optimal:
-                        model_objectives = m.component_data_objects(
-                            Objective, active=True
-                        )
-                        mobj = next(model_objectives)
-                        obj_val = value(mobj.expr)
-                        objectives.append(obj_val)
-                        if obj_val * obj_sign < obj_sign * best_objective:
-                            # objective has improved
-                            best_objective = obj_val
-                            best_model = m
-                            best_result = result
-
-                # if result.solution_status in {
-                #     SolutionStatus.feasible,
-                #     SolutionStatus.optimal,
-                # }:
-                #     result.solution_loader.load_vars()
-                #     logger.info(
-                #         f'solved NLP: {result.solution_status}, {result.termination_condition}'
-                #     )
-                #     if config.break_on_solution:
-                #         best_model = m
-                #         best_result = result
-                #         break
-
-                # if result.termination_condition is tc.optimal:
-                #     model_objectives = m.component_data_objects(Objective, active=True)
-                #     mobj = next(model_objectives)
-                #     obj_val = value(mobj.expr)
-                #     objectives.append(obj_val)
-                #     if obj_val * obj_sign < obj_sign * best_objective:
-                #         # objective has improved
-                #         best_objective = obj_val
-                #         best_model = m
-                #         best_result = result
+                        best_model = m
+                        best_result = result
+
                 if num_iter == 1:
                     # if it's the first iteration, set the best_model and
                     # best_result regardless of solution status in case the
diff --git a/pyomo/contrib/multistart/tests/test_multi.py b/pyomo/contrib/multistart/tests/test_multi.py
index 1c34138fbb3..cf58d8b14bd 100644
--- a/pyomo/contrib/multistart/tests/test_multi.py
+++ b/pyomo/contrib/multistart/tests/test_multi.py
@@ -134,18 +134,19 @@ def test_multiple_obj(self):
         with self.assertRaisesRegex(RuntimeError, "multiple active objectives"):
             SolverFactory('multistart').solve(m)
 
-    def test_no_obj(self):
-        m = ConcreteModel()
-        m.x = Var()
-        with self.assertRaisesRegex(RuntimeError, "no active objective"):
-            SolverFactory('multistart').solve(m)
-
-    def test_const_obj(self):
-        m = ConcreteModel()
-        m.x = Var()
-        m.o = Objective(expr=5)
-        with self.assertRaisesRegex(RuntimeError, "constant objective"):
-            SolverFactory('multistart').solve(m)
+    # Would like to remove these tests to allow for square model solves.
+    # def test_no_obj(self):
+    #     m = ConcreteModel()
+    #     m.x = Var()
+    #     with self.assertRaisesRegex(RuntimeError, "no active objective"):
+    #         SolverFactory('multistart').solve(m)
+
+    # def test_const_obj(self):
+    #     m = ConcreteModel()
+    #     m.x = Var()
+    #     m.o = Objective(expr=5)
+    #     with self.assertRaisesRegex(RuntimeError, "constant objective"):
+    #         SolverFactory('multistart').solve(m)
 
 
 def build_model():

From 331d8c8301fd81499a900990997210df84a3c50e Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Mon, 20 Jul 2026 14:50:34 -0600
Subject: [PATCH 26/26] Added new initialize bool. Ran black

---
 pyomo/contrib/multistart/multi.py | 58 ++++++++++++++++++++-----------
 1 file changed, 38 insertions(+), 20 deletions(-)

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index 55ac5906b18..3c05f7d6049 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -63,7 +63,11 @@ class MultiStart:
     )
     CONFIG.declare(
         "solver",
-        ConfigValue(default="ipopt", description="solver to use, defaults to ipopt"),
+        ConfigValue(
+            default="ipopt",
+            description="solver to use, defaults to ipopt"
+            "Should also be able to accept solver objects. In progress",
+        ),
     )
     CONFIG.declare(
         "solver_args",
@@ -155,6 +159,14 @@ class MultiStart:
         ),
     )
 
+    CONFIG.declare(
+        "initialize",
+        ConfigValue(
+            default=False,
+            description="Boolean for whether solver is being used to initialize model. Default is False.",
+        ),
+    )
+
     def available(self, exception_flag=True):
         """Check if solver is available.
 
@@ -177,6 +189,10 @@ def solve(self, model, **kwds):
             method=config.sampling_method, rng=config.rng, seed=config.seed
         )
 
+        if config.initialize == True:
+            config.solver_args["load_solutions"] = False
+            config.solver_args["raise_exception_on_nonoptimal_result"] = False
+
         solver = NewSolverFactory(config.solver)
 
         # Model sense
@@ -215,14 +231,15 @@ def solve(self, model, **kwds):
 
             best_result = result = solver.solve(model, **config.solver_args)
             # Check the solution status before loading variables into the model.
-            if result.solution_status in {
-                SolutionStatus.feasible,
-                SolutionStatus.optimal,
-            }:
-                result.solution_loader.load_vars()
-                logger.info(
-                    f'solved NLP: {result.solution_status}, {result.termination_condition}'
-                )
+            if config.initialize:
+                if result.solution_status in {
+                    SolutionStatus.feasible,
+                    SolutionStatus.optimal,
+                }:
+                    result.solution_loader.load_vars()
+                    logger.info(
+                        f'solved NLP: {result.solution_status}, {result.termination_condition}'
+                    )
 
             if best_result.solution_status is SolutionStatus.optimal:
                 obj_val = value(obj.expr)
@@ -258,17 +275,18 @@ def solve(self, model, **kwds):
                 result = solver.solve(m, **config.solver_args)  # , tee=True)
 
                 # Check the solution status before loading variables into the model.
-                if result.solution_status in {
-                    SolutionStatus.feasible,
-                    SolutionStatus.optimal,
-                }:
-                    result.solution_loader.load_vars()
-                    logger.info(
-                        f'solved NLP: {result.solution_status}, {result.termination_condition}'
-                    )
-                    # If we are looking for the first feasible solution, then return immediately
-                    if config.break_on_solution:
-                        return best_result
+                if config.initialize:
+                    if result.solution_status in {
+                        SolutionStatus.feasible,
+                        SolutionStatus.optimal,
+                    }:
+                        result.solution_loader.load_vars()
+                        logger.info(
+                            f'solved NLP: {result.solution_status}, {result.termination_condition}'
+                        )
+                        # If we are looking for the first feasible solution, then return immediately
+                        if config.break_on_solution:
+                            return best_result
 
                 if best_result.solution_status is SolutionStatus.optimal:
                     model_objectives = m.component_data_objects(Objective, active=True)