Merge pull request Pyomo#1476 from ZedongPeng/new-lp-nlp

Fix warmstart and add cycling check in MindtPy

doc/OnlineDocs/contributed_packages/mindtpy.rst

@@ -33,7 +33,7 @@ An example which includes the modeling approach may be found below.
   >>> model.x = Var(bounds=(1.0,10.0),initialize=5.0)
   >>> model.y = Var(within=Binary)
 
-  >>> model.c1 = Constraint(expr=(model.x-3.0)**2 <= 50.0*(1-model.y))
+  >>> model.c1 = Constraint(expr=(model.x-4.0)**2 - model.x <= 50.0*(1-model.y))
   >>> model.c2 = Constraint(expr=model.x*log(model.x)+5.0 <= 50.0*(model.y))
 
   >>> model.objective = Objective(expr=model.x, sense=minimize)
@@ -87,7 +87,7 @@ A usage example for single tree is as follows:
   >>> model.x = pyo.Var(bounds=(1.0, 10.0), initialize=5.0)
   >>> model.y = pyo.Var(within=Binary)
 
-  >>> model.c1 = pyo.Constraint(expr=(model.x-3.0)**2 <= 50.0*(1-model.y))
+  >>> model.c1 = Constraint(expr=(model.x-4.0)**2 - model.x <= 50.0*(1-model.y))
   >>> model.c2 = pyo.Constraint(expr=model.x*log(model.x)+5.0 <= 50.0*(model.y))
   
   >>> model.objective = pyo.Objective(expr=model.x, sense=pyo.minimize)

pyomo/contrib/gdpopt/util.py

@@ -105,7 +105,7 @@ def presolve_lp_nlp(solve_data, config):
     return False, None
 
 
-def process_objective(solve_data, config, move_linear_objective=False):
+def process_objective(solve_data, config, move_linear_objective=False, use_mcpp=True):
     """Process model objective function.
 
     Check that the model has only 1 valid objective.
@@ -144,10 +144,11 @@ def process_objective(solve_data, config, move_linear_objective=False):
         if move_linear_objective:
             config.logger.info("Moving objective to constraint set.")
         else:
-            config.logger.info("Objective is nonlinear. Moving it to constraint set.")
+            config.logger.info(
+                "Objective is nonlinear. Moving it to constraint set.")
 
         util_blk.objective_value = Var(domain=Reals, initialize=0)
-        if mcpp_available():
+        if mcpp_available() and use_mcpp:
             mc_obj = McCormick(main_obj.expr)
             util_blk.objective_value.setub(mc_obj.upper())
             util_blk.objective_value.setlb(mc_obj.lower())
@@ -206,8 +207,8 @@ def copy_var_list_values(from_list, to_list, config,
             # Check to see if this is just a tolerance issue
             if ignore_integrality \
                 and ('is not in domain Binary' in err_msg
-                or 'is not in domain Integers' in err_msg):
-               v_to.value = value(v_from, exception=False)
+                     or 'is not in domain Integers' in err_msg):
+                v_to.value = value(v_from, exception=False)
             elif 'is not in domain Binary' in err_msg and (
                     fabs(var_val - 1) <= config.integer_tolerance or
                     fabs(var_val) <= config.integer_tolerance):

pyomo/contrib/mindtpy/MindtPy.py

@@ -118,7 +118,7 @@ class MindtPySolver(object):
     ))
     CONFIG.declare("nlp_solver", ConfigValue(
         default="ipopt",
-        domain=In(["ipopt"]),
+        domain=In(["ipopt", "gams"]),
         description="NLP subsolver name",
         doc="Which NLP subsolver is going to be used for solving the nonlinear"
             "subproblems"
@@ -191,7 +191,7 @@ class MindtPySolver(object):
         description="Tolerance on variable bounds."
     ))
     CONFIG.declare("zero_tolerance", ConfigValue(
-        default=1E-10,
+        default=1E-8,
         description="Tolerance on variable equal to zero."
     ))
     CONFIG.declare("initial_feas", ConfigValue(
@@ -210,7 +210,7 @@ class MindtPySolver(object):
         domain=bool
     ))
     CONFIG.declare("add_integer_cuts", ConfigValue(
-        default=True,
+        default=False,
         description="Add integer cuts (no-good cuts) to binary variables to disallow same integer solution again."
                     "Note that 'integer_to_binary' flag needs to be used to apply it to actual integers and not just binaries.",
         domain=bool
@@ -231,6 +231,21 @@ class MindtPySolver(object):
                     "slack variables here are used to deal with nonconvex MINLP",
         domain=bool
     ))
+    CONFIG.declare("continuous_var_bound", ConfigValue(
+        default=1e10,
+        description="default bound added to unbounded continuous variables in nonlinear constraint if single tree is activated.",
+        domain=PositiveFloat
+    ))
+    CONFIG.declare("integer_var_bound", ConfigValue(
+        default=1e9,
+        description="default bound added to unbounded integral variables in nonlinear constraint if single tree is activated.",
+        domain=PositiveFloat
+    ))
+    CONFIG.declare("cycling_check", ConfigValue(
+        default=True,
+        description="check if OA algorithm is stalled in a cycle and terminate.",
+        domain=bool
+    ))
 
     def available(self, exception_flag=True):
         """Check if solver is available.
