diff --git a/pyomo/contrib/mindtpy/MindtPy.py b/pyomo/contrib/mindtpy/MindtPy.py
index 6eb27c4c649..bd873d950fd 100644
--- a/pyomo/contrib/mindtpy/MindtPy.py
+++ b/pyomo/contrib/mindtpy/MindtPy.py
@@ -50,6 +50,14 @@
 - Add single-tree implementation.
 - Add support for cplex_persistent solver.
 - Fix bug in OA cut expression in cut_generation.py.
+
+24.1.11 changes:
+- fix gurobi single tree termination check bug
+- fix Gurobi single tree cycle handling
+- fix bug in feasibility pump method
+- add special handling for infeasible relaxed NLP
+- update the log format of infeasible fixed NLP subproblems
+- create a new copy_var_list_values function
 """
 
 from pyomo.contrib.mindtpy import __version__
diff --git a/pyomo/contrib/mindtpy/__init__.py b/pyomo/contrib/mindtpy/__init__.py
index 8e2c2d9eaa4..8dcd085211f 100644
--- a/pyomo/contrib/mindtpy/__init__.py
+++ b/pyomo/contrib/mindtpy/__init__.py
@@ -1 +1 @@
-__version__ = (0, 1, 0)
+__version__ = (1, 0, 0)
diff --git a/pyomo/contrib/mindtpy/algorithm_base_class.py b/pyomo/contrib/mindtpy/algorithm_base_class.py
index 570e7c0a27d..3d5a7ebad03 100644
--- a/pyomo/contrib/mindtpy/algorithm_base_class.py
+++ b/pyomo/contrib/mindtpy/algorithm_base_class.py
@@ -55,7 +55,6 @@
     SuppressInfeasibleWarning,
     _DoNothing,
     lower_logger_level_to,
-    copy_var_list_values,
     get_main_elapsed_time,
     time_code,
 )
@@ -80,6 +79,7 @@
     set_solver_mipgap,
     set_solver_constraint_violation_tolerance,
     update_solver_timelimit,
+    copy_var_list_values,
 )
 
 single_tree, single_tree_available = attempt_import('pyomo.contrib.mindtpy.single_tree')
@@ -102,12 +102,14 @@ def __init__(self, **kwds):
         self.fixed_nlp = None
 
         # We store bounds, timing info, iteration count, incumbent, and the
-        # expression of the original (possibly nonlinear) objective function.
+        # Expression of the original (possibly nonlinear) objective function.
         self.results = SolverResults()
         self.timing = Bunch()
         self.curr_int_sol = []
         self.should_terminate = False
         self.integer_list = []
+        # Dictionary {integer solution (tuple): [cuts begin index, cuts end index] (list)}
+        self.integer_solution_to_cuts_index = dict()
 
         # Set up iteration counters
         self.nlp_iter = 0
@@ -123,9 +125,15 @@ def __init__(self, **kwds):
         self.log_formatter = (
             ' {:>9}   {:>15}   {:>15g}   {:>12g}   {:>12g}   {:>7.2%}   {:>7.2f}'
         )
+        self.termination_condition_log_formatter = (
+            ' {:>9}   {:>15}   {:>15}   {:>12g}   {:>12g}   {:>7.2%}   {:>7.2f}'
+        )
         self.fixed_nlp_log_formatter = (
             '{:1}{:>9}   {:>15}   {:>15g}   {:>12g}   {:>12g}   {:>7.2%}   {:>7.2f}'
         )
+        self.infeasible_fixed_nlp_log_formatter = (
+            '{:1}{:>9}   {:>15}   {:>15}   {:>12g}   {:>12g}   {:>7.2%}   {:>7.2f}'
+        )
         self.log_note_formatter = ' {:>9}   {:>15}   {:>15}'
 
         # Flag indicating whether the solution improved in the past
@@ -511,9 +519,9 @@ def get_primal_integral(self):
         return primal_integral
 
     def get_integral_info(self):
-        '''
+        """
         Obtain primal integral, dual integral and primal dual gap integral.
