Add Primal/Dual Integral in MindtPy

pyomo/contrib/mindtpy/MindtPy.py

@@ -12,21 +12,43 @@
 
 """Implementation of the MindtPy solver.
 
-The MindtPy (MINLP Decomposition Toolkit) solver applies a variety of
-decomposition-based approaches to solve nonlinear continuous-discrete problems.
-These approaches include:
+22.2.10 changes:
+- Add support for partitioning nonlinear-sum objective.
 
-- Outer approximation
-- Benders decomposition [pending]
-- Partial surrogate cuts [pending]
-- Extended cutting plane [pending]
+22.1.12 changes:
+- Improve the log.
 
-This solver implementation was developed by Carnegie Mellon University in the
-research group of Ignacio Grossmann.
+21.12.15 changes:
+- Improve the online doc.
 
-For nonconvex problems, the bounds solve_data.LB and solve_data.UB may not be
-rigorous. Questions: Please make a post at StackOverflow and/or David Bernal
-<https://github.com/bernalde>
+21.11.10 changes:
+- Add support for solution pool of MIP solvers.
+
+21.8.21 changes:
+- Add support for gurobi_persistent solver in (Regularized) LP/NLP-based B&B algorithm.
+
+21.5.19 changes:
+- Add Feasibility Pump strategy.
+- Add Regularized Outer Approximation method.
+- Restructure and simplify the MindyPy code.
+
+20.10.15 changes:
+- Add Extended Cutting Plane and Global Outer Approximation strategy.
+- Update online doc.
+
+20.6.30 changes:
+- Add support for different norms (L1, L2, L-infinity) of the objective function in the feasibility subproblem.
+- Add support for different differentiate_mode to calculate Jacobian.
+
+20.6.9 changes:
+- Add cycling check in Outer Approximation method.
+- Add support for GAMS solvers interface.
+- Fix warmstart for both OA and LP/NLP method.
+
+20.5.9 changes:
+- Add single-tree implementation.
+- Add support for cplex_persistent solver.
+- Fix bug in OA cut expression in cut_generation.py.
 
 """
 from __future__ import division
@@ -35,9 +57,9 @@
                                        time_code, setup_results_object, process_objective, lower_logger_level_to)
 from pyomo.contrib.mindtpy.initialization import MindtPy_initialize_main
 from pyomo.contrib.mindtpy.iterate import MindtPy_iteration_loop
-from pyomo.contrib.mindtpy.util import model_is_valid, set_up_solve_data, set_up_logger
+from pyomo.contrib.mindtpy.util import model_is_valid, set_up_solve_data, set_up_logger, get_primal_integral, get_dual_integral
 from pyomo.core import (Block, ConstraintList, NonNegativeReals,
-                        Set, Suffix, Var, VarList, TransformationFactory, Objective, RangeSet)
+                        Var, VarList, TransformationFactory, RangeSet, minimize)
 from pyomo.opt import SolverFactory
 from pyomo.contrib.mindtpy.config_options import _get_MindtPy_config, check_config
 from pyomo.common.config import add_docstring_list
@@ -52,7 +74,9 @@ class MindtPySolver(object):
     """
     Decomposition solver for Mixed-Integer Nonlinear Programming (MINLP) problems.
 
-    The MindtPy (Mixed-Integer Nonlinear Decomposition Toolbox in Pyomo) solver applies a variety of decomposition-based approaches to solve Mixed-Integer Nonlinear Programming (MINLP) problems. 
+    The MindtPy (Mixed-Integer Nonlinear Decomposition Toolbox in Pyomo) solver 
+    applies a variety of decomposition-based approaches to solve Mixed-Integer 
+    Nonlinear Programming (MINLP) problems. 
     These approaches include:
 
     - Outer approximation (OA)
@@ -63,7 +87,10 @@ class MindtPySolver(object):
     - Regularized LP/NLP based branch-and-bound (RLP/NLP)
     - Feasibility pump (FP)
 
-    This solver implementation has been developed by David Bernal (@bernalde) and Zedong Peng (@ZedongPeng) as part of research efforts at the Grossmann Research Group (http://egon.cheme.cmu.edu/) at the Department of Chemical Engineering at Carnegie Mellon University.
+    This solver implementation has been developed by David Bernal <https://github.com/bernalde>
+    and Zedong Peng <https://github.com/ZedongPeng> as part of research efforts at the Grossmann
+    Research Group (http://egon.cheme.cmu.edu/) at the Department of Chemical Engineering at 
+    Carnegie Mellon University.
     """
     CONFIG = _get_MindtPy_config()
 
