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/55] 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/55] 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/55] 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/55] 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/55] 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/55] 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/55] 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/55] 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/55] 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/55] 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/55] 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/55] 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/55] 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/55] 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/55] 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/55] 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/55] 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/55] 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/55] 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/55] 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/55] 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/55] 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/55] 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/55] 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/55] 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/55] 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)

From 99d738a744b754b5eee3166b05daee9257019dd8 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Tue, 21 Jul 2026 08:34:37 -0600
Subject: [PATCH 27/55] Fixed SolverFactory issue

---
 pyomo/contrib/multistart/multi.py            | 5 ++---
 pyomo/contrib/multistart/tests/test_multi.py | 2 +-
 2 files changed, 3 insertions(+), 4 deletions(-)

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index 3c05f7d6049..189cf29cce8 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -21,8 +21,7 @@
 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
-from pyomo.opt import SolverFactory, SolverStatus
-from pyomo.contrib.solver.common.factory import SolverFactory as NewSolverFactory
+from pyomo.contrib.solver.common.factory import SolverFactory
 from pyomo.contrib.solver.common.results import SolutionStatus
 from pyomo.common.dependencies.scipy import stats
 from pyomo.common.dependencies import numpy as np
@@ -193,7 +192,7 @@ def solve(self, model, **kwds):
             config.solver_args["load_solutions"] = False
             config.solver_args["raise_exception_on_nonoptimal_result"] = False
 
-        solver = NewSolverFactory(config.solver)
+        solver = SolverFactory(config.solver)
 
         # Model sense
         objectives = model.component_data_objects(Objective, active=True)
diff --git a/pyomo/contrib/multistart/tests/test_multi.py b/pyomo/contrib/multistart/tests/test_multi.py
index cf58d8b14bd..a96562ae195 100644
--- a/pyomo/contrib/multistart/tests/test_multi.py
+++ b/pyomo/contrib/multistart/tests/test_multi.py
@@ -21,12 +21,12 @@
     Constraint,
     NonNegativeReals,
     Objective,
-    SolverFactory,
     Var,
     maximize,
     sin,
     value,
 )
+from pyomo.contrib.solver.common.factory import SolverFactory
 
 
 @unittest.skipIf(not SolverFactory('ipopt').available(), "IPOPT not available")

From db10cfb77144ecada5959eafc75ef5bde42ba78a Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Tue, 21 Jul 2026 11:25:22 -0600
Subject: [PATCH 28/55] Added vectorized strategy, adjusted defaults, ran
 black.

---
 pyomo/contrib/multistart/multi.py  | 69 +++++++++++-------------------
 pyomo/contrib/multistart/reinit.py | 41 +++++++++++-------
 2 files changed, 52 insertions(+), 58 deletions(-)

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index 189cf29cce8..27c30025d2a 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -53,6 +53,7 @@ class MultiStart:
             doc="""Specify the restart strategy.
 
         - "rand": random choice between variable bounds
+        - "rand_vector": random choice, vectorized approach with sampler
         - "midpoint_guess_and_bound": midpoint between current value and farthest bound
         - "rand_guess_and_bound": random choice between current value and farthest bound
         - "rand_distributed": random choice among evenly distributed values
@@ -138,8 +139,9 @@ class MultiStart:
         "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'",
+            description="Method for sampling random starting points for reinitialization step. "
+            "Supported options are 'random_uniform', 'latin_hypercube', and 'sobol_sampling'. "
+            "Only utilized when config.strategy is 'rand_vector'.",
         ),
     )
     CONFIG.declare(
@@ -158,14 +160,6 @@ 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.
 
@@ -188,9 +182,9 @@ 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
+        # Set options so infeasible solve does not interrupt runs
+        config.solver_args["load_solutions"] = False
+        config.solver_args["raise_exception_on_nonoptimal_result"] = False
 
         solver = SolverFactory(config.solver)
 
@@ -230,15 +224,14 @@ 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.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 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)
@@ -274,18 +267,17 @@ 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.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 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)
@@ -298,13 +290,6 @@ def solve(self, model, **kwds):
                         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
-                    # model is infeasible.
-                    best_model = m
-                    best_result = result
-
             if using_HCS and not HCS_completed:
                 logger.warning(
                     "High confidence stopping rule was unable to complete "
@@ -354,8 +339,6 @@ def __init__(self, method="uniform", rng=None, seed=None):
         self.seed = seed
 
         # Define or create a random number generator
-        # All
-
         if rng is not None:
             self.rng = rng
         else:
diff --git a/pyomo/contrib/multistart/reinit.py b/pyomo/contrib/multistart/reinit.py
index 9eaee477cc5..856c3948deb 100644
--- a/pyomo/contrib/multistart/reinit.py
+++ b/pyomo/contrib/multistart/reinit.py
@@ -25,6 +25,14 @@ def rand(val, lb, ub, rng):
     return sample
 
 
+def rand_vector(lbs, ubs, sampler):
+    lowers = lbs
+    uppers = ubs
+    # Generate vector of samples using sampler
+    samples = sampler.sample_vector(lowers, uppers)
+    return samples
+
+
 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
@@ -34,7 +42,10 @@ def midpoint_guess_and_bound(val, lb, ub, rng=None):
 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 rng.uniform(val, far_bound)
+    if far_bound == ub:
+        return rng.uniform(val, far_bound)
+    else:
+        return rng.uniform(far_bound, val)
 
 
 def rand_distributed(val, lb, ub, rng, divisions=9):
@@ -54,6 +65,7 @@ def linspace(lower, upper, n):
 
 strategies = {
     "rand": rand,
+    "rand_vector": rand_vector,
     "midpoint_guess_and_bound": midpoint_guess_and_bound,
     "rand_guess_and_bound": rand_guess_and_bound,
     "rand_distributed": rand_distributed,
@@ -85,19 +97,17 @@ def reinitialize_variables(model, config, sampler):
 
         eligible_vars.append(var)
 
-    # Sample for new methods as a vector
-    if sampler.method in {"uniform", "lhs", "sobol"}:
+    if config.strategy == "rand_vector":
+
         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)
+        samples = rand_vector(lowers, uppers, sampler)
 
         # assign samples to variables
         for var, sample in zip(eligible_vars, samples):
@@ -105,12 +115,13 @@ def reinitialize_variables(model, config, sampler):
 
         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,
-        )
+    # Otherwise
+    else:
+        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,
+            )

From 93e1e7e008b0f9dc8ea471c6c973e8318eb27469 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Tue, 21 Jul 2026 13:35:13 -0600
Subject: [PATCH 29/55] Changed configs

---
 pyomo/contrib/multistart/multi.py                      |  2 --
 .../initialization/examples/init_polynomial_ex.py      | 10 ++++++++++
 pyomo/devel/initialization/initialize.py               |  6 ++++--
 pyomo/devel/initialization/multistart_init.py          |  2 +-
 4 files changed, 15 insertions(+), 5 deletions(-)

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index 27c30025d2a..5bd81d5420a 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -322,8 +322,6 @@ def __exit__(self, t, v, traceback):
 
 
 # Sampling class to organize and configure random samplers
-
-
 class SamplingManager:
     def __init__(self, method="uniform", rng=None, seed=None):
         aliases = {
diff --git a/pyomo/devel/initialization/examples/init_polynomial_ex.py b/pyomo/devel/initialization/examples/init_polynomial_ex.py
index 643c2faf180..45010f00e7b 100644
--- a/pyomo/devel/initialization/examples/init_polynomial_ex.py
+++ b/pyomo/devel/initialization/examples/init_polynomial_ex.py
@@ -53,6 +53,16 @@ def global_init_ex():
     return results.solution_status, m.x.value
 
 
+def multistart_init_ex():
+    m = build_model()
+    nlp_solver = SolverFactory('ipopt')
+    global_solver = SolverFactory('scip_direct')
+    results = ini.initialize_with_global_opt(
+        nlp=m, nlp_solver=nlp_solver, global_solver=global_solver
+    )
+
+    return results.solution_status, m.x.value
+
 if __name__ == '__main__':
     # stat, x = lp_init_ex()
     # stat, x = pwl_init_ex()
diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py
index e3ace56edc6..ca82acb7c81 100644
--- a/pyomo/devel/initialization/initialize.py
+++ b/pyomo/devel/initialization/initialize.py
@@ -362,7 +362,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.sampling_method = "lhs"
+        multistart_solver.CONFIG.break_on_solution=True
 
     if not skip_initial_nlp_solve:
         res = _try_nlp_solve(nlp, nlp_solver)
@@ -375,8 +376,9 @@ def initialize_with_multistart_opt(
         res = _initialize_with_multistart_solver(
             nlp=nlp,
             multistart_solver=multistart_solver,
-            default_bound=default_bound,
             seed=seed,
+            default_bound=default_bound,
+
         )
     finally:
         _cleanup(orig_var_data)
diff --git a/pyomo/devel/initialization/multistart_init.py b/pyomo/devel/initialization/multistart_init.py
index 46ea2ff7e80..f3b5624c9a2 100644
--- a/pyomo/devel/initialization/multistart_init.py
+++ b/pyomo/devel/initialization/multistart_init.py
@@ -22,7 +22,7 @@
 
 
 def _initialize_with_multistart_solver(
-    nlp: BlockData, multistart_solver, default_bound=1.0e8, seed=None
+    nlp: BlockData, multistart_solver, seed, default_bound=1.0e8, 
 ):
 
     # Make a shallow clone

From e4f564c4317987c5ebf8978333cd703111efd631 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Thu, 23 Jul 2026 14:04:52 -0600
Subject: [PATCH 30/55] Updated solverwrapper to meet new requirements,  ran
 black

---
 pyomo/contrib/multistart/multi.py      | 270 ++++++++++++++-----------
 pyomo/devel/initialization/__init__.py |   1 +
 2 files changed, 153 insertions(+), 118 deletions(-)

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index 5bd81d5420a..8f3c6957ea0 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -12,26 +12,169 @@
 
 from pyomo.common.config import (
     ConfigBlock,
+    ConfigDict,
     ConfigValue,
     In,
     document_kwargs_from_configdict,
+    document_class_CONFIG,
+    document_configdict,
+    ADVANCED_OPTION,
 )
+
 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
+from pyomo.contrib.solver.common.base import SolverBase
+from pyomo.contrib.solver.common.config import SolverConfig
 from pyomo.contrib.solver.common.factory import SolverFactory
 from pyomo.contrib.solver.common.results import SolutionStatus
+
 from pyomo.common.dependencies.scipy import stats
 from pyomo.common.dependencies import numpy as np
 
 logger = logging.getLogger('pyomo.contrib.multistart')
 
 
+@document_configdict()
+class MultistartConfig(SolverConfig):
+    def __init__(
+        self,
+        description=None,
+        doc=None,
+        implicit=False,
+        implicit_domain=None,
+        visibility=0,
+    ):
+        super().__init__(
+            description=description,
+            doc=doc,
+            implicit=implicit,
+            implicit_domain=implicit_domain,
+            visibility=visibility,
+        )
+
+        self.strategy = self.declare(
+            "strategy",
+            ConfigValue(
+                default="rand",
+                domain=In(strategies.keys()),
+                description="Specify the restart strategy. Defaults to rand.",
+                doc="""Specify the restart strategy.
+
+            - "rand": random choice between variable bounds
+            - "rand_vector": random choice, vectorized approach with sampler
+            - "midpoint_guess_and_bound": midpoint between current value and farthest bound
+            - "rand_guess_and_bound": random choice between current value and farthest bound
+            - "rand_distributed": random choice among evenly distributed values
+            - "midpoint": exact midpoint between the bounds. If using this option, multiple iterations are useless.
+            """,
+            ),
+        )
+        self.solver = self.declare(
+            "solver",
+            ConfigValue(
+                default="ipopt",
+                description="solver to use, defaults to ipopt"
+                "Should also be able to accept solver objects. In progress",
+            ),
+        )
+        self.solver_args = self.declare(
+            "solver_args",
+            ConfigValue(
+                default={},
+                description="Dictionary of keyword arguments to pass to the solver.",
+            ),
+        )
+        self.iterations = self.declare(
+            "iterations",
+            ConfigValue(
+                default=10,
+                description="Specify the number of iterations, defaults to 10. "
+                "If -1 is specified, the high confidence stopping rule will be used",
+            ),
+        )
+        self.stopping_mass = self.declare(
+            "stopping_mass",
+            ConfigValue(
+                default=0.5,
+                description="Maximum allowable estimated missing mass of optima.",
+                doc="""Maximum allowable estimated missing mass of optima for the
+            high confidence stopping rule, only used with the random strategy.
+            The lower the parameter, the stricter the rule.
+            Value bounded in (0, 1].""",
+            ),
+        )
+        self.stopping_delta = self.declare(
+            "stopping_delta",
+            ConfigValue(
+                default=0.5,
+                description="1 minus the confidence level required for the stopping rule.",
+                doc="""1 minus the confidence level required for the stopping rule for the
+            high confidence stopping rule, only used with the random strategy.
+            The lower the parameter, the stricter the rule.
+            Value bounded in (0, 1].""",
+            ),
+        )
+        # self.surpress_unbounded_warning = self.declare(
+        #     "suppress_unbounded_warning",
+        #     ConfigValue(
+        #         default=False,
+        #         domain=bool,
+        #         description="True to suppress warning for skipping unbounded variables.",
+        #     ),
+        # )
+        self.HCS_max_iterations = self.declare(
+            "HCS_max_iterations",
+            ConfigValue(
+                default=1000,
+                description="Maximum number of iterations before interrupting the high confidence stopping rule.",
+            ),
+        )
+        self.HCS_tolerance = self.declare(
+            "HCS_tolerance",
+            ConfigValue(
+                default=0,
+                description="Tolerance on HCS objective value equality. Defaults to Python float equality precision.",
+            ),
+        )
+        self.break_on_solution = self.declare(
+            "break_on_solution",
+            ConfigValue(
+                default=False,
+                description="Condition to break if a feasible or optimal solution is found. Defaults to False.",
+            ),
+        )
+        self.sampling_method = self.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'. "
+                "Only utilized when config.strategy is 'rand_vector'.",
+            ),
+        )
+        self.seed = self.declare(
+            "seed",
+            ConfigValue(
+                default=None,
+                description="Seed for reproducibility in random sampling methods.",
+            ),
+        )
+        self.rng = self.declare(
+            "rng",
+            ConfigValue(
+                default=None,
+                description="Random number generator for reproducibility in random sampling methods. \
+                    Preferred over seed.",
+            ),
+        )
+
+
 @SolverFactory.register('multistart', doc='MultiStart solver for NLPs')
-@document_kwargs_from_configdict('CONFIG')
-class MultiStart:
+@document_class_CONFIG(methods=['solve'])
+class MultiStart(SolverBase):
     """Solver wrapper that initializes at multiple starting points.
 
     # TODO: also return appropriate duals