-        '''
+        """
         self.primal_integral = self.get_primal_integral()
         self.dual_integral = self.get_dual_integral()
         self.primal_dual_gap_integral = self.primal_integral + self.dual_integral
@@ -797,7 +805,7 @@ def MindtPy_initialization(self):
             try:
                 self.curr_int_sol = get_integer_solution(self.working_model)
             except TypeError as e:
-                config.logger.error(e)
+                config.logger.error(e, exc_info=True)
                 raise ValueError(
                     'The initial integer combination is not provided or not complete. '
                     'Please provide the complete integer combination or use other initialization strategy.'
@@ -805,6 +813,10 @@ def MindtPy_initialization(self):
             self.integer_list.append(self.curr_int_sol)
             fixed_nlp, fixed_nlp_result = self.solve_subproblem()
             self.handle_nlp_subproblem_tc(fixed_nlp, fixed_nlp_result)
+            self.integer_solution_to_cuts_index[self.curr_int_sol] = [
+                1,
+                len(self.mip.MindtPy_utils.cuts.oa_cuts),
+            ]
         elif config.init_strategy == 'FP':
             self.init_rNLP()
             self.fp_loop()
@@ -840,6 +852,29 @@ def init_rNLP(self, add_oa_cuts=True):
             if len(results.solution) > 0:
                 self.rnlp.solutions.load_from(results)
         subprob_terminate_cond = results.solver.termination_condition
+
+        # Sometimes, the NLP solver might be trapped in a infeasible solution if the objective function is nonlinear and partition_obj_nonlinear_terms is True. If this happens, we will use the original objective function instead.
+        if (
+            subprob_terminate_cond == tc.infeasible
+            and config.partition_obj_nonlinear_terms
+            and self.rnlp.MindtPy_utils.objective_list[0].expr.polynomial_degree()
+            not in self.mip_objective_polynomial_degree
+        ):
+            config.logger.info(
+                'Initial relaxed NLP problem is infeasible. This might be related to partition_obj_nonlinear_terms. Trying to solve it again without partitioning nonlinear objective function.'
+            )
+            self.rnlp.MindtPy_utils.objective.deactivate()
+            self.rnlp.MindtPy_utils.objective_list[0].activate()
+            results = self.nlp_opt.solve(
+                self.rnlp,
+                tee=config.nlp_solver_tee,
+                load_solutions=self.load_solutions,
+                **nlp_args,
+            )
+            if len(results.solution) > 0:
+                self.rnlp.solutions.load_from(results)
+            subprob_terminate_cond = results.solver.termination_condition
+
         if subprob_terminate_cond in {tc.optimal, tc.feasible, tc.locallyOptimal}:
             main_objective = MindtPy.objective_list[-1]
             if subprob_terminate_cond == tc.optimal:
@@ -880,12 +915,14 @@ def init_rNLP(self, add_oa_cuts=True):
                     self.rnlp.MindtPy_utils.variable_list,
                     self.mip.MindtPy_utils.variable_list,
                     config,
+                    ignore_integrality=True,
                 )
                 if config.init_strategy == 'FP':
                     copy_var_list_values(
                         self.rnlp.MindtPy_utils.variable_list,
                         self.working_model.MindtPy_utils.variable_list,
                         config,
+                        ignore_integrality=True,
                     )
                 self.add_cuts(
                     dual_values=dual_values,
@@ -1050,7 +1087,7 @@ def solve_subproblem(self):
                         0, c_geq * (rhs - value(c.body))
                     )
                 except (ValueError, OverflowError) as e:
-                    config.logger.error(e)
+                    config.logger.error(e, exc_info=True)
                     self.fixed_nlp.tmp_duals[c] = None
                     evaluation_error = True
             if evaluation_error:
@@ -1067,8 +1104,9 @@ def solve_subproblem(self):
                 tolerance=config.constraint_tolerance,
             )
         except InfeasibleConstraintException as e:
+            config.logger.error(e, exc_info=True)
             config.logger.error(
-                str(e) + '\nInfeasibility detected in deactivate_trivial_constraints.'
+                'Infeasibility detected in deactivate_trivial_constraints.'
             )
             results = SolverResults()
             results.solver.termination_condition = tc.infeasible
@@ -1219,7 +1257,18 @@ def handle_subproblem_infeasible(self, fixed_nlp, cb_opt=None):
         # TODO try something else? Reinitialize with different initial
         # value?
         config = self.config
-        config.logger.info('NLP subproblem was locally infeasible.')
+        config.logger.info(
+            self.infeasible_fixed_nlp_log_formatter.format(
+                ' ',
+                self.nlp_iter,
+                'Fixed NLP',
+                'Infeasible',
+                self.primal_bound,
+                self.dual_bound,
+                self.rel_gap,
+                get_main_elapsed_time(self.timing),
+            )
+        )
         self.nlp_infeasible_counter += 1
         if config.calculate_dual_at_solution:
             for c in fixed_nlp.MindtPy_utils.constraint_list:
@@ -1241,7 +1290,6 @@ def handle_subproblem_infeasible(self, fixed_nlp, cb_opt=None):
         #         elif var.has_lb() and abs(value(var) - var.lb) < config.absolute_bound_tolerance:
         #             fixed_nlp.ipopt_zU_out[var] = -1
 
-        config.logger.info('Solving feasibility problem')
         feas_subproblem, feas_subproblem_results = self.solve_feasibility_subproblem()
         # TODO: do we really need this?
         if self.should_terminate:
@@ -1357,7 +1405,7 @@ def solve_feasibility_subproblem(self):
                     if len(feas_soln.solution) > 0:
                         feas_subproblem.solutions.load_from(feas_soln)
             except (ValueError, OverflowError) as e:
-                config.logger.error(e)
+                config.logger.error(e, exc_info=True)
                 for nlp_var, orig_val in zip(
                     MindtPy.variable_list, self.initial_var_values
                 ):
@@ -1375,6 +1423,18 @@ def solve_feasibility_subproblem(self):
         self.handle_feasibility_subproblem_tc(
             feas_soln.solver.termination_condition, MindtPy
         )
+        config.logger.info(
+            self.fixed_nlp_log_formatter.format(
+                ' ',
+                self.nlp_iter,
+                'Feasibility NLP',
+                value(feas_subproblem.MindtPy_utils.feas_obj),
+                self.primal_bound,
+                self.dual_bound,
+                self.rel_gap,
+                get_main_elapsed_time(self.timing),
+            )
+        )
         MindtPy.feas_opt.deactivate()
         for constr in MindtPy.nonlinear_constraint_list:
             constr.activate()
@@ -1486,9 +1546,8 @@ def fix_dual_bound(self, last_iter_cuts):
             try:
                 self.dual_bound = self.stored_bound[self.primal_bound]
             except KeyError as e:
-                config.logger.error(
-                    str(e) + '\nNo stored bound found. Bound fix failed.'
-                )
+                config.logger.error(e, exc_info=True)
+                config.logger.error('No stored bound found. Bound fix failed.')
         else:
             config.logger.info(
                 'Solve the main problem without the last no_good cut to fix the bound.'
@@ -1502,7 +1561,7 @@ def fix_dual_bound(self, last_iter_cuts):
                 self.handle_nlp_subproblem_tc(fixed_nlp, fixed_nlp_result)
 
             MindtPy = self.mip.MindtPy_utils
-            # deactivate the integer cuts generated after the best solution was found.
+            # Deactivate the integer cuts generated after the best solution was found.
             self.deactivate_no_good_cuts_when_fixing_bound(MindtPy.cuts.no_good_cuts)
             if (
                 config.add_regularization is not None
@@ -1601,6 +1660,7 @@ def solve_main(self):
         # setup main problem
         self.setup_main()
         mip_args = self.set_up_mip_solver()
+        update_solver_timelimit(self.mip_opt, config.mip_solver, self.timing, config)
 
         try:
             main_mip_results = self.mip_opt.solve(
@@ -1613,7 +1673,7 @@ def solve_main(self):
             if len(main_mip_results.solution) > 0:
                 self.mip.solutions.load_from(main_mip_results)
         except (ValueError, AttributeError, RuntimeError) as e:
-            config.logger.error(e)
+            config.logger.error(e, exc_info=True)
             if config.single_tree:
                 config.logger.warning('Single tree terminate.')
                 if get_main_elapsed_time(self.timing) >= config.time_limit:
@@ -1626,7 +1686,11 @@ def solve_main(self):
                         "No-good cuts are added and GOA algorithm doesn't converge within the time limit. "
                         'No integer solution is found, so the CPLEX solver will report an error status. '
                     )
-            return None, None
+            # Value error will be raised if the MIP problem is unbounded and appsi solver is used when loading solutions. Although the problem is unbounded, a valid result is provided and we do not return None to let the algorithm continue.
+            if 'main_mip_results' in locals():
+                return self.mip, main_mip_results
+            else:
+                return None, None
         if config.solution_pool:
             main_mip_results._solver_model = self.mip_opt._solver_model
             main_mip_results._pyomo_var_to_solver_var_map = (
@@ -1658,6 +1722,7 @@ def solve_fp_main(self):
         config = self.config
         self.setup_fp_main()
         mip_args = self.set_up_mip_solver()
+        update_solver_timelimit(self.mip_opt, config.mip_solver, self.timing, config)
 
         main_mip_results = self.mip_opt.solve(
             self.mip,
@@ -1791,7 +1856,6 @@ def handle_main_optimal(self, main_mip, update_bound=True):
                     f"Integer variable {var.name} not initialized.  "
                     "Setting it to its lower bound"
                 )
-                # nlp_var.bounds[0]
                 var.set_value(var.lb, skip_validation=True)
         # warm start for the nlp subproblem
         copy_var_list_values(
@@ -1857,11 +1921,6 @@ def handle_main_max_timelimit(self, main_mip, main_mip_results):
         """
         # If we have found a valid feasible solution, we take that. If not, we can at least use the dual bound.
         MindtPy = main_mip.MindtPy_utils
-        self.config.logger.info(
-            'Unable to optimize MILP main problem '
-            'within time limit. '
-            'Using current solver feasible solution.'
-        )
         copy_var_list_values(
             main_mip.MindtPy_utils.variable_list,
             self.fixed_nlp.MindtPy_utils.variable_list,
@@ -1870,10 +1929,10 @@ def handle_main_max_timelimit(self, main_mip, main_mip_results):
         )
         self.update_suboptimal_dual_bound(main_mip_results)
         self.config.logger.info(
-            self.log_formatter.format(
+            self.termination_condition_log_formatter.format(
                 self.mip_iter,
                 'MILP',
-                value(MindtPy.mip_obj.expr),
+                'maxTimeLimit',
                 self.primal_bound,
                 self.dual_bound,
                 self.rel_gap,
@@ -1900,8 +1959,18 @@ def handle_main_unbounded(self, main_mip):
         # to the constraints, and deactivated for the linear main problem.
         config = self.config
         MindtPy = main_mip.MindtPy_utils
+        config.logger.info(
+            self.termination_condition_log_formatter.format(
+                self.mip_iter,
+                'MILP',
+                'Unbounded',
+                self.primal_bound,
+                self.dual_bound,
+                self.rel_gap,
+                get_main_elapsed_time(self.timing),
+            )
+        )
         config.logger.warning(
-            'main MILP was unbounded. '
             'Resolving with arbitrary bound values of (-{0:.10g}, {0:.10g}) on the objective. '
             'You can change this bound with the option obj_bound.'.format(
                 config.obj_bound
@@ -2308,8 +2377,9 @@ def solve_fp_subproblem(self):
                 tolerance=config.constraint_tolerance,
             )
         except InfeasibleConstraintException as e:
+            config.logger.error(e, exc_info=True)
             config.logger.error(
-                str(e) + '\nInfeasibility detected in deactivate_trivial_constraints.'
+                'Infeasibility detected in deactivate_trivial_constraints.'
             )
             results = SolverResults()
             results.solver.termination_condition = tc.infeasible
@@ -2342,6 +2412,7 @@ def handle_fp_subproblem_optimal(self, fp_nlp):
             fp_nlp.MindtPy_utils.variable_list,
             self.working_model.MindtPy_utils.variable_list,
             self.config,
+            ignore_integrality=True,
         )
         add_orthogonality_cuts(self.working_model, self.mip, self.config)
 
@@ -2526,7 +2597,7 @@ def fp_loop(self):
             self.working_model.MindtPy_utils.cuts.del_component('fp_orthogonality_cuts')
 
     def initialize_mip_problem(self):
-        '''Deactivate the nonlinear constraints to create the MIP problem.'''
+        """Deactivate the nonlinear constraints to create the MIP problem."""
         # if single tree is activated, we need to add bounds for unbounded variables in nonlinear constraints to avoid unbounded main problem.
         config = self.config
         if config.single_tree:
@@ -2557,7 +2628,7 @@ def initialize_mip_problem(self):
 
         self.fixed_nlp = self.working_model.clone()
         TransformationFactory('core.fix_integer_vars').apply_to(self.fixed_nlp)
-        initialize_feas_subproblem(self.fixed_nlp, config)
+        initialize_feas_subproblem(self.fixed_nlp, config.feasibility_norm)
 
     def initialize_subsolvers(self):
         """Initialize and set options for MIP and NLP subsolvers."""
@@ -2585,7 +2656,7 @@ def initialize_subsolvers(self):
             self.nlp_opt, config.nlp_solver, config
         )
         set_solver_constraint_violation_tolerance(
-            self.feasibility_nlp_opt, config.nlp_solver, config
+            self.feasibility_nlp_opt, config.nlp_solver, config, warm_start=False
         )
 
         self.set_appsi_solver_update_config()
@@ -2942,10 +3013,12 @@ def MindtPy_iteration_loop(self):
 
         # if add_no_good_cuts is True, the bound obtained in the last iteration is no reliable.
         # we correct it after the iteration.
+        # There is no need to fix the dual bound if no feasible solution has been found.
         if (
             (config.add_no_good_cuts or config.use_tabu_list)
             and not self.should_terminate
             and config.add_regularization is None
+            and self.best_solution_found is not None
         ):
             self.fix_dual_bound(self.last_iter_cuts)
         config.logger.info(
diff --git a/pyomo/contrib/mindtpy/cut_generation.py b/pyomo/contrib/mindtpy/cut_generation.py
index 28d302104a3..4ee7a6ff07b 100644
--- a/pyomo/contrib/mindtpy/cut_generation.py
+++ b/pyomo/contrib/mindtpy/cut_generation.py
@@ -210,8 +210,8 @@ def add_oa_cuts_for_grey_box(
                             target_model_grey_box.inputs.values()
                         )
                     )
+                    - (output - value(output))
                 )
-                - (output - value(output))
                 - (slack_var if config.add_slack else 0)
                 <= 0
             )
@@ -271,8 +271,9 @@ def add_ecp_cuts(
                 try:
                     upper_slack = constr.uslack()
                 except (ValueError, OverflowError) as e:
+                    config.logger.error(e, exc_info=True)
                     config.logger.error(
-                        str(e) + '\nConstraint {} has caused either a '
+                        'Constraint {} has caused either a '
                         'ValueError or OverflowError.'
                         '\n'.format(constr)
                     )
@@ -300,8 +301,9 @@ def add_ecp_cuts(
                 try:
                     lower_slack = constr.lslack()
                 except (ValueError, OverflowError) as e:
+                    config.logger.error(e, exc_info=True)
                     config.logger.error(
-                        str(e) + '\nConstraint {} has caused either a '
+                        'Constraint {} has caused either a '
                         'ValueError or OverflowError.'
                         '\n'.format(constr)
                     )
@@ -424,9 +426,9 @@ def add_affine_cuts(target_model, config, timing):
             try:
                 mc_eqn = mc(constr.body)
             except MCPP_Error as e:
+                config.logger.error(e, exc_info=True)
                 config.logger.error(
-                    '\nSkipping constraint %s due to MCPP error %s'
-                    % (constr.name, str(e))
+                    'Skipping constraint %s due to MCPP error' % (constr.name)
                 )
                 continue  # skip to the next constraint
 
diff --git a/pyomo/contrib/mindtpy/extended_cutting_plane.py b/pyomo/contrib/mindtpy/extended_cutting_plane.py
index 446304b1361..0a98f88ed3f 100644
--- a/pyomo/contrib/mindtpy/extended_cutting_plane.py
+++ b/pyomo/contrib/mindtpy/extended_cutting_plane.py
@@ -66,12 +66,6 @@ def MindtPy_iteration_loop(self):
 
             add_ecp_cuts(self.mip, self.jacobians, self.config, self.timing)
 
-        # if add_no_good_cuts is True, the bound obtained in the last iteration is no reliable.
-        # we correct it after the iteration.
-        if (
-            self.config.add_no_good_cuts or self.config.use_tabu_list
-        ) and not self.should_terminate:
-            self.fix_dual_bound(self.last_iter_cuts)
         self.config.logger.info(
             ' ==============================================================================================='
         )
@@ -84,9 +78,12 @@ def check_config(self):
         super().check_config()
 
     def initialize_mip_problem(self):
-        '''Deactivate the nonlinear constraints to create the MIP problem.'''
+        """Deactivate the nonlinear constraints to create the MIP problem."""
         super().initialize_mip_problem()
-        self.jacobians = calc_jacobians(self.mip, self.config)  # preload jacobians
+        self.jacobians = calc_jacobians(
+            self.mip.MindtPy_utils.nonlinear_constraint_list,
+            self.config.differentiate_mode,
+        )  # preload jacobians
         self.mip.MindtPy_utils.cuts.ecp_cuts = ConstraintList(
             doc='Extended Cutting Planes'
         )
@@ -140,7 +137,7 @@ def all_nonlinear_constraint_satisfied(self):
                     lower_slack = nlc.lslack()
                 except (ValueError, OverflowError) as e:
                     # Set lower_slack (upper_slack below) less than -config.ecp_tolerance in this case.
-                    config.logger.error(e)
+                    config.logger.error(e, exc_info=True)
                     lower_slack = -10 * config.ecp_tolerance
                 if lower_slack < -config.ecp_tolerance:
                     config.logger.debug(
@@ -153,7 +150,7 @@ def all_nonlinear_constraint_satisfied(self):
                 try:
                     upper_slack = nlc.uslack()
                 except (ValueError, OverflowError) as e:
-                    config.logger.error(e)
+                    config.logger.error(e, exc_info=True)
                     upper_slack = -10 * config.ecp_tolerance
                 if upper_slack < -config.ecp_tolerance:
                     config.logger.debug(
diff --git a/pyomo/contrib/mindtpy/feasibility_pump.py b/pyomo/contrib/mindtpy/feasibility_pump.py
index 990f56b8f93..a34cceb014c 100644
--- a/pyomo/contrib/mindtpy/feasibility_pump.py
+++ b/pyomo/contrib/mindtpy/feasibility_pump.py
@@ -44,9 +44,12 @@ def check_config(self):
         super().check_config()
 
     def initialize_mip_problem(self):
-        '''Deactivate the nonlinear constraints to create the MIP problem.'''
+        """Deactivate the nonlinear constraints to create the MIP problem."""
         super().initialize_mip_problem()
-        self.jacobians = calc_jacobians(self.mip, self.config)  # preload jacobians
+        self.jacobians = calc_jacobians(
+            self.mip.MindtPy_utils.nonlinear_constraint_list,
+            self.config.differentiate_mode,
+        )  # preload jacobians
         self.mip.MindtPy_utils.cuts.oa_cuts = ConstraintList(
             doc='Outer approximation cuts'
         )
diff --git a/pyomo/contrib/mindtpy/global_outer_approximation.py b/pyomo/contrib/mindtpy/global_outer_approximation.py
index dfb7ef54630..70fc4cffb90 100644
--- a/pyomo/contrib/mindtpy/global_outer_approximation.py
+++ b/pyomo/contrib/mindtpy/global_outer_approximation.py
@@ -67,7 +67,7 @@ def check_config(self):
         super().check_config()
 
     def initialize_mip_problem(self):
-        '''Deactivate the nonlinear constraints to create the MIP problem.'''
+        """Deactivate the nonlinear constraints to create the MIP problem."""
         super().initialize_mip_problem()
         self.mip.MindtPy_utils.cuts.aff_cuts = ConstraintList(doc='Affine cuts')
 
@@ -108,4 +108,5 @@ def deactivate_no_good_cuts_when_fixing_bound(self, no_good_cuts):
             if self.config.use_tabu_list:
                 self.integer_list = self.integer_list[:valid_no_good_cuts_num]
         except KeyError as e:
-            self.config.logger.error(str(e) + '\nDeactivating no-good cuts failed.')
+            self.config.logger.error(e, exc_info=True)
+            self.config.logger.error('Deactivating no-good cuts failed.')
diff --git a/pyomo/contrib/mindtpy/outer_approximation.py b/pyomo/contrib/mindtpy/outer_approximation.py
index 6cf0b26cb37..f6e6147724e 100644
--- a/pyomo/contrib/mindtpy/outer_approximation.py
+++ b/pyomo/contrib/mindtpy/outer_approximation.py
@@ -94,9 +94,12 @@ def check_config(self):
         _MindtPyAlgorithm.check_config(self)
 
     def initialize_mip_problem(self):
-        '''Deactivate the nonlinear constraints to create the MIP problem.'''
+        """Deactivate the nonlinear constraints to create the MIP problem."""
         super().initialize_mip_problem()
-        self.jacobians = calc_jacobians(self.mip, self.config)  # preload jacobians
+        self.jacobians = calc_jacobians(
+            self.mip.MindtPy_utils.nonlinear_constraint_list,
+            self.config.differentiate_mode,
+        )  # preload jacobians
         self.mip.MindtPy_utils.cuts.oa_cuts = ConstraintList(
             doc='Outer approximation cuts'
         )
diff --git a/pyomo/contrib/mindtpy/single_tree.py b/pyomo/contrib/mindtpy/single_tree.py
index 5383624b6aa..c1e52ed72d3 100644
--- a/pyomo/contrib/mindtpy/single_tree.py
+++ b/pyomo/contrib/mindtpy/single_tree.py
@@ -16,12 +16,12 @@
 from pyomo.repn import generate_standard_repn
 import pyomo.core.expr as EXPR
 from math import copysign
-from pyomo.contrib.mindtpy.util import get_integer_solution
-from pyomo.contrib.gdpopt.util import (
+from pyomo.contrib.mindtpy.util import (
+    get_integer_solution,
     copy_var_list_values,
-    get_main_elapsed_time,
-    time_code,
+    set_var_valid_value,
 )
+from pyomo.contrib.gdpopt.util import get_main_elapsed_time, time_code
 from pyomo.opt import TerminationCondition as tc
 from pyomo.core import minimize, value
 from pyomo.core.expr import identify_variables
@@ -35,17 +35,9 @@ class LazyOACallback_cplex(
     """Inherent class in CPLEX to call Lazy callback."""
 
     def copy_lazy_var_list_values(
-        self,
-        opt,
-        from_list,
-        to_list,
-        config,
-        skip_stale=False,
-        skip_fixed=True,
-        ignore_integrality=False,
+        self, opt, from_list, to_list, config, skip_stale=False, skip_fixed=True
     ):
         """This function copies variable values from one list to another.
-
         Rounds to Binary/Integer if necessary.
         Sets to zero for NonNegativeReals if necessary.
 
@@ -54,17 +46,15 @@ def copy_lazy_var_list_values(
         opt : SolverFactory
             The cplex_persistent solver.
         from_list : list
-            The variables that provides the values to copy from.
+            The variable list that provides the values to copy from.
         to_list : list
-            The variables that need to set value.
+            The variable list that needs to set value.
         config : ConfigBlock
             The specific configurations for MindtPy.
         skip_stale : bool, optional
             Whether to skip the stale variables, by default False.
         skip_fixed : bool, optional
             Whether to skip the fixed variables, by default True.
-        ignore_integrality : bool, optional
-            Whether to ignore the integrality of integer variables, by default False.
         """
         for v_from, v_to in zip(from_list, to_list):
             if skip_stale and v_from.stale:
@@ -72,43 +62,13 @@ def copy_lazy_var_list_values(
             if skip_fixed and v_to.is_fixed():
                 continue  # Skip fixed variables.
             v_val = self.get_values(opt._pyomo_var_to_solver_var_map[v_from])
-            try:
-                # We don't want to trigger the reset of the global stale
-                # indicator, so we will set this variable to be "stale",
-                # knowing that set_value will switch it back to "not
-                # stale"
-                v_to.stale = True
-                # NOTE: PEP 2180 changes the var behavior so that domain
-                # / bounds violations no longer generate exceptions (and
-                # instead log warnings).  This means that the following
-                # will always succeed and the ValueError should never be
-                # raised.
-                v_to.set_value(v_val, skip_validation=True)
-            except ValueError as e:
-                # Snap the value to the bounds
-                config.logger.error(e)
-                if (
-                    v_to.has_lb()
-                    and v_val < v_to.lb
-                    and v_to.lb - v_val <= config.variable_tolerance
-                ):
-                    v_to.set_value(v_to.lb, skip_validation=True)
-                elif (
-                    v_to.has_ub()
-                    and v_val > v_to.ub
-                    and v_val - v_to.ub <= config.variable_tolerance
-                ):
-                    v_to.set_value(v_to.ub, skip_validation=True)
-                # ... or the nearest integer
-                elif v_to.is_integer():
-                    rounded_val = int(round(v_val))
-                    if (
-                        ignore_integrality
-                        or abs(v_val - rounded_val) <= config.integer_tolerance
-                    ) and rounded_val in v_to.domain:
-                        v_to.set_value(rounded_val, skip_validation=True)
-                else:
-                    raise
+            set_var_valid_value(
+                v_to,
+                v_val,
+                config.integer_tolerance,
+                config.zero_tolerance,
+                ignore_integrality=False,
+            )
 
     def add_lazy_oa_cuts(
         self,
@@ -309,12 +269,11 @@ def add_lazy_affine_cuts(self, mindtpy_solver, config, opt):
                 try:
                     mc_eqn = mc(constr.body)
                 except MCPP_Error as e:
+                    config.logger.error(e, exc_info=True)
                     config.logger.debug(
-                        'Skipping constraint %s due to MCPP error %s'
-                        % (constr.name, str(e))
+                        'Skipping constraint %s due to MCPP error' % (constr.name)
                     )
                     continue  # skip to the next constraint
-                # TODO: check if the value of ccSlope and cvSlope is not Nan or inf. If so, we skip this.
                 ccSlope = mc_eqn.subcc()
                 cvSlope = mc_eqn.subcv()
                 ccStart = mc_eqn.concave()
@@ -705,10 +664,11 @@ def __call__(self):
         main_mip = self.main_mip
         mindtpy_solver = self.mindtpy_solver
 
+        # The lazy constraint callback may be invoked during MIP start processing. In that case get_solution_source returns mip_start_solution.
         # Reference: https://www.ibm.com/docs/en/icos/22.1.1?topic=SSSA5P_22.1.1/ilog.odms.cplex.help/refpythoncplex/html/cplex.callbacks.SolutionSource-class.htm
         # Another solution source is user_solution = 118, but it will not be encountered in LazyConstraintCallback.
-        config.logger.debug(
-            "Solution source: %s (111 node_solution, 117 heuristic_solution, 119 mipstart_solution)".format(
+        config.logger.info(
+            "Solution source: {} (111 node_solution, 117 heuristic_solution, 119 mipstart_solution)".format(
                 self.get_solution_source()
             )
         )
@@ -717,6 +677,7 @@ def __call__(self):
         # Lazy constraints separated when processing a MIP start will be discarded after that MIP start has been processed.
         # This means that the callback may have to separate the same constraint again for the next MIP start or for a solution that is found later in the solution process.
         # https://www.ibm.com/docs/en/icos/22.1.1?topic=SSSA5P_22.1.1/ilog.odms.cplex.help/refpythoncplex/html/cplex.callbacks.LazyConstraintCallback-class.htm
+        # For the MINLP3_simple example, all the solutions are obtained from mip_start (solution source). Therefore, it will not go to a branch and bound process.Cause an error output.
         if (
             self.get_solution_source()
             != cplex.callbacks.SolutionSource.mipstart_solution
@@ -727,6 +688,7 @@ def __call__(self):
             mindtpy_solver.mip_start_lazy_oa_cuts = []
 
         if mindtpy_solver.should_terminate:
+            # TODO: check the performance difference if we don't use self.abort() and let cplex terminate by itself.
             self.abort()
             return
         self.handle_lazy_main_feasible_solution(main_mip, mindtpy_solver, config, opt)
@@ -744,9 +706,9 @@ def __call__(self):
                         mindtpy_solver.mip, None, mindtpy_solver, config, opt
                     )
                 except ValueError as e:
+                    config.logger.error(e, exc_info=True)
                     config.logger.error(
-                        str(e)
-                        + "\nUsually this error is caused by the MIP start solution causing a math domain error. "
+                        "Usually this error is caused by the MIP start solution causing a math domain error. "
                         "We will skip it."
                     )
                     return
@@ -782,6 +744,7 @@ def __call__(self):
                 )
             )
             mindtpy_solver.results.solver.termination_condition = tc.optimal
+            # TODO: check the performance difference if we don't use self.abort() and let cplex terminate by itself.
             self.abort()
             return
 
@@ -909,19 +872,7 @@ def LazyOACallback_gurobi(cb_m, cb_opt, cb_where, mindtpy_solver, config):
             if mindtpy_solver.dual_bound != mindtpy_solver.dual_bound_progress[0]:
                 mindtpy_solver.add_regularization()
 
-        if (
-            abs(mindtpy_solver.primal_bound - mindtpy_solver.dual_bound)
-            <= config.absolute_bound_tolerance
-        ):
-            config.logger.info(
-                'MindtPy exiting on bound convergence. '
-                '|Primal Bound: {} - Dual Bound: {}| <= (absolute tolerance {})  \n'.format(
-                    mindtpy_solver.primal_bound,
-                    mindtpy_solver.dual_bound,
-                    config.absolute_bound_tolerance,
-                )
-            )
-            mindtpy_solver.results.solver.termination_condition = tc.optimal
+        if mindtpy_solver.bounds_converged() or mindtpy_solver.reached_time_limit():
             cb_opt._solver_model.terminate()
             return
 
@@ -952,15 +903,31 @@ def LazyOACallback_gurobi(cb_m, cb_opt, cb_where, mindtpy_solver, config):
                     )
                 return
             elif config.strategy == 'OA':
+                # Refer to the official document of GUROBI.
+                # Your callback should be prepared to cut off solutions that violate any of your lazy constraints, including those that have already been added. Node solutions will usually respect previously added lazy constraints, but not always.
+                # https://www.gurobi.com/documentation/current/refman/cs_cb_addlazy.html
+                # If this happens, MindtPy will look for the index of corresponding cuts, instead of solving the fixed-NLP again.
+                begin_index, end_index = mindtpy_solver.integer_solution_to_cuts_index[
+                    mindtpy_solver.curr_int_sol
+                ]
+                for ind in range(begin_index, end_index + 1):
+                    cb_opt.cbLazy(mindtpy_solver.mip.MindtPy_utils.cuts.oa_cuts[ind])
                 return
         else:
             mindtpy_solver.integer_list.append(mindtpy_solver.curr_int_sol)
+            if config.strategy == 'OA':
+                cut_ind = len(mindtpy_solver.mip.MindtPy_utils.cuts.oa_cuts)
 
         # solve subproblem
         # The constraint linearization happens in the handlers
         fixed_nlp, fixed_nlp_result = mindtpy_solver.solve_subproblem()
 
         mindtpy_solver.handle_nlp_subproblem_tc(fixed_nlp, fixed_nlp_result, cb_opt)
+        if config.strategy == 'OA':
+            # store the cut index corresponding to current integer solution.
+            mindtpy_solver.integer_solution_to_cuts_index[
+                mindtpy_solver.curr_int_sol
+            ] = [cut_ind + 1, len(mindtpy_solver.mip.MindtPy_utils.cuts.oa_cuts)]
 
 
 def handle_lazy_main_feasible_solution_gurobi(cb_m, cb_opt, mindtpy_solver, config):
diff --git a/pyomo/contrib/mindtpy/tests/MINLP_simple_grey_box.py b/pyomo/contrib/mindtpy/tests/MINLP_simple_grey_box.py
index 547efc0a74c..9c1f33e80cc 100644
--- a/pyomo/contrib/mindtpy/tests/MINLP_simple_grey_box.py
+++ b/pyomo/contrib/mindtpy/tests/MINLP_simple_grey_box.py
@@ -114,7 +114,7 @@ def evaluate_jacobian_equality_constraints(self):
             """Evaluate the Jacobian of the equality constraints."""
             return None
 
-        '''
+        """
         def _extract_and_assemble_fim(self):
             M = np.zeros((self.n_parameters, self.n_parameters))
             for i in range(self.n_parameters):
@@ -122,7 +122,7 @@ def _extract_and_assemble_fim(self):
                     M[i,k] = self._input_values[self.ele_to_order[(i,k)]]
 
             return M
-        '''
+        """
 
         def evaluate_jacobian_outputs(self):
             """Evaluate the Jacobian of the outputs."""
