diff --git a/doc/OnlineDocs/contributed_packages/gdpopt.rst b/doc/OnlineDocs/contributed_packages/gdpopt.rst
index d550b0ced76..670d7633f6d 100644
--- a/doc/OnlineDocs/contributed_packages/gdpopt.rst
+++ b/doc/OnlineDocs/contributed_packages/gdpopt.rst
@@ -93,10 +93,10 @@ An example that includes the modeling approach may be found below.
     Variables:
       x : Size=1, Index=None
           Key  : Lower : Value : Upper : Fixed : Stale : Domain
-          None :  -1.2 :   0.0 :     2 : False : False :  Reals
+          None :  -1.2 :     0 :     2 : False : False :  Reals
       y : Size=1, Index=None
           Key  : Lower : Value : Upper : Fixed : Stale : Domain
-          None :   -10 :   1.0 :    10 : False : False :  Reals
+          None :   -10 :     1 :    10 : False : False :  Reals
   <BLANKLINE>
     Objectives:
       objective : Size=1, Index=None, Active=True
@@ -106,7 +106,7 @@ An example that includes the modeling approach may be found below.
     Constraints:
       c : Size=1
           Key  : Lower : Body : Upper
-          None :   1.0 :  1.0 :   1.0
+          None :   1.0 :    1 :   1.0
 
 .. note:: 
 
diff --git a/pyomo/contrib/appsi/base.py b/pyomo/contrib/appsi/base.py
index 6d2b5ccfcd4..9c7da1eb60b 100644
--- a/pyomo/contrib/appsi/base.py
+++ b/pyomo/contrib/appsi/base.py
@@ -1007,7 +1007,7 @@ def add_constraints(self, cons: List[ConstraintData]):
                 raise ValueError(
                     'constraint {name} has already been added'.format(name=con.name)
                 )
-            self._active_constraints[con] = (con.lower, con.body, con.upper)
+            self._active_constraints[con] = con.expr
             if self.use_extensions and cmodel_available:
                 tmp = cmodel.prep_for_repn(con.body, self._expr_types)
             else:
@@ -1363,40 +1363,13 @@ def update(self, timer: HierarchicalTimer = None):
         cons_to_remove_and_add = dict()
         need_to_set_objective = False
         if config.update_constraints:
-            cons_to_update = list()
-            sos_to_update = list()
             for c in current_cons_dict.keys():
-                if c not in new_cons_set:
-                    cons_to_update.append(c)
+                if c not in new_cons_set and c.expr is not self._active_constraints[c]:
+                    cons_to_remove_and_add[c] = None
+            sos_to_update = []
             for c in current_sos_dict.keys():
                 if c not in new_sos_set:
                     sos_to_update.append(c)
-            for c in cons_to_update:
-                lower, body, upper = self._active_constraints[c]
-                new_lower, new_body, new_upper = c.lower, c.body, c.upper
-                if new_body is not body:
-                    cons_to_remove_and_add[c] = None
-                    continue
-                if new_lower is not lower:
-                    if (
-                        type(new_lower) is NumericConstant
-                        and type(lower) is NumericConstant
-                        and new_lower.value == lower.value
-                    ):
-                        pass
-                    else:
-                        cons_to_remove_and_add[c] = None
-                        continue
-                if new_upper is not upper:
-                    if (
-                        type(new_upper) is NumericConstant
-                        and type(upper) is NumericConstant
-                        and new_upper.value == upper.value
-                    ):
-                        pass
-                    else:
-                        cons_to_remove_and_add[c] = None
-                        continue
             self.remove_sos_constraints(sos_to_update)
             self.add_sos_constraints(sos_to_update)
         timer.stop('cons')