@@ -43,122 +186,7 @@ class MultiStart:
 
     """
 
-    CONFIG = ConfigBlock("MultiStart")
-    CONFIG.declare(
-        "strategy",
-        ConfigValue(
-            default="rand",
-            domain=In(strategies.keys()),
-            description="Specify the restart strategy. Defaults to rand.",
-            doc="""Specify the restart strategy.
-
-        - "rand": random choice between variable bounds
-        - "rand_vector": random choice, vectorized approach with sampler
-        - "midpoint_guess_and_bound": midpoint between current value and farthest bound
-        - "rand_guess_and_bound": random choice between current value and farthest bound
-        - "rand_distributed": random choice among evenly distributed values
-        - "midpoint": exact midpoint between the bounds. If using this option, multiple iterations are useless.
-        """,
-        ),
-    )
-    CONFIG.declare(
-        "solver",
-        ConfigValue(
-            default="ipopt",
-            description="solver to use, defaults to ipopt"
-            "Should also be able to accept solver objects. In progress",
-        ),
-    )
-    CONFIG.declare(
-        "solver_args",
-        ConfigValue(
-            default={},
-            description="Dictionary of keyword arguments to pass to the solver.",
-        ),
-    )
-    CONFIG.declare(
-        "iterations",
-        ConfigValue(
-            default=10,
-            description="Specify the number of iterations, defaults to 10. "
-            "If -1 is specified, the high confidence stopping rule will be used",
-        ),
-    )
-    CONFIG.declare(
-        "stopping_mass",
-        ConfigValue(
-            default=0.5,
-            description="Maximum allowable estimated missing mass of optima.",
-            doc="""Maximum allowable estimated missing mass of optima for the
-        high confidence stopping rule, only used with the random strategy.
-        The lower the parameter, the stricter the rule.
-        Value bounded in (0, 1].""",
-        ),
-    )
-    CONFIG.declare(
-        "stopping_delta",
-        ConfigValue(
-            default=0.5,
-            description="1 minus the confidence level required for the stopping rule.",
-            doc="""1 minus the confidence level required for the stopping rule for the
-        high confidence stopping rule, only used with the random strategy.
-        The lower the parameter, the stricter the rule.
-        Value bounded in (0, 1].""",
-        ),
-    )
-    CONFIG.declare(
-        "suppress_unbounded_warning",
-        ConfigValue(
-            default=False,
-            domain=bool,
-            description="True to suppress warning for skipping unbounded variables.",
-        ),
-    )
-    CONFIG.declare(
-        "HCS_max_iterations",
-        ConfigValue(
-            default=1000,
-            description="Maximum number of iterations before interrupting the high confidence stopping rule.",
-        ),
-    )
-    CONFIG.declare(
-        "HCS_tolerance",
-        ConfigValue(
-            default=0,
-            description="Tolerance on HCS objective value equality. Defaults to Python float equality precision.",
-        ),
-    )
-    CONFIG.declare(
-        "break_on_solution",
-        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'. "
-            "Only utilized when config.strategy is 'rand_vector'.",
-        ),
-    )
-    CONFIG.declare(
-        "seed",
-        ConfigValue(
-            default=None,
-            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.",
-        ),
-    )
+    CONFIG = MultistartConfig()
 
     def available(self, exception_flag=True):
         """Check if solver is available.
@@ -169,6 +197,12 @@ def available(self, exception_flag=True):
         """
         return True
 
+    def version(self):
+        """Get solver version
+        TODO: This is a solver wrapper, unsure how to define version in this case."""
+
+        return
+
     def license_is_valid(self):
         return True
 
diff --git a/pyomo/devel/initialization/__init__.py b/pyomo/devel/initialization/__init__.py
index 49da58c9d5e..d3796e90c3f 100644
--- a/pyomo/devel/initialization/__init__.py
+++ b/pyomo/devel/initialization/__init__.py
@@ -11,4 +11,5 @@
     initialize_with_LP_approximation,
     initialize_with_piecewise_linear_approximation,
     initialize_with_global_opt,
+    initialize_with_multistart_opt,
 )

From 4a80136b7aed13a5514443ccbf748a1351f048b5 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Thu, 23 Jul 2026 14:07:13 -0600
Subject: [PATCH 31/55] Add back suppression warning

---
 pyomo/contrib/multistart/multi.py | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index 8f3c6957ea0..1c14baa1b58 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -117,14 +117,14 @@ def __init__(
             Value bounded in (0, 1].""",
             ),
         )
