Fix bug in initialization routines to return solve with bounded nlp

pyomo/devel/initialization/global_init.py

@@ -31,24 +31,24 @@ 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,
+        load_solutions=False,
         raise_exception_on_nonoptimal_result=False,
         solver_options=opts,
     )
     logger.info(
         f'solved NLP with {global_solver.name}: {res.solution_status}, {res.termination_condition}'
     )
-    res = nlp_solver.solve(
-        nlp, load_solutions=False, raise_exception_on_nonoptimal_result=False
-    )
-    logger.info(
-        f'solved NLP with {nlp_solver.name}: {res.solution_status}, {res.termination_condition}'
-    )
     if res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
         res.solution_loader.load_vars()
-    else:
-        logger.warning('initialization was not successful via global optimization')
 
     return res

pyomo/devel/initialization/initialize.py

@@ -77,6 +77,21 @@ def _try_nlp_solve(nlp: BlockData, nlp_solver: SolverBase):
     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're-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,
     nlp_solver: SolverBase | None = None,
@@ -138,7 +153,7 @@ def initialize_with_piecewise_linear_approximation(
             return res
 
     if mip_solver is None:
-        mip_solver = _get_solver('gurobi_persistent', 'MILP solver')
+        mip_solver = _get_solver('highs', 'MILP solver')
 
     orig_var_data = _setup(nlp)
 
@@ -156,7 +171,9 @@ def initialize_with_piecewise_linear_approximation(
     finally:
         _cleanup(orig_var_data)
 
-    return res
+    nlp_res = _retry_nlp_solve(nlp, nlp_solver)
+
+    return nlp_res
 
 
 def initialize_with_LP_approximation(
@@ -223,7 +240,7 @@ def initialize_with_LP_approximation(
     orig_var_data = _setup(nlp)
 
     if lp_solver is None:
-        lp_solver = _get_solver('gurobi_persistent', 'LP solver')
+        lp_solver = _get_solver('highs', 'LP solver')
 
     try:
         res = _initialize_with_LP_approximation(
@@ -239,7 +256,9 @@ def initialize_with_LP_approximation(
     finally:
         _cleanup(orig_var_data)
 
-    return res
+    nlp_res = _retry_nlp_solve(nlp, nlp_solver)
+
+    return nlp_res
 
 
 def initialize_with_global_opt(
@@ -284,7 +303,7 @@ def initialize_with_global_opt(
     orig_var_data = _setup(nlp)
 
     if global_solver is None:
-        global_solver = _get_solver('gurobi_direct_minlp', 'global NLP solver')
+        global_solver = _get_solver('scip_direct', 'global NLP solver')
 
     try:
         res = _initialize_with_global_solver(
@@ -293,4 +312,6 @@ def initialize_with_global_opt(
     finally:
         _cleanup(orig_var_data)
 
-    return res
+    nlp_res = _retry_nlp_solve(nlp, nlp_solver)
+
+    return nlp_res

pyomo/devel/initialization/lp_approx_init.py

@@ -214,21 +214,10 @@ def _initialize_with_LP_approximation(
 
     # solve the LP
     lp_res = lp_solver.solve(
-        lp, load_solutions=True, raise_exception_on_nonoptimal_result=False
+        lp, load_solutions=False, raise_exception_on_nonoptimal_result=False
     )
     logger.info(f'solved LP: {lp_res.solution_status}, {lp_res.termination_condition}')
+    if lp_res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
+        lp_res.solution_loader.load_vars()
 
-    # 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

pyomo/devel/initialization/pwl_init.py

@@ -326,9 +326,11 @@ def _initialize_with_piecewise_linear_approximation(
 
         # solve the MILP
         res = mip_solver.solve(
-            _pwl, load_solutions=True, raise_exception_on_nonoptimal_result=False
+            _pwl, load_solutions=False, raise_exception_on_nonoptimal_result=False
         )
         logger.info(f'solved MILP: {res.solution_status}, {res.termination_condition}')
+        if res.solution_status in {SolutionStatus.feasible, SolutionStatus.optimal}:
+            res.solution_loader.load_vars()
 
         # load the variable values back into orig_vars
         for ov, nv in zip(orig_vars, new_vars):

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')
@@ -231,4 +233,3 @@ def test_pwl_ineq(self):
 
     logging.basicConfig(level=logging.INFO)
     t = TestInit()
-    t.test_pwl_init()