@@ -273,6 +288,8 @@ def solve(self, model, **kwds):
         solve_data = MindtPySolveData()
         solve_data.results = SolverResults()
         solve_data.timing = Container()
+        solve_data.curr_int_sol = []
+        solve_data.prev_int_sol = []
 
         solve_data.original_model = model
         solve_data.working_model = model.clone()
@@ -288,7 +305,7 @@ def solve(self, model, **kwds):
 
             MindtPy = solve_data.working_model.MindtPy_utils
             setup_results_object(solve_data, config)
-            process_objective(solve_data, config)
+            process_objective(solve_data, config, use_mcpp=False)
 
             # Save model initial values.
             solve_data.initial_var_values = list(
@@ -416,6 +433,10 @@ def solve(self, model, **kwds):
 
         solve_data.results.solver.iterations = solve_data.mip_iter
 
+        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/initialization.py

@@ -15,13 +15,18 @@
 from pyomo.contrib.mindtpy.nlp_solve import (solve_NLP_subproblem,
                                              handle_NLP_subproblem_optimal, handle_NLP_subproblem_infeasible,
                                              handle_NLP_subproblem_other_termination)
+from pyomo.contrib.mindtpy.util import var_bound_add
 
 
 def MindtPy_initialize_master(solve_data, config):
     """Initialize the decomposition algorithm.
     This includes generating the initial cuts require to build the master
     problem.
     """
+    # if single tree is activated, we need to add bounds for unbounded variables in nonlinear constraints to avoid unbounded master problem.
+    if config.single_tree:
+        var_bound_add(solve_data, config)
+
     m = solve_data.mip = solve_data.working_model.clone()
     MindtPy = m.MindtPy_utils
     m.dual.deactivate()
@@ -58,7 +63,7 @@ def MindtPy_initialize_master(solve_data, config):
         # else:
 
         fixed_nlp, fixed_nlp_result = solve_NLP_subproblem(solve_data, config)
-        if fixed_nlp_result.solver.termination_condition is tc.optimal:
+        if fixed_nlp_result.solver.termination_condition is tc.optimal or fixed_nlp_result.solver.termination_condition is tc.locallyOptimal:
             handle_NLP_subproblem_optimal(fixed_nlp, solve_data, config)
         elif fixed_nlp_result.solver.termination_condition is tc.infeasible:
             handle_NLP_subproblem_infeasible(fixed_nlp, solve_data, config)
@@ -79,7 +84,7 @@ def init_rNLP(solve_data, config):
         results = SolverFactory(config.nlp_solver).solve(
             m, **config.nlp_solver_args)
     subprob_terminate_cond = results.solver.termination_condition
-    if subprob_terminate_cond is tc.optimal:
+    if subprob_terminate_cond is tc.optimal or subprob_terminate_cond is tc.locallyOptimal:
         main_objective = next(m.component_data_objects(Objective, active=True))
         nlp_solution_values = list(v.value for v in MindtPy.variable_list)
         dual_values = list(m.dual[c] for c in MindtPy.constraint_list)
@@ -144,7 +149,11 @@ def init_max_binaries(solve_data, config):
     opt = SolverFactory(config.mip_solver)
     if isinstance(opt, PersistentSolver):
         opt.set_instance(m)
-    results = opt.solve(m, options=config.mip_solver_args)
+    mip_args = dict(config.mip_solver_args)
+    if config.mip_solver == 'gams':
+        mip_args['add_options'] = mip_args.get('add_options', [])
+        mip_args['add_options'].append('option optcr=0.0;')
+    results = opt.solve(m, **mip_args)
 
     solve_terminate_cond = results.solver.termination_condition
     if solve_terminate_cond is tc.optimal:

pyomo/contrib/mindtpy/iterate.py

@@ -6,7 +6,7 @@
 from pyomo.contrib.mindtpy.nlp_solve import (solve_NLP_subproblem,
                                              handle_NLP_subproblem_optimal, handle_NLP_subproblem_infeasible,
                                              handle_NLP_subproblem_other_termination)
-from pyomo.core import minimize, Objective
+from pyomo.core import minimize, Objective, Var
 from pyomo.opt import TerminationCondition as tc
 from pyomo.contrib.gdpopt.util import get_main_elapsed_time
 
@@ -21,7 +21,7 @@ def MindtPy_iteration_loop(solve_data, config):
             '---MindtPy Master Iteration %s---'
             % solve_data.mip_iter)
 