diff --git a/pyomo/contrib/mindtpy/tests/test_mindtpy.py b/pyomo/contrib/mindtpy/tests/test_mindtpy.py
index e872eccc670..ae531f9bd84 100644
--- a/pyomo/contrib/mindtpy/tests/test_mindtpy.py
+++ b/pyomo/contrib/mindtpy/tests/test_mindtpy.py
@@ -327,6 +327,7 @@ def test_OA_APPSI_ipopt(self):
                     value(model.objective.expr), model.optimal_value, places=1
                 )
 
+    # CYIPOPT will raise WARNING (W1002) during loading solution.
     @unittest.skipUnless(
         SolverFactory('cyipopt').available(exception_flag=False),
         "APPSI_IPOPT not available.",
diff --git a/pyomo/contrib/mindtpy/tests/test_mindtpy_feas_pump.py b/pyomo/contrib/mindtpy/tests/test_mindtpy_feas_pump.py
index 697a63d17c8..dcb5c4bce75 100644
--- a/pyomo/contrib/mindtpy/tests/test_mindtpy_feas_pump.py
+++ b/pyomo/contrib/mindtpy/tests/test_mindtpy_feas_pump.py
@@ -17,7 +17,7 @@
 from pyomo.contrib.mindtpy.tests.feasibility_pump1 import FeasPump1
 from pyomo.contrib.mindtpy.tests.feasibility_pump2 import FeasPump2
 
-required_solvers = ('ipopt', 'cplex')
+required_solvers = ('ipopt', 'glpk')
 # TODO: 'appsi_highs' will fail here.
 if all(SolverFactory(s).available(exception_flag=False) for s in required_solvers):
     subsolvers_available = True
@@ -69,6 +69,22 @@ def test_FP(self):
                 log_infeasible_constraints(model)
                 self.assertTrue(is_feasible(model, self.get_config(opt)))
 
+    def test_FP_L1_norm(self):
+        """Test the feasibility pump algorithm."""
+        with SolverFactory('mindtpy') as opt:
+            for model in model_list:
+                model = model.clone()
+                results = opt.solve(
+                    model,
+                    strategy='FP',
+                    mip_solver=required_solvers[1],
+                    nlp_solver=required_solvers[0],
+                    absolute_bound_tolerance=1e-5,
+                    fp_main_norm='L1',
+                )
+                log_infeasible_constraints(model)
+                self.assertTrue(is_feasible(model, self.get_config(opt)))
+
     def test_FP_OA_8PP(self):
         """Test the FP-OA algorithm."""
         with SolverFactory('mindtpy') as opt:
diff --git a/pyomo/contrib/mindtpy/tests/unit_test.py b/pyomo/contrib/mindtpy/tests/unit_test.py
new file mode 100644
index 00000000000..a1ceadda41e
--- /dev/null
+++ b/pyomo/contrib/mindtpy/tests/unit_test.py
@@ -0,0 +1,101 @@
+#  ___________________________________________________________________________
+#
+#  Pyomo: Python Optimization Modeling Objects
+#  Copyright (c) 2008-2022
+#  National Technology and Engineering Solutions of Sandia, LLC
+#  Under the terms of Contract DE-NA0003525 with National Technology and
+#  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
+#  rights in this software.
+#  This software is distributed under the 3-clause BSD License.
+#  ___________________________________________________________________________
+
+import pyomo.common.unittest as unittest
+from pyomo.contrib.mindtpy.util import set_var_valid_value
+
+from pyomo.environ import Var, Integers, ConcreteModel, Integers
+from pyomo.contrib.mindtpy.algorithm_base_class import _MindtPyAlgorithm
+from pyomo.contrib.mindtpy.config_options import _get_MindtPy_OA_config
+from pyomo.contrib.mindtpy.tests.MINLP5_simple import SimpleMINLP5
+from pyomo.contrib.mindtpy.util import add_var_bound
+
+
+class UnitTestMindtPy(unittest.TestCase):
+    def test_set_var_valid_value(self):
+        m = ConcreteModel()
+        m.x1 = Var(within=Integers, bounds=(-1, 4), initialize=0)
+
+        set_var_valid_value(
+            m.x1,
+            var_val=5,
+            integer_tolerance=1e-6,
+            zero_tolerance=1e-6,
+            ignore_integrality=False,
+        )
+        self.assertEqual(m.x1.value, 4)
+
+        set_var_valid_value(
+            m.x1,
+            var_val=-2,
+            integer_tolerance=1e-6,
+            zero_tolerance=1e-6,
+            ignore_integrality=False,
+        )
+        self.assertEqual(m.x1.value, -1)
+
+        set_var_valid_value(
+            m.x1,
+            var_val=1.1,
+            integer_tolerance=1e-6,
+            zero_tolerance=1e-6,
+            ignore_integrality=True,
+        )
+        self.assertEqual(m.x1.value, 1.1)
+
+        set_var_valid_value(
+            m.x1,
+            var_val=2.00000001,
+            integer_tolerance=1e-6,
+            zero_tolerance=1e-6,
+            ignore_integrality=False,
+        )
+        self.assertEqual(m.x1.value, 2)
+
+        set_var_valid_value(
+            m.x1,
+            var_val=0.0000001,
+            integer_tolerance=1e-9,
+            zero_tolerance=1e-6,
+            ignore_integrality=False,
+        )
+        self.assertEqual(m.x1.value, 0)
+
+    def test_add_var_bound(self):
+        m = SimpleMINLP5().clone()
+        m.x.lb = None
+        m.x.ub = None
+        m.y.lb = None
+        m.y.ub = None
+        solver_object = _MindtPyAlgorithm()
+        solver_object.config = _get_MindtPy_OA_config()
+        solver_object.set_up_solve_data(m)
+        solver_object.create_utility_block(solver_object.working_model, 'MindtPy_utils')
+        add_var_bound(solver_object.working_model, solver_object.config)
+        self.assertEqual(
+            solver_object.working_model.x.lower,
+            -solver_object.config.continuous_var_bound - 1,
+        )
+        self.assertEqual(
+            solver_object.working_model.x.upper,
+            solver_object.config.continuous_var_bound,
+        )
+        self.assertEqual(
+            solver_object.working_model.y.lower,
+            -solver_object.config.integer_var_bound - 1,
+        )
+        self.assertEqual(
+            solver_object.working_model.y.upper, solver_object.config.integer_var_bound
+        )
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/pyomo/contrib/mindtpy/util.py b/pyomo/contrib/mindtpy/util.py
index cd2b31e5954..69c7ca5030a 100644
--- a/pyomo/contrib/mindtpy/util.py
+++ b/pyomo/contrib/mindtpy/util.py
@@ -23,6 +23,7 @@
     RangeSet,
     ConstraintList,
     TransformationFactory,
+    value,
 )
 from pyomo.repn import generate_standard_repn
 from pyomo.contrib.mcpp.pyomo_mcpp import mcpp_available, McCormick
