Document and Test Dual Sign Convention

doc/OnlineDocs/explanation/experimental/solvers.rst

@@ -357,3 +357,49 @@ implemented. For example, for ``ipopt``:
    :members:
    :show-inheritance:
    :inherited-members:
+
+
+Dual Sign Convention
+--------------------
+For all future solver interfaces, 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.
+

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.

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,177 @@ 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)