@@ -92,7 +119,8 @@ def solve(self, model, **kwds):
         results : SolverResults
             Results from solving the MINLP problem by MindtPy.
         """
-        config = self.CONFIG(kwds.pop('options', {}), preserve_implicit=True)  # TODO: do we need to set preserve_implicit=True?
+        config = self.CONFIG(kwds.pop('options', {
+        }), preserve_implicit=True)  # TODO: do we need to set preserve_implicit=True?
         config.set_value(kwds)
         set_up_logger(config)
         check_config(config)
@@ -115,13 +143,13 @@ def solve(self, model, **kwds):
 
             MindtPy = solve_data.working_model.MindtPy_utils
             setup_results_object(solve_data, config)
-            # In the process_objective function, as long as the objective function is nonlinear, it will be reformulated and the variable/constraint/objective lists will be updated. 
+            # In the process_objective function, as long as the objective function is nonlinear, it will be reformulated and the variable/constraint/objective lists will be updated.
             # For OA/GOA/LP-NLP algorithm, if the objective funtion is linear, it will not be reformulated as epigraph constraint.
             # If the objective function is linear, it will be reformulated as epigraph constraint only if the Feasibility Pump or ROA/RLP-NLP algorithm is activated. (move_linear_objective = True)
-            # In some cases, the variable/constraint/objective lists will not be updated even if the objective is epigraph-reformulated.  
+            # In some cases, the variable/constraint/objective lists will not be updated even if the objective is epigraph-reformulated.
             # In Feasibility Pump, since the distance calculation only includes discrete variables and the epigraph slack variables are continuous variables, the Feasibility Pump algorithm will not affected even if the variable list are updated.
             # In ROA and RLP/NLP, since the distance calculation does not include these epigraph slack variables, they should not be added to the variable list. (update_var_con_list = False)
-            # In the process_objective function, once the objective function has been reformulated as epigraph constraint, the variable/constraint/objective lists will not be updated only if the MINLP has a linear objective function and regularization is activated at the same time. 
+            # In the process_objective function, once the objective function has been reformulated as epigraph constraint, the variable/constraint/objective lists will not be updated only if the MINLP has a linear objective function and regularization is activated at the same time.
             # This is because the epigraph constraint is very "flat" for branching rules. The original objective function will be used for the main problem and epigraph reformulation will be used for the projection problem.
             # TODO: The logic here is too complicated, can we simplify it?
             process_objective(solve_data, config,
@@ -213,20 +241,29 @@ def solve(self, model, **kwds):
                         Var) if not i.fixed],
                     config)
                 # exclude fixed variables here. This is consistent with the definition of variable_list in GDPopt.util
-
-            solve_data.results.problem.lower_bound = solve_data.LB
-            solve_data.results.problem.upper_bound = solve_data.UB
-
-        solve_data.results.solver.timing = solve_data.timing
-        solve_data.results.solver.user_time = solve_data.timing.total
-        solve_data.results.solver.wallclock_time = solve_data.timing.total
-        solve_data.results.solver.iterations = solve_data.mip_iter
-        solve_data.results.solver.num_infeasible_nlp_subproblem = solve_data.nlp_infeasible_counter
-        solve_data.results.solver.best_solution_found_time = solve_data.best_solution_found_time
-
-        if config.single_tree:
-            solve_data.results.solver.num_nodes = solve_data.nlp_iter - \
-                (1 if config.init_strategy == 'rNLP' else 0)
+            if solve_data.objective_sense == minimize:
+                solve_data.results.problem.lower_bound = solve_data.dual_bound
+                solve_data.results.problem.upper_bound = solve_data.primal_bound
+            else:
+                solve_data.results.problem.lower_bound = solve_data.primal_bound
+                solve_data.results.problem.upper_bound = solve_data.dual_bound
+
+            solve_data.results.solver.timing = solve_data.timing
+            solve_data.results.solver.user_time = solve_data.timing.total
+            solve_data.results.solver.wallclock_time = solve_data.timing.total
+            solve_data.results.solver.iterations = solve_data.mip_iter
+            solve_data.results.solver.num_infeasible_nlp_subproblem = solve_data.nlp_infeasible_counter
+            solve_data.results.solver.best_solution_found_time = solve_data.best_solution_found_time
+            solve_data.results.solver.primal_integral = get_primal_integral(solve_data, config)
+            solve_data.results.solver.dual_integral = get_dual_integral(solve_data, config)
+            solve_data.results.solver.primal_dual_gap_integral = solve_data.results.solver.primal_integral + \
+                solve_data.results.solver.dual_integral
+            config.logger.info(' {:<25}:   {:>7.4f} '.format(
+                'Primal-dual gap integral', solve_data.results.solver.primal_dual_gap_integral))
+
+            if config.single_tree:
+                solve_data.results.solver.num_nodes = solve_data.nlp_iter - \
+                    (1 if config.init_strategy == 'rNLP' else 0)
 
         return solve_data.results
 

pyomo/contrib/mindtpy/config_options.py

@@ -301,7 +301,7 @@ def _add_tolerance_configs(CONFIG):
     CONFIG : ConfigBlock
         The specific configurations for MindtPy.
     """