@@ -40,27 +41,24 @@
 numpy = attempt_import('numpy')[0]
 
 
-def calc_jacobians(model, config):
+def calc_jacobians(constraint_list, differentiate_mode):
     """Generates a map of jacobians for the variables in the model.
 
     This function generates a map of jacobians corresponding to the variables in the
-    model.
+    constraint list.
 
     Parameters
     ----------
-    model : Pyomo model
-        Target model to calculate jacobian.
-    config : ConfigBlock
-        The specific configurations for MindtPy.
+    constraint_list : List
+        The list of constraints to calculate Jacobians.
+    differentiate_mode : String
+        The differentiate mode to calculate Jacobians.
     """
     # Map nonlinear_constraint --> Map(
     #     variable --> jacobian of constraint w.r.t. variable)
     jacobians = ComponentMap()
-    if config.differentiate_mode == 'reverse_symbolic':
-        mode = EXPR.differentiate.Modes.reverse_symbolic
-    elif config.differentiate_mode == 'sympy':
-        mode = EXPR.differentiate.Modes.sympy
-    for c in model.MindtPy_utils.nonlinear_constraint_list:
+    mode = EXPR.differentiate.Modes(differentiate_mode)
+    for c in constraint_list:
         vars_in_constr = list(EXPR.identify_variables(c.body))
         jac_list = EXPR.differentiate(c.body, wrt_list=vars_in_constr, mode=mode)
         jacobians[c] = ComponentMap(
@@ -69,7 +67,7 @@ def calc_jacobians(model, config):
     return jacobians
 
 
-def initialize_feas_subproblem(m, config):
+def initialize_feas_subproblem(m, feasibility_norm):
     """Adds feasibility slack variables according to config.feasibility_norm (given an infeasible problem).
        Defines the objective function of the feasibility subproblem.
 
@@ -77,14 +75,14 @@ def initialize_feas_subproblem(m, config):
     ----------
     m : Pyomo model
         The feasbility NLP subproblem.
-    config : ConfigBlock
-        The specific configurations for MindtPy.
+    feasibility_norm : String
+        The norm used to generate the objective function.
     """
     MindtPy = m.MindtPy_utils
     # generate new constraints
     for i, constr in enumerate(MindtPy.nonlinear_constraint_list, 1):
         if constr.has_ub():
-            if config.feasibility_norm in {'L1', 'L2'}:
+            if feasibility_norm in {'L1', 'L2'}:
                 MindtPy.feas_opt.feas_constraints.add(
                     constr.body - constr.upper <= MindtPy.feas_opt.slack_var[i]
                 )
@@ -93,7 +91,7 @@ def initialize_feas_subproblem(m, config):
                     constr.body - constr.upper <= MindtPy.feas_opt.slack_var
                 )
         if constr.has_lb():