diff --git a/pyomo/contrib/appsi/cmodel/src/fbbt_model.cpp b/pyomo/contrib/appsi/cmodel/src/fbbt_model.cpp
index bd8d7dbf854..ca865d429e2 100644
--- a/pyomo/contrib/appsi/cmodel/src/fbbt_model.cpp
+++ b/pyomo/contrib/appsi/cmodel/src/fbbt_model.cpp
@@ -205,7 +205,7 @@ void process_fbbt_constraints(FBBTModel *model, PyomoExprTypes &expr_types,
   py::handle con_body;
 
   for (py::handle c : cons) {
-    lower_body_upper = active_constraints[c];
+    lower_body_upper = c.attr("to_bounded_expression")();
     con_lb = lower_body_upper[0];
     con_body = lower_body_upper[1];
     con_ub = lower_body_upper[2];
diff --git a/pyomo/contrib/appsi/cmodel/src/lp_writer.cpp b/pyomo/contrib/appsi/cmodel/src/lp_writer.cpp
index 68baf2b8ae8..f33060ee523 100644
--- a/pyomo/contrib/appsi/cmodel/src/lp_writer.cpp
+++ b/pyomo/contrib/appsi/cmodel/src/lp_writer.cpp
@@ -289,7 +289,7 @@ void process_lp_constraints(py::list cons, py::object writer) {
   py::object nonlinear_expr;
   PyomoExprTypes expr_types = PyomoExprTypes();
   for (py::handle c : cons) {
-    lower_body_upper = active_constraints[c];
+    lower_body_upper = c.attr("to_bounded_expression")();
     cname = getSymbol(c, labeler);
     repn = generate_standard_repn(
         lower_body_upper[1], "compute_values"_a = false, "quadratic"_a = true);
diff --git a/pyomo/contrib/appsi/cmodel/src/nl_writer.cpp b/pyomo/contrib/appsi/cmodel/src/nl_writer.cpp
index 8de6cc74ab4..854262496ea 100644
--- a/pyomo/contrib/appsi/cmodel/src/nl_writer.cpp
+++ b/pyomo/contrib/appsi/cmodel/src/nl_writer.cpp
@@ -527,7 +527,7 @@ void process_nl_constraints(NLWriter *nl_writer, PyomoExprTypes &expr_types,
   py::handle repn_nonlinear_expr;
 
   for (py::handle c : cons) {
-    lower_body_upper = active_constraints[c];
+    lower_body_upper = c.attr("to_bounded_expression")();
     repn = generate_standard_repn(
         lower_body_upper[1], "compute_values"_a = false, "quadratic"_a = false);
     _const = appsi_expr_from_pyomo_expr(repn.attr("constant"), var_map,
diff --git a/pyomo/contrib/community_detection/community_graph.py b/pyomo/contrib/community_detection/community_graph.py
index 889940b5996..c67a8cd6690 100644
--- a/pyomo/contrib/community_detection/community_graph.py
+++ b/pyomo/contrib/community_detection/community_graph.py
@@ -123,7 +123,7 @@ def generate_model_graph(
         # Create a list of the variable numbers that occur in the given constraint equation
         numbered_variables_in_constraint_equation = [
             component_number_map[constraint_variable]
-            for constraint_variable in identify_variables(model_constraint.body)
+            for constraint_variable in identify_variables(model_constraint.expr)
         ]
 
         # Update constraint_variable_map
diff --git a/pyomo/contrib/fbbt/expression_bounds_walker.py b/pyomo/contrib/fbbt/expression_bounds_walker.py
index cb287d54df5..3cb32fcbf29 100644
--- a/pyomo/contrib/fbbt/expression_bounds_walker.py
+++ b/pyomo/contrib/fbbt/expression_bounds_walker.py
@@ -232,15 +232,15 @@ def _handle_unknowable_bounds(visitor, node, arg):
 
 
 def _handle_equality(visitor, node, arg1, arg2):
-    return eq(*arg1, *arg2)
+    return eq(*arg1, *arg2, feasibility_tol=visitor.feasibility_tol)
 
 
 def _handle_inequality(visitor, node, arg1, arg2):
-    return ineq(*arg1, *arg2)
+    return ineq(*arg1, *arg2, feasibility_tol=visitor.feasibility_tol)
 
 
 def _handle_ranged(visitor, node, arg1, arg2, arg3):
-    return ranged(*arg1, *arg2, *arg3)
+    return ranged(*arg1, *arg2, *arg3, feasibility_tol=visitor.feasibility_tol)
 
 
 def _handle_expr_if(visitor, node, arg1, arg2, arg3):
diff --git a/pyomo/contrib/fbbt/fbbt.py b/pyomo/contrib/fbbt/fbbt.py
index 1507c4a3cc5..4bd0e4552a1 100644
--- a/pyomo/contrib/fbbt/fbbt.py
+++ b/pyomo/contrib/fbbt/fbbt.py
@@ -12,6 +12,7 @@
 from collections import defaultdict
 from pyomo.common.collections import ComponentMap, ComponentSet
 from pyomo.contrib.fbbt.expression_bounds_walker import ExpressionBoundsVisitor
+import pyomo.core.expr.relational_expr as relational_expr
 import pyomo.core.expr.numeric_expr as numeric_expr
 from pyomo.core.expr.visitor import (
     ExpressionValueVisitor,
@@ -80,6 +81,27 @@ class FBBTException(PyomoException):
     pass
 
 
+def _prop_bnds_leaf_to_root_equality(visitor, node, arg1, arg2):
+    bnds_dict = visitor.bnds_dict
+    bnds_dict[node] = interval.eq(
+        *bnds_dict[arg1], *bnds_dict[arg2], visitor.feasibility_tol
+    )
+
+
+def _prop_bnds_leaf_to_root_inequality(visitor, node, arg1, arg2):
+    bnds_dict = visitor.bnds_dict
+    bnds_dict[node] = interval.ineq(
+        *bnds_dict[arg1], *bnds_dict[arg2], visitor.feasibility_tol
+    )
+
+
+def _prop_bnds_leaf_to_root_ranged(visitor, node, arg1, arg2, arg3):
+    bnds_dict = visitor.bnds_dict
+    bnds_dict[node] = interval.ranged(
+        *bnds_dict[arg1], *bnds_dict[arg2], *bnds_dict[arg3], visitor.feasibility_tol
+    )
+
+
 def _prop_bnds_leaf_to_root_ProductExpression(visitor, node, arg1, arg2):
     """
 
@@ -367,6 +389,9 @@ def _prop_bnds_leaf_to_root_NamedExpression(visitor, node, expr):
         numeric_expr.UnaryFunctionExpression: _prop_bnds_leaf_to_root_UnaryFunctionExpression,
         numeric_expr.LinearExpression: _prop_bnds_leaf_to_root_SumExpression,
         numeric_expr.AbsExpression: _prop_bnds_leaf_to_root_abs,
+        relational_expr.EqualityExpression: _prop_bnds_leaf_to_root_equality,
+        relational_expr.InequalityExpression: _prop_bnds_leaf_to_root_inequality,
+        relational_expr.RangedExpression: _prop_bnds_leaf_to_root_ranged,
         ExpressionData: _prop_bnds_leaf_to_root_NamedExpression,
         ScalarExpression: _prop_bnds_leaf_to_root_NamedExpression,
         ObjectiveData: _prop_bnds_leaf_to_root_NamedExpression,
@@ -375,6 +400,43 @@ def _prop_bnds_leaf_to_root_NamedExpression(visitor, node, expr):
 )
 
 
+def _prop_bnds_root_to_leaf_equality(node, bnds_dict, feasibility_tol):
+    assert bnds_dict[node][1]  # This expression is feasible
+    arg1, arg2 = node.args
+    lb1, ub1 = bnds_dict[arg1]
+    lb2, ub2 = bnds_dict[arg2]
+    bnds_dict[arg1] = bnds_dict[arg2] = max(lb1, lb2), min(ub1, ub2)
+
+
+def _prop_bnds_root_to_leaf_inequality(node, bnds_dict, feasibility_tol):
+    assert bnds_dict[node][1]  # This expression is feasible
+    arg1, arg2 = node.args
+    lb1, ub1 = bnds_dict[arg1]
+    lb2, ub2 = bnds_dict[arg2]
+    if lb1 > lb2:
+        bnds_dict[arg2] = lb1, ub2
+    if ub1 > ub2:
+        bnds_dict[arg1] = lb1, ub2
+
+
+def _prop_bnds_root_to_leaf_ranged(node, bnds_dict, feasibility_tol):
+    assert bnds_dict[node][1]  # This expression is feasible
+    arg1, arg2, arg3 = node.args
+    lb1, ub1 = bnds_dict[arg1]
+    lb2, ub2 = bnds_dict[arg2]
+    lb3, ub3 = bnds_dict[arg3]
+    if lb1 > lb2:
+        bnds_dict[arg2] = lb1, ub2
+        lb2 = lb1
+    if lb2 > lb3:
+        bnds_dict[arg3] = lb2, ub3
+    if ub2 > ub3:
+        bnds_dict[arg2] = lb2, ub3
+        ub2 = ub3
+    if ub1 > ub2:
+        bnds_dict[arg1] = lb1, ub2
+
+
 def _prop_bnds_root_to_leaf_ProductExpression(node, bnds_dict, feasibility_tol):
     """
 
@@ -953,6 +1015,16 @@ def _prop_bnds_root_to_leaf_NamedExpression(node, bnds_dict, feasibility_tol):
 _prop_bnds_root_to_leaf_map[ObjectiveData] = _prop_bnds_root_to_leaf_NamedExpression
 _prop_bnds_root_to_leaf_map[ScalarObjective] = _prop_bnds_root_to_leaf_NamedExpression
 
+_prop_bnds_root_to_leaf_map[relational_expr.EqualityExpression] = (
+    _prop_bnds_root_to_leaf_equality
+)
+_prop_bnds_root_to_leaf_map[relational_expr.InequalityExpression] = (
+    _prop_bnds_root_to_leaf_inequality
+)
+_prop_bnds_root_to_leaf_map[relational_expr.RangedExpression] = (
+    _prop_bnds_root_to_leaf_ranged
+)
+
 
 def _check_and_reset_bounds(var, lb, ub):
     """
@@ -1250,36 +1322,19 @@ def _fbbt_con(con, config):
 
     # a walker to propagate bounds from the variables to the root
     visitorA = _FBBTVisitorLeafToRoot(bnds_dict, feasibility_tol=config.feasibility_tol)
-    visitorA.walk_expression(con.body)
+    visitorA.walk_expression(con.expr)
 
-    # Now we need to replace the bounds in bnds_dict for the root
-    # node with the bounds on the constraint (if those bounds are
-    # better).
-    _lb = value(con.lower)
-    _ub = value(con.upper)
-    if _lb is None:
-        _lb = -interval.inf
-    if _ub is None:
-        _ub = interval.inf
-
-    lb, ub = bnds_dict[con.body]
+    always_feasible, possibly_feasible = bnds_dict[con.expr]
 
     # check if the constraint is infeasible
-    if lb > _ub + config.feasibility_tol or ub < _lb - config.feasibility_tol:
+    if not possibly_feasible:
         raise InfeasibleConstraintException(
             'Detected an infeasible constraint during FBBT: {0}'.format(str(con))
         )
 
     # check if the constraint is always satisfied
-    if config.deactivate_satisfied_constraints:
-        if lb >= _lb - config.feasibility_tol and ub <= _ub + config.feasibility_tol:
-            con.deactivate()
-
-    if _lb > lb:
-        lb = _lb
-    if _ub < ub:
-        ub = _ub
-    bnds_dict[con.body] = (lb, ub)
+    if config.deactivate_satisfied_constraints and always_feasible:
+        con.deactivate()
 
     # Now, propagate bounds back from the root to the variables
     visitorB = _FBBTVisitorRootToLeaf(
@@ -1287,7 +1342,7 @@ def _fbbt_con(con, config):
         integer_tol=config.integer_tol,
         feasibility_tol=config.feasibility_tol,
     )
-    visitorB.dfs_postorder_stack(con.body)
+    visitorB.dfs_postorder_stack(con.expr)
 
     new_var_bounds = ComponentMap()
     for _node, _bnds in bnds_dict.items():
@@ -1334,7 +1389,7 @@ def _fbbt_block(m, config):
     for c in m.component_data_objects(
         ctype=Constraint, active=True, descend_into=config.descend_into, sort=True
     ):
-        for v in identify_variables(c.body):
+        for v in identify_variables(c.expr):
             if v not in var_to_con_map:
                 var_to_con_map[v] = list()
             if v.lb is None:
@@ -1521,14 +1576,14 @@ def __init__(self, comp):
         if comp.ctype == Constraint:
             if comp.is_indexed():
                 for c in comp.values():
-                    self._vars.update(identify_variables(c.body))
+                    self._vars.update(identify_variables(c.expr))
             else:
-                self._vars.update(identify_variables(comp.body))
+                self._vars.update(identify_variables(comp.expr))
         else:
             for c in comp.component_data_objects(
                 Constraint, descend_into=True, active=True, sort=True
             ):
-                self._vars.update(identify_variables(c.body))
+                self._vars.update(identify_variables(c.expr))
 
     def save_bounds(self):
         bnds = ComponentMap()
diff --git a/pyomo/contrib/fbbt/interval.py b/pyomo/contrib/fbbt/interval.py
index a12d1a4529f..4b93d6e3f31 100644
--- a/pyomo/contrib/fbbt/interval.py
+++ b/pyomo/contrib/fbbt/interval.py
@@ -57,7 +57,7 @@ def BoolFlag(val):
     return _true if val else _false
 
 
-def ineq(xl, xu, yl, yu):
+def ineq(xl, xu, yl, yu, feasibility_tol):
     """Compute the "bounds" on an InequalityExpression
 
     Note this is *not* performing interval arithmetic: we are
@@ -67,9 +67,9 @@ def ineq(xl, xu, yl, yu):
 
     """
     ans = []
-    if yl < xu:
+    if yl < xu - feasibility_tol:
         ans.append(_false)
-    if xl <= yu:
+    if xl <= yu + feasibility_tol:
         ans.append(_true)
     assert ans
     if len(ans) == 1:
@@ -77,7 +77,7 @@ def ineq(xl, xu, yl, yu):
     return tuple(ans)
 
 
-def eq(xl, xu, yl, yu):
+def eq(xl, xu, yl, yu, feasibility_tol):
     """Compute the "bounds" on an EqualityExpression
 
     Note this is *not* performing interval arithmetic: we are
@@ -87,9 +87,13 @@ def eq(xl, xu, yl, yu):
 
     """
     ans = []
-    if xl != xu or yl != yu or xl != yl:
+    if (
+        abs(xl - xu) > feasibility_tol
+        or abs(yl - yu) > feasibility_tol
+        or abs(xl - yl) > feasibility_tol
+    ):
         ans.append(_false)
-    if xl <= yu and yl <= xu:
+    if xl <= yu + feasibility_tol and yl <= xu + feasibility_tol:
         ans.append(_true)
     assert ans
     if len(ans) == 1:
@@ -97,7 +101,7 @@ def eq(xl, xu, yl, yu):
     return tuple(ans)
 
 
-def ranged(xl, xu, yl, yu, zl, zu):
+def ranged(xl, xu, yl, yu, zl, zu, feasibility_tol):
     """Compute the "bounds" on a RangedExpression
 
     Note this is *not* performing interval arithmetic: we are
@@ -106,8 +110,8 @@ def ranged(xl, xu, yl, yu, zl, zu):
     `z` and `z`, `y` can be outside the range `x` and `z`, or both.
 
     """
-    lb = ineq(xl, xu, yl, yu)
-    ub = ineq(yl, yu, zl, zu)
+    lb = ineq(xl, xu, yl, yu, feasibility_tol)
+    ub = ineq(yl, yu, zl, zu, feasibility_tol)
     ans = []
     if not lb[0] or not ub[0]:
         ans.append(_false)
diff --git a/pyomo/contrib/fbbt/tests/test_fbbt.py b/pyomo/contrib/fbbt/tests/test_fbbt.py
index f7d08d11215..83e69233bb5 100644
--- a/pyomo/contrib/fbbt/tests/test_fbbt.py
+++ b/pyomo/contrib/fbbt/tests/test_fbbt.py
@@ -1335,3 +1335,31 @@ def test_named_expr(self):
 class TestFBBT(FbbtTestBase, unittest.TestCase):
     def setUp(self) -> None:
         self.tightener = fbbt
+
+    def test_ranged_expression(self):
+        # The python version of FBBT is slightly more flexible than
+        # APPSI's cmodel (it allows - and correctly handles -
+        # RangedExpressions with variable lower / upper bounds).  If we
+        # ever port that functionality into APPSI, then this test can be
+        # moved into the base class.
+        m = pyo.ConcreteModel()
+        m.l = pyo.Var(bounds=(2, None))
+        m.x = pyo.Var()
+        m.u = pyo.Var(bounds=(None, 8))
+        m.c = pyo.Constraint(expr=pyo.inequality(m.l, m.x, m.u))
+        self.tightener(m)
+        self.tightener(m)
+        self.assertEqual(m.l.bounds, (2, 8))
+        self.assertEqual(m.x.bounds, (2, 8))
+        self.assertEqual(m.u.bounds, (2, 8))
+
+        m = pyo.ConcreteModel()
+        m.l = pyo.Var(bounds=(2, None))
+        m.x = pyo.Var(bounds=(3, 7))
+        m.u = pyo.Var(bounds=(None, 8))
+        m.c = pyo.Constraint(expr=pyo.inequality(m.l, m.x, m.u))
+        self.tightener(m)
+        self.tightener(m)
+        self.assertEqual(m.l.bounds, (2, 7))
+        self.assertEqual(m.x.bounds, (3, 7))
+        self.assertEqual(m.u.bounds, (3, 8))
diff --git a/pyomo/contrib/sensitivity_toolbox/sens.py b/pyomo/contrib/sensitivity_toolbox/sens.py
index a3d69b2c7b1..34fbb92327a 100644
--- a/pyomo/contrib/sensitivity_toolbox/sens.py
+++ b/pyomo/contrib/sensitivity_toolbox/sens.py
@@ -9,16 +9,6 @@
 #  This software is distributed under the 3-clause BSD License.
 #  ___________________________________________________________________________
 
-# ______________________________________________________________________________
-#
-# Pyomo: Python Optimization Modeling Objects
-# Copyright (c) 2008-2024
-#  National Technology and Engineering Solutions of Sandia, LLC
-# Under the terms of Contract DE-NA0003525 with National Technology and
-# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
-# rights in this software.
-# This software is distributed under the 3-clause BSD License
-# ______________________________________________________________________________
 from pyomo.environ import (
     Param,
     Var,
@@ -34,8 +24,10 @@
 
 from pyomo.common.sorting import sorted_robust
 from pyomo.core.expr import ExpressionReplacementVisitor
+from pyomo.core.expr.numvalue import is_potentially_variable
 
 from pyomo.common.modeling import unique_component_name
+from pyomo.common.dependencies import numpy as np, scipy
 from pyomo.common.deprecation import deprecated
 from pyomo.common.tempfiles import TempfileManager
 from pyomo.opt import SolverFactory, SolverStatus
@@ -44,8 +36,6 @@
 import os
 import io
 import shutil
-from pyomo.common.dependencies import numpy as np, numpy_available
-from pyomo.common.dependencies import scipy, scipy_available
 
 logger = logging.getLogger('pyomo.contrib.sensitivity_toolbox')
 
@@ -684,25 +674,29 @@ def _replace_parameters_in_constraints(self, variableSubMap):
         )
         last_idx = 0
         for con in old_con_list:
-            if con.equality or con.lower is None or con.upper is None:
-                new_expr = param_replacer.walk_expression(con.expr)
-                block.constList.add(expr=new_expr)
+            new_expr = param_replacer.walk_expression(con.expr)
+            # TODO: We could only create new constraints for expressions
+            # where substitution actually happened, but that breaks some
+            # current tests:
+            #
+            # if new_expr is con.expr:
+            #     # No params were substituted.  We can ignore this constraint
+            #     continue
+            if new_expr.nargs() == 3 and (
+                is_potentially_variable(new_expr.arg(0))
+                or is_potentially_variable(new_expr.arg(2))
+            ):
+                # This is a potentially "invalid" range constraint: it
+                # may now have variables in the bounds.  For safety, we
+                # will split it into two simple inequalities.
+                block.constList.add(expr=(new_expr.arg(0) <= new_expr.arg(1)))
                 last_idx += 1
                 new_old_comp_map[block.constList[last_idx]] = con
-            else:
-                # Constraint must be a ranged inequality, break into
-                # separate constraints
-                new_body = param_replacer.walk_expression(con.body)
-                new_lower = param_replacer.walk_expression(con.lower)
-                new_upper = param_replacer.walk_expression(con.upper)
-
-                # Add constraint for lower bound
-                block.constList.add(expr=(new_lower <= new_body))
+                block.constList.add(expr=(new_expr.arg(1) <= new_expr.arg(2)))
                 last_idx += 1
                 new_old_comp_map[block.constList[last_idx]] = con
-
-                # Add constraint for upper bound
-                block.constList.add(expr=(new_body <= new_upper))
+            else:
+                block.constList.add(expr=new_expr)
                 last_idx += 1
                 new_old_comp_map[block.constList[last_idx]] = con
             con.deactivate()
diff --git a/pyomo/contrib/solver/persistent.py b/pyomo/contrib/solver/persistent.py
index 71322b7043e..65da81a0c08 100644
--- a/pyomo/contrib/solver/persistent.py
+++ b/pyomo/contrib/solver/persistent.py
@@ -111,8 +111,8 @@ def add_constraints(self, cons: List[ConstraintData]):
                 raise ValueError(
                     'constraint {name} has already been added'.format(name=con.name)
                 )
-            self._active_constraints[con] = (con.lower, con.body, con.upper)
-            tmp = collect_vars_and_named_exprs(con.body)
+            self._active_constraints[con] = con.expr
+            tmp = collect_vars_and_named_exprs(con.expr)
             named_exprs, variables, fixed_vars, external_functions = tmp
             self._check_for_new_vars(variables)
             self._named_expressions[con] = [(e, e.expr) for e in named_exprs]
@@ -417,40 +417,13 @@ def update(self, timer: HierarchicalTimer = None):
         cons_to_remove_and_add = {}
         need_to_set_objective = False
         if config.update_constraints:
-            cons_to_update = []
-            sos_to_update = []
             for c in current_cons_dict.keys():
-                if c not in new_cons_set:
-                    cons_to_update.append(c)
+                if c not in new_cons_set and c.expr is not self._active_constraints[c]:
+                    cons_to_remove_and_add[c] = None
+            sos_to_update = []
             for c in current_sos_dict.keys():
                 if c not in new_sos_set:
                     sos_to_update.append(c)
-            for c in cons_to_update:
-                lower, body, upper = self._active_constraints[c]
-                new_lower, new_body, new_upper = c.lower, c.body, c.upper
-                if new_body is not body:
-                    cons_to_remove_and_add[c] = None
-                    continue
-                if new_lower is not lower:
-                    if (
-                        type(new_lower) is NumericConstant
-                        and type(lower) is NumericConstant
-                        and new_lower.value == lower.value
-                    ):
-                        pass
-                    else:
-                        cons_to_remove_and_add[c] = None
-                        continue
-                if new_upper is not upper:
-                    if (
-                        type(new_upper) is NumericConstant
-                        and type(upper) is NumericConstant
-                        and new_upper.value == upper.value
-                    ):
-                        pass
-                    else:
-                        cons_to_remove_and_add[c] = None
-                        continue
             self.remove_sos_constraints(sos_to_update)
             self.add_sos_constraints(sos_to_update)
         timer.stop('cons')
diff --git a/pyomo/core/base/constraint.py b/pyomo/core/base/constraint.py
index e12860991c2..bc9a32f5404 100644
--- a/pyomo/core/base/constraint.py
+++ b/pyomo/core/base/constraint.py
@@ -64,6 +64,7 @@
     InequalityExpression,
     RangedExpression,
 }
+_strict_relational_exprs = {True, (False, True), (True, False), (True, True)}
 _rule_returned_none_error = """Constraint '%s': rule returned None.
 
 Constraint rules must return either a valid expression, a 2- or 3-member
@@ -151,7 +152,7 @@ class ConstraintData(ActiveComponentData):
         _active         A boolean that indicates whether this data is active
     """
 
-    __slots__ = ('_body', '_lower', '_upper', '_expr')
+    __slots__ = ('_expr',)
 
     # Set to true when a constraint class stores its expression
     # in linear canonical form
@@ -167,126 +168,172 @@ def __init__(self, expr=None, component=None):
         self._component = weakref_ref(component) if (component is not None) else None
         self._active = True
 
-        self._body = None
-        self._lower = None
-        self._upper = None
         self._expr = None
         if expr is not None:
             self.set_value(expr)
 
     def __call__(self, exception=True):
         """Compute the value of the body of this constraint."""
-        return value(self.body, exception=exception)
+        body = self.to_bounded_expression()[1]
+        if body.__class__ not in native_numeric_types:
+            body = value(self.body, exception=exception)
+        return body
+
+    def to_bounded_expression(self):
+        """Convert this constraint to a tuple of 3 expressions (lb, body, ub)
+
+        This method "standardizes" the expression into a 3-tuple of
+        expressions: (`lower_bound`, `body`, `upper_bound`).  Upon
+        conversion, `lower_bound` and `upper_bound` are guaranteed to be
+        `None`, numeric constants, or fixed (not necessarily constant)
+        expressions.
+
+        Note
+        ----
+        As this method operates on the *current state* of the
+        expression, any required expression manipulations (and by
+        extension, the result) can change after fixing / unfixing
+        :py:class:`Var` objects.
+
+        Raises
+        ------
+
+        ValueError: Raised if the expression cannot be mapped to this
+            form (i.e., :py:class:`RangedExpression` constraints with
+            variable lower or upper bounds.
+
+        """
+        expr = self._expr
+        if expr.__class__ is RangedExpression:
+            lb, body, ub = ans = expr.args
+            if (
+                lb.__class__ not in native_types
+                and lb.is_potentially_variable()
+                and not lb.is_fixed()
+            ):
+                raise ValueError(
+                    f"Constraint '{self.name}' is a Ranged Inequality with a "
+                    "variable lower bound.  Cannot normalize the "
+                    "constraint or send it to a solver."
+                )
+            if (
+                ub.__class__ not in native_types
+                and ub.is_potentially_variable()
+                and not ub.is_fixed()
+            ):
+                raise ValueError(
+                    f"Constraint '{self.name}' is a Ranged Inequality with a "
+                    "variable upper bound.  Cannot normalize the "
+                    "constraint or send it to a solver."
+                )
+            return ans
+        elif expr is not None:
+            lhs, rhs = expr.args
+            if rhs.__class__ in native_types or not rhs.is_potentially_variable():
+                return rhs if expr.__class__ is EqualityExpression else None, lhs, rhs
+            if lhs.__class__ in native_types or not lhs.is_potentially_variable():
+                return lhs, rhs, lhs if expr.__class__ is EqualityExpression else None
+            return 0 if expr.__class__ is EqualityExpression else None, lhs - rhs, 0
+        return None, None, None
 
     @property
     def body(self):
         """Access the body of a constraint expression."""
-        if self._body is not None:
-            return self._body
-        # The incoming RangedInequality had a potentially variable
-        # bound.  The "body" is fine, but the bounds may not be
-        # (although the responsibility for those checks lies with the
-        # lower/upper properties)
-        body = self._expr.arg(1)
-        if body.__class__ in native_types and body is not None:
-            return as_numeric(body)
-        return body
-
-    def _get_range_bound(self, range_arg):
-        # Equalities and simple inequalities can always be (directly)
-        # reformulated at construction time to force constant bounds.
-        # The only time we need to defer the determination of bounds is
-        # for ranged inequalities that contain non-constant bounds (so
-        # we *know* that the expr will have 3 args)
-        #
-        # It is possible that there is no expression at all (so catch that)
-        if self._expr is None:
-            return None
-        bound = self._expr.arg(range_arg)
-        if not is_fixed(bound):
-            raise ValueError(
-                "Constraint '%s' is a Ranged Inequality with a "
-                "variable %s bound.  Cannot normalize the "
-                "constraint or send it to a solver."
-                % (self.name, {0: 'lower', 2: 'upper'}[range_arg])
-            )
-        return bound
+        try:
+            ans = self.to_bounded_expression()[1]
+        except ValueError:
+            # It is possible that the expression is not currently valid
+            # (i.e., a ranged expression with a non-fixed bound).  We
+            # will catch that exception here and - if this actually *is*
+            # a RangedExpression - return the body.
+            if self._expr.__class__ is RangedExpression:
+                _, ans, _ = self._expr.args
+            else:
+                raise
+        if ans.__class__ in native_types and ans is not None:
+            # Historically, constraint.lower was guaranteed to return a type
+            # derived from Pyomo NumericValue (or None).  Replicate that.
+            #
+            # [JDS 6/2024: it would be nice to remove this behavior,
+            # although possibly unnecessary, as people should use
+            # to_bounded_expression() instead]
+            return as_numeric(ans)
+        return ans
 
     @property
     def lower(self):
         """Access the lower bound of a constraint expression."""
-        bound = self._lower if self._body is not None else self._get_range_bound(0)
-        # Historically, constraint.lower was guaranteed to return a type
-        # derived from Pyomo NumericValue (or None).  Replicate that
-        # functionality, although clients should in almost all cases
-        # move to using ConstraintData.lb instead of accessing
-        # lower/body/upper to avoid the unnecessary creation (and
-        # inevitable destruction) of the NumericConstant wrappers.
-        if bound is None:
-            return None
-        return as_numeric(bound)
+        ans = self.to_bounded_expression()[0]
+        if ans.__class__ in native_types and ans is not None:
+            # Historically, constraint.lower was guaranteed to return a type
+            # derived from Pyomo NumericValue (or None).  Replicate that
+            # functionality, although clients should in almost all cases
+            # move to using ConstraintData.lb instead of accessing
+            # lower/body/upper to avoid the unnecessary creation (and
+            # inevitable destruction) of the NumericConstant wrappers.
+            return as_numeric(ans)
+        return ans
 
     @property
     def upper(self):
         """Access the upper bound of a constraint expression."""
-        bound = self._upper if self._body is not None else self._get_range_bound(2)
-        # Historically, constraint.upper was guaranteed to return a type
-        # derived from Pyomo NumericValue (or None).  Replicate that
-        # functionality, although clients should in almost all cases
-        # move to using ConstraintData.ub instead of accessing
-        # lower/body/upper to avoid the unnecessary creation (and
-        # inevitable destruction) of the NumericConstant wrappers.
-        if bound is None:
-            return None
-        return as_numeric(bound)
+        ans = self.to_bounded_expression()[2]
+        if ans.__class__ in native_types and ans is not None:
+            # Historically, constraint.upper was guaranteed to return a type
+            # derived from Pyomo NumericValue (or None).  Replicate that
+            # functionality, although clients should in almost all cases
+            # move to using ConstraintData.lb instead of accessing
+            # lower/body/upper to avoid the unnecessary creation (and
+            # inevitable destruction) of the NumericConstant wrappers.
+            return as_numeric(ans)
+        return ans
 
     @property
     def lb(self):
         """Access the value of the lower bound of a constraint expression."""
-        bound = self._lower if self._body is not None else self._get_range_bound(0)
+        bound = self.to_bounded_expression()[0]
+        if bound is None:
+            return None
         if bound.__class__ not in native_numeric_types:
-            if bound is None:
-                return None
             bound = float(value(bound))
+        # Note that "bound != bound" catches float('nan')
         if bound in _nonfinite_values or bound != bound:
-            # Note that "bound != bound" catches float('nan')
             if bound == -_inf:
                 return None
-            else:
-                raise ValueError(
-                    "Constraint '%s' created with an invalid non-finite "
-                    "lower bound (%s)." % (self.name, bound)
-                )
+            raise ValueError(
+                f"Constraint '{self.name}' created with an invalid non-finite "
+                f"lower bound ({bound})."
+            )
         return bound
 
     @property
     def ub(self):
         """Access the value of the upper bound of a constraint expression."""
-        bound = self._upper if self._body is not None else self._get_range_bound(2)
+        bound = self.to_bounded_expression()[2]
+        if bound is None:
+            return None
         if bound.__class__ not in native_numeric_types:
-            if bound is None:
-                return None
             bound = float(value(bound))
+        # Note that "bound != bound" catches float('nan')
         if bound in _nonfinite_values or bound != bound:
-            # Note that "bound != bound" catches float('nan')
             if bound == _inf:
                 return None
-            else:
-                raise ValueError(
-                    "Constraint '%s' created with an invalid non-finite "
-                    "upper bound (%s)." % (self.name, bound)
-                )
+            raise ValueError(
+                f"Constraint '{self.name}' created with an invalid non-finite "
+                f"upper bound ({bound})."
+            )
         return bound
 
     @property
     def equality(self):
         """A boolean indicating whether this is an equality constraint."""
-        if self._expr.__class__ is EqualityExpression:
+        expr = self.expr
+        if expr.__class__ is EqualityExpression:
             return True
-        elif self._expr.__class__ is RangedExpression:
+        elif expr.__class__ is RangedExpression:
             # TODO: this is a very restrictive form of structural equality.
-            lb = self._expr.arg(0)
-            if lb is not None and lb is self._expr.arg(2):
+            lb = expr.arg(0)
+            if lb is not None and lb is expr.arg(2):
                 return True
         return False
 
@@ -317,15 +364,22 @@ def expr(self):
 
     def get_value(self):
         """Get the expression on this constraint."""
-        return self._expr
+        return self.expr
 
     def set_value(self, expr):
         """Set the expression on this constraint."""
         # Clear any previously-cached normalized constraint
-        self._lower = self._upper = self._body = self._expr = None
-
+        self._expr = None
         if expr.__class__ in _known_relational_expressions:
+            if getattr(expr, 'strict', False) in _strict_relational_exprs:
+                raise ValueError(
+                    "Constraint '%s' encountered a strict "
+                    "inequality expression ('>' or '<').  All "
+                    "constraints must be formulated using "
+                    "using '<=', '>=', or '=='." % (self.name,)
+                )
             self._expr = expr
+
         elif expr.__class__ is tuple:  # or expr_type is list:
             for arg in expr:
                 if (
@@ -422,120 +476,6 @@ def set_value(self, expr):
                     "\n   (0, model.price[item], 50)" % (self.name, str(expr))
                 )
                 raise ValueError(msg)
-        #
-        # Normalize the incoming expressions, if we can
-        #
-        args = self._expr.args
-        if self._expr.__class__ is InequalityExpression:
-            if self._expr.strict:
-                raise ValueError(
-                    "Constraint '%s' encountered a strict "
-                    "inequality expression ('>' or '< '). All"
-                    " constraints must be formulated using "
-                    "using '<=', '>=', or '=='." % (self.name,)
-                )
-            if (
-                args[1] is None
-                or args[1].__class__ in native_numeric_types
-                or not args[1].is_potentially_variable()
-            ):
-                self._body = args[0]
-                self._upper = args[1]
-            elif (
-                args[0] is None
-                or args[0].__class__ in native_numeric_types
-                or not args[0].is_potentially_variable()
-            ):
-                self._lower = args[0]
-                self._body = args[1]
-            else:
-                self._body = args[0] - args[1]
-                self._upper = 0
-        elif self._expr.__class__ is EqualityExpression:
-            if args[0] is None or args[1] is None:
-                # Error check: ensure equality does not have infinite RHS
-                raise ValueError(
-                    "Equality constraint '%s' defined with "
-                    "non-finite term (%sHS == None)."
-                    % (self.name, 'L' if args[0] is None else 'R')
-                )
-            if (
-                args[0].__class__ in native_numeric_types
-                or not args[0].is_potentially_variable()
-            ):
-                self._lower = self._upper = args[0]
-                self._body = args[1]
-            elif (
-                args[1].__class__ in native_numeric_types
-                or not args[1].is_potentially_variable()
-            ):
-                self._lower = self._upper = args[1]
-                self._body = args[0]
-            else:
-                self._lower = self._upper = 0
-                self._body = args[0] - args[1]
-            # The following logic is caught below when checking for
-            # invalid non-finite bounds:
-            #
-            # if self._lower.__class__ in native_numeric_types and \
-            #    not math.isfinite(self._lower):
-            #     raise ValueError(
-            #         "Equality constraint '%s' defined with "
-            #         "non-finite term." % (self.name))
-        elif self._expr.__class__ is RangedExpression:
-            if any(self._expr.strict):
-                raise ValueError(
-                    "Constraint '%s' encountered a strict "
-                    "inequality expression ('>' or '< '). All"
-                    " constraints must be formulated using "
-                    "using '<=', '>=', or '=='." % (self.name,)
-                )
-            if all(
-                (
-                    arg is None
-                    or arg.__class__ in native_numeric_types
-                    or not arg.is_potentially_variable()
-                )
-                for arg in (args[0], args[2])
-            ):
-                self._lower, self._body, self._upper = args
-        else:
-            # Defensive programming: we currently only support three
-            # relational expression types.  This will only be hit if
-            # someone defines a fourth...
-            raise DeveloperError(
-                "Unrecognized relational expression type: %s"
-                % (self._expr.__class__.__name__,)
-            )
-
-        # We have historically forced the body to be a numeric expression.
-        # TODO: remove this requirement
-        if self._body.__class__ in native_types and self._body is not None:
-            self._body = as_numeric(self._body)
-
-        # We have historically mapped incoming inf to None
-        if self._lower.__class__ in native_numeric_types:
-            bound = self._lower
-            if bound in _nonfinite_values or bound != bound:
-                # Note that "bound != bound" catches float('nan')
-                if bound == -_inf:
-                    self._lower = None
-                else:
-                    raise ValueError(
-                        "Constraint '%s' created with an invalid non-finite "
-                        "lower bound (%s)." % (self.name, self._lower)
-                    )
-        if self._upper.__class__ in native_numeric_types:
-            bound = self._upper
-            if bound in _nonfinite_values or bound != bound:
-                # Note that "bound != bound" catches float('nan')
-                if bound == _inf:
-                    self._upper = None
-                else:
-                    raise ValueError(
-                        "Constraint '%s' created with an invalid non-finite "
-                        "upper bound (%s)." % (self.name, self._upper)
-                    )
 
     def lslack(self):
         """
@@ -911,6 +851,7 @@ class SimpleConstraint(metaclass=RenamedClass):
     {
         'add',
         'set_value',
+        'to_bounded_expression',
         'body',
         'lower',
         'upper',
diff --git a/pyomo/core/kernel/constraint.py b/pyomo/core/kernel/constraint.py
index 6aa4abc4bfe..fe8eb8b2c1f 100644
--- a/pyomo/core/kernel/constraint.py
+++ b/pyomo/core/kernel/constraint.py
@@ -177,6 +177,9 @@ class _MutableBoundsConstraintMixin(object):
     # Define some of the IConstraint abstract methods
     #
 
+    def to_bounded_expression(self):
+        return self.lower, self.body, self.upper
+
     @property
     def lower(self):
         """The expression for the lower bound of the constraint"""
diff --git a/pyomo/core/plugins/transform/add_slack_vars.py b/pyomo/core/plugins/transform/add_slack_vars.py
index 39903384729..31c1107d692 100644
--- a/pyomo/core/plugins/transform/add_slack_vars.py
+++ b/pyomo/core/plugins/transform/add_slack_vars.py
@@ -150,26 +150,29 @@ def _apply_to_impl(self, instance, **kwds):
             if not cons.active:
                 continue
             cons_name = cons.getname(fully_qualified=True)
-            if cons.lower is not None:
+            lower = cons.lower
+            body = cons.body
+            upper = cons.upper
+            if lower is not None:
                 # we add positive slack variable to body:
                 # declare positive slack
                 varName = "_slack_plus_" + cons_name
                 posSlack = Var(within=NonNegativeReals)
                 xblock.add_component(varName, posSlack)
                 # add positive slack to body expression
-                cons._body += posSlack
+                body += posSlack
                 # penalize slack in objective
                 obj_expr += posSlack
-            if cons.upper is not None:
+            if upper is not None:
                 # we subtract a positive slack variable from the body:
                 # declare slack
                 varName = "_slack_minus_" + cons_name
                 negSlack = Var(within=NonNegativeReals)
                 xblock.add_component(varName, negSlack)
                 # add negative slack to body expression
-                cons._body -= negSlack
+                body -= negSlack
                 # add slack to objective
                 obj_expr += negSlack
-
+            cons.set_value((lower, body, upper))
         # make a new objective that minimizes sum of slack variables
         xblock._slack_objective = Objective(expr=obj_expr)
diff --git a/pyomo/core/tests/unit/test_con.py b/pyomo/core/tests/unit/test_con.py
index 15f190e281e..07c7eb3af8e 100644
--- a/pyomo/core/tests/unit/test_con.py
+++ b/pyomo/core/tests/unit/test_con.py
@@ -84,21 +84,55 @@ def rule(model):
         self.assertEqual(model.c.upper, 0)
 
     def test_tuple_construct_inf_equality(self):
-        model = self.create_model(abstract=True)
-
-        def rule(model):
-            return (model.x, float('inf'))
-
-        model.c = Constraint(rule=rule)
-        self.assertRaises(ValueError, model.create_instance)
-
-        model = self.create_model(abstract=True)
-
-        def rule(model):
-            return (float('inf'), model.x)
+        model = self.create_model(abstract=True).create_instance()
 
-        model.c = Constraint(rule=rule)
-        self.assertRaises(ValueError, model.create_instance)
+        model.c = Constraint(expr=(model.x, float('inf')))
+        self.assertEqual(model.c.equality, True)
+        self.assertEqual(model.c.lower, float('inf'))
+        self.assertIs(model.c.body, model.x)
+        self.assertEqual(model.c.upper, float('inf'))
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'c' created with an invalid non-finite lower bound \(inf\).",
+        ):
+            model.c.lb
+        self.assertEqual(model.c.ub, None)
+
+        model.d = Constraint(expr=(float('inf'), model.x))
+        self.assertEqual(model.d.equality, True)
+        self.assertEqual(model.d.lower, float('inf'))
+        self.assertIs(model.d.body, model.x)
+        self.assertEqual(model.d.upper, float('inf'))
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'd' created with an invalid non-finite lower bound \(inf\).",
+        ):
+            model.d.lb
+        self.assertEqual(model.d.ub, None)
+
+        model.e = Constraint(expr=(model.x, float('-inf')))
+        self.assertEqual(model.e.equality, True)
+        self.assertEqual(model.e.lower, float('-inf'))
+        self.assertIs(model.e.body, model.x)
+        self.assertEqual(model.e.upper, float('-inf'))
+        self.assertEqual(model.e.lb, None)
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'e' created with an invalid non-finite upper bound \(-inf\).",
+        ):
+            model.e.ub
+
+        model.f = Constraint(expr=(float('-inf'), model.x))
+        self.assertEqual(model.f.equality, True)
+        self.assertEqual(model.f.lower, float('-inf'))
+        self.assertIs(model.f.body, model.x)
+        self.assertEqual(model.f.upper, float('-inf'))
+        self.assertEqual(model.f.lb, None)
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'f' created with an invalid non-finite upper bound \(-inf\).",
+        ):
+            model.f.ub
 
     def test_tuple_construct_1sided_inequality(self):
         model = self.create_model()
@@ -134,9 +168,11 @@ def rule(model):
         model.c = Constraint(rule=rule)
 
         self.assertEqual(model.c.equality, False)
-        self.assertEqual(model.c.lower, None)
+        self.assertEqual(model.c.lower, float('-inf'))
         self.assertIs(model.c.body, model.y)
         self.assertEqual(model.c.upper, 1)
+        self.assertEqual(model.c.lb, None)
+        self.assertEqual(model.c.ub, 1)
 
         model = self.create_model()
 
@@ -148,7 +184,9 @@ def rule(model):
         self.assertEqual(model.c.equality, False)
         self.assertEqual(model.c.lower, 0)
         self.assertIs(model.c.body, model.y)
-        self.assertEqual(model.c.upper, None)
+        self.assertEqual(model.c.upper, float('inf'))
+        self.assertEqual(model.c.lb, 0)
+        self.assertEqual(model.c.ub, None)
 
     def test_tuple_construct_unbounded_inequality(self):
         model = self.create_model()
@@ -171,9 +209,11 @@ def rule(model):
         model.c = Constraint(rule=rule)
 
         self.assertEqual(model.c.equality, False)
-        self.assertEqual(model.c.lower, None)
+        self.assertEqual(model.c.lower, float('-inf'))
         self.assertIs(model.c.body, model.y)
-        self.assertEqual(model.c.upper, None)
+        self.assertEqual(model.c.upper, float('inf'))
+        self.assertEqual(model.c.lb, None)
+        self.assertEqual(model.c.ub, None)
 
     def test_tuple_construct_invalid_1sided_inequality(self):
         model = self.create_model(abstract=True)
@@ -229,7 +269,11 @@ def rule(model):
         ):
             instance.c.lower
         self.assertIs(instance.c.body, instance.y)
-        self.assertEqual(instance.c.upper, 1)
+        with self.assertRaisesRegex(
+            ValueError,
+            "Constraint 'c' is a Ranged Inequality with a variable lower bound",
+        ):
+            instance.c.upper
         instance.x.fix(3)
         self.assertEqual(value(instance.c.lower), 3)
 
@@ -240,7 +284,11 @@ def rule(model):
 
         model.c = Constraint(rule=rule)
         instance = model.create_instance()
-        self.assertEqual(instance.c.lower, 0)
+        with self.assertRaisesRegex(
+            ValueError,
+            "Constraint 'c' is a Ranged Inequality with a variable upper bound",
+        ):
+            instance.c.lower
         self.assertIs(instance.c.body, instance.y)
         with self.assertRaisesRegex(
             ValueError,
@@ -276,21 +324,23 @@ def rule(model):
         self.assertEqual(model.c.upper, 0)
 
     def test_expr_construct_inf_equality(self):
-        model = self.create_model(abstract=True)
-
-        def rule(model):
-            return model.x == float('inf')
-
-        model.c = Constraint(rule=rule)
-        self.assertRaises(ValueError, model.create_instance)
+        model = self.create_model(abstract=True).create_instance()
 
-        model = self.create_model(abstract=True)
-
-        def rule(model):
-            return float('inf') == model.x
+        model.c = Constraint(expr=model.x == float('inf'))
+        self.assertEqual(model.c.ub, None)
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'c' created with an invalid non-finite lower bound \(inf\).",
+        ):
+            model.c.lb
 
-        model.c = Constraint(rule=rule)
-        self.assertRaises(ValueError, model.create_instance)
+        model.d = Constraint(expr=model.x == float('-inf'))
+        self.assertEqual(model.d.lb, None)
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'd' created with an invalid non-finite upper bound \(-inf\).",
+        ):
+            model.d.ub
 
     def test_expr_construct_1sided_inequality(self):
         model = self.create_model()
@@ -350,9 +400,11 @@ def rule(model):
         model.c = Constraint(rule=rule)
 
         self.assertEqual(model.c.equality, False)
-        self.assertEqual(model.c.lower, None)
+        self.assertIs(model.c.lower, None)
         self.assertIs(model.c.body, model.y)
-        self.assertEqual(model.c.upper, None)
+        self.assertEqual(model.c.upper, float('inf'))
+        self.assertIs(model.c.ub, None)
+        self.assertIs(model.c.lb, None)
 
         model = self.create_model()
 
@@ -362,9 +414,11 @@ def rule(model):
         model.c = Constraint(rule=rule)
 
         self.assertEqual(model.c.equality, False)
-        self.assertEqual(model.c.lower, None)
+        self.assertEqual(model.c.lower, float('-inf'))
         self.assertIs(model.c.body, model.y)
         self.assertEqual(model.c.upper, None)
+        self.assertIs(model.c.ub, None)
+        self.assertIs(model.c.lb, None)
 
         model = self.create_model()
 
@@ -374,9 +428,11 @@ def rule(model):
         model.c = Constraint(rule=rule)
 
         self.assertEqual(model.c.equality, False)
-        self.assertEqual(model.c.lower, None)
+        self.assertEqual(model.c.lower, float('-inf'))
         self.assertIs(model.c.body, model.y)
         self.assertEqual(model.c.upper, None)
+        self.assertIs(model.c.ub, None)
+        self.assertIs(model.c.lb, None)
 
         model = self.create_model()
 
@@ -388,40 +444,40 @@ def rule(model):
         self.assertEqual(model.c.equality, False)
         self.assertEqual(model.c.lower, None)
         self.assertIs(model.c.body, model.y)
-        self.assertEqual(model.c.upper, None)
+        self.assertEqual(model.c.upper, float('inf'))
+        self.assertIs(model.c.ub, None)
+        self.assertIs(model.c.lb, None)
 
     def test_expr_construct_invalid_unbounded_inequality(self):
-        model = self.create_model(abstract=True)
-
-        def rule(model):
-            return model.y <= float('-inf')
-
-        model.c = Constraint(rule=rule)
-        self.assertRaises(ValueError, model.create_instance)
-
-        model = self.create_model(abstract=True)
-
-        def rule(model):
-            return float('inf') <= model.y
-
-        model.c = Constraint(rule=rule)
-        self.assertRaises(ValueError, model.create_instance)
-
-        model = self.create_model(abstract=True)
-
-        def rule(model):
-            return model.y >= float('inf')
+        model = self.create_model(abstract=True).create_instance()
 
-        model.c = Constraint(rule=rule)
-        self.assertRaises(ValueError, model.create_instance)
+        model.c = Constraint(expr=model.y <= float('-inf'))
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'c' created with an invalid non-finite upper bound \(-inf\).",
+        ):
+            model.c.ub
 
-        model = self.create_model(abstract=True)
+        model.d = Constraint(expr=float('inf') <= model.y)
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'd' created with an invalid non-finite lower bound \(inf\).",
+        ):
+            model.d.lb
 
-        def rule(model):
-            return float('-inf') >= model.y
+        model.e = Constraint(expr=model.y >= float('inf'))
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'e' created with an invalid non-finite lower bound \(inf\).",
+        ):
+            model.e.lb
 
-        model.c = Constraint(rule=rule)
-        self.assertRaises(ValueError, model.create_instance)
+        model.f = Constraint(expr=float('-inf') >= model.y)
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'f' created with an invalid non-finite upper bound \(-inf\).",
+        ):
+            model.f.ub
 
     def test_expr_construct_invalid(self):
         m = ConcreteModel()
@@ -484,9 +540,6 @@ def test_nondata_bounds(self):
         model.e2 = Expression()
         model.e3 = Expression()
         model.c.set_value((model.e1, model.e2, model.e3))
-        self.assertIsNone(model.c._lower)
-        self.assertIsNone(model.c._body)
-        self.assertIsNone(model.c._upper)
         self.assertIs(model.c.lower, model.e1)
         self.assertIs(model.c.body, model.e2)
         self.assertIs(model.c.upper, model.e3)
@@ -507,7 +560,7 @@ def test_nondata_bounds(self):
         self.assertIs(model.c.body.expr, model.v[2])
         with self.assertRaisesRegex(
             ValueError,
-            "Constraint 'c' is a Ranged Inequality with a variable upper bound",
+            "Constraint 'c' is a Ranged Inequality with a variable lower bound",
         ):
             model.c.upper
 
@@ -1574,10 +1627,30 @@ def rule1(model):
         self.assertIs(instance.c.body, instance.x)
         with self.assertRaisesRegex(
             ValueError,
-            "Constraint 'c' is a Ranged Inequality with a variable upper bound",
+            "Constraint 'c' is a Ranged Inequality with a variable lower bound",
         ):
             instance.c.upper
+
         #
+        def rule1(model):
+            return (0, model.x, model.z)
+
+        model = AbstractModel()
+        model.x = Var()
+        model.z = Var()
+        model.c = Constraint(rule=rule1)
+        instance = model.create_instance()
+        with self.assertRaisesRegex(
+            ValueError,
+            "Constraint 'c' is a Ranged Inequality with a variable upper bound",
+        ):
+            instance.c.lower
+        self.assertIs(instance.c.body, instance.x)
+        with self.assertRaisesRegex(
+            ValueError,
+            "Constraint 'c' is a Ranged Inequality with a variable upper bound",
+        ):
+            instance.c.upper
 
     def test_expression_constructor_coverage(self):
         def rule1(model):
@@ -1807,23 +1880,39 @@ def test_potentially_variable_bounds(self):
             r"Constraint 'c' is a Ranged Inequality with a variable lower bound",
         ):
             m.c.lower
