dual sign convention

Author: michaelbynum

doc/OnlineDocs/explanation/experimental/solvers.rst

@@ -437,3 +437,12 @@ As a result, the signs of the duals change. The KKT conditions are
    \nu_i g_i(x) = 0 \\
    \delta_i h_i(x) = 0
    \end{align}
+
+
+Pyomo also supports "range constraints" which are inequalities with both upper and lower bounds, where the bounds are not equal. For example,
+
+.. math::
+
+   -1 \leq x + y \leq 1
+
+These are handled very similarly to variable bounds in terms of dual sign conventions. For these, at most one "side" of the inequality can be active at a time. If neither side is active, then the dual will be zero. If the dual is nonzero, then the dual corresponds to the side of the constraint that is active. The dual for the other side will be implicitly zero. When accessing duals, the keys are the constraints. As a result, there is only one key for a range constraint, even though it is really two constraints. Therefore, the dual for the inactive side will not be reported explicitly. Again, the sign convention is based on the ``(lower, body, upper)`` representation of the constraint. Therefore, the left side of this inequality belongs to :math:`\mathcal{L}` and the right side belongs to :math:`\mathcal{U}`.

pyomo/contrib/solver/tests/solvers/test_solvers.py

@@ -246,6 +246,58 @@ def test_bounds(self, name: str, opt_class: Type[SolverBase], use_presolve: bool
         rc = res.solution_loader.get_reduced_costs()
         self.assertAlmostEqual(rc[m.x], -1)
 
+    @parameterized.expand(input=_load_tests(all_solvers))
+    def test_range(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 = pyo.ConcreteModel()
+        m.x = pyo.Var()
+        m.y = pyo.Var()
+        m.c1 = pyo.Constraint(expr=(-1, m.x + m.y, 1))
+        m.c2 = pyo.Constraint(expr=m.y - m.x >= -1)
+        m.obj = pyo.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)
+        assertExpressionsEqual(self, m.c2.body, m.y - m.x)
+        self.assertAlmostEqual(duals[m.c1], 0.5)
+        self.assertAlmostEqual(duals[m.c2], 0.5)
+
+        # now test the other side of the range constraint
+        m = pyo.ConcreteModel()
+        m.x = pyo.Var()
+        m.y = pyo.Var()
+        m.c1 = pyo.Constraint(expr=(-1, m.x + m.y, 1))
+        m.c2 = pyo.Constraint(expr=m.y - m.x <= 1)
+        m.obj = pyo.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)
+        assertExpressionsEqual(self, m.c2.body, m.y - m.x)
+        self.assertAlmostEqual(duals[m.c1], -0.5)
+        self.assertAlmostEqual(duals[m.c2], -0.5)
+
     @parameterized.expand(input=_load_tests(all_solvers))
     def test_equality_max(self, name: str, opt_class: Type[SolverBase], use_presolve: bool):
         opt: SolverBase = opt_class()
@@ -415,6 +467,58 @@ def test_bounds_max(self, name: str, opt_class: Type[SolverBase], use_presolve:
         rc = res.solution_loader.get_reduced_costs()
         self.assertAlmostEqual(rc[m.x], 1)
 
+    @parameterized.expand(input=_load_tests(all_solvers))
+    def test_range_max(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 = pyo.ConcreteModel()
+        m.x = pyo.Var()
+        m.y = pyo.Var()
+        m.c1 = pyo.Constraint(expr=(-1, m.x + m.y, 1))
+        m.c2 = pyo.Constraint(expr=m.y - m.x >= -1)
+        m.obj = pyo.Objective(expr=-m.y, sense=pyo.maximize)
+        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)
+        assertExpressionsEqual(self, m.c2.body, m.y - m.x)
+        self.assertAlmostEqual(duals[m.c1], -0.5)
+        self.assertAlmostEqual(duals[m.c2], -0.5)
+
+        # now test the other side of the range constraint
+        m = pyo.ConcreteModel()
+        m.x = pyo.Var()
+        m.y = pyo.Var()
+        m.c1 = pyo.Constraint(expr=(-1, m.x + m.y, 1))
+        m.c2 = pyo.Constraint(expr=m.y - m.x <= 1)
+        m.obj = pyo.Objective(expr=m.y, sense=pyo.maximize)
+        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)
+        assertExpressionsEqual(self, m.c2.body, m.y - m.x)
+        self.assertAlmostEqual(duals[m.c1], 0.5)
+        self.assertAlmostEqual(duals[m.c2], 0.5)
+
 
 @unittest.skipUnless(numpy_available, 'numpy is not available')
 class TestSolvers(unittest.TestCase):