-        # self.surpress_unbounded_warning = self.declare(
-        #     "suppress_unbounded_warning",
-        #     ConfigValue(
-        #         default=False,
-        #         domain=bool,
-        #         description="True to suppress warning for skipping unbounded variables.",
-        #     ),
-        # )
+        self.suppress_unbounded_warning = self.declare(
+            "suppress_unbounded_warning",
+            ConfigValue(
+                default=False,
+                domain=bool,
+                description="True to suppress warning for skipping unbounded variables.",
+            ),
+        )
         self.HCS_max_iterations = self.declare(
             "HCS_max_iterations",
             ConfigValue(

From 99a51f1df41afa02bf7bc561a2bb609467980fe4 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Thu, 23 Jul 2026 14:20:39 -0600
Subject: [PATCH 32/55] Ran black in devel

---
 pyomo/devel/initialization/examples/init_polynomial_ex.py | 1 +
 pyomo/devel/initialization/initialize.py                  | 3 +--
 pyomo/devel/initialization/multistart_init.py             | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/pyomo/devel/initialization/examples/init_polynomial_ex.py b/pyomo/devel/initialization/examples/init_polynomial_ex.py
index 45010f00e7b..0cf6b178ae6 100644
--- a/pyomo/devel/initialization/examples/init_polynomial_ex.py
+++ b/pyomo/devel/initialization/examples/init_polynomial_ex.py
@@ -63,6 +63,7 @@ def multistart_init_ex():
 
     return results.solution_status, m.x.value
 
+
 if __name__ == '__main__':
     # stat, x = lp_init_ex()
     # stat, x = pwl_init_ex()
diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py
index ca82acb7c81..6e24d5281ed 100644
--- a/pyomo/devel/initialization/initialize.py
+++ b/pyomo/devel/initialization/initialize.py
@@ -363,7 +363,7 @@ def initialize_with_multistart_opt(
         multistart_solver = pyo.SolverFactory("multistart")
         multistart_solver.CONFIG.seed = seed
         multistart_solver.CONFIG.sampling_method = "lhs"
-        multistart_solver.CONFIG.break_on_solution=True
+        multistart_solver.CONFIG.break_on_solution = True
 
     if not skip_initial_nlp_solve:
         res = _try_nlp_solve(nlp, nlp_solver)
@@ -378,7 +378,6 @@ def initialize_with_multistart_opt(
             multistart_solver=multistart_solver,
             seed=seed,
             default_bound=default_bound,
-
         )
     finally:
         _cleanup(orig_var_data)
diff --git a/pyomo/devel/initialization/multistart_init.py b/pyomo/devel/initialization/multistart_init.py
index f3b5624c9a2..f8addbe52b2 100644
--- a/pyomo/devel/initialization/multistart_init.py
+++ b/pyomo/devel/initialization/multistart_init.py
@@ -22,7 +22,7 @@
 
 
 def _initialize_with_multistart_solver(
-    nlp: BlockData, multistart_solver, seed, default_bound=1.0e8, 
+    nlp: BlockData, multistart_solver, seed, default_bound=1.0e8
 ):
 
     # Make a shallow clone

From 003801f72f3bafa1ad44b7104e8905ae4d59b0fd Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Thu, 23 Jul 2026 14:38:28 -0600
Subject: [PATCH 33/55] Added multistart to solver list in contrib.solvers

---
 pyomo/contrib/solver/tests/solvers/test_solvers.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/pyomo/contrib/solver/tests/solvers/test_solvers.py b/pyomo/contrib/solver/tests/solvers/test_solvers.py
index 0b30a6eb923..a7fddb6356d 100644
--- a/pyomo/contrib/solver/tests/solvers/test_solvers.py
+++ b/pyomo/contrib/solver/tests/solvers/test_solvers.py
@@ -48,6 +48,8 @@
 from pyomo.core.expr.compare import assertExpressionsEqual
 from pyomo.core.expr.numeric_expr import LinearExpression
 
+from pyomo.contrib.multistart.multi import MultiStart
+
 np, numpy_available = attempt_import('numpy')
 parameterized, param_available = attempt_import('parameterized')
 parameterized = parameterized.parameterized
@@ -80,6 +82,8 @@ def param_as_standalone_func(cls, p, func, name):
     ('scip_persistent', ScipPersistent),
     ('gams', GAMS),
     ('knitro_direct', KnitroDirectSolver),
+    ('multistart', MultiStart)
+
 ]
 mip_solvers = [
     ('gurobi_persistent', GurobiPersistent),

From c5e2109213f7571b457072d9de6708e1984adf49 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Thu, 23 Jul 2026 14:44:18 -0600
Subject: [PATCH 34/55] Ran black

---
 pyomo/contrib/solver/tests/solvers/test_solvers.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/pyomo/contrib/solver/tests/solvers/test_solvers.py b/pyomo/contrib/solver/tests/solvers/test_solvers.py
index a7fddb6356d..b687d6ff6a4 100644
--- a/pyomo/contrib/solver/tests/solvers/test_solvers.py
+++ b/pyomo/contrib/solver/tests/solvers/test_solvers.py
@@ -82,8 +82,7 @@ def param_as_standalone_func(cls, p, func, name):
     ('scip_persistent', ScipPersistent),
     ('gams', GAMS),
     ('knitro_direct', KnitroDirectSolver),
-    ('multistart', MultiStart)
-
+    ('multistart', MultiStart),
 ]
 mip_solvers = [
     ('gurobi_persistent', GurobiPersistent),

From 281ff011fda1fb7c3c5c095027fd6d7e632b3e4c Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Thu, 23 Jul 2026 17:17:13 -0600
Subject: [PATCH 35/55] Finished example, ran black

---
 pyomo/contrib/multistart/multi.py              |  6 +++---
 .../examples/init_polynomial_ex.py             | 18 +++++++++++++-----
 2 files changed, 16 insertions(+), 8 deletions(-)

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index 1c14baa1b58..dfaba5805b5 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -120,7 +120,7 @@ def __init__(
         self.suppress_unbounded_warning = self.declare(
             "suppress_unbounded_warning",
             ConfigValue(
-                default=False,
+                default=True,
                 domain=bool,
                 description="True to suppress warning for skipping unbounded variables.",
             ),
@@ -217,8 +217,8 @@ def solve(self, model, **kwds):
         )
 
         # Set options so infeasible solve does not interrupt runs
-        config.solver_args["load_solutions"] = False
-        config.solver_args["raise_exception_on_nonoptimal_result"] = False
+        # config.solver_args["load_solutions"] = False
+        # config.solver_args["raise_exception_on_nonoptimal_result"] = False
 
         solver = SolverFactory(config.solver)
 
diff --git a/pyomo/devel/initialization/examples/init_polynomial_ex.py b/pyomo/devel/initialization/examples/init_polynomial_ex.py
index 0cf6b178ae6..bcda8527450 100644
--- a/pyomo/devel/initialization/examples/init_polynomial_ex.py
+++ b/pyomo/devel/initialization/examples/init_polynomial_ex.py
@@ -55,17 +55,25 @@ def global_init_ex():
 
 def multistart_init_ex():
     m = build_model()
+    m.obj = pyo.Objective(expr=0)
     nlp_solver = SolverFactory('ipopt')
-    global_solver = SolverFactory('scip_direct')
-    results = ini.initialize_with_global_opt(
-        nlp=m, nlp_solver=nlp_solver, global_solver=global_solver
-    )
+    multistart_solver = pyo.SolverFactory('multistart')
 
-    return results.solution_status, m.x.value
+    multistart_solver.CONFIG.iterations = 50
+
+    opts = {"load_solutions": False, "raise_exception_on_nonoptimal_result": False}
+    multistart_solver.CONFIG.strategy = "rand"
+    multistart_solver.CONFIG.solver_args = opts
+
+    results = ini.initialize_with_multistart_opt(
+        nlp=m, nlp_solver=nlp_solver, multistart_solver=multistart_solver, seed=145
+    )
+    return results, m.x.value
 
 
 if __name__ == '__main__':
     # stat, x = lp_init_ex()
     # stat, x = pwl_init_ex()
     stat, x = global_init_ex()
+    # stat, x = multistart_init_ex()
     print(stat, round(x, 4))

From 127a973e67c16a45d2b1c067d85b133a3cb3300c Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Mon, 27 Jul 2026 08:11:37 -0600
Subject: [PATCH 36/55] Used utility to get vars from shallow clone

---
 pyomo/contrib/multistart/multi.py  | 26 ++++++++++++++++++++------
 pyomo/contrib/multistart/reinit.py |  6 +++---
 2 files changed, 23 insertions(+), 9 deletions(-)

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index dfaba5805b5..cbb67f4e869 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -25,11 +25,12 @@
 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
+from pyomo.core import Objective, Constraint, Var, minimize, value
 from pyomo.contrib.solver.common.base import SolverBase
 from pyomo.contrib.solver.common.config import SolverConfig
 from pyomo.contrib.solver.common.factory import SolverFactory
 from pyomo.contrib.solver.common.results import SolutionStatus
+from pyomo.util.vars_from_expressions import get_vars_from_components
 
 from pyomo.common.dependencies.scipy import stats
 from pyomo.common.dependencies import numpy as np
@@ -120,7 +121,7 @@ def __init__(
         self.suppress_unbounded_warning = self.declare(
             "suppress_unbounded_warning",
             ConfigValue(
-                default=True,
+                default=False,
                 domain=bool,
                 description="True to suppress warning for skipping unbounded variables.",
             ),
@@ -208,7 +209,7 @@ def license_is_valid(self):
 
     def solve(self, model, **kwds):
         # initialize keyword args
-        config = self.CONFIG(kwds.pop('options', {}))
+        config = self.config(kwds.pop('options', {}))
         config.set_value(kwds)
 
         # Create centralized sampler once
@@ -253,9 +254,22 @@ def solve(self, model, **kwds):
             setattr(
                 model,
                 tmp_var_list_name,
-                list(model.component_data_objects(ctype=Var, descend_into=True)),
+                list(model.component_data_objects(Var, descend_into=True)),
             )
-
+            # If the list has nothing in it, check components
+            print(len(model._vars_list))
+            print(model._vars_list)
+            if len(model._vars_list) == 0:
+                setattr(
+                    model,
+                    tmp_var_list_name,
+                    list(get_vars_from_components(model, 
+                                                  ctype=(Constraint, Objective),
+                                                  active=True,
+                                                  ))
+                    
+                )
+            print(model._vars_list)
             best_result = result = solver.solve(model, **config.solver_args)
             # Check the solution status before loading variables into the model.
             if result.solution_status in {
@@ -357,7 +371,7 @@ def __exit__(self, t, v, traceback):
 
 # Sampling class to organize and configure random samplers
 class SamplingManager:
-    def __init__(self, method="uniform", rng=None, seed=None):
+    def __init__(self, method="lhs", rng=None, seed=None):
         aliases = {
             "random_uniform": "uniform",
             "uniform": "uniform",
diff --git a/pyomo/contrib/multistart/reinit.py b/pyomo/contrib/multistart/reinit.py
index 856c3948deb..5465a9429b3 100644
--- a/pyomo/contrib/multistart/reinit.py
+++ b/pyomo/contrib/multistart/reinit.py
@@ -13,7 +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.util.vars_from_expressions import get_vars_from_components
 
 from pyomo.core import Var
 
@@ -81,8 +81,8 @@ def reinitialize_variables(model, config, sampler):
     """
 
     eligible_vars = []
-
-    for var in model.component_data_objects(ctype=Var, descend_into=True):
+    print(model._vars_list[0].value)
+    for var in model._vars_list:
         if var.is_fixed() or not var.is_continuous():
             continue
         if var.lb is None or var.ub is None:

From e087e99814b08bdd45a7f680d8ab665b8b132af6 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Mon, 27 Jul 2026 08:16:56 -0600
Subject: [PATCH 37/55] Remove extra print statements

---
 pyomo/contrib/multistart/multi.py  | 3 ---
 pyomo/contrib/multistart/reinit.py | 1 -
 2 files changed, 4 deletions(-)

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index cbb67f4e869..0ec586f7852 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -257,8 +257,6 @@ def solve(self, model, **kwds):
                 list(model.component_data_objects(Var, descend_into=True)),
             )
             # If the list has nothing in it, check components
-            print(len(model._vars_list))
-            print(model._vars_list)
             if len(model._vars_list) == 0:
                 setattr(
                     model,
@@ -269,7 +267,6 @@ def solve(self, model, **kwds):
                                                   ))
                     
                 )
-            print(model._vars_list)
             best_result = result = solver.solve(model, **config.solver_args)
             # Check the solution status before loading variables into the model.
             if result.solution_status in {
diff --git a/pyomo/contrib/multistart/reinit.py b/pyomo/contrib/multistart/reinit.py
index 5465a9429b3..e688c039901 100644
--- a/pyomo/contrib/multistart/reinit.py
+++ b/pyomo/contrib/multistart/reinit.py
@@ -81,7 +81,6 @@ def reinitialize_variables(model, config, sampler):
     """
 
     eligible_vars = []
-    print(model._vars_list[0].value)
     for var in model._vars_list:
         if var.is_fixed() or not var.is_continuous():
             continue

From 78f0869d439203eacbb7d165aa4d1e6e716b6122 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Mon, 27 Jul 2026 08:17:49 -0600
Subject: [PATCH 38/55] Modify config

---
 .../initialization/examples/init_polynomial_ex.py   | 13 ++++++-------
 pyomo/devel/initialization/initialize.py            |  6 +++---
 pyomo/devel/initialization/multistart_init.py       |  1 +
 3 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/pyomo/devel/initialization/examples/init_polynomial_ex.py b/pyomo/devel/initialization/examples/init_polynomial_ex.py
index bcda8527450..e48eb7fa649 100644
--- a/pyomo/devel/initialization/examples/init_polynomial_ex.py
+++ b/pyomo/devel/initialization/examples/init_polynomial_ex.py
@@ -57,13 +57,12 @@ def multistart_init_ex():
     m = build_model()
     m.obj = pyo.Objective(expr=0)
     nlp_solver = SolverFactory('ipopt')
-    multistart_solver = pyo.SolverFactory('multistart')
-
-    multistart_solver.CONFIG.iterations = 50
+    multistart_solver = SolverFactory('multistart')
 
     opts = {"load_solutions": False, "raise_exception_on_nonoptimal_result": False}
-    multistart_solver.CONFIG.strategy = "rand"
-    multistart_solver.CONFIG.solver_args = opts
+    multistart_solver.config.strategy = "rand_vector"
+    multistart_solver.config.solver_args = opts
+    multistart_solver.config.iterations = 10
 
     results = ini.initialize_with_multistart_opt(
         nlp=m, nlp_solver=nlp_solver, multistart_solver=multistart_solver, seed=145
@@ -74,6 +73,6 @@ def multistart_init_ex():
 if __name__ == '__main__':
     # stat, x = lp_init_ex()
     # stat, x = pwl_init_ex()
-    stat, x = global_init_ex()
-    # stat, x = multistart_init_ex()
+    # stat, x = global_init_ex()
+    stat, x = multistart_init_ex()
     print(stat, round(x, 4))
diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py
index 6e24d5281ed..55862e79402 100644
--- a/pyomo/devel/initialization/initialize.py
+++ b/pyomo/devel/initialization/initialize.py
@@ -361,9 +361,9 @@ def initialize_with_multistart_opt(
 
     if multistart_solver is None:
         multistart_solver = pyo.SolverFactory("multistart")
-        multistart_solver.CONFIG.seed = seed
-        multistart_solver.CONFIG.sampling_method = "lhs"
-        multistart_solver.CONFIG.break_on_solution = True
+        multistart_solver.config.seed = seed
+        multistart_solver.config.sampling_method = "lhs"
+        multistart_solver.config.break_on_solution = 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 f8addbe52b2..46489b32c3d 100644
--- a/pyomo/devel/initialization/multistart_init.py
+++ b/pyomo/devel/initialization/multistart_init.py
@@ -30,6 +30,7 @@ def _initialize_with_multistart_solver(
     # bounds on the nonlinear variables
     bound_all_nonlinear_variables(nlp, default_bound=default_bound)
 
+    multistart_solver.config.seed=seed
     res = multistart_solver.solve(nlp)
     logger.info('Finished multistart optimization iterations.')
 

From 45dd8c5f4f86dfb28f5b8eea4df663de11fa9a06 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Mon, 27 Jul 2026 08:18:46 -0600
Subject: [PATCH 39/55] Ran black

---
 pyomo/contrib/multistart/multi.py             | 10 +++++-----
 pyomo/devel/initialization/multistart_init.py |  2 +-
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index 0ec586f7852..12920420d84 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -261,11 +261,11 @@ def solve(self, model, **kwds):
                 setattr(
                     model,
                     tmp_var_list_name,
-                    list(get_vars_from_components(model, 
-                                                  ctype=(Constraint, Objective),
-                                                  active=True,
-                                                  ))
-                    
+                    list(
+                        get_vars_from_components(
+                            model, ctype=(Constraint, Objective), active=True
+                        )
+                    ),
                 )
             best_result = result = solver.solve(model, **config.solver_args)
             # Check the solution status before loading variables into the model.
diff --git a/pyomo/devel/initialization/multistart_init.py b/pyomo/devel/initialization/multistart_init.py
index 46489b32c3d..b9e5a7f55f3 100644
--- a/pyomo/devel/initialization/multistart_init.py
+++ b/pyomo/devel/initialization/multistart_init.py
@@ -30,7 +30,7 @@ def _initialize_with_multistart_solver(
     # bounds on the nonlinear variables
     bound_all_nonlinear_variables(nlp, default_bound=default_bound)
 
-    multistart_solver.config.seed=seed
+    multistart_solver.config.seed = seed
     res = multistart_solver.solve(nlp)
     logger.info('Finished multistart optimization iterations.')
 

From d58d937a555efce44782c0d62c759202c3f5cfdc Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Mon, 27 Jul 2026 08:37:41 -0600
Subject: [PATCH 40/55] Added break on solution

---
 pyomo/devel/initialization/examples/init_polynomial_ex.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/pyomo/devel/initialization/examples/init_polynomial_ex.py b/pyomo/devel/initialization/examples/init_polynomial_ex.py
index e48eb7fa649..ee15911d70d 100644
--- a/pyomo/devel/initialization/examples/init_polynomial_ex.py
+++ b/pyomo/devel/initialization/examples/init_polynomial_ex.py
@@ -60,9 +60,11 @@ def multistart_init_ex():
     multistart_solver = SolverFactory('multistart')
 
     opts = {"load_solutions": False, "raise_exception_on_nonoptimal_result": False}
+
     multistart_solver.config.strategy = "rand_vector"
     multistart_solver.config.solver_args = opts
     multistart_solver.config.iterations = 10
+    multistart_solver.config.break_on_solution = True
 
     results = ini.initialize_with_multistart_opt(
         nlp=m, nlp_solver=nlp_solver, multistart_solver=multistart_solver, seed=145

From b7262fbb5322f4a65ca77a42767dba0c9415923c Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Mon, 27 Jul 2026 12:52:26 -0600
Subject: [PATCH 41/55] Fixing equal bound variables

---
 pyomo/devel/initialization/multistart_init.py | 10 +++++++++-
 1 file changed, 9 insertions(+), 1 deletion(-)

diff --git a/pyomo/devel/initialization/multistart_init.py b/pyomo/devel/initialization/multistart_init.py
index b9e5a7f55f3..6eda5d16e14 100644
--- a/pyomo/devel/initialization/multistart_init.py
+++ b/pyomo/devel/initialization/multistart_init.py
@@ -15,7 +15,10 @@
 from pyomo.devel.initialization.bounds.bound_variables import (
     bound_all_nonlinear_variables,
 )
-from pyomo.devel.initialization.utils import shallow_clone
+from pyomo.devel.initialization.utils import (
+    shallow_clone, 
+    fix_vars_with_equal_bounds,
+)
 import logging
 
 logger = logging.getLogger(__name__)
@@ -29,6 +32,11 @@ def _initialize_with_multistart_solver(
     nlp = shallow_clone(nlp)
     # bounds on the nonlinear variables
     bound_all_nonlinear_variables(nlp, default_bound=default_bound)
+    logger.info('bounded nonlinear variables')
+
+    # fix variables with equal bounds for sampler
+    fix_vars_with_equal_bounds(nlp)
+    logger.info('fixed variables with equal bounds')
 
     multistart_solver.config.seed = seed
     res = multistart_solver.solve(nlp)

From aa9ba5b19348f8dd4c51f502d96bcde37a58cf39 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Tue, 28 Jul 2026 14:28:16 -0600
Subject: [PATCH 42/55] Starting to fix issues with test_solvers

---
 pyomo/contrib/multistart/multi.py             | 97 ++++++++++++-------
 pyomo/contrib/multistart/tests/test_multi.py  | 15 ---
 .../solver/tests/solvers/test_solvers.py      |  2 +-
 3 files changed, 64 insertions(+), 50 deletions(-)

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index 12920420d84..3eaceb6041e 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -21,6 +21,8 @@
     ADVANCED_OPTION,
 )
 
+from typing import Any
+from pyomo.common.timing import HierarchicalTimer, default_timer
 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
@@ -30,6 +32,7 @@
 from pyomo.contrib.solver.common.config import SolverConfig
 from pyomo.contrib.solver.common.factory import SolverFactory
 from pyomo.contrib.solver.common.results import SolutionStatus
+from pyomo.contrib.solver.common.util import NoOptimalSolutionError, NoSolutionError, NoFeasibleSolutionError
 from pyomo.util.vars_from_expressions import get_vars_from_components
 
 from pyomo.common.dependencies.scipy import stats
@@ -150,7 +153,7 @@ def __init__(
         self.sampling_method = self.declare(
             "sampling_method",
             ConfigValue(
-                default="random_uniform",
+                default="latin_hypercube",
                 description="Method for sampling random starting points for reinitialization step. "
                 "Supported options are 'random_uniform', 'latin_hypercube', and 'sobol_sampling'. "
                 "Only utilized when config.strategy is 'rand_vector'.",
@@ -189,62 +192,80 @@ class MultiStart(SolverBase):
 
     CONFIG = MultistartConfig()
 
+    def __init__(self, **kwds: Any) -> None:
+        super().__init__(**kwds)
+
+        #: Instance configuration;
+        self.config = self.config
+
+
     def available(self, exception_flag=True):
         """Check if solver is available.
 
-        TODO: For now, it is always available. However, sub-solvers may not
-        always be available, and so this should reflect that possibility.
-
+        The multistart solver wrapper should always be available,
+        but it is not guaranteed the subsolvers will be. 
+        Check if the selected subsolver is available, which by default is ipopt.
         """
-        return True
 
-    def version(self):
-        """Get solver version
-        TODO: This is a solver wrapper, unsure how to define version in this case."""
+        subsolver = SolverFactory(self.config.solver)
+        return subsolver.available()
 
-        return
+    def version(self):
+        """Get solver version."""
+        """
+            Original implementation: 0.1.0,
+            Current implementation: 0.2.0,
+        """
+        return (0, 2, 0)
 
     def license_is_valid(self):
         return True
 
     def solve(self, model, **kwds):
+        start_time = default_timer()
+
         # initialize keyword args
         config = self.config(kwds.pop('options', {}))
         config.set_value(kwds)
 
+        timer = config.timer
+        if timer is None:
+            timer = config.timer = HierarchicalTimer()
+
         # Create centralized sampler once
         sampler = SamplingManager(
             method=config.sampling_method, rng=config.rng, seed=config.seed
         )
 
-        # Set options so infeasible solve does not interrupt runs
-        # config.solver_args["load_solutions"] = False
-        # config.solver_args["raise_exception_on_nonoptimal_result"] = False
+        # Set sub-solver options so infeasible solve does not interrupt runs
+        config.solver_args["load_solutions"] = False
+        config.solver_args["raise_exception_on_nonoptimal_result"] = False
+
+        # config.raise_exception_on_nonoptimal_result = False
 
         solver = SolverFactory(config.solver)
 
         # Model sense
-        objectives = model.component_data_objects(Objective, active=True)
-        obj = next(objectives, None)
-        # Check model validity
-        if next(objectives, None) is not None:
+        objectives = list(model.component_data_objects(Objective, active=True))
+        # Check length
+        print(len(objectives))
+        if len(objectives) > 1:
             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"
-        #     )
+        elif len(objectives) == 1:
+            obj = objectives[0]
+            obj.sign = 1 if obj.sense == minimize else -1
+            obj_sign = obj.sign
+
+        else:
+            obj_sign = 1      
+
+        best_objective = float('inf') * obj_sign
 
         # store objective values and objective/result information for best
         # solution obtained
         objectives = []
-        obj_sign = 1 if obj.sense == minimize else -1
-        best_objective = float('inf') * obj_sign
         best_model = model
         best_result = None
 
@@ -309,7 +330,7 @@ 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)
 
                 # Check the solution status before loading variables into the model.
                 if result.solution_status in {
@@ -326,6 +347,7 @@ def solve(self, model, **kwds):
 
                 if best_result.solution_status is SolutionStatus.optimal:
                     model_objectives = m.component_data_objects(Objective, active=True)
+                    print("Model objs", model_objectives)
                     mobj = next(model_objectives)
                     obj_val = value(mobj.expr)
                     objectives.append(obj_val)
@@ -335,16 +357,23 @@ def solve(self, model, **kwds):
                         best_model = m
                         best_result = result
 
-            if using_HCS and not HCS_completed:
-                logger.warning(
-                    "High confidence stopping rule was unable to complete "
-                    "after %s iterations. To increase this limit, change the "
-                    "HCS_max_iterations flag." % num_iter
-                )
+            if using_HCS:
+                if not HCS_completed:
+                    logger.warning(
+                        "High confidence stopping rule was unable to complete "
+                        "after %s iterations. To increase this limit, change the "
+                        "HCS_max_iterations flag." % num_iter
+                    )
 
+            else:
+                if config.raise_exception_on_nonoptimal_result and not using_HCS:
+                    if best_result.solution_status != SolutionStatus.optimal:
+                        raise NoOptimalSolutionError
+                    
+            
             # if no better result was found than initial solve, then return
             # that without needing to copy variables.
-            if best_model is model:
+            if best_model is model:                             
                 return best_result
 
             # reassign the given models vars to the new models vars
diff --git a/pyomo/contrib/multistart/tests/test_multi.py b/pyomo/contrib/multistart/tests/test_multi.py
index a96562ae195..b273ab72177 100644
--- a/pyomo/contrib/multistart/tests/test_multi.py
+++ b/pyomo/contrib/multistart/tests/test_multi.py
@@ -134,21 +134,6 @@ def test_multiple_obj(self):
         with self.assertRaisesRegex(RuntimeError, "multiple active objectives"):
             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():
     """Simple non-convex model with many local minima"""
     model = ConcreteModel()
diff --git a/pyomo/contrib/solver/tests/solvers/test_solvers.py b/pyomo/contrib/solver/tests/solvers/test_solvers.py
index b687d6ff6a4..957658a8a56 100644
--- a/pyomo/contrib/solver/tests/solvers/test_solvers.py
+++ b/pyomo/contrib/solver/tests/solvers/test_solvers.py
@@ -1232,7 +1232,7 @@ def test_trivial_constraints(
         opt.config.load_solutions = False
         res = opt.solve(m)
         self.assertNotEqual(res.solution_status, SolutionStatus.optimal)
-        if isinstance(opt, Ipopt):
+        if isinstance(opt, Ipopt) or isinstance(opt, MultiStart):
             acceptable_termination_conditions = {
                 TerminationCondition.locallyInfeasible,
                 TerminationCondition.unbounded,

From a312c49e499573132b926a7512020c2c2f3a637e Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Tue, 28 Jul 2026 14:29:53 -0600
Subject: [PATCH 43/55] Ran black

---
 pyomo/contrib/multistart/multi.py             | 16 +++++++++-------
 pyomo/contrib/multistart/tests/test_multi.py  |  1 +
 pyomo/devel/initialization/multistart_init.py |  5 +----
 3 files changed, 11 insertions(+), 11 deletions(-)

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index 3eaceb6041e..ab27d97a60c 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -32,7 +32,11 @@
 from pyomo.contrib.solver.common.config import SolverConfig
 from pyomo.contrib.solver.common.factory import SolverFactory
 from pyomo.contrib.solver.common.results import SolutionStatus
-from pyomo.contrib.solver.common.util import NoOptimalSolutionError, NoSolutionError, NoFeasibleSolutionError
+from pyomo.contrib.solver.common.util import (
+    NoOptimalSolutionError,
+    NoSolutionError,
+    NoFeasibleSolutionError,
+)
 from pyomo.util.vars_from_expressions import get_vars_from_components
 
 from pyomo.common.dependencies.scipy import stats
@@ -198,12 +202,11 @@ def __init__(self, **kwds: Any) -> None:
         #: Instance configuration;
         self.config = self.config
 
-
     def available(self, exception_flag=True):
         """Check if solver is available.
 
         The multistart solver wrapper should always be available,
-        but it is not guaranteed the subsolvers will be. 
+        but it is not guaranteed the subsolvers will be.
         Check if the selected subsolver is available, which by default is ipopt.
         """
 
@@ -259,7 +262,7 @@ def solve(self, model, **kwds):
             obj_sign = obj.sign
 
         else:
-            obj_sign = 1      
+            obj_sign = 1
 
         best_objective = float('inf') * obj_sign
 
@@ -369,11 +372,10 @@ def solve(self, model, **kwds):
                 if config.raise_exception_on_nonoptimal_result and not using_HCS:
                     if best_result.solution_status != SolutionStatus.optimal:
                         raise NoOptimalSolutionError
-                    
-            
+
             # if no better result was found than initial solve, then return
             # that without needing to copy variables.
-            if best_model is model:                             
+            if best_model is model:
                 return best_result
 
             # reassign the given models vars to the new models vars
diff --git a/pyomo/contrib/multistart/tests/test_multi.py b/pyomo/contrib/multistart/tests/test_multi.py
index b273ab72177..ea2edb6e236 100644
--- a/pyomo/contrib/multistart/tests/test_multi.py
+++ b/pyomo/contrib/multistart/tests/test_multi.py
@@ -134,6 +134,7 @@ def test_multiple_obj(self):
         with self.assertRaisesRegex(RuntimeError, "multiple active objectives"):
             SolverFactory('multistart').solve(m)
 
+
 def build_model():
     """Simple non-convex model with many local minima"""
     model = ConcreteModel()
diff --git a/pyomo/devel/initialization/multistart_init.py b/pyomo/devel/initialization/multistart_init.py
index 6eda5d16e14..07eafe1b4ce 100644
--- a/pyomo/devel/initialization/multistart_init.py
+++ b/pyomo/devel/initialization/multistart_init.py
@@ -15,10 +15,7 @@
 from pyomo.devel.initialization.bounds.bound_variables import (
     bound_all_nonlinear_variables,
 )
-from pyomo.devel.initialization.utils import (
-    shallow_clone, 
-    fix_vars_with_equal_bounds,
-)
+from pyomo.devel.initialization.utils import shallow_clone, fix_vars_with_equal_bounds
 import logging
 
 logger = logging.getLogger(__name__)

From 836fa9c2b1e740e91a99a17f1decfad2c5fc8568 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Wed, 29 Jul 2026 09:53:48 -0600
Subject: [PATCH 44/55] Mod test for new expected error

---
 pyomo/contrib/multistart/tests/test_multi.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/pyomo/contrib/multistart/tests/test_multi.py b/pyomo/contrib/multistart/tests/test_multi.py
index ea2edb6e236..bcc12e4e6ee 100644
--- a/pyomo/contrib/multistart/tests/test_multi.py
+++ b/pyomo/contrib/multistart/tests/test_multi.py
@@ -27,6 +27,7 @@
     value,
 )
 from pyomo.contrib.solver.common.factory import SolverFactory
+from pyomo.contrib.solver.common.util import NoOptimalSolutionError
 
 
 @unittest.skipIf(not SolverFactory('ipopt').available(), "IPOPT not available")
@@ -106,7 +107,9 @@ def test_model_infeasible(self):
         m.x = Var(bounds=(0, 1))
         m.c = Constraint(expr=m.x >= 2)
         m.o = Objective(expr=m.x)
-        SolverFactory('multistart').solve(m, iterations=2)
+
+        with self.assertRaises(NoOptimalSolutionError):
+            SolverFactory('multistart').solve(m, iterations=2)
         output = StringIO()
         with LoggingIntercept(output, 'pyomo.contrib.multistart', logging.WARNING):
             SolverFactory('multistart').solve(m, iterations=-1, HCS_max_iterations=3)

From 0c805cd1605413eb00ff26971d3c09fe7ccf266a Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Thu, 30 Jul 2026 10:42:27 -0600
Subject: [PATCH 45/55] Work progress on testing for all solvers

---
 pyomo/contrib/multistart/multi.py             | 147 ++++++++++++------
 pyomo/contrib/multistart/reinit.py            |  22 +--
 pyomo/contrib/multistart/tests/test_multi.py  |   1 +
 .../contrib/solver/common/solution_loader.py  |   2 +-
 .../solver/tests/solvers/test_solvers.py      |   6 +-
 5 files changed, 119 insertions(+), 59 deletions(-)

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index ab27d97a60c..ad2caafd9d8 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -22,6 +22,7 @@
 )
 
 from typing import Any
+import datetime
 from pyomo.common.timing import HierarchicalTimer, default_timer
 from pyomo.common.modeling import unique_component_name
 from pyomo.common.dependencies import numpy as np
@@ -31,7 +32,12 @@
 from pyomo.contrib.solver.common.base import SolverBase
 from pyomo.contrib.solver.common.config import SolverConfig
 from pyomo.contrib.solver.common.factory import SolverFactory
-from pyomo.contrib.solver.common.results import SolutionStatus
+from pyomo.contrib.solver.common.results import (
+    Results,
+    TerminationCondition,
+    SolutionStatus,
+)
+
 from pyomo.contrib.solver.common.util import (
     NoOptimalSolutionError,
     NoSolutionError,
@@ -80,16 +86,16 @@ def __init__(
             """,
             ),
         )
-        self.solver = self.declare(
-            "solver",
+        self.subsolver = self.declare(
+            "subsolver",
             ConfigValue(
                 default="ipopt",
                 description="solver to use, defaults to ipopt"
                 "Should also be able to accept solver objects. In progress",
             ),
         )
-        self.solver_args = self.declare(
-            "solver_args",
+        self.subsolver_args = self.declare(
+            "subsolver_args",
             ConfigValue(
                 default={},
                 description="Dictionary of keyword arguments to pass to the solver.",
@@ -210,7 +216,7 @@ def available(self, exception_flag=True):
         Check if the selected subsolver is available, which by default is ipopt.
         """
 
-        subsolver = SolverFactory(self.config.solver)
+        subsolver = SolverFactory(self.config.subsolver)
         return subsolver.available()
 
     def version(self):
@@ -235,18 +241,26 @@ def solve(self, model, **kwds):
         if timer is None:
             timer = config.timer = HierarchicalTimer()
 
+        # Allocate the results object so we can populate it as we go
+        results = Results()
+        results.timing_info.start_timestamp = datetime.datetime.now(
+            datetime.timezone.utc
+        )
+        results.solver_name = self.name
+
         # Create centralized sampler once
         sampler = SamplingManager(
             method=config.sampling_method, rng=config.rng, seed=config.seed
         )
 
-        # Set sub-solver options so infeasible solve does not interrupt runs
-        config.solver_args["load_solutions"] = False
-        config.solver_args["raise_exception_on_nonoptimal_result"] = False
+        # Set sub-solver options
+        config.subsolver_args["load_solutions"] = False
+        config.subsolver_args["raise_exception_on_nonoptimal_result"] = False
 
-        # config.raise_exception_on_nonoptimal_result = False
+        config.subsolver_args["time_limit"] = config.time_limit
+        config.subsolver_args["tee"] = config.tee
 
-        solver = SolverFactory(config.solver)
+        solver = SolverFactory(config.subsolver)
 
         # Model sense
         objectives = list(model.component_data_objects(Objective, active=True))
@@ -262,7 +276,9 @@ def solve(self, model, **kwds):
             obj_sign = obj.sign
 
         else:
+            obj = None
             obj_sign = 1
+            config.break_on_solution = True
 
         best_objective = float('inf') * obj_sign
 
@@ -273,6 +289,7 @@ def solve(self, model, **kwds):
         best_result = None
 
         try:
+            timer.start('initial_solve')
             # create temporary variable list for value transfer
             tmp_var_list_name = unique_component_name(model, "_vars_list")
             setattr(
@@ -291,21 +308,25 @@ def solve(self, model, **kwds):
                         )
                     ),
                 )
-            best_result = result = solver.solve(model, **config.solver_args)
+            best_result = result = solver.solve(model, **config.subsolver_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.load_solutions:
+                result.solution_loader.load_solution()
 
-            if best_result.solution_status is SolutionStatus.optimal:
-                obj_val = value(obj.expr)
-                best_objective = obj_val
-                objectives.append(obj_val)
+            if result.solution_status is SolutionStatus.optimal:
+                if obj is not None:
+                    obj_val = value(obj.expr)
+                    best_objective = obj_val
+                    objectives.append(obj_val)
+            timer.stop('initial_solve')
 
             num_iter = 0
             max_iter = config.iterations
@@ -319,7 +340,9 @@ def solve(self, model, **kwds):
                 ), "High confidence stopping rule requires rand strategy."
                 max_iter = config.HCS_max_iterations
 
+            timer.start('iterative_solves')
             while num_iter < max_iter:
+                # timer.start(f"timer_iter_{num_iter}")
                 if using_HCS and should_stop(
                     objectives,
                     config.stopping_mass,
@@ -327,38 +350,50 @@ def solve(self, model, **kwds):
                     config.HCS_tolerance,
                 ):
                     HCS_completed = True
+                    # timer.stop(f"timer_iter_{num_iter}")
                     break
-                logger.info(f"num_iter: {num_iter}\n")
+
                 num_iter += 1
+                logger.info(f"num_iter: {num_iter}\n")
+
                 # 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)
-
+                result = solver.solve(m, **config.subsolver_args)
+                # if config.load_solutions:
                 # 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.load_solutions:
+                    result.solution_loader.load_solution()
                     # 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)
-                    print("Model objs", model_objectives)
-                    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
+                        # timer.stop(f"timer_iter_{num_iter}")
+                        break
+
+                if result.solution_status is SolutionStatus.optimal:
+                    if obj is not None:
+                        model_objectives = m.component_data_objects(Objective, active=True)
+                        # print("Model objs", model_objectives)
+                        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
+                # timer.stop(f"timer_iter_{num_iter}")
+
+            timer.stop('iterative_solves')
+            print(num_iter)
 
             if using_HCS:
                 if not HCS_completed:
@@ -368,24 +403,46 @@ def solve(self, model, **kwds):
                         "HCS_max_iterations flag." % num_iter
                     )
 
-            else:
-                if config.raise_exception_on_nonoptimal_result and not using_HCS:
-                    if best_result.solution_status != SolutionStatus.optimal:
-                        raise NoOptimalSolutionError
+            # if config.raise_exception_on_nonoptimal_result:
+            #     if best_result.solution_status != SolutionStatus.optimal:
+            #         raise NoOptimalSolutionError
+
+            if (
+                config.raise_exception_on_nonoptimal_result
+                and best_result.solution_status != SolutionStatus.optimal
+            ):
+                raise NoOptimalSolutionError()
 
             # if no better result was found than initial solve, then return
             # that without needing to copy variables.
-            if best_model is model:
-                return best_result
-
-            # reassign the given models vars to the new models vars
             orig_var_list = getattr(model, tmp_var_list_name)
             best_soln_var_list = getattr(best_model, tmp_var_list_name)
-            for orig_var, new_var in zip(orig_var_list, best_soln_var_list):
-                if not orig_var.is_fixed():
-                    orig_var.set_value(new_var.value, skip_validation=True)
+            if config.load_solutions:
+                if best_result.solution_status == SolutionStatus.noSolution:
+                    raise NoSolutionError()
+
+                if best_model is not model:
+                    # reassign the given models vars to the new models vars
+                    for orig_var, new_var in zip(orig_var_list, best_soln_var_list):
+                        if not orig_var.is_fixed():
+                            orig_var.set_value(new_var.value, skip_validation=True)
 
+
+
+            # # # results.subsolver_results = best_result
+            # results.solution_loader = best_result.solution_loader
+            # results.termination_condition = best_result.termination_condition
+            # results.solution_status = best_result.solution_status
+            # results.solution_loader = best_result.solution_loader
+
+            # results.timing_info.timer = timer
+            # results.timing_info.wall_time = default_timer() - start_time
+            # return results
+
+            best_result.timing_info.timer = timer
+            best_result.timing_info.wall_time = default_timer() - start_time
             return best_result
+
         finally:
             # Remove temporary variable list
             delattr(model, tmp_var_list_name)
@@ -425,9 +482,9 @@ def _ensure_qmc(self, dim):
             return
 
         if self.method == "lhs":
-            self.qmc_sampler = stats.qmc.LatinHypercube(d=dim, seed=self.seed)
+            self.qmc_sampler = stats.qmc.LatinHypercube(d=dim, rng=self.rng)
         elif self.method == "sobol":
-            self.qmc_sampler = stats.qmc.Sobol(d=dim, scramble=True, seed=self.seed)
+            self.qmc_sampler = stats.qmc.Sobol(d=dim, scramble=True, seed=self.rng)
         else:
             raise ValueError(f"QMC sampler not valid for method '{self.method}'")
 
diff --git a/pyomo/contrib/multistart/reinit.py b/pyomo/contrib/multistart/reinit.py
index e688c039901..bb7f0ee81c8 100644
--- a/pyomo/contrib/multistart/reinit.py
+++ b/pyomo/contrib/multistart/reinit.py
@@ -20,8 +20,10 @@
 logger = logging.getLogger('pyomo.contrib.multistart')
 
 
-def rand(val, lb, ub, rng):
-    sample = rng.uniform(lb, ub)  # uniform distribution between lb and ub
+def rand(val, lb, ub, sampler):
+    # if sampler.method == "uniform":
+    sample = sampler.rng.uniform(lb, ub)  # uniform distribution between lb and ub
+    
     return sample
 
 
@@ -33,28 +35,28 @@ def rand_vector(lbs, ubs, sampler):
     return samples
 
 
-def midpoint_guess_and_bound(val, lb, ub, rng=None):
+def midpoint_guess_and_bound(val, lb, ub, sampler=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
 
 
-def rand_guess_and_bound(val, lb, ub, rng):
+def rand_guess_and_bound(val, lb, ub, sampler):
     """Random choice between current value and farthest bound."""
     far_bound = ub if ((ub - val) >= (val - lb)) else lb  # farther bound
     if far_bound == ub:
-        return rng.uniform(val, far_bound)
+        return sampler.rng.uniform(val, far_bound)
     else:
-        return rng.uniform(far_bound, val)
+        return sampler.rng.uniform(far_bound, val)
 
 
-def rand_distributed(val, lb, ub, rng, divisions=9):
+def rand_distributed(val, lb, ub, sampler, divisions=9):
     """Random choice among evenly distributed set of values between bounds."""
     set_distributed_vals = linspace(lb, ub, divisions)
-    return rng.choice(set_distributed_vals)
+    return sampler.rng.choice(set_distributed_vals)
 
 
-def simple_midpoint(val, lb, ub, rng=None):
+def simple_midpoint(val, lb, ub, sampler=None):
     return (lb + ub) * 0.5
 
 
@@ -121,6 +123,6 @@ def reinitialize_variables(model, config, sampler):
             # print(f"val = {val}\n")
             # apply reinitialization strategy to variable
             var.set_value(
-                strategies[config.strategy](val, var.lb, var.ub, sampler.rng),
+                strategies[config.strategy](val, var.lb, var.ub, sampler),
                 skip_validation=True,
             )
diff --git a/pyomo/contrib/multistart/tests/test_multi.py b/pyomo/contrib/multistart/tests/test_multi.py
index bcc12e4e6ee..e6f5795d1d6 100644
--- a/pyomo/contrib/multistart/tests/test_multi.py
+++ b/pyomo/contrib/multistart/tests/test_multi.py
@@ -61,6 +61,7 @@ def test_as_good_as_standard(self):
                 clone_objective_value, standard_objective_value
             )  # assumes maximization
 
+
     def test_as_good_with_HCS_rule(self):
         """test that the high confidence stopping rule with very lenient
         parameters does no worse.
diff --git a/pyomo/contrib/solver/common/solution_loader.py b/pyomo/contrib/solver/common/solution_loader.py
index faaf7c75685..2fda01e5c32 100644
--- a/pyomo/contrib/solver/common/solution_loader.py
+++ b/pyomo/contrib/solver/common/solution_loader.py
@@ -63,7 +63,7 @@ def solution(self, solution_id: Any) -> "SolutionLoaderView":
 
            results = solver.solve(model)
            results.solution(2).load_vars()
-           results.solution(2).load_import_suffixes()
+           results.solution(l2).load_import_suffixes()
 
         Parameters
         ----------
diff --git a/pyomo/contrib/solver/tests/solvers/test_solvers.py b/pyomo/contrib/solver/tests/solvers/test_solvers.py
index 957658a8a56..9fc95e32d91 100644
--- a/pyomo/contrib/solver/tests/solvers/test_solvers.py
+++ b/pyomo/contrib/solver/tests/solvers/test_solvers.py
@@ -1153,7 +1153,7 @@ def test_results_infeasible(
         opt.config.load_solutions = False
         res = opt.solve(m)
         self.assertNotEqual(res.solution_status, SolutionStatus.optimal)
-        if isinstance(opt, Ipopt):
+        if isinstance(opt, Ipopt) or isinstance(opt, MultiStart):
             acceptable_termination_conditions = {
                 TerminationCondition.locallyInfeasible,
                 TerminationCondition.unbounded,
@@ -1168,7 +1168,7 @@ def test_results_infeasible(
         self.assertAlmostEqual(m.y.value, None)
         self.assertTrue(res.incumbent_objective is None)
 
-        if not isinstance(opt, Ipopt):
+        if not isinstance(opt, Ipopt) and not isinstance(opt, MultiStart):
             # ipopt can return the values of the variables/duals at the last iterate
             # even if it did not converge; raise_exception_on_nonoptimal_result
             # is set to False, so we are free to load infeasible solutions
@@ -1912,7 +1912,7 @@ def test_time_limit(
                 constant=0,
             )
             m.c2[t] = expr == 1
-        if isinstance(opt, Ipopt):
+        if isinstance(opt, Ipopt) or isinstance(opt, MultiStart):
             opt.config.time_limit = 1e-6
         else:
             opt.config.time_limit = 0

From 746b9969763c606abd9c88bf4de8d1db137ebf6d Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Tue, 4 Aug 2026 09:23:11 -0600
Subject: [PATCH 46/55] Swapping to new design with custom results

---
 pyomo/contrib/multistart/multi.py | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index ad2caafd9d8..85c4230aeec 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -185,6 +185,10 @@ def __init__(
             ),
         )
 
+# class MultiStartResults(Results):
+
+
+
 
 @SolverFactory.register('multistart', doc='MultiStart solver for NLPs')
 @document_class_CONFIG(methods=['solve'])
@@ -323,7 +327,7 @@ def solve(self, model, **kwds):
 
             if result.solution_status is SolutionStatus.optimal:
                 if obj is not None:
-                    obj_val = value(obj.expr)
+                    obj_val = result.incumbent_objective
                     best_objective = obj_val
                     objectives.append(obj_val)
             timer.stop('initial_solve')

From d5c517d7adc5c52c32faa48a619ab0cd62b6960c Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Wed, 5 Aug 2026 12:16:34 -0600
Subject: [PATCH 47/55] Progress on multistart redesign

---
 pyomo/contrib/multistart/multi.py             | 312 +++++++++---------
 .../examples/init_polynomial_ex.py            |   6 +-
 2 files changed, 168 insertions(+), 150 deletions(-)

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index 85c4230aeec..979b1376ca9 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -21,7 +21,7 @@
     ADVANCED_OPTION,
 )
 
-from typing import Any
+from typing import Any, Optional
 import datetime
 from pyomo.common.timing import HierarchicalTimer, default_timer
 from pyomo.common.modeling import unique_component_name
@@ -47,10 +47,13 @@
 
 from pyomo.common.dependencies.scipy import stats
 from pyomo.common.dependencies import numpy as np
+from pyomo.core.staleflag import StaleFlagManager
+
 
 logger = logging.getLogger('pyomo.contrib.multistart')
 
 
+
 @document_configdict()
 class MultistartConfig(SolverConfig):
     def __init__(
@@ -127,7 +130,7 @@ def __init__(
                 description="1 minus the confidence level required for the stopping rule.",
                 doc="""1 minus the confidence level required for the stopping rule for the
             high confidence stopping rule, only used with the random strategy.
-            The lower the parameter, the stricter the rule.
+            The lower the parameter, the stricter the rule.               visibility=DEVELOPER_OPTION,
             Value bounded in (0, 1].""",
             ),
         )
@@ -185,9 +188,28 @@ def __init__(
             ),
         )
 
-# class MultiStartResults(Results):
-
-
+class MultiStartResults(Results):
+    def __init__(
+    self,
+    description=None,
+    doc=None,
+    implicit=False,
+    implicit_domain=None,
+    visibility=0,
+):
+        super().__init__(
+            description=description,
+            doc=doc,
+            implicit=implicit,
+            implicit_domain=implicit_domain,
+            visibility=visibility,
+        )
+        self.feasible_solution_list: Optional[list] = self.declare(
+            'feasible_solution_list',
+            ConfigValue(
+                description="Object for loading the solution back into the model.",
+            ),
+        )
 
 
 @SolverFactory.register('multistart', doc='MultiStart solver for NLPs')
@@ -246,11 +268,12 @@ def solve(self, model, **kwds):
             timer = config.timer = HierarchicalTimer()
 
         # Allocate the results object so we can populate it as we go
-        results = Results()
+        results = MultiStartResults()
         results.timing_info.start_timestamp = datetime.datetime.now(
             datetime.timezone.utc
         )
-        results.solver_name = self.name
+        # As we are about to run a solver, update the stale flag
+        StaleFlagManager.mark_all_as_stale()
 
         # Create centralized sampler once
         sampler = SamplingManager(
@@ -260,7 +283,6 @@ def solve(self, model, **kwds):
         # Set sub-solver options
         config.subsolver_args["load_solutions"] = False
         config.subsolver_args["raise_exception_on_nonoptimal_result"] = False
-
         config.subsolver_args["time_limit"] = config.time_limit
         config.subsolver_args["tee"] = config.tee
 
@@ -283,173 +305,169 @@ def solve(self, model, **kwds):
             obj = None
             obj_sign = 1
             config.break_on_solution = True
-
-        best_objective = float('inf') * obj_sign
-
+     
         # store objective values and objective/result information for best
         # solution obtained
         objectives = []
         best_model = model
         best_result = None
+        best_objective = float('inf') * obj_sign
+        results.feasible_solution_list = []
 
-        try:
-            timer.start('initial_solve')
-            # create temporary variable list for value transfer
-            tmp_var_list_name = unique_component_name(model, "_vars_list")
+
+
+        timer.start('initial_solve')
+        # create temporary variable list for value transfer
+        tmp_var_list_name = unique_component_name(model, "_vars_list")
+        setattr(
+            model,
+            tmp_var_list_name,
+            list(model.component_data_objects(Var, descend_into=True)),
+        )
+        # If the list has nothing in it, check components
+        if len(model._vars_list) == 0:
             setattr(
                 model,
                 tmp_var_list_name,
-                list(model.component_data_objects(Var, descend_into=True)),
+                list(
+                    get_vars_from_components(
+                        model, ctype=(Constraint, Objective), active=True
+                    )
+                ),
             )
-            # If the list has nothing in it, check components
-            if len(model._vars_list) == 0:
-                setattr(
-                    model,
-                    tmp_var_list_name,
-                    list(
-                        get_vars_from_components(
-                            model, ctype=(Constraint, Objective), active=True
-                        )
-                    ),
-                )
-            best_result = result = solver.solve(model, **config.subsolver_args)
+        best_result = result = solver.solve(model, **config.subsolver_args)
+        # Check the solution status before loading variables into the model.
+        if result.solution_status in {
+            SolutionStatus.feasible,
+            SolutionStatus.optimal,
+        }:
+
+            logger.info(
+                f'solved NLP: {result.solution_status}, {result.termination_condition}'
+            )
+            # if config.load_solutions:
+            result.solution_loader.load_solution()
+
+        if result.solution_status is SolutionStatus.optimal:
+            if obj is not None:
+                obj_val = result.incumbent_objective
+                best_objective = obj_val
+                objectives.append(obj_val)
+        timer.stop('initial_solve')
+
+        num_iter = 0
+        max_iter = config.iterations
+        # if HCS rule is specified, reinitialize completely randomly until
+        # rule specifies stopping
+        using_HCS = config.iterations == -1
+        HCS_completed = False
+        if using_HCS:
+            assert (
+                config.strategy == "rand"
+            ), "High confidence stopping rule requires rand strategy."
+            max_iter = config.HCS_max_iterations
+
+        timer.start('iterative_solves')
+        while num_iter < max_iter:
+            # timer.start(f"timer_iter_{num_iter}")
+            if using_HCS and should_stop(
+                objectives,
+                config.stopping_mass,
+                config.stopping_delta,
+                config.HCS_tolerance,
+            ):
+                HCS_completed = True
+                # timer.stop(f"timer_iter_{num_iter}")
+                break
+
+            num_iter += 1
+            logger.info(f"num_iter: {num_iter}\n")
+
+            # at first iteration, solve the originally passed model
+            m = model
+            reinitialize_variables(m, config, sampler)
+            result = solver.solve(m, **config.subsolver_args)
+            # if config.load_solutions:
             # Check the solution status before loading variables into the model.
             if result.solution_status in {
                 SolutionStatus.feasible,
                 SolutionStatus.optimal,
             }:
-
                 logger.info(
                     f'solved NLP: {result.solution_status}, {result.termination_condition}'
                 )
                 # if config.load_solutions:
-                result.solution_loader.load_solution()
+                results.feasible_solution_list.append(result)
+                # If we are looking for the first feasible solution, then return immediately
+                if config.break_on_solution:
+                    best_model = m
+                    best_result = result
+                    # timer.stop(f"timer_iter_{num_iter}")
+                    break
 
             if result.solution_status is SolutionStatus.optimal:
                 if obj is not None:
-                    obj_val = result.incumbent_objective
-                    best_objective = obj_val
+                    model_objectives = m.component_data_objects(Objective, active=True)
+                    # print("Model objs", model_objectives)
+                    mobj = next(model_objectives)
+                    obj_val = value(mobj.expr)
                     objectives.append(obj_val)
-            timer.stop('initial_solve')
-
-            num_iter = 0
-            max_iter = config.iterations
-            # if HCS rule is specified, reinitialize completely randomly until
-            # rule specifies stopping
-            using_HCS = config.iterations == -1
-            HCS_completed = False
-            if using_HCS:
-                assert (
-                    config.strategy == "rand"
-                ), "High confidence stopping rule requires rand strategy."
-                max_iter = config.HCS_max_iterations
-
-            timer.start('iterative_solves')
-            while num_iter < max_iter:
-                # timer.start(f"timer_iter_{num_iter}")
-                if using_HCS and should_stop(
-                    objectives,
-                    config.stopping_mass,
-                    config.stopping_delta,
-                    config.HCS_tolerance,
-                ):
-                    HCS_completed = True
-                    # timer.stop(f"timer_iter_{num_iter}")
-                    break
-
-                num_iter += 1
-                logger.info(f"num_iter: {num_iter}\n")
-
-                # 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.subsolver_args)
-                # if config.load_solutions:
-                # Check the solution status before loading variables into the model.
-                if result.solution_status in {
-                    SolutionStatus.feasible,
-                    SolutionStatus.optimal,
-                }:
-                    logger.info(
-                        f'solved NLP: {result.solution_status}, {result.termination_condition}'
-                    )
-                    # if config.load_solutions:
-                    result.solution_loader.load_solution()
-                    # If we are looking for the first feasible solution, then return immediately
-                    if config.break_on_solution:
+                    if obj_val * obj_sign < obj_sign * best_objective:
+                        # objective has improved
+                        best_objective = obj_val
                         best_model = m
                         best_result = result
-                        # timer.stop(f"timer_iter_{num_iter}")
-                        break
-
-                if result.solution_status is SolutionStatus.optimal:
-                    if obj is not None:
-                        model_objectives = m.component_data_objects(Objective, active=True)
-                        # print("Model objs", model_objectives)
-                        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
-                # timer.stop(f"timer_iter_{num_iter}")
-
-            timer.stop('iterative_solves')
-            print(num_iter)
-
-            if using_HCS:
-                if not HCS_completed:
-                    logger.warning(
-                        "High confidence stopping rule was unable to complete "
-                        "after %s iterations. To increase this limit, change the "
-                        "HCS_max_iterations flag." % num_iter
-                    )
-
-            # if config.raise_exception_on_nonoptimal_result:
-            #     if best_result.solution_status != SolutionStatus.optimal:
-            #         raise NoOptimalSolutionError
-
-            if (
-                config.raise_exception_on_nonoptimal_result
-                and best_result.solution_status != SolutionStatus.optimal
-            ):
-                raise NoOptimalSolutionError()
+            # timer.stop(f"timer_iter_{num_iter}")
 
-            # if no better result was found than initial solve, then return
-            # that without needing to copy variables.
-            orig_var_list = getattr(model, tmp_var_list_name)
-            best_soln_var_list = getattr(best_model, tmp_var_list_name)
-            if config.load_solutions:
-                if best_result.solution_status == SolutionStatus.noSolution:
-                    raise NoSolutionError()
+        timer.stop('iterative_solves')
+        print(num_iter)
 
-                if best_model is not model:
-                    # reassign the given models vars to the new models vars
-                    for orig_var, new_var in zip(orig_var_list, best_soln_var_list):
-                        if not orig_var.is_fixed():
-                            orig_var.set_value(new_var.value, skip_validation=True)
-
-
-
-            # # # results.subsolver_results = best_result
-            # results.solution_loader = best_result.solution_loader
-            # results.termination_condition = best_result.termination_condition
-            # results.solution_status = best_result.solution_status
-            # results.solution_loader = best_result.solution_loader
-
-            # results.timing_info.timer = timer
-            # results.timing_info.wall_time = default_timer() - start_time
-            # return results
-
-            best_result.timing_info.timer = timer
-            best_result.timing_info.wall_time = default_timer() - start_time
-            return best_result
+        if using_HCS:
+            if not HCS_completed:
+                logger.warning(
+                    "High confidence stopping rule was unable to complete "
+                    "after %s iterations. To increase this limit, change the "
+                    "HCS_max_iterations flag." % num_iter
+                )
 
-        finally:
-            # Remove temporary variable list
-            delattr(model, tmp_var_list_name)
+        # if config.raise_exception_on_nonoptimal_result:
+        #     if best_result.solution_status != SolutionStatus.optimal:
+        #         raise NoOptimalSolutionError
+
+        if (
+            config.raise_exception_on_nonoptimal_result
+            and best_result.solution_status != SolutionStatus.optimal
+        ):
+            raise NoOptimalSolutionError()
+
+        # if no better result was found than initial solve, then return
+        # that without needing to copy variables.
+        orig_var_list = getattr(model, tmp_var_list_name)
+        best_soln_var_list = getattr(best_model, tmp_var_list_name)
+        if config.load_solutions:
+            if best_result.solution_status == SolutionStatus.noSolution:
+                raise NoSolutionError()
+
+            if best_model is not model:
+                # reassign the given models vars to the new models vars
+                for orig_var, new_var in zip(orig_var_list, best_soln_var_list):
+                    if not orig_var.is_fixed():
+                        orig_var.set_value(new_var.value, skip_validation=True)
+
+
+        best_result.timing_info.timer = timer
+        best_result.timing_info.wall_time = default_timer() - start_time
+        return best_result
+
+        # # results.subsolver_results = best_result
+        results.solution_loader = best_result.solution_loader
+        results.termination_condition = best_result.termination_condition
+        results.solution_status = best_result.solution_status
+        results.solution_loader = best_result.solution_loader
+
+        results.timing_info.timer = timer
+        results.timing_info.wall_time = default_timer() - start_time
+        return results
 
     def __enter__(self):
         return self
diff --git a/pyomo/devel/initialization/examples/init_polynomial_ex.py b/pyomo/devel/initialization/examples/init_polynomial_ex.py
index ee15911d70d..8afd438e007 100644
--- a/pyomo/devel/initialization/examples/init_polynomial_ex.py
+++ b/pyomo/devel/initialization/examples/init_polynomial_ex.py
@@ -55,14 +55,14 @@ def global_init_ex():
 
 def multistart_init_ex():
     m = build_model()
-    m.obj = pyo.Objective(expr=0)
+    # m.obj = pyo.Objective(expr=0)
     nlp_solver = SolverFactory('ipopt')
     multistart_solver = SolverFactory('multistart')
 
-    opts = {"load_solutions": False, "raise_exception_on_nonoptimal_result": False}
+    # opts = {"load_solutions": False, "raise_exception_on_nonoptimal_result": False}
 
     multistart_solver.config.strategy = "rand_vector"
-    multistart_solver.config.solver_args = opts
+    # multistart_solver.config.solver_args = opts
     multistart_solver.config.iterations = 10
     multistart_solver.config.break_on_solution = True
 

From bcfbc05ddbe3de52673e64e29ed2e53591402a33 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Wed, 5 Aug 2026 13:16:33 -0600
Subject: [PATCH 48/55] All test_solver tests passing!

---
 pyomo/contrib/multistart/multi.py | 53 ++++++++++---------------------
 1 file changed, 17 insertions(+), 36 deletions(-)

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index 979b1376ca9..c0fdd800042 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -251,7 +251,8 @@ def version(self):
             Original implementation: 0.1.0,
             Current implementation: 0.2.0,
         """
-        return (0, 2, 0)
+        current = (0,2,0)
+        return current
 
     def license_is_valid(self):
         return True
@@ -309,13 +310,10 @@ def solve(self, model, **kwds):
         # store objective values and objective/result information for best
         # solution obtained
         objectives = []
-        best_model = model
         best_result = None
         best_objective = float('inf') * obj_sign
         results.feasible_solution_list = []
 
-
-
         timer.start('initial_solve')
         # create temporary variable list for value transfer
         tmp_var_list_name = unique_component_name(model, "_vars_list")
@@ -341,12 +339,10 @@ def solve(self, model, **kwds):
             SolutionStatus.feasible,
             SolutionStatus.optimal,
         }:
-
+            results.feasible_solution_list.append(result)
             logger.info(
                 f'solved NLP: {result.solution_status}, {result.termination_condition}'
             )
-            # if config.load_solutions:
-            result.solution_loader.load_solution()
 
         if result.solution_status is SolutionStatus.optimal:
             if obj is not None:
@@ -396,11 +392,9 @@ def solve(self, model, **kwds):
                 logger.info(
                     f'solved NLP: {result.solution_status}, {result.termination_condition}'
                 )
-                # if config.load_solutions:
                 results.feasible_solution_list.append(result)
                 # If we are looking for the first feasible solution, then return immediately
                 if config.break_on_solution:
-                    best_model = m
                     best_result = result
                     # timer.stop(f"timer_iter_{num_iter}")
                     break
@@ -410,7 +404,7 @@ def solve(self, model, **kwds):
                     model_objectives = m.component_data_objects(Objective, active=True)
                     # print("Model objs", model_objectives)
                     mobj = next(model_objectives)
-                    obj_val = value(mobj.expr)
+                    obj_val = result.incumbent_objective
                     objectives.append(obj_val)
                     if obj_val * obj_sign < obj_sign * best_objective:
                         # objective has improved
@@ -422,52 +416,39 @@ def solve(self, model, **kwds):
         timer.stop('iterative_solves')
         print(num_iter)
 
-        if using_HCS:
-            if not HCS_completed:
+        if using_HCS and not HCS_completed:
                 logger.warning(
                     "High confidence stopping rule was unable to complete "
                     "after %s iterations. To increase this limit, change the "
                     "HCS_max_iterations flag." % num_iter
                 )
 
-        # if config.raise_exception_on_nonoptimal_result:
-        #     if best_result.solution_status != SolutionStatus.optimal:
-        #         raise NoOptimalSolutionError
-
         if (
             config.raise_exception_on_nonoptimal_result
             and best_result.solution_status != SolutionStatus.optimal
         ):
             raise NoOptimalSolutionError()
 
-        # if no better result was found than initial solve, then return
-        # that without needing to copy variables.
-        orig_var_list = getattr(model, tmp_var_list_name)
-        best_soln_var_list = getattr(best_model, tmp_var_list_name)
-        if config.load_solutions:
-            if best_result.solution_status == SolutionStatus.noSolution:
-                raise NoSolutionError()
-
-            if best_model is not model:
-                # reassign the given models vars to the new models vars
-                for orig_var, new_var in zip(orig_var_list, best_soln_var_list):
-                    if not orig_var.is_fixed():
-                        orig_var.set_value(new_var.value, skip_validation=True)
-
 
-        best_result.timing_info.timer = timer
-        best_result.timing_info.wall_time = default_timer() - start_time
-        return best_result
-
-        # # results.subsolver_results = best_result
         results.solution_loader = best_result.solution_loader
         results.termination_condition = best_result.termination_condition
         results.solution_status = best_result.solution_status
-        results.solution_loader = best_result.solution_loader
+        results.incumbent_objective = best_result.incumbent_objective
+        results.solver_log = best_result.solver_log
+
+        if config.load_solutions:
+            if results.solution_status == SolutionStatus.noSolution:
+                raise NoSolutionError()
+
+            results.solution_loader.load_solution()
 
+        results.solver_name = self.name
+        results.solver_version = self.version()
+        results.solver_config = config
         results.timing_info.timer = timer
         results.timing_info.wall_time = default_timer() - start_time
         return results
+    
 
     def __enter__(self):
         return self

From 1a5c276a2a74795424512fb60d6f949f66973124 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Wed, 5 Aug 2026 13:31:35 -0600
Subject: [PATCH 49/55] All og multistart tests also pass!

---
 pyomo/contrib/multistart/multi.py            | 34 +++++++++-----------
 pyomo/contrib/multistart/reinit.py           |  2 +-
 pyomo/contrib/multistart/tests/test_multi.py |  8 +++--
 3 files changed, 22 insertions(+), 22 deletions(-)

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index c0fdd800042..aee2b7a51fb 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -49,11 +49,9 @@
 from pyomo.common.dependencies import numpy as np
 from pyomo.core.staleflag import StaleFlagManager
 
-
 logger = logging.getLogger('pyomo.contrib.multistart')
 
 
-
 @document_configdict()
 class MultistartConfig(SolverConfig):
     def __init__(
@@ -188,15 +186,16 @@ def __init__(
             ),
         )
 
+
 class MultiStartResults(Results):
     def __init__(
-    self,
-    description=None,
-    doc=None,
-    implicit=False,
-    implicit_domain=None,
-    visibility=0,
-):
+        self,
+        description=None,
+        doc=None,
+        implicit=False,
+        implicit_domain=None,
+        visibility=0,
+    ):
         super().__init__(
             description=description,
             doc=doc,
@@ -207,7 +206,7 @@ def __init__(
         self.feasible_solution_list: Optional[list] = self.declare(
             'feasible_solution_list',
             ConfigValue(
-                description="Object for loading the solution back into the model.",
+                description="Object for loading the solution back into the model."
             ),
         )
 
@@ -251,7 +250,7 @@ def version(self):
             Original implementation: 0.1.0,
             Current implementation: 0.2.0,
         """
-        current = (0,2,0)
+        current = (0, 2, 0)
         return current
 
     def license_is_valid(self):
@@ -306,7 +305,7 @@ def solve(self, model, **kwds):
             obj = None
             obj_sign = 1
             config.break_on_solution = True
-     
+
         # store objective values and objective/result information for best
         # solution obtained
         objectives = []
@@ -335,10 +334,7 @@ def solve(self, model, **kwds):
             )
         best_result = result = solver.solve(model, **config.subsolver_args)
         # Check the solution status before loading variables into the model.
-        if result.solution_status in {
-            SolutionStatus.feasible,
-            SolutionStatus.optimal,
-        }:
+        if result.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
             results.feasible_solution_list.append(result)
             logger.info(
                 f'solved NLP: {result.solution_status}, {result.termination_condition}'
@@ -414,9 +410,11 @@ def solve(self, model, **kwds):
             # timer.stop(f"timer_iter_{num_iter}")
 
         timer.stop('iterative_solves')
+        delattr(model, tmp_var_list_name)
         print(num_iter)
 
-        if using_HCS and not HCS_completed:
+        if using_HCS:
+            if not HCS_completed:
                 logger.warning(
                     "High confidence stopping rule was unable to complete "
                     "after %s iterations. To increase this limit, change the "
@@ -429,7 +427,6 @@ def solve(self, model, **kwds):
         ):
             raise NoOptimalSolutionError()
 
-
         results.solution_loader = best_result.solution_loader
         results.termination_condition = best_result.termination_condition
         results.solution_status = best_result.solution_status
@@ -448,7 +445,6 @@ def solve(self, model, **kwds):
         results.timing_info.timer = timer
         results.timing_info.wall_time = default_timer() - start_time
         return results
-    
 
     def __enter__(self):
         return self
diff --git a/pyomo/contrib/multistart/reinit.py b/pyomo/contrib/multistart/reinit.py
index bb7f0ee81c8..4bdb3568bba 100644
--- a/pyomo/contrib/multistart/reinit.py
+++ b/pyomo/contrib/multistart/reinit.py
@@ -23,7 +23,7 @@
 def rand(val, lb, ub, sampler):
     # if sampler.method == "uniform":
     sample = sampler.rng.uniform(lb, ub)  # uniform distribution between lb and ub
-    
+
     return sample
 
 
diff --git a/pyomo/contrib/multistart/tests/test_multi.py b/pyomo/contrib/multistart/tests/test_multi.py
index e6f5795d1d6..cd6ca6878c2 100644
--- a/pyomo/contrib/multistart/tests/test_multi.py
+++ b/pyomo/contrib/multistart/tests/test_multi.py
@@ -61,7 +61,6 @@ def test_as_good_as_standard(self):
                 clone_objective_value, standard_objective_value
             )  # assumes maximization
 
-
     def test_as_good_with_HCS_rule(self):
         """test that the high confidence stopping rule with very lenient
         parameters does no worse.
@@ -113,7 +112,12 @@ def test_model_infeasible(self):
             SolverFactory('multistart').solve(m, iterations=2)
         output = StringIO()
         with LoggingIntercept(output, 'pyomo.contrib.multistart', logging.WARNING):
-            SolverFactory('multistart').solve(m, iterations=-1, HCS_max_iterations=3)
+            SolverFactory('multistart').solve(
+                m,
+                iterations=-1,
+                HCS_max_iterations=3,
+                raise_exception_on_nonoptimal_result=False,
+            )
             self.assertIn(
                 "High confidence stopping rule was unable to "
                 "complete after 3 iterations.",

From 2ace3cebd4511c2d48708b79e05a759d41094da2 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Wed, 5 Aug 2026 14:25:25 -0600
Subject: [PATCH 50/55] Add support for all sampling methods in all strategies

---
 pyomo/contrib/multistart/multi.py  | 17 +++++++++++------
 pyomo/contrib/multistart/reinit.py | 11 ++++-------
 2 files changed, 15 insertions(+), 13 deletions(-)

diff --git a/pyomo/contrib/multistart/multi.py b/pyomo/contrib/multistart/multi.py
index aee2b7a51fb..92e0f95d8cc 100644
--- a/pyomo/contrib/multistart/multi.py
+++ b/pyomo/contrib/multistart/multi.py
@@ -291,7 +291,6 @@ def solve(self, model, **kwds):
         # Model sense
         objectives = list(model.component_data_objects(Objective, active=True))
         # Check length
-        print(len(objectives))
         if len(objectives) > 1:
             raise RuntimeError(
                 "Multistart solver is unable to handle model with multiple active objectives."
@@ -397,21 +396,16 @@ def solve(self, model, **kwds):
 
             if result.solution_status is SolutionStatus.optimal:
                 if obj is not None:
-                    model_objectives = m.component_data_objects(Objective, active=True)
-                    # print("Model objs", model_objectives)
-                    mobj = next(model_objectives)
                     obj_val = result.incumbent_objective
                     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
             # timer.stop(f"timer_iter_{num_iter}")
 
         timer.stop('iterative_solves')
         delattr(model, tmp_var_list_name)
-        print(num_iter)
 
         if using_HCS:
             if not HCS_completed:
@@ -487,6 +481,17 @@ def _ensure_qmc(self, dim):
         else:
             raise ValueError(f"QMC sampler not valid for method '{self.method}'")
 
+    def sample_scalar(self, lower, upper):
+        if self.method == "uniform":
+            return self.rng.uniform(lower, upper)
+
+        if self.method in ("lhs", "sobol"):
+            self._ensure_qmc(dim=1)
+            x = self.qmc_sampler.random(n=1)  # shape (1, d)
+            return stats.qmc.scale(x, lower, upper).item()
+
+        raise ValueError(f"Unknown sampling 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)
diff --git a/pyomo/contrib/multistart/reinit.py b/pyomo/contrib/multistart/reinit.py
index 4bdb3568bba..04aa40cf4a1 100644
--- a/pyomo/contrib/multistart/reinit.py
+++ b/pyomo/contrib/multistart/reinit.py
@@ -15,15 +15,12 @@
 from pyomo.common.dependencies.scipy import stats
 from pyomo.util.vars_from_expressions import get_vars_from_components
 
-from pyomo.core import Var
-
 logger = logging.getLogger('pyomo.contrib.multistart')
 
 
 def rand(val, lb, ub, sampler):
-    # if sampler.method == "uniform":
-    sample = sampler.rng.uniform(lb, ub)  # uniform distribution between lb and ub
-
+    # sample = sampler.rng.uniform(lb, ub)
+    sample = sampler.sample_scalar(lb, ub)  # uniform distribution between lb and ub
     return sample
 
 
@@ -45,9 +42,9 @@ def rand_guess_and_bound(val, lb, ub, sampler):
     """Random choice between current value and farthest bound."""
     far_bound = ub if ((ub - val) >= (val - lb)) else lb  # farther bound
     if far_bound == ub:
-        return sampler.rng.uniform(val, far_bound)
+        return sampler.sample_scalar(val, far_bound)
     else:
-        return sampler.rng.uniform(far_bound, val)
+        return sampler.sample_scalar(far_bound, val)
 
 
 def rand_distributed(val, lb, ub, sampler, divisions=9):

From 961580f969f5516f88c2bc9767b021a2f55e7f1d Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Wed, 5 Aug 2026 14:32:51 -0600
Subject: [PATCH 51/55] Update init example

---
 .../initialization/examples/init_polynomial_ex.py      | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/pyomo/devel/initialization/examples/init_polynomial_ex.py b/pyomo/devel/initialization/examples/init_polynomial_ex.py
index 8afd438e007..7402ab55321 100644
--- a/pyomo/devel/initialization/examples/init_polynomial_ex.py
+++ b/pyomo/devel/initialization/examples/init_polynomial_ex.py
@@ -55,14 +55,14 @@ def global_init_ex():
 
 def multistart_init_ex():
     m = build_model()
-    # m.obj = pyo.Objective(expr=0)
     nlp_solver = SolverFactory('ipopt')
     multistart_solver = SolverFactory('multistart')
 
-    # opts = {"load_solutions": False, "raise_exception_on_nonoptimal_result": False}
-
-    multistart_solver.config.strategy = "rand_vector"
-    # multistart_solver.config.solver_args = opts
+    # multistart_solver.config.strategy = "rand"
+    # multistart_solver.config.strategy = "rand_vector"
+    # multistart_solver.config.sampling_method = "uniform"
+    # multistart_solver.config.sampling_method = "lhs"
+    # multistart_solver.config.sampling_method = "sobol"
     multistart_solver.config.iterations = 10
     multistart_solver.config.break_on_solution = True
 

From 01bd9e04a4d7b980649464ff8a4f722945edf2da Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Wed, 5 Aug 2026 16:07:50 -0600
Subject: [PATCH 52/55] Removing gams from list to see if other checks pass

---
 pyomo/contrib/solver/tests/solvers/test_solvers.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyomo/contrib/solver/tests/solvers/test_solvers.py b/pyomo/contrib/solver/tests/solvers/test_solvers.py
index 9fc95e32d91..5f07455ea6f 100644
--- a/pyomo/contrib/solver/tests/solvers/test_solvers.py
+++ b/pyomo/contrib/solver/tests/solvers/test_solvers.py
@@ -80,7 +80,7 @@ def param_as_standalone_func(cls, p, func, name):
     ('highs', Highs),
     ('scip_direct', ScipDirect),
     ('scip_persistent', ScipPersistent),
-    ('gams', GAMS),
+    # ('gams', GAMS),
     ('knitro_direct', KnitroDirectSolver),
     ('multistart', MultiStart),
 ]

From 00c827ffba8f41857f5fc4c89b3f0049f0624da1 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Wed, 5 Aug 2026 16:16:18 -0600
Subject: [PATCH 53/55] Add back gams, misunderstood

---
 pyomo/contrib/solver/tests/solvers/test_solvers.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyomo/contrib/solver/tests/solvers/test_solvers.py b/pyomo/contrib/solver/tests/solvers/test_solvers.py
index 5f07455ea6f..9fc95e32d91 100644
--- a/pyomo/contrib/solver/tests/solvers/test_solvers.py
+++ b/pyomo/contrib/solver/tests/solvers/test_solvers.py
@@ -80,7 +80,7 @@ def param_as_standalone_func(cls, p, func, name):
     ('highs', Highs),
     ('scip_direct', ScipDirect),
     ('scip_persistent', ScipPersistent),
-    # ('gams', GAMS),
+    ('gams', GAMS),
     ('knitro_direct', KnitroDirectSolver),
     ('multistart', MultiStart),
 ]

From 34504f7dae1338914bf9e21679fdb23668eb5941 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Thu, 6 Aug 2026 08:42:46 -0600
Subject: [PATCH 54/55] Add multistart init ex to testing

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

diff --git a/pyomo/devel/initialization/tests/test_initialization.py b/pyomo/devel/initialization/tests/test_initialization.py
index 042bb890828..7ab8d778123 100644
--- a/pyomo/devel/initialization/tests/test_initialization.py
+++ b/pyomo/devel/initialization/tests/test_initialization.py
@@ -16,6 +16,7 @@
     lp_init_ex,
     pwl_init_ex,
     global_init_ex,
+    multistart_init_ex
 )
 from pyomo.common import unittest
 from pyomo.common.dependencies import scipy_available
@@ -85,6 +86,11 @@ def test_poly_lp(self):
         self.assertEqual(stat, SolutionStatus.optimal)
         self.assertAlmostEqual(x, -9.920159607881597)
 
+    def test_poly_multistart(self):
+        stat, x = lp_init_ex()
+        self.assertEqual(stat, SolutionStatus.optimal)
+        self.assertAlmostEqual(x, -9.920159607881597)
+        
 
 class TestInit(unittest.TestCase):
     @unittest.skipUnless(highs.available(), 'highs is not available')
@@ -229,6 +235,8 @@ def test_pwl_ineq(self):
 
 
 if __name__ == '__main__':
+
+    
     import logging
 
     logging.basicConfig(level=logging.INFO)

From fa26f40ebe2a60d094e52fdb3d034e42cb102d86 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Thu, 6 Aug 2026 09:15:12 -0600
Subject: [PATCH 55/55] Add default False for solution loading, ran black

---
 pyomo/devel/initialization/multistart_init.py           | 5 +++++
 pyomo/devel/initialization/tests/test_initialization.py | 5 ++---
 2 files changed, 7 insertions(+), 3 deletions(-)

diff --git a/pyomo/devel/initialization/multistart_init.py b/pyomo/devel/initialization/multistart_init.py
index 07eafe1b4ce..0f33c3d7e45 100644
--- a/pyomo/devel/initialization/multistart_init.py
+++ b/pyomo/devel/initialization/multistart_init.py
@@ -36,7 +36,12 @@ def _initialize_with_multistart_solver(
     logger.info('fixed variables with equal bounds')
 
     multistart_solver.config.seed = seed
+    multistart_solver.config.load_solutions = False
+    multistart_solver.config.raise_exception_on_nonoptimal_result = False
+
     res = multistart_solver.solve(nlp)
     logger.info('Finished multistart optimization iterations.')
+    if res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
+        res.solution_loader.load_vars()
 
     return res
diff --git a/pyomo/devel/initialization/tests/test_initialization.py b/pyomo/devel/initialization/tests/test_initialization.py
index 7ab8d778123..1d39a9c3dae 100644
--- a/pyomo/devel/initialization/tests/test_initialization.py
+++ b/pyomo/devel/initialization/tests/test_initialization.py
@@ -16,7 +16,7 @@
     lp_init_ex,
     pwl_init_ex,
     global_init_ex,
-    multistart_init_ex
+    multistart_init_ex,
 )
 from pyomo.common import unittest
 from pyomo.common.dependencies import scipy_available
@@ -90,7 +90,7 @@ def test_poly_multistart(self):
         stat, x = lp_init_ex()
         self.assertEqual(stat, SolutionStatus.optimal)
         self.assertAlmostEqual(x, -9.920159607881597)
-        
+
 
 class TestInit(unittest.TestCase):
     @unittest.skipUnless(highs.available(), 'highs is not available')
@@ -236,7 +236,6 @@ def test_pwl_ineq(self):
 
 if __name__ == '__main__':
 
-    
     import logging
 
     logging.basicConfig(level=logging.INFO)