-        self.assertIs(m.c.upper, m.u)
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'c' is a Ranged Inequality with a variable lower bound",
+        ):
+            self.assertIs(m.c.upper, m.u)
         with self.assertRaisesRegex(
             ValueError,
             r"Constraint 'c' is a Ranged Inequality with a variable lower bound",
         ):
             m.c.lb
-        self.assertEqual(m.c.ub, 10)
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'c' is a Ranged Inequality with a variable lower bound",
+        ):
+            self.assertEqual(m.c.ub, 10)
 
         m.l = 15
         m.u.expr = m.x
-        self.assertIs(m.c.lower, m.l)
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'c' is a Ranged Inequality with a variable upper bound",
+        ):
+            self.assertIs(m.c.lower, m.l)
         with self.assertRaisesRegex(
             ValueError,
             r"Constraint 'c' is a Ranged Inequality with a variable upper bound",
         ):
             m.c.upper
-        self.assertEqual(m.c.lb, 15)
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'c' is a Ranged Inequality with a variable upper bound",
+        ):
+            self.assertEqual(m.c.lb, 15)
         with self.assertRaisesRegex(
             ValueError,
             r"Constraint 'c' is a Ranged Inequality with a variable upper bound",
@@ -1890,17 +1979,16 @@ def test_tuple_expression(self):
         ):
             m.c = (m.x, None)
 