-            if config.feasibility_norm in {'L1', 'L2'}:
+            if feasibility_norm in {'L1', 'L2'}:
                 MindtPy.feas_opt.feas_constraints.add(
                     constr.body - constr.lower >= -MindtPy.feas_opt.slack_var[i]
                 )
@@ -102,11 +100,11 @@ def initialize_feas_subproblem(m, config):
                     constr.body - constr.lower >= -MindtPy.feas_opt.slack_var
                 )
     # Setup objective function for the feasibility subproblem.
-    if config.feasibility_norm == 'L1':
+    if feasibility_norm == 'L1':
         MindtPy.feas_obj = Objective(
             expr=sum(s for s in MindtPy.feas_opt.slack_var.values()), sense=minimize
         )
-    elif config.feasibility_norm == 'L2':
+    elif feasibility_norm == 'L2':
         MindtPy.feas_obj = Objective(
             expr=sum(s * s for s in MindtPy.feas_opt.slack_var.values()), sense=minimize
         )
@@ -133,12 +131,12 @@ def add_var_bound(model, config):
         for var in EXPR.identify_variables(c.body):
             if var.has_lb() and var.has_ub():
                 continue
-            elif not var.has_lb():
+            if not var.has_lb():
                 if var.is_integer():
                     var.setlb(-config.integer_var_bound - 1)
                 else:
                     var.setlb(-config.continuous_var_bound - 1)