-        if algorithm_should_terminate(solve_data, config):
+        if algorithm_should_terminate(solve_data, config, check_cycling=False):
             break
 
         solve_data.mip_subiter = 0
@@ -39,15 +39,15 @@ def MindtPy_iteration_loop(solve_data, config):
         else:
             raise NotImplementedError()
 
-        if algorithm_should_terminate(solve_data, config):
+        if algorithm_should_terminate(solve_data, config, check_cycling=True):
             break
 
         if config.single_tree is False:  # if we don't use lazy callback, i.e. LP_NLP
             # Solve NLP subproblem
             # The constraint linearization happens in the handlers
             fixed_nlp, fixed_nlp_result = solve_NLP_subproblem(
                 solve_data, config)
-            if fixed_nlp_result.solver.termination_condition is tc.optimal:
+            if fixed_nlp_result.solver.termination_condition is tc.optimal or fixed_nlp_result.solver.termination_condition is tc.locallyOptimal:
                 handle_NLP_subproblem_optimal(fixed_nlp, solve_data, config)
             elif fixed_nlp_result.solver.termination_condition is tc.infeasible:
                 handle_NLP_subproblem_infeasible(fixed_nlp, solve_data, config)
@@ -93,7 +93,7 @@ def MindtPy_iteration_loop(solve_data, config):
         #         config.strategy = 'OA'
 
 
