diff --git a/pyomo/contrib/solver/ipopt.py b/pyomo/contrib/solver/ipopt.py index 3ac1a5ac4a2..dc632adb184 100644 --- a/pyomo/contrib/solver/ipopt.py +++ b/pyomo/contrib/solver/ipopt.py @@ -17,7 +17,11 @@ from pyomo.common import Executable from pyomo.common.config import ConfigValue, document_kwargs_from_configdict, ConfigDict -from pyomo.common.errors import PyomoException, DeveloperError +from pyomo.common.errors import ( + PyomoException, + DeveloperError, + InfeasibleConstraintException, +) from pyomo.common.tempfiles import TempfileManager from pyomo.common.timing import HierarchicalTimer from pyomo.core.base.var import _GeneralVarData @@ -72,11 +76,7 @@ def __init__( ), ) self.writer_config: ConfigDict = self.declare( - 'writer_config', - ConfigValue( - default=NLWriter.CONFIG(), - description="Configuration that controls options in the NL writer.", - ), + 'writer_config', NLWriter.CONFIG() ) @@ -314,15 +314,19 @@ def solve(self, model, **kwds): ) as row_file, open(basename + '.col', 'w') as col_file: timer.start('write_nl_file') self._writer.config.set_value(config.writer_config) - nl_info = self._writer.write( - model, - nl_file, - row_file, - col_file, - symbolic_solver_labels=config.symbolic_solver_labels, - ) + try: + nl_info = self._writer.write( + model, + nl_file, + row_file, + col_file, + symbolic_solver_labels=config.symbolic_solver_labels, + ) + proven_infeasible = False + except InfeasibleConstraintException: + proven_infeasible = True timer.stop('write_nl_file') - if len(nl_info.variables) > 0: + if not proven_infeasible and len(nl_info.variables) > 0: # Get a copy of the environment to pass to the subprocess env = os.environ.copy() if nl_info.external_function_libraries: @@ -361,11 +365,17 @@ def solve(self, model, **kwds): timer.stop('subprocess') # This is the stuff we need to parse to get the iterations # and time - iters, ipopt_time_nofunc, ipopt_time_func, ipopt_total_time = ( + (iters, ipopt_time_nofunc, ipopt_time_func, ipopt_total_time) = ( self._parse_ipopt_output(ostreams[0]) ) - if len(nl_info.variables) == 0: + if proven_infeasible: + results = Results() + results.termination_condition = TerminationCondition.provenInfeasible + results.solution_loader = SolSolutionLoader(None, None) + results.iteration_count = 0 + results.timing_info.total_seconds = 0 + elif len(nl_info.variables) == 0: if len(nl_info.eliminated_vars) == 0: results = Results() results.termination_condition = TerminationCondition.emptyModel @@ -457,7 +467,7 @@ def solve(self, model, **kwds): ) results.solver_configuration = config - if len(nl_info.variables) > 0: + if not proven_infeasible and len(nl_info.variables) > 0: results.solver_log = ostreams[0].getvalue() # Capture/record end-time / wall-time diff --git a/pyomo/contrib/solver/tests/solvers/test_solvers.py b/pyomo/contrib/solver/tests/solvers/test_solvers.py index cf5f6cf5c57..a4f4a3bc389 100644 --- a/pyomo/contrib/solver/tests/solvers/test_solvers.py +++ b/pyomo/contrib/solver/tests/solvers/test_solvers.py @@ -1508,6 +1508,71 @@ def test_bug_2(self, name: str, opt_class: Type[SolverBase], use_presolve: bool) res = opt.solve(m) self.assertAlmostEqual(res.incumbent_objective, -18, 5) + @parameterized.expand(input=_load_tests(nl_solvers)) + def test_presolve_with_zero_coef( + self, name: str, opt_class: Type[SolverBase], use_presolve: bool + ): + opt: SolverBase = opt_class() + if not opt.available(): + raise unittest.SkipTest(f'Solver {opt.name} not available.') + if use_presolve: + opt.config.writer_config.linear_presolve = True + else: + opt.config.writer_config.linear_presolve = False + + """ + when c2 gets presolved out, c1 becomes + x - y + y = 0 which becomes + x - 0*y == 0 which is the zero we are testing for + """ + m = pe.ConcreteModel() + m.x = pe.Var() + m.y = pe.Var() + m.z = pe.Var() + m.obj = pe.Objective(expr=m.x**2 + m.y**2 + m.z**2) + m.c1 = pe.Constraint(expr=m.x == m.y + m.z + 1.5) + m.c2 = pe.Constraint(expr=m.z == -m.y) + + res = opt.solve(m) + self.assertAlmostEqual(res.incumbent_objective, 2.25) + self.assertAlmostEqual(m.x.value, 1.5) + self.assertAlmostEqual(m.y.value, 0) + self.assertAlmostEqual(m.z.value, 0) + + m.x.setlb(2) + res = opt.solve( + m, load_solutions=False, raise_exception_on_nonoptimal_result=False + ) + if use_presolve: + exp = TerminationCondition.provenInfeasible + else: + exp = TerminationCondition.locallyInfeasible + self.assertEqual(res.termination_condition, exp) + + m = pe.ConcreteModel() + m.w = pe.Var() + m.x = pe.Var() + m.y = pe.Var() + m.z = pe.Var() + m.obj = pe.Objective(expr=m.x**2 + m.y**2 + m.z**2 + m.w**2) + m.c1 = pe.Constraint(expr=m.x + m.w == m.y + m.z) + m.c2 = pe.Constraint(expr=m.z == -m.y) + m.c3 = pe.Constraint(expr=m.x == -m.w) + + res = opt.solve(m) + self.assertAlmostEqual(res.incumbent_objective, 0) + self.assertAlmostEqual(m.w.value, 0) + self.assertAlmostEqual(m.x.value, 0) + self.assertAlmostEqual(m.y.value, 0) + self.assertAlmostEqual(m.z.value, 0) + + del m.c1 + m.c1 = pe.Constraint(expr=m.x + m.w == m.y + m.z + 1.5) + res = opt.solve( + m, load_solutions=False, raise_exception_on_nonoptimal_result=False + ) + self.assertEqual(res.termination_condition, exp) + @parameterized.expand(input=_load_tests(all_solvers)) def test_scaling(self, name: str, opt_class: Type[SolverBase], use_presolve: bool): opt: SolverBase = opt_class() diff --git a/pyomo/repn/plugins/nl_writer.py b/pyomo/repn/plugins/nl_writer.py index f3ff94ea8c9..a256cd1b900 100644 --- a/pyomo/repn/plugins/nl_writer.py +++ b/pyomo/repn/plugins/nl_writer.py @@ -1760,7 +1760,7 @@ def _linear_presolve(self, comp_by_linear_var, lcon_by_linear_nnz, var_bounds): id2_isdiscrete = var_map[id2].domain.isdiscrete() if var_map[_id].domain.isdiscrete() ^ id2_isdiscrete: # if only one variable is discrete, then we need to - # substiitute out the other + # substitute out the other if id2_isdiscrete: _id, id2 = id2, _id coef, coef2 = coef2, coef @@ -1820,10 +1820,15 @@ def _linear_presolve(self, comp_by_linear_var, lcon_by_linear_nnz, var_bounds): # appropriately (that expr_info is persisting in the # eliminated_vars dict - and we will use that to # update other linear expressions later.) + old_nnz = len(expr_info.linear) c = expr_info.linear.pop(_id, 0) + nnz = old_nnz - 1 expr_info.const += c * b if x in expr_info.linear: expr_info.linear[x] += c * a + if expr_info.linear[x] == 0: + nnz -= 1 + coef = expr_info.linear.pop(x) elif a: expr_info.linear[x] = c * a # replacing _id with x... NNZ is not changing, @@ -1831,10 +1836,17 @@ def _linear_presolve(self, comp_by_linear_var, lcon_by_linear_nnz, var_bounds): # this constraint comp_by_linear_var[x].append((con_id, expr_info)) continue - # NNZ has been reduced by 1 - nnz = len(expr_info.linear) - _old = lcon_by_linear_nnz[nnz + 1] + _old = lcon_by_linear_nnz[old_nnz] if con_id in _old: + if not nnz: + if abs(expr_info.const) > TOL: + # constraint is trivially infeasible + raise InfeasibleConstraintException( + "model contains a trivially infeasible constraint " + f"{expr_info.const} == {coef}*{var_map[x]}" + ) + # constraint is trivially feasible + eliminated_cons.add(con_id) lcon_by_linear_nnz[nnz][con_id] = _old.pop(con_id) # If variables were replaced by the variable that # we are currently eliminating, then we need to update diff --git a/pyomo/repn/tests/ampl/test_nlv2.py b/pyomo/repn/tests/ampl/test_nlv2.py index 8b95fc03bdb..215715dba10 100644 --- a/pyomo/repn/tests/ampl/test_nlv2.py +++ b/pyomo/repn/tests/ampl/test_nlv2.py @@ -1683,6 +1683,131 @@ def test_presolve_named_expressions(self): G0 2 #obj 0 0 1 0 +""", + OUT.getvalue(), + ) + ) + + def test_presolve_zero_coef(self): + m = ConcreteModel() + m.x = Var() + m.y = Var() + m.z = Var() + m.obj = Objective(expr=m.x**2 + m.y**2 + m.z**2) + m.c1 = Constraint(expr=m.x == m.y + m.z + 1.5) + m.c2 = Constraint(expr=m.z == -m.y) + + OUT = io.StringIO() + with LoggingIntercept() as LOG: + nlinfo = nl_writer.NLWriter().write( + m, OUT, symbolic_solver_labels=True, linear_presolve=True + ) + self.assertEqual(LOG.getvalue(), "") + + self.assertEqual(nlinfo.eliminated_vars[0], (m.x, 1.5)) + self.assertIs(nlinfo.eliminated_vars[1][0], m.y) + self.assertExpressionsEqual( + nlinfo.eliminated_vars[1][1], LinearExpression([-1.0 * m.z]) + ) + + self.assertEqual( + *nl_diff( + """g3 1 1 0 # problem unknown + 1 0 1 0 0 #vars, constraints, objectives, ranges, eqns + 0 1 0 0 0 0 #nonlinear constrs, objs; ccons: lin, nonlin, nd, nzlb + 0 0 #network constraints: nonlinear, linear + 0 1 0 #nonlinear vars in constraints, objectives, both + 0 0 0 1 #linear network variables; functions; arith, flags + 0 0 0 0 0 #discrete variables: binary, integer, nonlinear (b,c,o) + 0 1 #nonzeros in Jacobian, obj. gradient + 3 1 #max name lengths: constraints, variables + 0 0 0 0 0 #common exprs: b,c,o,c1,o1 +O0 0 #obj +o54 #sumlist +3 #(n) +o5 #^ +n1.5 +n2 +o5 #^ +o16 #- +v0 #z +n2 +o5 #^ +v0 #z +n2 +x0 #initial guess +r #0 ranges (rhs's) +b #1 bounds (on variables) +3 #z +k0 #intermediate Jacobian column lengths +G0 1 #obj +0 0 +""", + OUT.getvalue(), + ) + ) + + m.c3 = Constraint(expr=m.x == 2) + OUT = io.StringIO() + with LoggingIntercept() as LOG: + with self.assertRaisesRegex( + nl_writer.InfeasibleConstraintException, + r"model contains a trivially infeasible constraint 0.5 == 0.0\*y", + ): + nlinfo = nl_writer.NLWriter().write( + m, OUT, symbolic_solver_labels=True, linear_presolve=True + ) + self.assertEqual(LOG.getvalue(), "") + + m.c1.set_value(m.x >= m.y + m.z + 1.5) + OUT = io.StringIO() + with LoggingIntercept() as LOG: + nlinfo = nl_writer.NLWriter().write( + m, OUT, symbolic_solver_labels=True, linear_presolve=True + ) + self.assertEqual(LOG.getvalue(), "") + + self.assertIs(nlinfo.eliminated_vars[0][0], m.y) + self.assertExpressionsEqual( + nlinfo.eliminated_vars[0][1], LinearExpression([-1.0 * m.z]) + ) + self.assertEqual(nlinfo.eliminated_vars[1], (m.x, 2)) + + self.assertEqual( + *nl_diff( + """g3 1 1 0 # problem unknown + 1 1 1 0 0 #vars, constraints, objectives, ranges, eqns + 0 1 0 0 0 0 #nonlinear constrs, objs; ccons: lin, nonlin, nd, nzlb + 0 0 #network constraints: nonlinear, linear + 0 1 0 #nonlinear vars in constraints, objectives, both + 0 0 0 1 #linear network variables; functions; arith, flags + 0 0 0 0 0 #discrete variables: binary, integer, nonlinear (b,c,o) + 0 1 #nonzeros in Jacobian, obj. gradient + 3 1 #max name lengths: constraints, variables + 0 0 0 0 0 #common exprs: b,c,o,c1,o1 +C0 #c1 +n0 +O0 0 #obj +o54 #sumlist +3 #(n) +o5 #^ +n2 +n2 +o5 #^ +o16 #- +v0 #z +n2 +o5 #^ +v0 #z +n2 +x0 #initial guess +r #1 ranges (rhs's) +1 0.5 #c1 +b #1 bounds (on variables) +3 #z +k0 #intermediate Jacobian column lengths +G0 1 #obj +0 0 """, OUT.getvalue(), )