From 39a78631a7aecf88687d87c2b8e6226d8ed9a26d Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Sun, 16 Mar 2025 22:55:52 -0600 Subject: [PATCH 1/4] documenting dual sign convention for new solver interfaces --- .../explanation/experimental/solvers.rst | 46 +++++++++++++++++++ doc/OnlineDocs/reference/bibliography.rst | 3 ++ 2 files changed, 49 insertions(+) diff --git a/doc/OnlineDocs/explanation/experimental/solvers.rst b/doc/OnlineDocs/explanation/experimental/solvers.rst index 43d6f0a378a..217b4e35d84 100644 --- a/doc/OnlineDocs/explanation/experimental/solvers.rst +++ b/doc/OnlineDocs/explanation/experimental/solvers.rst @@ -354,3 +354,49 @@ implemented. For example, for ``ipopt``: :members: :show-inheritance: :inherited-members: + + +Dual Sign Convention +-------------------- +Regardless of the solver, Pyomo adopts the following sign convention. Given the problem + +.. math:: + + \begin{align} + \min & f(x) \\ + \text{s.t.} \\ + & c_i(x) = 0 && \forall i \in \mathcal{E} \\ + & g_i(x) \leq 0 && \forall i \in \mathcal{U} \\ + & h_i(x) \geq 0 && \forall i \in \mathcal{L} + \end{align} + +We define the Lagrangian as + +.. math:: + + \begin{align} + L(x, \lambda, \nu, \delta) = f(x) - \sum_{i \in \mathcal{E}} \lambda_i c_i(x) - \sum_{i \in \mathcal{U}} \nu_i g_i(x) - \sum_{i \in \mathcal{L}} \delta_i h_i(x) + \end{align} + +Then, the KKT conditions are [NW99]_ + +.. math:: + + \begin{align} + \nabla_x L(x, \lambda, \nu, \delta) = 0 \\ + c(x) = 0 \\ + g(x) \leq 0 \\ + h(x) \geq 0 \\ + \nu \leq 0 \\ + \delta \geq 0 \\ + \nu_i g_i(x) = 0 \\ + \delta_i h_i(x) = 0 + \end{align} + +Note that this sign convention is based on the ``(lower, body, upper)`` +representation of constraints rather than the expression provided by a +user. Users can specify constraints with variables on both the left- and +right-hand sides of equalities and inequalities. However, the +``(lower, body, upper)`` representation ensures that all variables +appear in the ``body``, matching the form of the problem above. + diff --git a/doc/OnlineDocs/reference/bibliography.rst b/doc/OnlineDocs/reference/bibliography.rst index be64644a0f2..4a4ea37f3e5 100644 --- a/doc/OnlineDocs/reference/bibliography.rst +++ b/doc/OnlineDocs/reference/bibliography.rst @@ -149,3 +149,6 @@ Bibliography .. [VAN10] J. P. Vielma, S. Ahmed, and G. Nemhauser. "Mixed-Integer Models for Non-separable Piecewise Linear Optimization: Unifying framework and Extensions", *Operations Research* 58(2), 303-315. 2010. + +.. [NW99] Nocedal, Jorge, and Stephen J. Wright, eds. Numerical + optimization. New York, NY: Springer New York, 1999. From 3cceefd7469b652e05cb6ed4e684ddbb1069d28d Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Tue, 18 Mar 2025 08:57:48 -0600 Subject: [PATCH 2/4] tests for dual sign convention --- .../solver/tests/solvers/test_solvers.py | 170 ++++++++++++++++++ 1 file changed, 170 insertions(+) diff --git a/pyomo/contrib/solver/tests/solvers/test_solvers.py b/pyomo/contrib/solver/tests/solvers/test_solvers.py index 6de5cf5e30a..de378fcea56 100644 --- a/pyomo/contrib/solver/tests/solvers/test_solvers.py +++ b/pyomo/contrib/solver/tests/solvers/test_solvers.py @@ -33,6 +33,7 @@ from pyomo.contrib.solver.solvers.gurobi_direct import GurobiDirect from pyomo.contrib.solver.solvers.highs import Highs from pyomo.core.expr.numeric_expr import LinearExpression +from pyomo.core.expr.compare import assertExpressionsEqual np, numpy_available = attempt_import('numpy') @@ -75,6 +76,175 @@ def _load_tests(solver_list): return res +class TestDualSignConvention(unittest.TestCase): + @parameterized.expand(input=_load_tests(all_solvers)) + def test_equality( + 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.') + + # for now, we don't support getting duals if linear_presolve = True + if any(name.startswith(i) for i in nl_solvers_set): + if use_presolve: + raise unittest.SkipTest(f'cannot yet get duals if NLWriter presolve is on') + else: + opt.config.writer_config.linear_presolve = False + + m = pe.ConcreteModel() + m.x = pe.Var() + m.y = pe.Var() + m.c1 = pe.Constraint(expr=m.y - m.x - 1 == 0) + m.c2 = pe.Constraint(expr=m.y + m.x - 1 == 0) + m.obj = pe.Objective(expr=m.x + m.y) + res = opt.solve(m) + self.assertEqual(res.solution_status, SolutionStatus.optimal) + self.assertAlmostEqual(m.x.value, 0) + self.assertAlmostEqual(m.y.value, 1) + duals = res.solution_loader.get_duals() + # the sign convention is based on the (lower, body, upper) representation of the constraint, + # so we need to make sure the constraint body is what we expect + assertExpressionsEqual(self, m.c1.body, m.y - m.x - 1) + assertExpressionsEqual(self, m.c2.body, m.y + m.x - 1) + self.assertAlmostEqual(duals[m.c1], 0) + self.assertAlmostEqual(duals[m.c2], 1) + + # multiply the constraints by -1 and make sure the sign of the dual flips + m = pe.ConcreteModel() + m.x = pe.Var() + m.y = pe.Var() + m.c1 = pe.Constraint(expr=-m.y + m.x + 1 == 0) + m.c2 = pe.Constraint(expr=-m.y - m.x + 1 == 0) + m.obj = pe.Objective(expr=m.x + m.y) + res = opt.solve(m) + self.assertEqual(res.solution_status, SolutionStatus.optimal) + self.assertAlmostEqual(m.x.value, 0) + self.assertAlmostEqual(m.y.value, 1) + duals = res.solution_loader.get_duals() + # the sign convention is based on the (lower, body, upper) representation of the constraint, + # so we need to make sure the constraint body is what we expect + assertExpressionsEqual(self, m.c1.body, -m.y + m.x + 1) + assertExpressionsEqual(self, m.c2.body, -m.y - m.x + 1) + self.assertAlmostEqual(duals[m.c1], 0) + self.assertAlmostEqual(duals[m.c2], -1) + + @parameterized.expand(input=_load_tests(all_solvers)) + def test_inequality( + 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.') + + # for now, we don't support getting duals if linear_presolve = True + if any(name.startswith(i) for i in nl_solvers_set): + if use_presolve: + raise unittest.SkipTest(f'cannot yet get duals if NLWriter presolve is on') + else: + opt.config.writer_config.linear_presolve = False + + m = pe.ConcreteModel() + m.x = pe.Var() + m.y = pe.Var() + m.c1 = pe.Constraint(expr=m.x - m.y + 1 <= 0) + m.c2 = pe.Constraint(expr=-m.x - m.y + 1 <= 0) + m.obj = pe.Objective(expr=m.y) + res = opt.solve(m) + self.assertEqual(res.solution_status, SolutionStatus.optimal) + self.assertAlmostEqual(m.x.value, 0) + self.assertAlmostEqual(m.y.value, 1) + duals = res.solution_loader.get_duals() + # the sign convention is based on the (lower, body, upper) representation of the constraint, + # so we need to make sure the constraint body is what we expect + assertExpressionsEqual(self, m.c1.body, m.x - m.y + 1) + assertExpressionsEqual(self, m.c2.body, -m.x - m.y + 1) + self.assertAlmostEqual(m.c1.ub, 0) + self.assertIsNone(m.c1.lb) + self.assertAlmostEqual(m.c2.ub, 0) + self.assertIsNone(m.c2.lb) + self.assertAlmostEqual(duals[m.c1], -0.5) + self.assertAlmostEqual(duals[m.c2], -0.5) + + # multiply the constraints by -1 and make sure the sign of the dual flips + m = pe.ConcreteModel() + m.x = pe.Var() + m.y = pe.Var() + m.c1 = pe.Constraint(expr=-m.x + m.y - 1 >= 0) + m.c2 = pe.Constraint(expr=m.x + m.y - 1 >= 0) + m.obj = pe.Objective(expr=m.y) + res = opt.solve(m) + self.assertEqual(res.solution_status, SolutionStatus.optimal) + self.assertAlmostEqual(m.x.value, 0) + self.assertAlmostEqual(m.y.value, 1) + duals = res.solution_loader.get_duals() + # the sign convention is based on the (lower, body, upper) representation of the constraint, + # so we need to make sure the constraint body is what we expect + assertExpressionsEqual(self, m.c1.body, -m.x + m.y - 1) + assertExpressionsEqual(self, m.c2.body, m.x + m.y - 1) + self.assertAlmostEqual(m.c1.lb, 0) + self.assertIsNone(m.c1.ub) + self.assertAlmostEqual(m.c2.lb, 0) + self.assertIsNone(m.c2.ub) + self.assertAlmostEqual(duals[m.c1], 0.5) + self.assertAlmostEqual(duals[m.c2], 0.5) + + @parameterized.expand(input=_load_tests(all_solvers)) + def test_bounds( + 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.') + + # for now, we don't support getting duals if linear_presolve = True + if any(name.startswith(i) for i in nl_solvers_set): + if use_presolve: + raise unittest.SkipTest(f'cannot yet get duals if NLWriter presolve is on') + else: + opt.config.writer_config.linear_presolve = False + + # first check the lower bound + m = pe.ConcreteModel() + m.x = pe.Var(bounds=(0, None)) + m.y = pe.Var() + m.c1 = pe.Constraint(expr=m.x - m.y + 1 <= 0) + m.obj = pe.Objective(expr=m.y) + res = opt.solve(m) + self.assertEqual(res.solution_status, SolutionStatus.optimal) + self.assertAlmostEqual(m.x.value, 0) + self.assertAlmostEqual(m.y.value, 1) + duals = res.solution_loader.get_duals() + # the sign convention is based on the (lower, body, upper) representation of the constraint, + # so we need to make sure the constraint body is what we expect + assertExpressionsEqual(self, m.c1.body, m.x - m.y + 1) + self.assertAlmostEqual(m.c1.ub, 0) + self.assertIsNone(m.c1.lb) + self.assertAlmostEqual(duals[m.c1], -1) + rc = res.solution_loader.get_reduced_costs() + self.assertAlmostEqual(rc[m.x], 1) + + # now check the upper bound + m = pe.ConcreteModel() + m.x = pe.Var(bounds=(None, 0)) + m.y = pe.Var() + m.c1 = pe.Constraint(expr=-m.x - m.y + 1 <= 0) + m.obj = pe.Objective(expr=m.y) + res = opt.solve(m) + self.assertEqual(res.solution_status, SolutionStatus.optimal) + self.assertAlmostEqual(m.x.value, 0) + self.assertAlmostEqual(m.y.value, 1) + duals = res.solution_loader.get_duals() + # the sign convention is based on the (lower, body, upper) representation of the constraint, + # so we need to make sure the constraint body is what we expect + assertExpressionsEqual(self, m.c1.body, -m.x - m.y + 1) + self.assertAlmostEqual(m.c1.ub, 0) + self.assertIsNone(m.c1.lb) + self.assertAlmostEqual(duals[m.c1], -1) + rc = res.solution_loader.get_reduced_costs() + self.assertAlmostEqual(rc[m.x], -1) + + @unittest.skipUnless(numpy_available, 'numpy is not available') class TestSolvers(unittest.TestCase): @parameterized.expand(input=all_solvers) From 31985df235152161e329580cde068944e93fbd3c Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Thu, 20 Mar 2025 09:06:43 -0600 Subject: [PATCH 3/4] run black --- .../solver/tests/solvers/test_solvers.py | 20 ++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/pyomo/contrib/solver/tests/solvers/test_solvers.py b/pyomo/contrib/solver/tests/solvers/test_solvers.py index de378fcea56..ba268115b7d 100644 --- a/pyomo/contrib/solver/tests/solvers/test_solvers.py +++ b/pyomo/contrib/solver/tests/solvers/test_solvers.py @@ -78,9 +78,7 @@ def _load_tests(solver_list): class TestDualSignConvention(unittest.TestCase): @parameterized.expand(input=_load_tests(all_solvers)) - def test_equality( - self, name: str, opt_class: Type[SolverBase], use_presolve: bool - ): + def test_equality(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.') @@ -88,7 +86,9 @@ def test_equality( # for now, we don't support getting duals if linear_presolve = True if any(name.startswith(i) for i in nl_solvers_set): if use_presolve: - raise unittest.SkipTest(f'cannot yet get duals if NLWriter presolve is on') + raise unittest.SkipTest( + f'cannot yet get duals if NLWriter presolve is on' + ) else: opt.config.writer_config.linear_presolve = False @@ -140,7 +140,9 @@ def test_inequality( # for now, we don't support getting duals if linear_presolve = True if any(name.startswith(i) for i in nl_solvers_set): if use_presolve: - raise unittest.SkipTest(f'cannot yet get duals if NLWriter presolve is on') + raise unittest.SkipTest( + f'cannot yet get duals if NLWriter presolve is on' + ) else: opt.config.writer_config.linear_presolve = False @@ -190,9 +192,7 @@ def test_inequality( self.assertAlmostEqual(duals[m.c2], 0.5) @parameterized.expand(input=_load_tests(all_solvers)) - def test_bounds( - self, name: str, opt_class: Type[SolverBase], use_presolve: bool - ): + def test_bounds(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.') @@ -200,7 +200,9 @@ def test_bounds( # for now, we don't support getting duals if linear_presolve = True if any(name.startswith(i) for i in nl_solvers_set): if use_presolve: - raise unittest.SkipTest(f'cannot yet get duals if NLWriter presolve is on') + raise unittest.SkipTest( + f'cannot yet get duals if NLWriter presolve is on' + ) else: opt.config.writer_config.linear_presolve = False From f6380e36424b32e56368e4cfd6c0eed0a2130dde Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Thu, 20 Mar 2025 09:20:21 -0600 Subject: [PATCH 4/4] update solver docs --- doc/OnlineDocs/explanation/experimental/solvers.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/OnlineDocs/explanation/experimental/solvers.rst b/doc/OnlineDocs/explanation/experimental/solvers.rst index 950e01624f3..14ed53a0f9c 100644 --- a/doc/OnlineDocs/explanation/experimental/solvers.rst +++ b/doc/OnlineDocs/explanation/experimental/solvers.rst @@ -361,7 +361,7 @@ implemented. For example, for ``ipopt``: Dual Sign Convention -------------------- -Regardless of the solver, Pyomo adopts the following sign convention. Given the problem +For all future solver interfaces, Pyomo adopts the following sign convention. Given the problem .. math::