-def algorithm_should_terminate(solve_data, config):
+def algorithm_should_terminate(solve_data, config, check_cycling):
     """Check if the algorithm should terminate.
 
     Termination conditions based on solver options and progress.
@@ -133,6 +133,30 @@ def algorithm_should_terminate(solve_data, config):
             format(solve_data.LB, solve_data.UB))
         solve_data.results.solver.termination_condition = tc.maxTimeLimit
         return True
+
+    # Cycling check
+    if config.cycling_check == True and solve_data.mip_iter >= 1 and check_cycling:
+        temp = []
+        for var in solve_data.mip.component_data_objects(ctype=Var):
+            if var.is_integer():
+                temp.append(int(round(var.value)))
+        solve_data.curr_int_sol = temp
+
+        if solve_data.curr_int_sol == solve_data.prev_int_sol:
+            config.logger.info(
+                'Cycling happens after {} master iterations. '
+                'This issue happens when the NLP subproblem violates constraint qualification. '
+                'Convergence to optimal solution is not guaranteed.'
+                .format(solve_data.mip_iter))
+            config.logger.info(
+                'Final bound values: LB: {}  UB: {}'.
+                format(solve_data.LB, solve_data.UB))
+            # TODO determine solve_data.LB, solve_data.UB is inf or -inf.
+            solve_data.results.solver.termination_condition = tc.feasible
+            return True
+
+        solve_data.prev_int_sol = solve_data.curr_int_sol
+
     # if not algorithm_is_making_progress(solve_data, config):
     #     config.logger.debug(
     #         'Algorithm is not making enough progress. '

pyomo/contrib/mindtpy/mip_solve.py

@@ -83,16 +83,20 @@ def solve_OA_master(solve_data, config):
         masteropt._solver_model.set_log_stream(None)
         masteropt._solver_model.set_error_stream(None)
         masteropt.options['timelimit'] = config.time_limit
+    mip_args = dict(config.mip_solver_args)
+    if config.mip_solver == 'gams':
+        mip_args['add_options'] = mip_args.get('add_options', [])
+        mip_args['add_options'].append('option optcr=0.0;')
     master_mip_results = masteropt.solve(
-        solve_data.mip, **config.mip_solver_args)  # , tee=True)
+        solve_data.mip, **mip_args)  # , tee=True)
 
     if master_mip_results.solver.termination_condition is tc.optimal:
         if config.single_tree:
             if main_objective.sense == minimize:
                 solve_data.LB = max(
                     master_mip_results.problem.lower_bound, solve_data.LB)
                 solve_data.LB_progress.append(solve_data.LB)
-
+            else:
                 solve_data.UB = min(
                     master_mip_results.problem.upper_bound, solve_data.UB)
                 solve_data.UB_progress.append(solve_data.UB)
@@ -115,10 +119,10 @@ def handle_master_mip_optimal(master_mip, solve_data, config, copy=True):
         master_mip.component_data_objects(Objective, active=True))
     # check if the value of binary variable is valid
     for var in MindtPy.variable_list:
-        if var.value == None:
+        if var.value == None and var.is_integer():
             config.logger.warning(
-                "Variables {} not initialized are set to it's lower bound when using the initial_binary initialization method".format(var.name))
-            var.value = 0  # nlp_var.bounds[0]
+                "Integer variable {} not initialized. It is set to it's lower bound when using the initial_binary initialization method".format(var.name))
+            var.value = var.lb  # nlp_var.bounds[0]
     copy_var_list_values(
         master_mip.MindtPy_utils.variable_list,
         solve_data.working_model.MindtPy_utils.variable_list,
@@ -189,10 +193,8 @@ def handle_master_mip_infeasible(master_mip, solve_data, config):
     main_objective = next(
         master_mip.component_data_objects(Objective, active=True))
     if main_objective.sense == minimize:
-        solve_data.LB = float('inf')
-        solve_data.LB_progress.append(solve_data.UB)
+        solve_data.LB_progress.append(solve_data.LB)
     else:
-        solve_data.UB = float('-inf')
         solve_data.UB_progress.append(solve_data.UB)
 
 

pyomo/contrib/mindtpy/nlp_solve.py

@@ -35,13 +35,6 @@ def solve_NLP_subproblem(solve_data, config):
     # Set up NLP
     TransformationFactory('core.fix_integer_vars').apply_to(fixed_nlp)
 
-    # restore original variable values
-    for nlp_var, orig_val in zip(
-            MindtPy.variable_list,
-            solve_data.initial_var_values):
-        if not nlp_var.fixed and not nlp_var.is_binary():
-            nlp_var.value = orig_val
-
     MindtPy.MindtPy_linear_cuts.deactivate()
     fixed_nlp.tmp_duals = ComponentMap()
     # tmp_duals are the value of the dual variables stored before using deactivate trivial contraints
@@ -53,18 +46,28 @@ def solve_NLP_subproblem(solve_data, config):
     # | g(x) <= b  | -1    | g(x1) > b    | g(x1) - b            |
     # | g(x) >= b  | +1    | g(x1) >= b   | 0                    |
     # | g(x) >= b  | +1    | g(x1) < b    | b - g(x1)            |
-
+    evaluation_error = False
     for c in fixed_nlp.component_data_objects(ctype=Constraint, active=True,
                                               descend_into=True):
         # We prefer to include the upper bound as the right hand side since we are
         # considering c by default a (hopefully) convex function, which would make
         # c >= lb a nonconvex inequality which we wouldn't like to add linearizations
         # if we don't have to
-        rhs = c.upper if c. has_ub() else c.lower
+        rhs = c.upper if c.has_ub() else c.lower
         c_geq = -1 if c.has_ub() else 1
         # c_leq = 1 if c.has_ub else -1
-        fixed_nlp.tmp_duals[c] = c_geq * max(
-            0, c_geq*(rhs - value(c.body)))
+        try:
+            fixed_nlp.tmp_duals[c] = c_geq * max(
+                0, c_geq*(rhs - value(c.body)))
+        except (ValueError, OverflowError) as error:
+            fixed_nlp.tmp_duals[c] = None
+            evaluation_error = True
+    if evaluation_error:
+        for nlp_var, orig_val in zip(
+                MindtPy.variable_list,
+                solve_data.initial_var_values):
+            if not nlp_var.fixed and not nlp_var.is_binary():
+                nlp_var.value = orig_val
         # fixed_nlp.tmp_duals[c] = c_leq * max(
         #     0, c_leq*(value(c.body) - rhs))
         # TODO: change logic to c_leq based on benchmarking
@@ -217,7 +220,7 @@ def solve_NLP_feas(solve_data, config):
         feas_soln = SolverFactory(config.nlp_solver).solve(
             fixed_nlp, **config.nlp_solver_args)
     subprob_terminate_cond = feas_soln.solver.termination_condition
-    if subprob_terminate_cond is tc.optimal:
+    if subprob_terminate_cond is tc.optimal or subprob_terminate_cond is tc.locallyOptimal:
         copy_var_list_values(
             MindtPy.variable_list,
             solve_data.working_model.MindtPy_utils.variable_list,

pyomo/contrib/mindtpy/single_tree.py

@@ -124,19 +124,6 @@ def handle_lazy_master_mip_feasible_sol(self, master_mip, solve_data, config, op
                                        master_mip.MindtPy_utils.variable_list,
                                        solve_data.working_model.MindtPy_utils.variable_list,
                                        config)
-        # update the bound
-        if main_objective.sense == minimize:
-            solve_data.LB = max(
-                self.get_objective_value(),
-                # self.get_best_objective_value(),
-                solve_data.LB)
-            solve_data.LB_progress.append(solve_data.LB)
-        else:
-            solve_data.UB = min(
-                self.get_objective_value(),
-                # self.get_best_objective_value(),
-                solve_data.UB)
-            solve_data.UB_progress.append(solve_data.UB)
         config.logger.info(
             'MIP %s: OBJ: %s  LB: %s  UB: %s'
             % (solve_data.mip_iter, value(MindtPy.MindtPy_oa_obj.expr),
@@ -241,7 +228,7 @@ def __call__(self):
         fixed_nlp, fixed_nlp_result = solve_NLP_subproblem(solve_data, config)
 
         # add oa cuts
-        if fixed_nlp_result.solver.termination_condition is tc.optimal:
+        if fixed_nlp_result.solver.termination_condition is tc.optimal or fixed_nlp_result.solver.termination_condition is tc.locallyOptimal:
             self.handle_lazy_NLP_subproblem_optimal(
                 fixed_nlp, solve_data, config, opt)
         elif fixed_nlp_result.solver.termination_condition is tc.infeasible:
@@ -250,3 +237,4 @@ def __call__(self):
         else:
             self.handle_lazy_NLP_subproblem_other_termination(fixed_nlp, fixed_nlp_result.solver.termination_condition,
                                                               solve_data, config)
+