-            elif not var.has_ub():
+            if not var.has_ub():
                 if var.is_integer():
                     var.setub(config.integer_var_bound)
                 else:
@@ -568,7 +566,9 @@ def set_solver_mipgap(opt, solver_name, config):
         opt.options['add_options'].append('option optcr=%s;' % config.mip_solver_mipgap)
 
 
-def set_solver_constraint_violation_tolerance(opt, solver_name, config):
+def set_solver_constraint_violation_tolerance(
+    opt, solver_name, config, warm_start=True
+):
     """Set constraint violation tolerance for solvers.
 
     Parameters
@@ -602,15 +602,16 @@ def set_solver_constraint_violation_tolerance(opt, solver_name, config):
                 opt.options['add_options'].append(
                     'constr_viol_tol ' + str(config.zero_tolerance)
                 )
-                # Ipopt warmstart options
-                opt.options['add_options'].append(
-                    'warm_start_init_point       yes\n'
-                    'warm_start_bound_push       1e-9\n'
-                    'warm_start_bound_frac       1e-9\n'
-                    'warm_start_slack_bound_frac 1e-9\n'
-                    'warm_start_slack_bound_push 1e-9\n'
-                    'warm_start_mult_bound_push  1e-9\n'
-                )
+                if warm_start:
+                    # Ipopt warmstart options
+                    opt.options['add_options'].append(
+                        'warm_start_init_point       yes\n'
+                        'warm_start_bound_push       1e-9\n'
+                        'warm_start_bound_frac       1e-9\n'
+                        'warm_start_slack_bound_frac 1e-9\n'
+                        'warm_start_slack_bound_push 1e-9\n'
+                        'warm_start_mult_bound_push  1e-9\n'
+                    )
             elif config.nlp_solver_args['solver'] == 'conopt':
                 opt.options['add_options'].append(
                     'RTNWMA ' + str(config.zero_tolerance)
@@ -684,41 +685,20 @@ def copy_var_list_values_from_solution_pool(
         Whether to ignore the integrality of integer variables, by default False.
     """
     for v_from, v_to in zip(from_list, to_list):
