diff --git a/pyomo/common/dependencies.py b/pyomo/common/dependencies.py
index 266074201b0..4ddbe1c9ee8 100644
--- a/pyomo/common/dependencies.py
+++ b/pyomo/common/dependencies.py
@@ -778,6 +778,7 @@ def _finalize_matplotlib(module, available):
 def _finalize_numpy(np, available):
     if not available:
         return
+    numeric_types.native_types.add(np.ndarray)
     numeric_types.RegisterLogicalType(np.bool_)
     for t in (
         np.int_,
diff --git a/pyomo/common/unittest.py b/pyomo/common/unittest.py
index 5ab79d668f3..c53b7f50398 100644
--- a/pyomo/common/unittest.py
+++ b/pyomo/common/unittest.py
@@ -31,6 +31,7 @@
 import pytest as pytest
 
 from pyomo.common.collections import Mapping, Sequence
+from pyomo.common.errors import InvalidValueError
 from pyomo.common.tee import capture_output
 
 from unittest import mock
@@ -52,11 +53,11 @@ def _floatOrCall(val):
     """
     try:
         return float(val)
-    except TypeError:
+    except (TypeError, InvalidValueError):
         pass
     try:
         return float(val())
-    except TypeError:
+    except (TypeError, InvalidValueError):
         pass
     try:
         return val.value
diff --git a/pyomo/contrib/pyros/tests/test_grcs.py b/pyomo/contrib/pyros/tests/test_grcs.py
index 3214c987804..0d24b799b99 100644
--- a/pyomo/contrib/pyros/tests/test_grcs.py
+++ b/pyomo/contrib/pyros/tests/test_grcs.py
@@ -4380,9 +4380,13 @@ def test_separation_subsolver_error(self):
             ),
         )
 
+    # FIXME: This test is expected to fail now, as writing out invalid
+    # models generates an exception in the problem writer (and is never
+    # actually sent to the solver)
     @unittest.skipUnless(
         baron_license_is_valid, "Global NLP solver is not available and licensed."
     )
+    @unittest.expectedFailure
     def test_discrete_separation_subsolver_error(self):
         """
         Test PyROS for two-stage problem with discrete type set,
diff --git a/pyomo/core/base/constraint.py b/pyomo/core/base/constraint.py
index 8316088dbf9..53afa35c70c 100644
--- a/pyomo/core/base/constraint.py
+++ b/pyomo/core/base/constraint.py
@@ -187,12 +187,12 @@ def __call__(self, exception=True):
     def has_lb(self):
         """Returns :const:`False` when the lower bound is
         :const:`None` or negative infinity"""
-        return self.lower is not None
+        return self.lb is not None
 
     def has_ub(self):
         """Returns :const:`False` when the upper bound is
         :const:`None` or positive infinity"""
-        return self.upper is not None
+        return self.ub is not None
 
     def lslack(self):
         """
@@ -338,50 +338,39 @@ def body(self):
         """Access the body of a constraint expression."""
         if self._body is not None:
             return self._body
-        else:
-            # 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 _lb(self):
-        if self._body is not None:
-            bound = self._lower
-        elif self._expr is None:
-            return None
-        else:
-            bound = self._expr.arg(0)
-            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, 'lower')
-                )
-        return bound
-
-    def _ub(self):
-        if self._body is not None:
-            bound = self._upper
-        elif self._expr is None:
+        # 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
-        else:
-            bound = self._expr.arg(2)
-            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, 'upper')
-                )
+        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
 
     @property
     def lower(self):
         """Access the lower bound of a constraint expression."""
-        bound = self._lb()
+        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
@@ -395,7 +384,7 @@ def lower(self):
     @property
     def upper(self):
         """Access the upper bound of a constraint expression."""
-        bound = self._ub()
+        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
@@ -409,13 +398,15 @@ def upper(self):
     @property
     def lb(self):
         """Access the value of the lower bound of a constraint expression."""
-        bound = self._lb()
-        if bound.__class__ not in native_types:
-            bound = value(bound)
+        bound = self._lower if self._body is not None else self._get_range_bound(0)
+        if bound.__class__ not in native_numeric_types:
+            if bound is None:
+                return None
+            bound = float(value(bound))
         if bound in _nonfinite_values or bound != bound:
             # Note that "bound != bound" catches float('nan')
             if bound == -_inf:
-                bound = None
+                return None
             else:
                 raise ValueError(
                     "Constraint '%s' created with an invalid non-finite "
@@ -426,13 +417,15 @@ def lb(self):
     @property
     def ub(self):
         """Access the value of the upper bound of a constraint expression."""
-        bound = self._ub()
-        if bound.__class__ not in native_types:
-            bound = value(bound)
+        bound = self._upper if self._body is not None else self._get_range_bound(2)
+        if bound.__class__ not in native_numeric_types:
+            if bound is None:
+                return None
+            bound = float(value(bound))
         if bound in _nonfinite_values or bound != bound:
             # Note that "bound != bound" catches float('nan')
             if bound == _inf:
-                bound = None
+                return None
             else:
                 raise ValueError(
                     "Constraint '%s' created with an invalid non-finite "
diff --git a/pyomo/core/base/var.py b/pyomo/core/base/var.py
index 5dba9d12396..3e6b541ded4 100644
--- a/pyomo/core/base/var.py
+++ b/pyomo/core/base/var.py
@@ -60,8 +60,7 @@
 
 _inf = float('inf')
 _ninf = -_inf
-_no_lower_bound = {None, _ninf}
-_no_upper_bound = {None, _inf}
+_nonfinite_values = {_inf, _ninf}
 _known_global_real_domains = dict(
     [(_, True) for _ in real_global_set_ids]
     + [(_, False) for _ in integer_global_set_ids]
@@ -110,12 +109,12 @@ class _VarData(ComponentData, NumericValue):
     def has_lb(self):
         """Returns :const:`False` when the lower bound is
         :const:`None` or negative infinity"""
-        return self.lb not in _no_lower_bound
+        return self.lb is not None
 
     def has_ub(self):
         """Returns :const:`False` when the upper bound is
         :const:`None` or positive infinity"""
-        return self.ub not in _no_upper_bound
+        return self.ub is not None
 
     # TODO: deprecate this?  Properties are generally preferred over "set*()"
     def setlb(self, val):
@@ -450,54 +449,94 @@ def domain(self, domain):
     def bounds(self):
         # Custom implementation of _VarData.bounds to avoid unnecessary
         # expression generation and duplicate calls to domain.bounds()
-        domain_bounds = self.domain.bounds()
-        if self._lb is None:
-            lb = domain_bounds[0]
-        else:
-            lb = self._lb
-            if lb.__class__ not in native_types:
-                lb = lb()
-            if domain_bounds[0] is not None:
-                lb = max(lb, domain_bounds[0])
-        if self._ub is None:
-            ub = domain_bounds[1]
-        else:
-            ub = self._ub
-            if ub.__class__ not in native_types:
-                ub = ub()
-            if domain_bounds[1] is not None:
-                ub = min(ub, domain_bounds[1])
-        return None if lb == _ninf else lb, None if ub == _inf else ub
+        domain_lb, domain_ub = self.domain.bounds()
+        # lb is the tighter of the domain and bounds
+        lb = self._lb
+        if lb.__class__ not in native_numeric_types:
+            if lb is not None:
+                lb = float(value(lb))
+        if lb in _nonfinite_values or lb != lb:
+            if lb == _ninf:
+                lb = None
+            else:
+                raise ValueError(
+                    "Var '%s' created with an invalid non-finite "
+                    "lower bound (%s)." % (self.name, lb)
+                )
+        if domain_lb is not None:
+            if lb is None:
+                lb = domain_lb
+            else:
+                lb = max(lb, domain_lb)
+        # ub is the tighter of the domain and bounds
+        ub = self._ub
+        if ub.__class__ not in native_numeric_types:
+            if ub is not None:
+                ub = float(value(ub))
+        if ub in _nonfinite_values or ub != ub:
+            if ub == _inf:
+                ub = None
+            else:
+                raise ValueError(
+                    "Var '%s' created with an invalid non-finite "
+                    "upper bound (%s)." % (self.name, ub)
+                )
+        if domain_ub is not None:
+            if ub is None:
+                ub = domain_ub
+            else:
+                ub = min(ub, domain_ub)
+        return lb, ub
 
     @_VarData.lb.getter
     def lb(self):
         # Custom implementation of _VarData.lb to avoid unnecessary
         # expression generation
-        dlb, _ = self.domain.bounds()
-        if self._lb is None:
-            lb = dlb
-        else:
-            lb = self._lb
-            if lb.__class__ not in native_types:
-                lb = lb()
-            if dlb is not None:
-                lb = max(lb, dlb)
-        return None if lb == _ninf else lb
+        domain_lb, domain_ub = self.domain.bounds()
+        # lb is the tighter of the domain and bounds
+        lb = self._lb
+        if lb.__class__ not in native_numeric_types:
+            if lb is not None:
+                lb = float(value(lb))
+        if lb in _nonfinite_values or lb != lb:
+            if lb == _ninf:
+                lb = None
+            else:
+                raise ValueError(
+                    "Var '%s' created with an invalid non-finite "
+                    "lower bound (%s)." % (self.name, lb)
+                )
+        if domain_lb is not None:
+            if lb is None:
+                lb = domain_lb
+            else:
+                lb = max(lb, domain_lb)
+        return lb
 
     @_VarData.ub.getter
     def ub(self):
         # Custom implementation of _VarData.ub to avoid unnecessary
         # expression generation
-        _, dub = self.domain.bounds()
-        if self._ub is None:
-            ub = dub
-        else:
-            ub = self._ub
-            if ub.__class__ not in native_types:
-                ub = ub()
-            if dub is not None:
-                ub = min(ub, dub)
-        return None if ub == _inf else ub
+        domain_lb, domain_ub = self.domain.bounds()
+        # ub is the tighter of the domain and bounds
+        ub = self._ub
+        if ub.__class__ not in native_numeric_types:
+            if ub is not None:
+                ub = float(value(ub))
+        if ub in _nonfinite_values or ub != ub:
+            if ub == _inf:
+                ub = None
+            else:
+                raise ValueError(
+                    "Var '%s' created with an invalid non-finite "
+                    "upper bound (%s)." % (self.name, ub)
+                )
+        if domain_ub is not None:
+            if ub is None:
+                ub = domain_ub
+            else:
+                ub = min(ub, domain_ub)
+        return ub
 
     @property
     def lower(self):
diff --git a/pyomo/core/kernel/constraint.py b/pyomo/core/kernel/constraint.py
index 3c5b8164fe1..7c7969cb025 100644
--- a/pyomo/core/kernel/constraint.py
+++ b/pyomo/core/kernel/constraint.py
@@ -218,7 +218,10 @@ def upper(self, ub):
     @property
     def lb(self):
         """The value of the lower bound of the constraint"""
-        return value(self._lb)
+        lb = value(self.lower)
+        if lb == _neg_inf:
+            return None
+        return lb
 
     @lb.setter
     def lb(self, lb):
@@ -227,7 +230,10 @@ def lb(self, lb):
     @property
     def ub(self):
         """The value of the upper bound of the constraint"""
-        return value(self._ub)
+        ub = value(self.upper)
+        if ub == _pos_inf:
+            return None
+        return ub
 
     @ub.setter
     def ub(self, ub):
diff --git a/pyomo/core/kernel/matrix_constraint.py b/pyomo/core/kernel/matrix_constraint.py
index ab278eb59c2..1dc0fa7ddc3 100644
--- a/pyomo/core/kernel/matrix_constraint.py
+++ b/pyomo/core/kernel/matrix_constraint.py
@@ -19,6 +19,8 @@
 from pyomo.core.kernel.constraint import IConstraint, constraint_tuple
 
 _noarg = object()
+_pos_inf = float('inf')
+_neg_inf = float('-inf')
 
 #
 # Note: This class is experimental. The implementation may
@@ -131,7 +133,10 @@ def upper(self, ub):
     @property
     def lb(self):
         """The value of the lower bound of the constraint"""
-        return value(self.lower)
+        lb = value(self.lower)
+        if lb == _neg_inf:
+            return None
+        return lb
 
     @lb.setter
     def lb(self, lb):
@@ -140,7 +145,10 @@ def lb(self, lb):
     @property
     def ub(self):
         """The value of the upper bound of the constraint"""
-        return value(self.upper)
+        ub = value(self.upper)
+        if ub == _pos_inf:
+            return None
+        return ub
 
     @ub.setter
     def ub(self, ub):
diff --git a/pyomo/core/kernel/variable.py b/pyomo/core/kernel/variable.py
index 75f3118620c..ff54bcb2fca 100644
--- a/pyomo/core/kernel/variable.py
+++ b/pyomo/core/kernel/variable.py
@@ -102,7 +102,10 @@ def bounds(self, bounds_tuple):
     @property
     def lb(self):
         """Return the numeric value of the variable lower bound."""
-        return value(self.lower)
+        lb = value(self.lower)
+        if lb == _neg_inf:
+            return None
+        return lb
 
     @lb.setter
     def lb(self, val):
@@ -111,7 +114,10 @@ def lb(self, val):
     @property
     def ub(self):
         """Return the numeric value of the variable upper bound."""
-        return value(self.upper)
+        ub = value(self.upper)
+        if ub == _pos_inf:
+            return None
+        return ub
 
     @ub.setter
     def ub(self, val):
diff --git a/pyomo/core/tests/unit/kernel/test_constraint.py b/pyomo/core/tests/unit/kernel/test_constraint.py
index f3ddfe50697..f2f219cc66f 100644
--- a/pyomo/core/tests/unit/kernel/test_constraint.py
+++ b/pyomo/core/tests/unit/kernel/test_constraint.py
@@ -109,26 +109,22 @@ def test_has_lb_ub(self):
 
         c.lb = float('-inf')
         self.assertEqual(c.has_lb(), False)
-        self.assertEqual(c.lb, float('-inf'))
-        self.assertEqual(type(c.lb), float)
+        self.assertEqual(c.lb, None)
         self.assertEqual(c.has_ub(), False)
         self.assertIs(c.ub, None)
 
         c.ub = float('inf')
         self.assertEqual(c.has_lb(), False)
-        self.assertEqual(c.lb, float('-inf'))
-        self.assertEqual(type(c.lb), float)
+        self.assertEqual(c.lb, None)
         self.assertEqual(c.has_ub(), False)
-        self.assertEqual(c.ub, float('inf'))
-        self.assertEqual(type(c.ub), float)
+        self.assertEqual(c.ub, None)
 
         c.lb = 0
         self.assertEqual(c.has_lb(), True)
         self.assertEqual(c.lb, 0)
         self.assertEqual(type(c.lb), int)
         self.assertEqual(c.has_ub(), False)
-        self.assertEqual(c.ub, float('inf'))
-        self.assertEqual(type(c.ub), float)
+        self.assertEqual(c.ub, None)
 
         c.ub = 0
         self.assertEqual(c.has_lb(), True)
@@ -159,13 +155,11 @@ def test_has_lb_ub(self):
         self.assertEqual(c.lb, float('inf'))
         self.assertEqual(type(c.lb), float)
         self.assertEqual(c.has_ub(), False)
-        self.assertEqual(c.ub, float('inf'))
-        self.assertEqual(type(c.ub), float)
+        self.assertEqual(c.ub, None)
 
         c.rhs = float('-inf')
         self.assertEqual(c.has_lb(), False)
-        self.assertEqual(c.lb, float('-inf'))
-        self.assertEqual(type(c.lb), float)
+        self.assertEqual(c.lb, None)
         self.assertEqual(c.has_ub(), True)
         self.assertEqual(c.ub, float('-inf'))
         self.assertEqual(type(c.ub), float)
@@ -187,22 +181,22 @@ def test_has_lb_ub(self):
 
         pL.value = float('-inf')
         self.assertEqual(c.has_lb(), False)
-        self.assertEqual(c.lb, float('-inf'))
+        self.assertEqual(c.lb, None)
         with self.assertRaises(ValueError):
             self.assertEqual(c.has_ub(), False)
         self.assertIs(c.upper, pU)
 
         pU.value = float('inf')
         self.assertEqual(c.has_lb(), False)
-        self.assertEqual(c.lb, float('-inf'))
+        self.assertEqual(c.lb, None)
         self.assertEqual(c.has_ub(), False)
-        self.assertEqual(c.ub, float('inf'))
+        self.assertEqual(c.ub, None)
 
         pL.value = 0
         self.assertEqual(c.has_lb(), True)
         self.assertEqual(c.lb, 0)
         self.assertEqual(c.has_ub(), False)
-        self.assertEqual(c.ub, float('inf'))
+        self.assertEqual(c.ub, None)
 
         pU.value = 0
         self.assertEqual(c.has_lb(), True)
@@ -227,12 +221,12 @@ def test_has_lb_ub(self):
         self.assertEqual(c.has_lb(), True)
         self.assertEqual(c.lb, float('inf'))
         self.assertEqual(c.has_ub(), False)
-        self.assertEqual(c.ub, float('inf'))
+        self.assertEqual(c.ub, None)
 
         pL.value = float('-inf')
         c.rhs = pL
         self.assertEqual(c.has_lb(), False)
-        self.assertEqual(c.lb, float('-inf'))
+        self.assertEqual(c.lb, None)
         self.assertEqual(c.has_ub(), True)
         self.assertEqual(c.ub, float('-inf'))
 
@@ -784,18 +778,14 @@ def test_tuple_construct_inf_equality(self):
         c = constraint((x, float('inf')))
         self.assertEqual(c.equality, True)
         self.assertEqual(c.lb, float('inf'))
-        self.assertEqual(type(c.lb), float)
-        self.assertEqual(c.ub, float('inf'))
-        self.assertEqual(type(c.ub), float)
+        self.assertEqual(c.ub, None)
         self.assertEqual(c.rhs, float('inf'))
         self.assertEqual(type(c.rhs), float)
         self.assertIs(c.body, x)
         c = constraint((float('inf'), x))
         self.assertEqual(c.equality, True)
         self.assertEqual(c.lb, float('inf'))
-        self.assertEqual(type(c.lb), float)
-        self.assertEqual(c.ub, float('inf'))
-        self.assertEqual(type(c.ub), float)
+        self.assertEqual(c.ub, None)
         self.assertEqual(c.rhs, float('inf'))
         self.assertEqual(type(c.rhs), float)
         self.assertIs(c.body, x)
@@ -818,8 +808,7 @@ def test_tuple_construct_1sided_inf_inequality(self):
         y = variable()
         c = constraint((float('-inf'), y, 1))
         self.assertEqual(c.equality, False)
-        self.assertEqual(c.lb, float('-inf'))
-        self.assertEqual(type(c.lb), float)
+        self.assertEqual(c.lb, None)
         self.assertIs(c.body, y)
         self.assertEqual(c.ub, 1)
         self.assertEqual(type(c.ub), int)
@@ -829,8 +818,7 @@ def test_tuple_construct_1sided_inf_inequality(self):
         self.assertEqual(c.lb, 0)
         self.assertEqual(type(c.lb), int)
         self.assertIs(c.body, y)
-        self.assertEqual(c.ub, float('inf'))
-        self.assertEqual(type(c.ub), float)
+        self.assertEqual(c.ub, None)
 
     def test_tuple_construct_unbounded_inequality(self):
         y = variable()
@@ -842,11 +830,9 @@ def test_tuple_construct_unbounded_inequality(self):
 
         c = constraint((float('-inf'), y, float('inf')))
         self.assertEqual(c.equality, False)
-        self.assertEqual(c.lb, float('-inf'))
-        self.assertEqual(type(c.lb), float)
+        self.assertEqual(c.lb, None)
         self.assertIs(c.body, y)
-        self.assertEqual(c.ub, float('inf'))
-        self.assertEqual(type(c.ub), float)
+        self.assertEqual(c.ub, None)
 
     def test_tuple_construct_invalid_1sided_inequality(self):
         x = variable()
@@ -913,13 +899,13 @@ def test_expr_construct_equality(self):
         c.expr = x == float('inf')
         self.assertEqual(c.equality, True)
         self.assertEqual(c.lb, float('inf'))
-        self.assertEqual(c.ub, float('inf'))
+        self.assertEqual(c.ub, None)
         self.assertEqual(c.rhs, float('inf'))
         self.assertIs(c.body, x)
         c.expr = float('inf') == x
         self.assertEqual(c.equality, True)
         self.assertEqual(c.lb, float('inf'))
-        self.assertEqual(c.ub, float('inf'))
+        self.assertEqual(c.ub, None)
         self.assertEqual(c.rhs, float('inf'))
         self.assertIs(c.body, x)
 
@@ -964,13 +950,13 @@ def test_expr_construct_inf_equality(self):
         c = constraint(x == float('inf'))
         self.assertEqual(c.equality, True)
         self.assertEqual(c.lb, float('inf'))
-        self.assertEqual(c.ub, float('inf'))
+        self.assertEqual(c.ub, None)
         self.assertEqual(c.rhs, float('inf'))
         self.assertIs(c.body, x)
         c = constraint(float('inf') == x)
         self.assertEqual(c.equality, True)
         self.assertEqual(c.lb, float('inf'))
-        self.assertEqual(c.ub, float('inf'))
+        self.assertEqual(c.ub, None)
         self.assertEqual(c.rhs, float('inf'))
         self.assertIs(c.body, x)
 
@@ -1006,17 +992,17 @@ def test_expr_construct_unbounded_inequality(self):
         self.assertEqual(c.equality, False)
         self.assertIs(c.lb, None)
         self.assertIs(c.body, y)
-        self.assertEqual(c.ub, float('inf'))
+        self.assertEqual(c.ub, None)
 
         c = constraint(float('-inf') <= y)
         self.assertEqual(c.equality, False)
-        self.assertEqual(c.lb, float('-inf'))
+        self.assertEqual(c.lb, None)
         self.assertIs(c.body, y)
         self.assertIs(c.ub, None)
 
         c = constraint(y >= float('-inf'))
         self.assertEqual(c.equality, False)
-        self.assertEqual(c.lb, float('-inf'))
+        self.assertEqual(c.lb, None)
         self.assertIs(c.body, y)
         self.assertIs(c.ub, None)
 
@@ -1024,7 +1010,7 @@ def test_expr_construct_unbounded_inequality(self):
         self.assertEqual(c.equality, False)
         self.assertIs(c.lb, None)
         self.assertIs(c.body, y)
-        self.assertEqual(c.ub, float('inf'))
+        self.assertEqual(c.ub, None)
 
     def test_expr_construct_unbounded_inequality(self):
         y = variable()
@@ -1087,48 +1073,48 @@ def test_equality_infinite(self):
         c.expr = v == float('inf')
         self.assertEqual(c.equality, True)
         self.assertEqual(c.lb, float('inf'))
-        self.assertEqual(c.ub, float('inf'))
+        self.assertEqual(c.ub, None)
         self.assertEqual(c.rhs, float('inf'))
         self.assertIs(c.body, v)
         c.expr = (v, float('inf'))
         self.assertEqual(c.equality, True)
         self.assertEqual(c.lb, float('inf'))
-        self.assertEqual(c.ub, float('inf'))
+        self.assertEqual(c.ub, None)
         self.assertEqual(c.rhs, float('inf'))
         self.assertIs(c.body, v)
         c.expr = float('inf') == v
         self.assertEqual(c.equality, True)
         self.assertEqual(c.lb, float('inf'))
-        self.assertEqual(c.ub, float('inf'))
+        self.assertEqual(c.ub, None)
         self.assertEqual(c.rhs, float('inf'))
         self.assertIs(c.body, v)
         c.expr = (float('inf'), v)
         self.assertEqual(c.equality, True)
         self.assertEqual(c.lb, float('inf'))
-        self.assertEqual(c.ub, float('inf'))
+        self.assertEqual(c.ub, None)
         self.assertEqual(c.rhs, float('inf'))
         self.assertIs(c.body, v)
         c.expr = v == float('-inf')
         self.assertEqual(c.equality, True)
-        self.assertEqual(c.lb, float('-inf'))
+        self.assertEqual(c.lb, None)
         self.assertEqual(c.ub, float('-inf'))
         self.assertEqual(c.rhs, float('-inf'))
         self.assertIs(c.body, v)
         c.expr = (v, float('-inf'))
         self.assertEqual(c.equality, True)
-        self.assertEqual(c.lb, float('-inf'))
+        self.assertEqual(c.lb, None)
         self.assertEqual(c.ub, float('-inf'))
         self.assertEqual(c.rhs, float('-inf'))
         self.assertIs(c.body, v)
         c.expr = float('-inf') == v
         self.assertEqual(c.equality, True)
-        self.assertEqual(c.lb, float('-inf'))
+        self.assertEqual(c.lb, None)
         self.assertEqual(c.ub, float('-inf'))
         self.assertEqual(c.rhs, float('-inf'))
         self.assertIs(c.body, v)
         c.expr = (float('-inf'), v)
         self.assertEqual(c.equality, True)
-        self.assertEqual(c.lb, float('-inf'))
+        self.assertEqual(c.lb, None)
         self.assertEqual(c.ub, float('-inf'))
         self.assertEqual(c.rhs, float('-inf'))
         self.assertIs(c.body, v)
diff --git a/pyomo/core/tests/unit/kernel/test_matrix_constraint.py b/pyomo/core/tests/unit/kernel/test_matrix_constraint.py
index fcb29c9f09c..c986e5eda96 100644
--- a/pyomo/core/tests/unit/kernel/test_matrix_constraint.py
+++ b/pyomo/core/tests/unit/kernel/test_matrix_constraint.py
@@ -115,8 +115,8 @@ def test_init(self):
         self.assertTrue((ctuple.ub == numpy.inf).all())
         self.assertTrue((ctuple.equality == False).all())
         for c in ctuple:
-            self.assertEqual(c.lb, -numpy.inf)
-            self.assertEqual(c.ub, numpy.inf)
+            self.assertEqual(c.lb, None)
+            self.assertEqual(c.ub, None)
             self.assertEqual(c.equality, False)
             self.assertEqual(c(), 0)
             self.assertEqual(c.slack, float('inf'))
@@ -568,8 +568,8 @@ def test_data_bounds(self):
         with self.assertRaises(ValueError):
             ctuple.rhs
         for c in ctuple:
-            self.assertEqual(c.lb, -numpy.inf)
-            self.assertEqual(c.ub, numpy.inf)
+            self.assertEqual(c.lb, None)
+            self.assertEqual(c.ub, None)
             self.assertEqual(c.equality, False)
             with self.assertRaises(ValueError):
                 c.rhs
@@ -652,8 +652,8 @@ def test_data_bounds(self):
         with self.assertRaises(ValueError):
             ctuple.rhs
         for c in ctuple:
-            self.assertEqual(c.lb, -numpy.inf)
-            self.assertEqual(c.ub, numpy.inf)
+            self.assertEqual(c.lb, None)
+            self.assertEqual(c.ub, None)
             self.assertEqual(c.equality, False)
             with self.assertRaises(ValueError):
                 c.rhs
@@ -681,8 +681,8 @@ def test_data_bounds(self):
         with self.assertRaises(ValueError):
             ctuple.rhs
         for c in ctuple:
-            self.assertEqual(c.lb, -numpy.inf)
-            self.assertEqual(c.ub, numpy.inf)
+            self.assertEqual(c.lb, None)
+            self.assertEqual(c.ub, None)
             self.assertEqual(c.equality, False)
             with self.assertRaises(ValueError):
                 c.rhs
diff --git a/pyomo/core/tests/unit/kernel/test_variable.py b/pyomo/core/tests/unit/kernel/test_variable.py
index 4127bd1adfe..e360240f3b2 100644
--- a/pyomo/core/tests/unit/kernel/test_variable.py
+++ b/pyomo/core/tests/unit/kernel/test_variable.py
@@ -577,9 +577,9 @@ def test_binary_type(self):
         self.assertEqual(v.is_discrete(), False)
         self.assertEqual(v.is_binary(), False)
         self.assertEqual(v.is_integer(), False)
-        self.assertEqual(v.lb, float('-inf'))
-        self.assertEqual(v.ub, float('inf'))
-        self.assertEqual(v.bounds, (float('-inf'), float('inf')))
+        self.assertEqual(v.lb, None)
+        self.assertEqual(v.ub, None)
+        self.assertEqual(v.bounds, (None, None))
         self.assertEqual(v.has_lb(), False)
         self.assertEqual(v.has_ub(), False)
 
@@ -662,8 +662,8 @@ def test_binary_type(self):
         self.assertEqual(v.is_discrete(), True)
         self.assertEqual(v.is_binary(), False)
         self.assertEqual(v.is_integer(), True)
-        self.assertEqual(v.lb, float('-inf'))
-        self.assertEqual(v.ub, float('inf'))
+        self.assertEqual(v.lb, None)
+        self.assertEqual(v.ub, None)
         self.assertEqual(v.has_lb(), False)
         self.assertEqual(v.has_ub(), False)
 
@@ -695,21 +695,21 @@ def test_has_lb_ub(self):
 
         v.lb = float('-inf')
         self.assertEqual(v.has_lb(), False)
-        self.assertEqual(v.lb, float('-inf'))
+        self.assertEqual(v.lb, None)
         self.assertEqual(v.has_ub(), False)
         self.assertEqual(v.ub, None)
 
         v.ub = float('inf')
         self.assertEqual(v.has_lb(), False)
-        self.assertEqual(v.lb, float('-inf'))
+        self.assertEqual(v.lb, None)
         self.assertEqual(v.has_ub(), False)
-        self.assertEqual(v.ub, float('inf'))
+        self.assertEqual(v.ub, None)
 
         v.lb = 0
         self.assertEqual(v.has_lb(), True)
         self.assertEqual(v.lb, 0)
         self.assertEqual(v.has_ub(), False)
-        self.assertEqual(v.ub, float('inf'))
+        self.assertEqual(v.ub, None)
 
         v.ub = 0
         self.assertEqual(v.has_lb(), True)
diff --git a/pyomo/core/tests/unit/test_var.py b/pyomo/core/tests/unit/test_var.py
index bbd4bba030f..3df33f0c761 100644
--- a/pyomo/core/tests/unit/test_var.py
+++ b/pyomo/core/tests/unit/test_var.py
@@ -139,6 +139,77 @@ def test_upper_bound_setter(self):
         m.y[1].setub(3)
         self.assertEqual(m.y[1].ub, 3)
 
+    def test_lb(self):
+        m = ConcreteModel()
+        m.x = Var()
+        self.assertEqual(m.x.lb, None)
+        m.x.domain = NonNegativeReals
+        self.assertEqual(m.x.lb, 0)
+        m.x.lb = float('inf')
+        with self.assertRaisesRegex(
+            ValueError, r'invalid non-finite lower bound \(inf\)'
+        ):
+            m.x.lb
+        m.x.lb = float('nan')
+        with self.assertRaisesRegex(
+            ValueError, r'invalid non-finite lower bound \(nan\)'
+        ):
+            m.x.lb
+
+    def test_ub(self):
+        m = ConcreteModel()
+        m.x = Var()
+        self.assertEqual(m.x.ub, None)
+        m.x.domain = NonPositiveReals
+        self.assertEqual(m.x.ub, 0)
+        m.x.ub = float('-inf')
+        with self.assertRaisesRegex(
+            ValueError, r'invalid non-finite upper bound \(-inf\)'
+        ):
+            m.x.ub
+        m.x.ub = float('nan')
+        with self.assertRaisesRegex(
+            ValueError, r'invalid non-finite upper bound \(nan\)'
+        ):
+            m.x.ub
+
+    def test_bounds(self):
+        m = ConcreteModel()
+        m.x = Var()
+        lb, ub = m.x.bounds
+        self.assertEqual(lb, None)
+        self.assertEqual(ub, None)
+        m.x.domain = NonNegativeReals
+        lb, ub = m.x.bounds
+        self.assertEqual(lb, 0)
+        self.assertEqual(ub, None)
+        m.x.lb = float('inf')
+        with self.assertRaisesRegex(
+            ValueError, r'invalid non-finite lower bound \(inf\)'
+        ):
+            lb, ub = m.x.bounds
+        m.x.lb = float('nan')
+        with self.assertRaisesRegex(
+            ValueError, r'invalid non-finite lower bound \(nan\)'
+        ):
+            lb, ub = m.x.bounds
+
+        m.x.lb = None
+        m.x.domain = NonPositiveReals
+        lb, ub = m.x.bounds
+        self.assertEqual(lb, None)
+        self.assertEqual(ub, 0)
+        m.x.ub = float('-inf')
+        with self.assertRaisesRegex(
+            ValueError, r'invalid non-finite upper bound \(-inf\)'
+        ):
+            lb, ub = m.x.bounds
+        m.x.ub = float('nan')
+        with self.assertRaisesRegex(
+            ValueError, r'invalid non-finite upper bound \(nan\)'
+        ):
+            lb, ub = m.x.bounds
+
 
 class PyomoModel(unittest.TestCase):
     def setUp(self):
diff --git a/pyomo/repn/plugins/lp_writer.py b/pyomo/repn/plugins/lp_writer.py
index 8b04ebf1750..5f3972abce6 100644
--- a/pyomo/repn/plugins/lp_writer.py
+++ b/pyomo/repn/plugins/lp_writer.py
@@ -45,6 +45,7 @@
     FileDeterminism_to_SortComponents,
     categorize_valid_components,
     initialize_var_map_from_column_order,
+    int_float,
     ordered_active_constraints,
 )
 
@@ -382,8 +383,11 @@ def write(self, model):
             if with_debug_timing and con.parent_component() is not last_parent:
                 timer.toc('Constraint %s', last_parent, level=logging.DEBUG)
                 last_parent = con.parent_component()
+            # Note: Constraint.lb/ub guarantee a return value that is
+            # either a (finite) native_numeric_type, or None
             lb = con.lb
             ub = con.ub
+
             if lb is None and ub is None:
                 # Note: you *cannot* output trivial (unbounded)
                 # constraints in LP format.  I suppose we could add a
@@ -399,6 +403,8 @@ def write(self, model):
 
             # Pull out the constant: we will move it to the bounds
             offset = repn.constant
+            if offset.__class__ not in int_float:
+                offset = float(offset)
             repn.constant = 0
 
             if repn.linear or getattr(repn, 'quadratic', None):
@@ -423,14 +429,20 @@ def write(self, model):
                 repn.linear[id(ONE_VAR_CONSTANT)] = 0
 
             symbol = labeler(con)
-            if lb == ub and lb is not None:
-                label = f'c_e_{symbol}_'
-                addSymbol(con, label)
-                ostream.write(f'\n{label}:\n')
-                self.write_expression(ostream, repn, False)
-                ostream.write(f'= {(lb - offset)!r}\n')
-            elif lb is not None and lb != neg_inf:
-                if ub is not None and ub != inf:
+            if lb is not None:
+                if ub is None:
+                    label = f'c_l_{symbol}_'
+                    addSymbol(con, label)
+                    ostream.write(f'\n{label}:\n')
+                    self.write_expression(ostream, repn, False)
+                    ostream.write(f'>= {(lb - offset)!r}\n')
+                elif lb == ub:
+                    label = f'c_e_{symbol}_'
+                    addSymbol(con, label)
+                    ostream.write(f'\n{label}:\n')
+                    self.write_expression(ostream, repn, False)
+                    ostream.write(f'= {(lb - offset)!r}\n')
+                else:
                     # We will need the constraint body twice.  Generate
                     # in a buffer so we only have to do that once.
                     buf = StringIO()
@@ -447,13 +459,7 @@ def write(self, model):
                     ostream.write(f'\n{label}:\n')
                     ostream.write(buf)
                     ostream.write(f'<= {(ub - offset)!r}\n')
-                else:
-                    label = f'c_l_{symbol}_'
-                    addSymbol(con, label)
-                    ostream.write(f'\n{label}:\n')
-                    self.write_expression(ostream, repn, False)
-                    ostream.write(f'>= {(lb - offset)!r}\n')
-            elif ub is not None and ub != inf:
+            elif ub is not None:
                 label = f'c_u_{symbol}_'
                 addSymbol(con, label)
                 ostream.write(f'\n{label}:\n')
@@ -495,9 +501,11 @@ def write(self, model):
             elif v.is_integer():
                 integer_vars.append(v_symbol)
 
+            # Note: Var.bounds guarantees the values are either (finite)
+            # native_numeric_types or None
             lb, ub = v.bounds
             lb = '-inf' if lb is None else repr(lb)
-            ub = '+inf' if ub is None or ub == inf else repr(ub)
+            ub = '+inf' if ub is None else repr(ub)
             ostream.write(f"\n   {lb} <= {v_symbol} <= {ub}")
 
         if integer_vars:
@@ -532,6 +540,8 @@ def write(self, model):
             for soscon in sos:
                 ostream.write(f'\n{getSymbol(soscon)}: S{soscon.level}::\n')
                 for v, w in getattr(soscon, 'get_items', soscon.items)():
+                    if w.__class__ not in int_float:
+                        w = float(f)
                     ostream.write(f"  {getSymbol(v)}:{w!r}\n")
 
         ostream.write("\nend\n")
@@ -550,6 +560,8 @@ def write_expression(self, ostream, expr, is_objective):
             for vid, coef in sorted(
                 expr.linear.items(), key=lambda x: getVarOrder(x[0])
             ):
+                if coef.__class__ not in int_float:
+                    coef = float(coef)
                 if coef < 0:
                     ostream.write(f'{coef!r} {getSymbol(getVar(vid))}\n')
                 else:
@@ -571,6 +583,8 @@ def _normalize_constraint(data):
                 else:
                     col = c1, c2
                     sym = f' {getSymbol(getVar(vid1))} * {getSymbol(getVar(vid2))}\n'
+                if coef.__class__ not in int_float:
+                    coef = float(coef)
                 if coef < 0:
                     return col, repr(coef) + sym
                 else:
diff --git a/pyomo/repn/plugins/nl_writer.py b/pyomo/repn/plugins/nl_writer.py
index b94cec257f5..d3846580c52 100644
--- a/pyomo/repn/plugins/nl_writer.py
+++ b/pyomo/repn/plugins/nl_writer.py
@@ -77,6 +77,7 @@
     categorize_valid_components,
     complex_number_error,
     initialize_var_map_from_column_order,
+    int_float,
     ordered_active_constraints,
     nan,
     sum_like_expression_types,
@@ -97,7 +98,6 @@
 inf = float('inf')
 minus_inf = -inf
 
-
 _CONSTANT = ExprType.CONSTANT
 _MONOMIAL = ExprType.MONOMIAL
 _GENERAL = ExprType.GENERAL
@@ -357,8 +357,8 @@ def __init__(self, name, column_order, row_order, obj_order, model_id):
     def update(self, suffix):
         missing_component = missing_other = 0
         self.datatype.add(suffix.datatype)
-        for item in suffix.items():
-            missing = self._store(*item)
+        for obj, val in suffix.items():
+            missing = self._store(obj, val)
             if missing:
                 if missing > 0:
                     missing_component += missing
@@ -405,25 +405,34 @@ def store(self, obj, val):
             )
 
     def _store(self, obj, val):
-        missing_ct = 0
         _id = id(obj)
         if _id in self._column_order:
-            self.var[self._column_order[_id]] = val
+            obj = self.var
+            key = self._column_order[_id]
         elif _id in self._row_order:
-            self.con[self._row_order[_id]] = val
+            obj = self.con
+            key = self._row_order[_id]
         elif _id in self._obj_order:
-            self.obj[self._obj_order[_id]] = val
+            obj = self.obj
+            key = self._obj_order[_id]
         elif _id == self._model_id:
-            self.prob[0] = val
-        elif isinstance(obj, PyomoObject):
-            if obj.is_indexed():
-                for o in obj.values():
-                    missing_ct += self._store(o, val)
-            else:
-                missing_ct = 1
+            obj = self.prob
+            key = 0
         else:
-            missing_ct = -1
-        return missing_ct
+            missing_ct = 0
+            if isinstance(obj, PyomoObject):
+                if obj.is_indexed():
+                    for o in obj.values():
+                        missing_ct += self._store(o, val)
+                else:
+                    missing_ct = 1
+            else:
+                missing_ct = -1
+            return missing_ct
+        if val.__class__ not in int_float:
+            val = float(val)
+        obj[key] = val
+        return 0
 
 
 class _NLWriter_impl(object):
@@ -562,16 +571,17 @@ def write(self, model):
             expr = visitor.walk_expression((con.body, con, 0))
             if expr.named_exprs:
                 self._record_named_expression_usage(expr.named_exprs, con, 0)
+            # Note: Constraint.lb/ub guarantee a return value that is
+            # either a (finite) native_numeric_type, or None
+            const = expr.const
+            if const.__class__ not in int_float:
+                const = float(const)
             lb = con.lb
-            if lb == minus_inf:
-                lb = None
-            elif lb is not None:
-                lb = repr(lb - expr.const)
+            if lb is not None:
+                lb = repr(lb - const)
             ub = con.ub
-            if ub == inf:
-                ub = None
-            elif ub is not None:
-                ub = repr(ub - expr.const)
+            if ub is not None:
+                ub = repr(ub - const)
             _type = _RANGE_TYPE(lb, ub)
             if _type == 4:
                 n_equality += 1
@@ -809,14 +819,12 @@ def write(self, model):
         self.column_order = column_order = {_id: i for i, _id in enumerate(variables)}
         for idx, _id in enumerate(variables):
             v = var_map[_id]
+            # Note: Var.bounds guarantees the values are either (finite)
+            # native_numeric_types or None
             lb, ub = v.bounds
-            if lb == minus_inf:
-                lb = None
-            elif lb is not None:
+            if lb is not None:
                 lb = repr(lb)
-            if ub == inf:
-                ub = None
-            elif ub is not None:
+            if ub is not None:
                 ub = repr(ub)
             variables[idx] = (v, _id, _RANGE_TYPE(lb, ub), lb, ub)
         timer.toc("Computed variable bounds", level=logging.DEBUG)
@@ -1121,6 +1129,7 @@ def write(self, model):
                 if not _vals:
                     continue
                 ostream.write(f"S{_field|_float} {len(_vals)} {name}\n")
+                # Note: _SuffixData store/update guarantee the value is int/float
                 ostream.write(
                     ''.join(f"{_id} {_vals[_id]!r}\n" for _id in sorted(_vals))
                 )
@@ -1210,6 +1219,7 @@ def write(self, model):
                 logger.warning("ignoring 'dual' suffix for Model")
             if _data.con:
                 ostream.write(f"d{len(_data.con)}\n")
+                # Note: _SuffixData store/update guarantee the value is int/float
                 ostream.write(
                     ''.join(f"{_id} {_data.con[_id]!r}\n" for _id in sorted(_data.con))
                 )
@@ -1218,15 +1228,20 @@ def write(self, model):
         # "x" lines (variable initialization)
         #
         _init_lines = [
-            f'{var_idx} {info[0].value!r}{col_comments[var_idx]}\n'
-            for var_idx, info in enumerate(variables)
-            if info[0].value is not None
+            (var_idx, val if val.__class__ in int_float else float(val))
+            for var_idx, val in enumerate(v[0].value for v in variables)
+            if val is not None
         ]
         ostream.write(
             'x%d%s\n'
             % (len(_init_lines), "\t# initial guess" if symbolic_solver_labels else '')
         )
-        ostream.write(''.join(_init_lines))
+        ostream.write(
+            ''.join(
+                f'{var_idx} {val!r}{col_comments[var_idx]}\n'
+                for var_idx, val in _init_lines
+            )
+        )
 
         #
         # "r" lines (constraint bounds)
@@ -1311,7 +1326,10 @@ def write(self, model):
                 continue
             ostream.write(f'J{row_idx} {len(linear)}{row_comments[row_idx]}\n')
             for _id in sorted(linear.keys(), key=column_order.__getitem__):
-                ostream.write(f'{column_order[_id]} {linear[_id]!r}\n')
+                val = linear[_id]
+                if val.__class__ not in int_float:
+                    val = float(val)
+                ostream.write(f'{column_order[_id]} {val!r}\n')
 
         #
         # "G" lines (non-empty terms in the Objective)
@@ -1324,7 +1342,10 @@ def write(self, model):
                 continue
             ostream.write(f'G{obj_idx} {len(linear)}{row_comments[obj_idx + n_cons]}\n')
             for _id in sorted(linear.keys(), key=column_order.__getitem__):
-                ostream.write(f'{column_order[_id]} {linear[_id]!r}\n')
+                val = linear[_id]
+                if val.__class__ not in int_float:
+                    val = float(val)
+                ostream.write(f'{column_order[_id]} {val!r}\n')
 
         # Generate the return information
         info = NLWriterInfo(
@@ -1481,10 +1502,28 @@ def _write_nl_expression(self, repn, include_const):
             if include_const and repn.const:
                 # Add the constant to the NL expression.  AMPL adds the
                 # constant as the second argument, so we will too.
-                nl = self.template.binary_sum + nl + (self.template.const % repn.const)
+                nl = (
+                    self.template.binary_sum
+                    + nl
+                    + (
+                        self.template.const
+                        % (
+                            repn.const
+                            if repn.const.__class__ in int_float
+                            else float(repn.const)
+                        )
+                    )
+                )
             self.ostream.write(nl % tuple(map(self.var_id_to_nl.__getitem__, args)))
         elif include_const:
-            self.ostream.write(self.template.const % repn.const)
+            self.ostream.write(
+                self.template.const
+                % (
+                    repn.const
+                    if repn.const.__class__ in int_float
+                    else float(repn.const)
+                )
+            )
         else:
             self.ostream.write(self.template.const % 0)
 
@@ -1504,7 +1543,10 @@ def _write_v_line(self, expr_id, k):
         #
         ostream.write(f'V{self.next_V_line_id} {len(linear)} {k}{lbl}\n')
         for _id in sorted(linear, key=column_order.__getitem__):
-            ostream.write(f'{column_order[_id]} {linear[_id]!r}\n')
+            val = linear[_id]
+            if val.__class__ not in int_float:
+                val = float(val)
+            ostream.write(f'{column_order[_id]} {val!r}\n')
         self._write_nl_expression(info[1], True)
         self.next_V_line_id += 1
 
@@ -1629,7 +1671,9 @@ def compile_repn(self, visitor, prefix='', args=None, named_exprs=None):
                 args.extend(self.nonlinear[1])
         if self.const:
             nterms += 1
-            nl_sum += template.const % self.const
+            nl_sum += template.const % (
+                self.const if self.const.__class__ in int_float else float(self.const)
+            )
 
         if nterms > 2:
             return (prefix + (template.nary_sum % nterms) + nl_sum, args, named_exprs)
diff --git a/pyomo/repn/tests/ampl/test_nlv2.py b/pyomo/repn/tests/ampl/test_nlv2.py
index a3b2e262e6d..fe04847b6cb 100644
--- a/pyomo/repn/tests/ampl/test_nlv2.py
+++ b/pyomo/repn/tests/ampl/test_nlv2.py
@@ -21,6 +21,7 @@
 from pyomo.repn.util import InvalidNumber
 from pyomo.repn.tests.nl_diff import nl_diff
 
+from pyomo.common.dependencies import numpy, numpy_available
 from pyomo.common.log import LoggingIntercept
 from pyomo.common.tempfiles import TempfileManager
 from pyomo.core.expr import Expr_if, inequality, LinearExpression
@@ -1093,6 +1094,79 @@ def OBJ(m):
 1 2
 2 3
 3 1
+""",
+                OUT.getvalue(),
+            )
+        )
+
+    @unittest.skipUnless(numpy_available, "test requires numpy")
+    def test_nonfloat_constants(self):
+        import pyomo.environ as pyo
+
+        v = numpy.array([[8], [3], [6], [11]])
+        w = numpy.array([[5], [7], [4], [3]])
+        # Create model
+        m = pyo.ConcreteModel()
+        m.I = pyo.Set(initialize=range(4))
+        # Variables: note initialization with non-numeric value
+        m.zero = pyo.Param(initialize=numpy.array([0]), mutable=True)
+        m.one = pyo.Param(initialize=numpy.array([1]), mutable=True)
+        m.x = pyo.Var(m.I, bounds=(m.zero, m.one), domain=pyo.Integers, initialize=True)
+        # Params: initialize with 1-member ndarrays
+        m.limit = pyo.Param(initialize=numpy.array([14]), mutable=True)
+        m.v = pyo.Param(m.I, initialize=v, mutable=True)
+        m.w = pyo.Param(m.I, initialize=w, mutable=True)
+        # Objective: note use of numpy in coefficients
+        m.value = pyo.Objective(expr=pyo.sum_product(m.v, m.x), sense=pyo.maximize)
+        # Constraint: note use of numpy in coefficients and RHS
+        m.weight = pyo.Constraint(expr=pyo.sum_product(m.w, m.x) <= m.limit)
+
+        OUT = io.StringIO()
+        with LoggingIntercept() as LOG:
+            nl_writer.NLWriter().write(m, OUT, symbolic_solver_labels=True)
+        self.assertEqual(LOG.getvalue(), "")
+        self.assertEqual(
+            *nl_diff(
+                """g3 1 1 0       #problem unknown
+ 4 1 1 0 0     #vars, constraints, objectives, ranges, eqns
+ 0 0 0 0 0 0   #nonlinear constrs, objs; ccons: lin, nonlin, nd, nzlb
+ 0 0   #network constraints: nonlinear, linear
+ 0 0 0 #nonlinear vars in constraints, objectives, both
+ 0 0 0 1       #linear network variables; functions; arith, flags
+ 0 4 0 0 0     #discrete variables: binary, integer, nonlinear (b,c,o)
+ 4 4   #nonzeros in Jacobian, obj. gradient
+ 6 4   #max name lengths: constraints, variables
+ 0 0 0 0 0     #common exprs: b,c,o,c1,o1
+C0     #weight
+n0
+O0 1   #value
+n0
+x4     #initial guess
+0 1.0  #x[0]
+1 1.0  #x[1]
+2 1.0  #x[2]
+3 1.0  #x[3]
+r      #1 ranges (rhs's)
+1 14.0 #weight
+b      #4 bounds (on variables)
+0 0 1  #x[0]
+0 0 1  #x[1]
+0 0 1  #x[2]
+0 0 1  #x[3]
+k3     #intermediate Jacobian column lengths
+1
+2
+3
+J0 4   #weight
+0 5
+1 7
+2 4
+3 3
+G0 4   #value
+0 8
+1 3
+2 6
+3 11
 """,
                 OUT.getvalue(),
             )
diff --git a/pyomo/repn/tests/test_linear.py b/pyomo/repn/tests/test_linear.py
index 04b36d67058..cee887e9bab 100644
--- a/pyomo/repn/tests/test_linear.py
+++ b/pyomo/repn/tests/test_linear.py
@@ -1605,6 +1605,6 @@ def test_nonnumeric(self):
         self.assertEqual(cfg.var_map, {})
         self.assertEqual(cfg.var_order, {})
         self.assertEqual(repn.multiplier, 1)
-        self.assertEqual(str(repn.constant), 'InvalidNumber([3 4])')
+        self.assertEqual(str(repn.constant), 'InvalidNumber(array([3, 4]))')
         self.assertEqual(repn.linear, {})
         self.assertEqual(repn.nonlinear, None)
diff --git a/pyomo/repn/tests/test_util.py b/pyomo/repn/tests/test_util.py
index a26bf97ca31..58cbbe049cf 100644
--- a/pyomo/repn/tests/test_util.py
+++ b/pyomo/repn/tests/test_util.py
@@ -177,10 +177,26 @@ def test_InvalidNumber(self):
 
         # TODO: eventually these should raise exceptions
         d = InvalidNumber('abc')
-        self.assertEqual(repr(b), "5")
-        self.assertEqual(repr(d), "'abc'")
-        self.assertEqual(f'{b}', "5")
-        self.assertEqual(f'{d}', "abc")
+        with self.assertRaisesRegex(
+            InvalidValueError,
+            r"Cannot emit InvalidNumber\(5\) in compiled representation",
+        ):
+            repr(b)
+        with self.assertRaisesRegex(
+            InvalidValueError,
+            r"Cannot emit InvalidNumber\('abc'\) in compiled representation",
+        ):
+            repr(d)
+        with self.assertRaisesRegex(
+            InvalidValueError,
+            r"Cannot emit InvalidNumber\(5\) in compiled representation",
+        ):
+            f'{b}'
+        with self.assertRaisesRegex(
+            InvalidValueError,
+            r"Cannot emit InvalidNumber\('abc'\) in compiled representation",
+        ):
+            f'{d}'
 
     def test_apply_operation(self):
         m = ConcreteModel()
diff --git a/pyomo/repn/util.py b/pyomo/repn/util.py
index 73d7c1cc921..c660abdf7a1 100644
--- a/pyomo/repn/util.py
+++ b/pyomo/repn/util.py
@@ -46,6 +46,7 @@
 
 HALT_ON_EVALUATION_ERROR = False
 nan = float('nan')
+int_float = {int, float}
 
 
 class ExprType(enum.IntEnum):
@@ -154,24 +155,35 @@ def __le__(self, other):
     def __ge__(self, other):
         return self._cmp(operator.ge, other)
 
+    def _error(self, msg):
+        causes = list(filter(None, self.causes))
+        if causes:
+            msg += "\nThe InvalidNumber was generated by:\n\t"
+            msg += "\n\t".join(causes)
+        raise InvalidValueError(msg)
+
     def __str__(self):
-        return f'InvalidNumber({self.value})'
+        # We will support simple conversion of InvalidNumber to strings
+        # (for reporting purposes)
+        return f'InvalidNumber({self.value!r})'
 
     def __repr__(self):
-        # FIXME: We want to move to where converting InvalidNumber to
-        # string (with either repr() or f"") should raise a
-        # InvalidValueError.  However, at the moment, this breaks some
-        # tests in PyROS.
-        return repr(self.value)
-        # raise InvalidValueError(f'Cannot emit {str(self)} in compiled representation')
+        # We want attempts to convert InvalidNumber to a string
+        # representation to raise a InvalidValueError.
+        self._error(f'Cannot emit {str(self)} in compiled representation')
 
     def __format__(self, format_spec):
         # FIXME: We want to move to where converting InvalidNumber to
         # string (with either repr() or f"") should raise a
         # InvalidValueError.  However, at the moment, this breaks some
         # tests in PyROS.
-        return self.value.__format__(format_spec)
-        # raise InvalidValueError(f'Cannot emit {str(self)} in compiled representation')
+        # return self.value.__format__(format_spec)
+        self._error(f'Cannot emit {str(self)} in compiled representation')
+
+    def __float__(self):
+        # We want attempts to convert InvalidNumber to a float
+        # representation to raise a InvalidValueError.
+        self._error(f'Cannot convert {str(self)} to float')
 
     def __neg__(self):
         return self._op(operator.neg, self)