+        # You can create it with an infinite value, but then one of the
+        # bounds will fail:
+        m.c = (m.x, float('inf'))
+        self.assertIsNone(m.c.ub)
         with self.assertRaisesRegex(
             ValueError,
             r"Constraint 'c' created with an invalid "
             r"non-finite lower bound \(inf\)",
         ):
-            m.c = (m.x, float('inf'))
-
-        with self.assertRaisesRegex(
-            ValueError, r"Equality constraint 'c' defined with non-finite term"
-        ):
-            m.c = EqualityExpression((m.x, None))
+            m.c.lb
 
 
 if __name__ == "__main__":
diff --git a/pyomo/gdp/plugins/bilinear.py b/pyomo/gdp/plugins/bilinear.py
index 67390801348..bc91836ea9c 100644
--- a/pyomo/gdp/plugins/bilinear.py
+++ b/pyomo/gdp/plugins/bilinear.py
@@ -77,9 +77,10 @@ def _transformBlock(self, block, instance):
         for component in block.component_data_objects(
             Constraint, active=True, descend_into=False
         ):
-            expr = self._transformExpression(component.body, instance)
-            instance.bilinear_data_.c_body[id(component)] = component.body
-            component._body = expr
+            lb, body, ub = component.to_bounded_expression()
+            expr = self._transformExpression(body, instance)
+            instance.bilinear_data_.c_body[id(component)] = body
+            component.set_value((lb, expr, ub))
 
     def _transformExpression(self, expr, instance):
         if expr.polynomial_degree() > 2:
diff --git a/pyomo/gdp/plugins/cuttingplane.py b/pyomo/gdp/plugins/cuttingplane.py
index 6c77a582987..4cef098eba9 100644
--- a/pyomo/gdp/plugins/cuttingplane.py
+++ b/pyomo/gdp/plugins/cuttingplane.py
@@ -400,7 +400,8 @@ def back_off_constraint_with_calculated_cut_violation(
     val = value(transBlock_rHull.infeasibility_objective) - TOL
     if val <= 0:
         logger.info("\tBacking off cut by %s" % val)
-        cut._body += abs(val)
+        lb, body, ub = cut.to_bounded_expression()
+        cut.set_value((lb, body + abs(val), ub))
     # else there is nothing to do: restore the objective
     transBlock_rHull.del_component(transBlock_rHull.infeasibility_objective)
     transBlock_rHull.separation_objective.activate()
@@ -424,7 +425,8 @@ def back_off_constraint_by_fixed_tolerance(
                    this callback
     TOL: An absolute tolerance to be added to make cut more conservative.
     """
-    cut._body += TOL
+    lb, body, ub = cut.to_bounded_expression()
+    cut.set_value((lb, body + TOL, ub))
 
 
 @TransformationFactory.register(
diff --git a/pyomo/solvers/plugins/solvers/persistent_solver.py b/pyomo/solvers/plugins/solvers/persistent_solver.py
index 3c2a9e52eab..ef883fe5496 100644
--- a/pyomo/solvers/plugins/solvers/persistent_solver.py
+++ b/pyomo/solvers/plugins/solvers/persistent_solver.py
@@ -262,7 +262,9 @@ def _add_and_collect_column_data(self, var, obj_coef, constraints, coefficients)
         coeff_list = list()
         constr_list = list()
         for val, c in zip(coefficients, constraints):
-            c._body += val * var
+            lb, body, ub = c.to_bounded_expression()
+            body += val * var
+            c.set_value((lb, body, ub))
             self._vars_referenced_by_con[c].add(var)
 
             cval = _convert_to_const(val)
diff --git a/pyomo/util/infeasible.py b/pyomo/util/infeasible.py
index 961d5b35036..6a90a4c3773 100644
--- a/pyomo/util/infeasible.py
+++ b/pyomo/util/infeasible.py
@@ -159,7 +159,7 @@ def log_infeasible_constraints(
         if log_variables:
             line += ''.join(
                 f"\n  - VAR {v.name}: {v.value}"
-                for v in identify_variables(constr.body, include_fixed=True)
+                for v in identify_variables(constr.expr, include_fixed=True)
             )
 
         logger.info(line)