-    CONFIG.declare('bound_tolerance', ConfigValue(
+    CONFIG.declare('absolute_bound_tolerance', ConfigValue(
         default=1E-4,
         domain=PositiveFloat,
         description='Bound tolerance',
@@ -312,7 +312,7 @@ def _add_tolerance_configs(CONFIG):
         domain=PositiveFloat,
         description='Relative bound tolerance',
         doc='Relative tolerance for bound feasibility checks.'
-            '(UB - LB) / (1e-10+|bestinteger|) <= relative tolerance.'
+            '|Primal Bound - Dual Bound| / (1e-10 + |Primal Bound|) <= relative tolerance.'
     ))
     CONFIG.declare('small_dual_tolerance', ConfigValue(
         default=1E-8,
@@ -369,6 +369,11 @@ def _add_bound_configs(CONFIG):
         description='Default bound added to unbounded integral variables in nonlinear constraint if single tree is activated.',
         domain=PositiveFloat
     ))
+    CONFIG.declare('initial_bound_coef', ConfigValue(
+        default=1E-1,
+        domain=PositiveFloat,
+        description='The coefficient used to approximate the initial primal/dual bound.'
+    ))
 
 
 def _add_fp_configs(CONFIG):
@@ -536,7 +541,7 @@ def check_config(config):
             config.equality_relaxation = False
     # if ecp tolerance is not provided use bound tolerance
     if config.ecp_tolerance is None:
-        config.ecp_tolerance = config.bound_tolerance
+        config.ecp_tolerance = config.absolute_bound_tolerance
 
     if config.solver_tee:
         config.mip_solver_tee = True

pyomo/contrib/mindtpy/feasibility_pump.py

@@ -71,10 +71,10 @@ def solve_fp_subproblem(solve_data, config):
     fp_nlp.MindtPy_utils.objective_list[-1].deactivate()
     if solve_data.objective_sense == minimize:
         fp_nlp.improving_objective_cut = Constraint(
-            expr=sum(fp_nlp.MindtPy_utils.objective_value[:]) <= solve_data.UB)
+            expr=sum(fp_nlp.MindtPy_utils.objective_value[:]) <= solve_data.primal_bound)
     else:
         fp_nlp.improving_objective_cut = Constraint(
-            expr=sum(fp_nlp.MindtPy_utils.objective_value[:]) >= solve_data.LB)
+            expr=sum(fp_nlp.MindtPy_utils.objective_value[:]) >= solve_data.primal_bound)
 
     # Add norm_constraint, which guarantees the monotonicity of the norm objective value sequence of all iterations
     # Ref: Paper 'A storm of feasibility pumps for nonconvex MINLP'   https://doi.org/10.1007/s10107-012-0608-x