-        try:
-            if config.mip_solver == 'cplex_persistent':
-                var_val = solver_model.solution.pool.get_values(
-                    solution_name, var_map[v_from]
-                )
-            elif config.mip_solver == 'gurobi_persistent':
-                solver_model.setParam(gurobipy.GRB.Param.SolutionNumber, solution_name)
-                var_val = var_map[v_from].Xn
-            # We don't want to trigger the reset of the global stale
-            # indicator, so we will set this variable to be "stale",
-            # knowing that set_value will switch it back to "not
-            # stale"
-            v_to.stale = True
-            # NOTE: PEP 2180 changes the var behavior so that domain /
-            # bounds violations no longer generate exceptions (and
-            # instead log warnings).  This means that the following will
-            # always succeed and the ValueError should never be raised.
-            v_to.set_value(var_val, skip_validation=True)
-        except ValueError as e:
-            config.logger.error(e)
-            rounded_val = int(round(var_val))
-            # Check to see if this is just a tolerance issue
-            if ignore_integrality and v_to.is_integer():
-                v_to.set_value(var_val, skip_validation=True)
-            elif v_to.is_integer() and (
-                abs(var_val - rounded_val) <= config.integer_tolerance
-            ):
-                v_to.set_value(rounded_val, skip_validation=True)
-            elif abs(var_val) <= config.zero_tolerance and 0 in v_to.domain:
-                v_to.set_value(0, skip_validation=True)
-            else:
-                config.logger.error(
-                    'Unknown validation domain error setting variable %s' % (v_to.name,)
-                )
-                raise
+        if config.mip_solver == 'cplex_persistent':
+            var_val = solver_model.solution.pool.get_values(
+                solution_name, var_map[v_from]
+            )
+        elif config.mip_solver == 'gurobi_persistent':
+            solver_model.setParam(gurobipy.GRB.Param.SolutionNumber, solution_name)
+            var_val = var_map[v_from].Xn
+        set_var_valid_value(
+            v_to,
+            var_val,
+            config.integer_tolerance,
+            config.zero_tolerance,
+            ignore_integrality,
+        )
 
 
 class GurobiPersistent4MindtPy(GurobiPersistent):