@@ -181,7 +181,7 @@ def fp_loop(solve_data, config):
             config.logger.info(solve_data.log_formatter.format(
                 solve_data.fp_iter, 'FP-NLP', value(
                     fp_nlp.MindtPy_utils.fp_nlp_obj),
-                solve_data.LB, solve_data.UB, solve_data.rel_gap,
+                solve_data.primal_bound, solve_data.dual_bound, solve_data.rel_gap,
                 get_main_elapsed_time(solve_data.timing)))
             handle_fp_subproblem_optimal(fp_nlp, solve_data, config)
         elif fp_nlp_result.solver.termination_condition in {tc.infeasible, tc.noSolution}:
@@ -296,7 +296,7 @@ def handle_fp_main_tc(feas_main_results, solve_data, config):
         config.logger.info(solve_data.log_formatter.format(
             solve_data.fp_iter, 'FP-MIP', value(
                 solve_data.mip.MindtPy_utils.fp_mip_obj),
-            solve_data.LB, solve_data.UB, solve_data.rel_gap, get_main_elapsed_time(solve_data.timing)))
+            solve_data.primal_bound, solve_data.dual_bound, solve_data.rel_gap, get_main_elapsed_time(solve_data.timing)))
         return False
     elif feas_main_results.solver.termination_condition is tc.maxTimeLimit:
         config.logger.warning('FP-MIP reaches max TimeLimit')

pyomo/contrib/mindtpy/initialization.py

@@ -82,10 +82,10 @@ def MindtPy_initialize_main(solve_data, config):
         '{} is the initial strategy being used.'
         '\n'.format(config.init_strategy))
     config.logger.info(
-        ' =============================================================================================')
+        ' ===============================================================================================')
     config.logger.info(
-        ' {:>9} | {:>15} | {:>15} | {:>11} | {:>11} | {:^7} | {:>7}\n'.format('Iteration', 'Subproblem Type', 'Objective Value', 'Lower Bound',
-                                                                              'Upper Bound', ' Gap ', 'Time(s)'))
+        ' {:>9} | {:>15} | {:>15} | {:>12} | {:>12} | {:^7} | {:>7}\n'.format('Iteration', 'Subproblem Type', 'Objective Value', 'Primal Bound',
+                                                                              'Dual Bound', ' Gap ', 'Time(s)'))
     # Do the initialization
     if config.init_strategy == 'rNLP':
         init_rNLP(solve_data, config)
@@ -142,7 +142,7 @@ def init_rNLP(solve_data, config):
         dual_values = list(
             m.dual[c] for c in MindtPy.constraint_list) if config.calculate_dual else None
         config.logger.info(solve_data.log_formatter.format('-', 'Relaxed NLP', value(main_objective.expr),
-                                                           solve_data.LB, solve_data.UB, solve_data.rel_gap,
+                                                           solve_data.primal_bound, solve_data.dual_bound, solve_data.rel_gap,
                                                            get_main_elapsed_time(solve_data.timing)))
         # Add OA cut
         if config.strategy in {'OA', 'GOA', 'FP'}:
@@ -236,7 +236,7 @@ def init_max_binaries(solve_data, config):
             solve_data.working_model.MindtPy_utils.variable_list,
             config)
         config.logger.info(solve_data.log_formatter.format('-', 'Max binary MILP', value(MindtPy.max_binary_obj.expr),
-                                                           solve_data.LB, solve_data.UB, solve_data.rel_gap,
+                                                           solve_data.primal_bound, solve_data.dual_bound, solve_data.rel_gap,
                                                            get_main_elapsed_time(solve_data.timing)))
     elif solve_terminate_cond is tc.infeasible:
         raise ValueError(