@@ -743,25 +723,6 @@ def f(gurobi_model, where):
         return f
 
 
-def set_up_logger(config):
-    """Set up the formatter and handler for logger.
-
-    Parameters
-    ----------
-    config : ConfigBlock
-        The specific configurations for MindtPy.
-    """
-    config.logger.handlers.clear()
-    config.logger.propagate = False
-    ch = logging.StreamHandler()
-    ch.setLevel(config.logging_level)
-    # create formatter and add it to the handlers
-    formatter = logging.Formatter('%(message)s')
-    ch.setFormatter(formatter)
-    # add the handlers to logger
-    config.logger.addHandler(ch)
-
-
 def epigraph_reformulation(exp, slack_var_list, constraint_list, use_mcpp, sense):
     """Epigraph reformulation.
 
@@ -965,3 +926,101 @@ def generate_norm_constraint(fp_nlp_model, mip_model, config):
             mip_model.MindtPy_utils.discrete_variable_list,
         ):
             fp_nlp_model.norm_constraint.add(nlp_var - mip_var.value <= rhs)
+
+
+def copy_var_list_values(
+    from_list,
+    to_list,
+    config,
+    skip_stale=False,
+    skip_fixed=True,
+    ignore_integrality=False,
+):
+    """Copy variable values from one list to another.
+    Rounds to Binary/Integer if necessary
+    Sets to zero for NonNegativeReals if necessary
+
+    from_list : list
+        The variables that provide the values to copy from.
+    to_list : list
+        The variables that need to set value.
+    config : ConfigBlock
+        The specific configurations for MindtPy.
+    skip_stale : bool, optional
+        Whether to skip the stale variables, by default False.
+    skip_fixed : bool, optional
+        Whether to skip the fixed variables, by default True.
+    ignore_integrality : bool, optional
+        Whether to ignore the integrality of integer variables, by default False.
+    """
+    for v_from, v_to in zip(from_list, to_list):
+        if skip_stale and v_from.stale:
+            continue  # Skip stale variable values.
+        if skip_fixed and v_to.is_fixed():
+            continue  # Skip fixed variables.
+        var_val = value(v_from, exception=False)
+        set_var_valid_value(
+            v_to,
+            var_val,
+            config.integer_tolerance,
+            config.zero_tolerance,
+            ignore_integrality,
+        )
+
+
+def set_var_valid_value(
+    var, var_val, integer_tolerance, zero_tolerance, ignore_integrality
+):
+    """This function tries to set a valid value for variable with the given input.
+    Rounds to Binary/Integer if necessary.
+    Sets to zero for NonNegativeReals if necessary.
+
+    Parameters
+    ----------
+    var : Var
+        The variable that needs to set value.
+    var_val : float
+        The desired value to set for var.
+    integer_tolerance: float
+        Tolerance on integral values.
+    zero_tolerance: float
+        Tolerance on variable equal to zero.
+    ignore_integrality : bool, optional
+        Whether to ignore the integrality of integer variables, by default False.
+
+    Raises
+    ------
+    ValueError
+        Cannot successfully set the value to the variable.
+    """
+    # NOTE: PEP 2180 changes the var behavior so that domain
+    # bounds violations no longer generate exceptions (and
+    # instead log warnings).  This means that the set_value method
+    # will always succeed and the ValueError should never be raised.
+
+    # We don't want to trigger the reset of the global stale
+    # indicator, so we will set this variable to be "stale",
+    # knowing that set_value will switch it back to "not stale".
+    var.stale = True
+    rounded_val = int(round(var_val))
+    if (
+        var_val in var.domain
+        and not ((var.has_lb() and var_val < var.lb))
+        and not ((var.has_ub() and var_val > var.ub))
+    ):
+        var.set_value(var_val)
+    elif var.has_lb() and var_val < var.lb:
+        var.set_value(var.lb)
+    elif var.has_ub() and var_val > var.ub:
+        var.set_value(var.ub)
+    elif ignore_integrality and var.is_integer():
+        var.set_value(var_val, skip_validation=True)
+    elif var.is_integer() and (math.fabs(var_val - rounded_val) <= integer_tolerance):
+        var.set_value(rounded_val)
+    elif abs(var_val) <= zero_tolerance and 0 in var.domain:
+        var.set_value(0)
+    else:
+        raise ValueError(
+            "set_var_valid_value failed with variable {}, value = {} and rounded value = {}"
+            "".format(var.name, var_val, rounded_val)
+        )