diff --git a/doc/OnlineDocs/conf.py b/doc/OnlineDocs/conf.py
index 04fe458407b..1aab4cd76c2 100644
--- a/doc/OnlineDocs/conf.py
+++ b/doc/OnlineDocs/conf.py
@@ -271,7 +271,7 @@ def check_output(self, want, got, optionflags):
     yaml_available, networkx_available, matplotlib_available,
     pympler_available, dill_available,
 )
-pint_available = attempt_import('pint', defer_check=False)[1]
+pint_available = attempt_import('pint', defer_import=False)[1]
 from pyomo.contrib.parmest.parmest import parmest_available
 
 import pyomo.environ as _pe # (trigger all plugin registrations)
diff --git a/doc/OnlineDocs/contributed_packages/pyros.rst b/doc/OnlineDocs/contributed_packages/pyros.rst
index aad37a9685a..76a751dd994 100644
--- a/doc/OnlineDocs/contributed_packages/pyros.rst
+++ b/doc/OnlineDocs/contributed_packages/pyros.rst
@@ -689,7 +689,7 @@ could have been equivalently written as:
   ...     },
   ... )
   ==============================================================================
-  PyROS: The Pyomo Robust Optimization Solver.
+  PyROS: The Pyomo Robust Optimization Solver...
   ...
   ------------------------------------------------------------------------------
   Robust optimal solution identified.
diff --git a/pyomo/common/dependencies.py b/pyomo/common/dependencies.py
index 9e96fdd5860..472b0011edb 100644
--- a/pyomo/common/dependencies.py
+++ b/pyomo/common/dependencies.py
@@ -9,13 +9,16 @@
 #  This software is distributed under the 3-clause BSD License.
 #  ___________________________________________________________________________
 
-from collections.abc import Mapping
 import inspect
 import importlib
 import logging
 import sys
 import warnings
 
+from collections.abc import Mapping
+from types import ModuleType
+from typing import List
+
 from .deprecation import deprecated, deprecation_warning, in_testing_environment
 from .errors import DeferredImportError
 
@@ -127,7 +130,7 @@ class DeferredImportModule(object):
 
     This object is returned by :py:func:`attempt_import()` in lieu of
     the module when :py:func:`attempt_import()` is called with
-    ``defer_check=True``.  Any attempts to access attributes on this
+    ``defer_import=True``.  Any attempts to access attributes on this
     object will trigger the actual module import and return either the
     appropriate module attribute or else if the module import fails,
     raise a :py:class:`.DeferredImportError` exception.
@@ -312,6 +315,12 @@ def __init__(
         self._module = None
         self._available = None
         self._deferred_submodules = deferred_submodules
+        # If this import has a callback, then record this deferred
+        # import so that any direct imports of this module also trigger
+        # the resolution of this DeferredImportIndicator (and the
+        # corresponding callback)
+        if callback is not None:
+            DeferredImportCallbackFinder._callbacks.setdefault(name, []).append(self)
 
     def __bool__(self):
         self.resolve()
@@ -433,6 +442,83 @@ def check_min_version(module, min_version):
 check_min_version._parser = None
 
 
+#
+# Note that we are duck-typing the Loader and MetaPathFinder base
+# classes from importlib.abc.  This avoids a (surprisingly costly)
+# import of importlib.abc
+#
+class DeferredImportCallbackLoader:
+    """Custom Loader to resolve registered :py:class:`DeferredImportIndicator` objects
+
+    This :py:class:`importlib.abc.Loader` loader wraps a regular loader
+    and automatically resolves the registered
+    :py:class:`DeferredImportIndicator` objects after the module is
+    loaded.
+
+    """
+
+    def __init__(self, loader, deferred_indicators: List[DeferredImportIndicator]):
+        self._loader = loader
+        self._deferred_indicators = deferred_indicators
+
+    def module_repr(self, module: ModuleType) -> str:
+        return self._loader.module_repr(module)
+
+    def create_module(self, spec) -> ModuleType:
+        return self._loader.create_module(spec)
+
+    def exec_module(self, module: ModuleType) -> None:
+        self._loader.exec_module(module)
+        # Now that the module has been loaded, trigger the resolution of
+        # the deferred indicators (and their associated callbacks)
+        for deferred in self._deferred_indicators:
+            deferred.resolve()
+
+    def load_module(self, fullname) -> ModuleType:
+        return self._loader.load_module(fullname)
+
+
+class DeferredImportCallbackFinder:
+    """Custom Finder that will wrap the normal loader to trigger callbacks
+
+    This :py:class:`importlib.abc.MetaPathFinder` finder will wrap the
+    normal loader returned by ``PathFinder`` with a loader that will
+    trigger custom callbacks after the module is loaded.  We use this to
+    trigger the post import callbacks registered through
+    :py:func:`attempt_import` even when a user imports the target library
+    directly (and not through attribute access on the
+    :py:class:`DeferredImportModule`.
+
+    """
+
+    _callbacks = {}
+
+    def find_spec(self, fullname, path, target=None):
+        if fullname not in self._callbacks:
+            return None
+
+        spec = importlib.machinery.PathFinder.find_spec(fullname, path, target)
+        if spec is None:
+            # Module not found.  Returning None will proceed to the next
+            # finder (which is likely to raise a ModuleNotFoundError)
+            return None
+        spec.loader = DeferredImportCallbackLoader(
+            spec.loader, self._callbacks[fullname]
+        )
+        return spec
+
+    def invalidate_caches(self):
+        pass
+
+
+_DeferredImportCallbackFinder = DeferredImportCallbackFinder()
+# Insert the DeferredImportCallbackFinder at the beginning of the
+# sys.meta_path so that it is found before the standard finders (so that
+# we can correctly inject the resolution of the DeferredImportIndicators
+# -- which triggers the needed callbacks)
+sys.meta_path.insert(0, _DeferredImportCallbackFinder)
+
+
 def attempt_import(
     name,
     error_message=None,
@@ -441,7 +527,8 @@ def attempt_import(
     alt_names=None,
     callback=None,
     importer=None,
-    defer_check=True,
+    defer_check=None,
+    defer_import=None,
     deferred_submodules=None,
     catch_exceptions=None,
 ):
@@ -495,7 +582,8 @@ def attempt_import(
         The message for the exception raised by :py:class:`ModuleUnavailable`
 
     only_catch_importerror: bool, optional
-        DEPRECATED: use catch_exceptions instead or only_catch_importerror.
+        DEPRECATED: use ``catch_exceptions`` instead of ``only_catch_importerror``.
+
         If True (the default), exceptions other than ``ImportError`` raised
         during module import will be reraised.  If False, any exception
         will result in returning a :py:class:`ModuleUnavailable` object.
@@ -506,13 +594,14 @@ def attempt_import(
         ``module.__version__``)
 
     alt_names: list, optional
-        DEPRECATED: alt_names no longer needs to be specified and is ignored.
+        DEPRECATED: ``alt_names`` no longer needs to be specified and is ignored.
+
         A list of common alternate names by which to look for this
         module in the ``globals()`` namespaces.  For example, the alt_names
         for NumPy would be ``['np']``.  (deprecated in version 6.0)
 
-    callback: function, optional
-        A function with the signature "``fcn(module, available)``" that
+    callback: Callable[[ModuleType, bool], None], optional
+        A function with the signature ``fcn(module, available)`` that
         will be called after the import is first attempted.
 
     importer: function, optional
@@ -522,10 +611,16 @@ def attempt_import(
         want to import/return the first one that is available.
 
     defer_check: bool, optional
-        If True (the default), then the attempted import is deferred
-        until the first use of either the module or the availability
-        flag.  The method will return instances of :py:class:`DeferredImportModule`
-        and :py:class:`DeferredImportIndicator`.
+        DEPRECATED: renamed to ``defer_import`` (deprecated in version 6.7.2.dev0)
+
+    defer_import: bool, optional
+        If True, then the attempted import is deferred until the first
+        use of either the module or the availability flag.  The method
+        will return instances of :py:class:`DeferredImportModule` and
+        :py:class:`DeferredImportIndicator`.  If False, the import will
+        be attempted immediately.  If not set, then the import will be
+        deferred unless the ``name`` is already present in
+        ``sys.modules``.
 
     deferred_submodules: Iterable[str], optional
         If provided, an iterable of submodule names within this module
@@ -576,9 +671,26 @@ def attempt_import(
     if catch_exceptions is None:
         catch_exceptions = (ImportError,)
 
+    if defer_check is not None:
+        deprecation_warning(
+            'defer_check=%s is deprecated.  Please use defer_import' % (defer_check,),
+            version='6.7.2.dev0',
+        )
+        assert defer_import is None
+        defer_import = defer_check
+
+    # If the module has already been imported, there is no reason to
+    # further defer things: just import it.
+    if defer_import is None:
+        if name in sys.modules:
+            defer_import = False
+            deferred_submodules = None
+        else:
+            defer_import = True
+
     # If we are going to defer the check until later, return the
     # deferred import module object
-    if defer_check:
+    if defer_import:
         if deferred_submodules:
             if isinstance(deferred_submodules, Mapping):
                 deprecation_warning(
@@ -621,7 +733,7 @@ def attempt_import(
         return DeferredImportModule(indicator, deferred, None), indicator
 
     if deferred_submodules:
-        raise ValueError("deferred_submodules is only valid if defer_check==True")
+        raise ValueError("deferred_submodules is only valid if defer_import==True")
 
     return _perform_import(
         name=name,
@@ -672,6 +784,11 @@ def _perform_import(
     return module, False
 
 
+@deprecated(
+    "``declare_deferred_modules_as_importable()`` is deprecated.  "
+    "Use the :py:class:`declare_modules_as_importable` context manager.",
+    version='6.7.2.dev0',
+)
 def declare_deferred_modules_as_importable(globals_dict):
     """Make all :py:class:`DeferredImportModules` in ``globals_dict`` importable
 
@@ -698,6 +815,7 @@ def declare_deferred_modules_as_importable(globals_dict):
        ...     'scipy', callback=_finalize_scipy,
        ...     deferred_submodules=['stats', 'sparse', 'spatial', 'integrate'])
        >>> declare_deferred_modules_as_importable(globals())
+       WARNING: DEPRECATED: ...
 
     Which enables users to use:
 
@@ -712,20 +830,87 @@ def declare_deferred_modules_as_importable(globals_dict):
     :py:class:`ModuleUnavailable` instance.
 
     """
-    _global_name = globals_dict['__name__'] + '.'
-    deferred = list(
-        (k, v) for k, v in globals_dict.items() if type(v) is DeferredImportModule
-    )
-    while deferred:
-        name, mod = deferred.pop(0)
-        mod.__path__ = None
-        mod.__spec__ = None
-        sys.modules[_global_name + name] = mod
-        deferred.extend(
-            (name + '.' + k, v)
-            for k, v in mod.__dict__.items()
-            if type(v) is DeferredImportModule
-        )
+    return declare_modules_as_importable(globals_dict).__exit__(None, None, None)
+
+
+class declare_modules_as_importable(object):
+    """Make all :py:class:`ModuleType` and :py:class:`DeferredImportModules`
+    importable through the ``globals_dict`` context.
+
+    This context manager will detect all modules imported into the
+    specified ``globals_dict`` environment (either directly or through
+    :py:func:`attempt_import`) and will make those modules importable
+    from the specified ``globals_dict`` context.  It works by detecting
+    changes in the specified ``globals_dict`` dictionary and adding any new
+    modules or instances of :py:class:`DeferredImportModule` that it
+    finds (and any of their deferred submodules) to ``sys.modules`` so
+    that the modules can be imported through the ``globals_dict``
+    namespace.
+
+    For example, ``pyomo/common/dependencies.py`` declares:
+
+    .. doctest::
+       :hide:
+
+       >>> from pyomo.common.dependencies import (
+       ...     attempt_import, _finalize_scipy, __dict__ as dep_globals,
+       ...     declare_modules_as_importable, )
+       >>> # Sphinx does not provide a proper globals()
+       >>> def globals(): return dep_globals
+
+    .. doctest::
+
+       >>> with declare_modules_as_importable(globals()):
+       ...     scipy, scipy_available = attempt_import(
+       ...        'scipy', callback=_finalize_scipy,
+       ...        deferred_submodules=['stats', 'sparse', 'spatial', 'integrate'])
+
+    Which enables users to use:
+
+    .. doctest::
+
+       >>> import pyomo.common.dependencies.scipy.sparse as spa
+
+    If the deferred import has not yet been triggered, then the
+    :py:class:`DeferredImportModule` is returned and named ``spa``.
+    However, if the import has already been triggered, then ``spa`` will
+    either be the ``scipy.sparse`` module, or a
+    :py:class:`ModuleUnavailable` instance.
+
+    """
+
+    def __init__(self, globals_dict):
+        self.globals_dict = globals_dict
+        self.init_dict = {}
+        self.init_modules = None
+
+    def __enter__(self):
+        self.init_dict.update(self.globals_dict)
+        self.init_modules = set(sys.modules)
+
+    def __exit__(self, exc_type, exc_value, traceback):
+        _global_name = self.globals_dict['__name__'] + '.'
+        deferred = {
+            k: v
+            for k, v in self.globals_dict.items()
+            if k not in self.init_dict
+            and isinstance(v, (ModuleType, DeferredImportModule))
+        }
+        if self.init_modules:
+            for name in set(sys.modules) - self.init_modules:
+                if '.' in name and name.split('.', 1)[0] in deferred:
+                    sys.modules[_global_name + name] = sys.modules[name]
+        while deferred:
+            name, mod = deferred.popitem()
+            mod.__path__ = None
+            mod.__spec__ = None
+            sys.modules[_global_name + name] = mod
+            if isinstance(mod, DeferredImportModule):
+                deferred.update(
+                    (name + '.' + k, v)
+                    for k, v in mod.__dict__.items()
+                    if type(v) is DeferredImportModule
+                )
 
 
 #
@@ -778,11 +963,18 @@ def _finalize_matplotlib(module, available):
     if in_testing_environment():
         module.use('Agg')
     import matplotlib.pyplot
+    import matplotlib.pylab
+    import matplotlib.backends
 
 
 def _finalize_numpy(np, available):
     if not available:
         return
+    # scipy has a dependence on numpy.testing, and if we don't import it
+    # as part of resolving numpy, then certain deferred scipy imports
+    # fail when run under pytest.
+    import numpy.testing
+
     from . import numeric_types
 
     # Register ndarray as a native type to prevent 1-element ndarrays
@@ -842,41 +1034,42 @@ def _pyutilib_importer():
     return importlib.import_module('pyutilib')
 
 
-# Standard libraries that are slower to import and not strictly required
-# on all platforms / situations.
-ctypes, _ = attempt_import(
-    'ctypes', deferred_submodules=['util'], callback=_finalize_ctypes
-)
-random, _ = attempt_import('random')
-
-# Commonly-used optional dependencies
-dill, dill_available = attempt_import('dill')
-mpi4py, mpi4py_available = attempt_import('mpi4py')
-networkx, networkx_available = attempt_import('networkx')
-numpy, numpy_available = attempt_import('numpy', callback=_finalize_numpy)
-pandas, pandas_available = attempt_import('pandas')
-plotly, plotly_available = attempt_import('plotly')
-pympler, pympler_available = attempt_import('pympler', callback=_finalize_pympler)
-pyutilib, pyutilib_available = attempt_import('pyutilib', importer=_pyutilib_importer)
-scipy, scipy_available = attempt_import(
-    'scipy',
-    callback=_finalize_scipy,
-    deferred_submodules=['stats', 'sparse', 'spatial', 'integrate'],
-)
-yaml, yaml_available = attempt_import('yaml', callback=_finalize_yaml)
-
-# Note that matplotlib.pyplot can generate a runtime error on OSX when
-# not installed as a Framework (as is the case in the CI systems)
-matplotlib, matplotlib_available = attempt_import(
-    'matplotlib',
-    callback=_finalize_matplotlib,
-    deferred_submodules=['pyplot', 'pylab'],
-    catch_exceptions=(ImportError, RuntimeError),
-)
+with declare_modules_as_importable(globals()):
+    # Standard libraries that are slower to import and not strictly required
+    # on all platforms / situations.
+    ctypes, _ = attempt_import(
+        'ctypes', deferred_submodules=['util'], callback=_finalize_ctypes
+    )
+    random, _ = attempt_import('random')
+
+    # Commonly-used optional dependencies
+    dill, dill_available = attempt_import('dill')
+    mpi4py, mpi4py_available = attempt_import('mpi4py')
+    networkx, networkx_available = attempt_import('networkx')
+    numpy, numpy_available = attempt_import('numpy', callback=_finalize_numpy)
+    pandas, pandas_available = attempt_import('pandas')
+    plotly, plotly_available = attempt_import('plotly')
+    pympler, pympler_available = attempt_import('pympler', callback=_finalize_pympler)
+    pyutilib, pyutilib_available = attempt_import(
+        'pyutilib', importer=_pyutilib_importer
+    )
+    scipy, scipy_available = attempt_import(
+        'scipy',
+        callback=_finalize_scipy,
+        deferred_submodules=['stats', 'sparse', 'spatial', 'integrate'],
+    )
+    yaml, yaml_available = attempt_import('yaml', callback=_finalize_yaml)
+
+    # Note that matplotlib.pyplot can generate a runtime error on OSX when
+    # not installed as a Framework (as is the case in the CI systems)
+    matplotlib, matplotlib_available = attempt_import(
+        'matplotlib',
+        callback=_finalize_matplotlib,
+        deferred_submodules=['pyplot', 'pylab', 'backends'],
+        catch_exceptions=(ImportError, RuntimeError),
+    )
 
 try:
     import cPickle as pickle
 except ImportError:
     import pickle
-
-declare_deferred_modules_as_importable(globals())
diff --git a/pyomo/common/numeric_types.py b/pyomo/common/numeric_types.py
index ba104203667..ca2ce0f9c6c 100644
--- a/pyomo/common/numeric_types.py
+++ b/pyomo/common/numeric_types.py
@@ -12,7 +12,6 @@
 import logging
 import sys
 
-from pyomo.common.dependencies import numpy_available
 from pyomo.common.deprecation import deprecated, relocated_module_attribute
 from pyomo.common.errors import TemplateExpressionError
 
@@ -208,13 +207,6 @@ def check_if_numeric_type(obj):
     if obj_class in native_types:
         return obj_class in native_numeric_types
 
-    if 'numpy' in obj_class.__module__:
-        # trigger the resolution of numpy_available and check if this
-        # type was automatically registered
-        bool(numpy_available)
-        if obj_class in native_types:
-            return obj_class in native_numeric_types
-
     try:
         obj_plus_0 = obj + 0
         obj_p0_class = obj_plus_0.__class__
diff --git a/pyomo/common/tests/dep_mod.py b/pyomo/common/tests/dep_mod.py
index f6add596ed4..34c7219c6eb 100644
--- a/pyomo/common/tests/dep_mod.py
+++ b/pyomo/common/tests/dep_mod.py
@@ -13,8 +13,8 @@
 
 __version__ = '1.5'
 
-numpy, numpy_available = attempt_import('numpy', defer_check=True)
+numpy, numpy_available = attempt_import('numpy', defer_import=True)
 
 bogus_nonexisting_module, bogus_nonexisting_module_available = attempt_import(
-    'bogus_nonexisting_module', alt_names=['bogus_nem'], defer_check=True
+    'bogus_nonexisting_module', alt_names=['bogus_nem'], defer_import=True
 )
diff --git a/pyomo/common/tests/deps.py b/pyomo/common/tests/deps.py
index d00281553f4..5f8c1fffdf8 100644
--- a/pyomo/common/tests/deps.py
+++ b/pyomo/common/tests/deps.py
@@ -23,15 +23,16 @@
     bogus_nonexisting_module_available as has_bogus_nem,
 )
 
-bogus, bogus_available = attempt_import('nonexisting.module.bogus', defer_check=True)
+bogus, bogus_available = attempt_import('nonexisting.module.bogus', defer_import=True)
 
 pkl_test, pkl_available = attempt_import(
-    'nonexisting.module.pickle_test', deferred_submodules=['submod'], defer_check=True
+    'nonexisting.module.pickle_test', deferred_submodules=['submod'], defer_import=True
 )
 
 pyo, pyo_available = attempt_import(
     'pyomo',
     alt_names=['pyo'],
+    defer_import=True,
     deferred_submodules={'version': None, 'common.tests.dep_mod': ['dm']},
 )
 
diff --git a/pyomo/common/tests/test_config.py b/pyomo/common/tests/test_config.py
index 12657481764..a47f5e0d8af 100644
--- a/pyomo/common/tests/test_config.py
+++ b/pyomo/common/tests/test_config.py
@@ -2098,7 +2098,6 @@ def test_generate_custom_documentation(self):
                 "generate_documentation is deprecated.",
                 LOG,
             )
-        self.maxDiff = None
         # print(test)
         self.assertEqual(test, reference)
 
@@ -2113,7 +2112,6 @@ def test_generate_custom_documentation(self):
                 )
             )
         self.assertEqual(LOG.getvalue(), "")
-        self.maxDiff = None
         # print(test)
         self.assertEqual(test, reference)
 
@@ -2159,7 +2157,6 @@ def test_generate_custom_documentation(self):
                 "generate_documentation is deprecated.",
                 LOG,
             )
-        self.maxDiff = None
         # print(test)
         self.assertEqual(test, reference)
 
@@ -2577,7 +2574,6 @@ def test_argparse_help_implicit_disable(self):
         parser = argparse.ArgumentParser(prog='tester')
         self.config.initialize_argparse(parser)
         help = parser.format_help()
-        self.maxDiff = None
         self.assertIn(
             """
   -h, --help            show this help message and exit
@@ -3106,8 +3102,6 @@ def test_declare_from(self):
             cfg2.declare_from({})
 
     def test_docstring_decorator(self):
-        self.maxDiff = None
-
         @document_kwargs_from_configdict('CONFIG')
         class ExampleClass(object):
             CONFIG = ExampleConfig()
diff --git a/pyomo/common/tests/test_dependencies.py b/pyomo/common/tests/test_dependencies.py
index 30822a4f81f..31f9520b613 100644
--- a/pyomo/common/tests/test_dependencies.py
+++ b/pyomo/common/tests/test_dependencies.py
@@ -45,7 +45,7 @@ def test_import_error(self):
         module_obj, module_available = attempt_import(
             '__there_is_no_module_named_this__',
             'Testing import of a non-existent module',
-            defer_check=False,
+            defer_import=False,
         )
         self.assertFalse(module_available)
         with self.assertRaisesRegex(
@@ -85,7 +85,7 @@ def test_pickle(self):
 
     def test_import_success(self):
         module_obj, module_available = attempt_import(
-            'ply', 'Testing import of ply', defer_check=False
+            'ply', 'Testing import of ply', defer_import=False
         )
         self.assertTrue(module_available)
         import ply
@@ -123,7 +123,7 @@ def test_imported_deferred_import(self):
 
     def test_min_version(self):
         mod, avail = attempt_import(
-            'pyomo.common.tests.dep_mod', minimum_version='1.0', defer_check=False
+            'pyomo.common.tests.dep_mod', minimum_version='1.0', defer_import=False
         )
         self.assertTrue(avail)
         self.assertTrue(inspect.ismodule(mod))
@@ -131,7 +131,7 @@ def test_min_version(self):
         self.assertFalse(check_min_version(mod, '2.0'))
 
         mod, avail = attempt_import(
-            'pyomo.common.tests.dep_mod', minimum_version='2.0', defer_check=False
+            'pyomo.common.tests.dep_mod', minimum_version='2.0', defer_import=False
         )
         self.assertFalse(avail)
         self.assertIs(type(mod), ModuleUnavailable)
@@ -146,7 +146,7 @@ def test_min_version(self):
             'pyomo.common.tests.dep_mod',
             error_message="Failed import",
             minimum_version='2.0',
-            defer_check=False,
+            defer_import=False,
         )
         self.assertFalse(avail)
         self.assertIs(type(mod), ModuleUnavailable)
@@ -159,10 +159,10 @@ def test_min_version(self):
 
         # Verify check_min_version works with deferred imports
 
-        mod, avail = attempt_import('pyomo.common.tests.dep_mod', defer_check=True)
+        mod, avail = attempt_import('pyomo.common.tests.dep_mod', defer_import=True)
         self.assertTrue(check_min_version(mod, '1.0'))
 
-        mod, avail = attempt_import('pyomo.common.tests.dep_mod', defer_check=True)
+        mod, avail = attempt_import('pyomo.common.tests.dep_mod', defer_import=True)
         self.assertFalse(check_min_version(mod, '2.0'))
 
         # Verify check_min_version works when called directly
@@ -174,10 +174,10 @@ def test_min_version(self):
         self.assertFalse(check_min_version(mod, '1.0'))
 
     def test_and_or(self):
-        mod0, avail0 = attempt_import('ply', defer_check=True)
-        mod1, avail1 = attempt_import('pyomo.common.tests.dep_mod', defer_check=True)
+        mod0, avail0 = attempt_import('ply', defer_import=True)
+        mod1, avail1 = attempt_import('pyomo.common.tests.dep_mod', defer_import=True)
         mod2, avail2 = attempt_import(
-            'pyomo.common.tests.dep_mod', minimum_version='2.0', defer_check=True
+            'pyomo.common.tests.dep_mod', minimum_version='2.0', defer_import=True
         )
 
         _and = avail0 & avail1
@@ -233,11 +233,11 @@ def test_callbacks(self):
         def _record_avail(module, avail):
             ans.append(avail)
 
-        mod0, avail0 = attempt_import('ply', defer_check=True, callback=_record_avail)
+        mod0, avail0 = attempt_import('ply', defer_import=True, callback=_record_avail)
         mod1, avail1 = attempt_import(
             'pyomo.common.tests.dep_mod',
             minimum_version='2.0',
-            defer_check=True,
+            defer_import=True,
             callback=_record_avail,
         )
 
@@ -250,7 +250,7 @@ def _record_avail(module, avail):
     def test_import_exceptions(self):
         mod, avail = attempt_import(
             'pyomo.common.tests.dep_mod_except',
-            defer_check=True,
+            defer_import=True,
             only_catch_importerror=True,
         )
         with self.assertRaisesRegex(ValueError, "cannot import module"):
@@ -260,7 +260,7 @@ def test_import_exceptions(self):
 
         mod, avail = attempt_import(
             'pyomo.common.tests.dep_mod_except',
-            defer_check=True,
+            defer_import=True,
             only_catch_importerror=False,
         )
         self.assertFalse(avail)
@@ -268,7 +268,7 @@ def test_import_exceptions(self):
 
         mod, avail = attempt_import(
             'pyomo.common.tests.dep_mod_except',
-            defer_check=True,
+            defer_import=True,
             catch_exceptions=(ImportError, ValueError),
         )
         self.assertFalse(avail)
@@ -280,7 +280,7 @@ def test_import_exceptions(self):
         ):
             mod, avail = attempt_import(
                 'pyomo.common.tests.dep_mod_except',
-                defer_check=True,
+                defer_import=True,
                 only_catch_importerror=True,
                 catch_exceptions=(ImportError,),
             )
@@ -288,7 +288,7 @@ def test_import_exceptions(self):
     def test_generate_warning(self):
         mod, avail = attempt_import(
             'pyomo.common.tests.dep_mod_except',
-            defer_check=True,
+            defer_import=True,
             only_catch_importerror=False,
         )
 
@@ -324,7 +324,7 @@ def test_generate_warning(self):
     def test_log_warning(self):
         mod, avail = attempt_import(
             'pyomo.common.tests.dep_mod_except',
-            defer_check=True,
+            defer_import=True,
             only_catch_importerror=False,
         )
         log = StringIO()
@@ -366,9 +366,9 @@ def test_importer(self):
 
         def _importer():
             attempted_import.append(True)
-            return attempt_import('pyomo.common.tests.dep_mod', defer_check=False)[0]
+            return attempt_import('pyomo.common.tests.dep_mod', defer_import=False)[0]
 
-        mod, avail = attempt_import('foo', importer=_importer, defer_check=True)
+        mod, avail = attempt_import('foo', importer=_importer, defer_import=True)
 
         self.assertEqual(attempted_import, [])
         self.assertIsInstance(mod, DeferredImportModule)
@@ -401,17 +401,17 @@ def test_deferred_submodules(self):
         self.assertTrue(inspect.ismodule(deps.dm))
 
         with self.assertRaisesRegex(
-            ValueError, "deferred_submodules is only valid if defer_check==True"
+            ValueError, "deferred_submodules is only valid if defer_import==True"
         ):
             mod, mod_available = attempt_import(
                 'nonexisting.module',
-                defer_check=False,
+                defer_import=False,
                 deferred_submodules={'submod': None},
             )
 
         mod, mod_available = attempt_import(
             'nonexisting.module',
-            defer_check=True,
+            defer_import=True,
             deferred_submodules={'submod.subsubmod': None},
         )
         self.assertIs(type(mod), DeferredImportModule)
@@ -427,7 +427,7 @@ def test_UnavailableClass(self):
         module_obj, module_available = attempt_import(
             '__there_is_no_module_named_this__',
             'Testing import of a non-existent module',
-            defer_check=False,
+            defer_import=False,
         )
 
         class A_Class(UnavailableClass(module_obj)):
diff --git a/pyomo/common/tests/test_log.py b/pyomo/common/tests/test_log.py
index 64691c0015a..166e1e44cdb 100644
--- a/pyomo/common/tests/test_log.py
+++ b/pyomo/common/tests/test_log.py
@@ -511,7 +511,6 @@ def test_verbatim(self):
             "\n"
             "    quote block\n"
         )
-        self.maxDiff = None
         self.assertEqual(self.stream.getvalue(), ans)
 
 
diff --git a/pyomo/common/tests/test_timing.py b/pyomo/common/tests/test_timing.py
index 0a4224c5476..48288746882 100644
--- a/pyomo/common/tests/test_timing.py
+++ b/pyomo/common/tests/test_timing.py
@@ -107,7 +107,6 @@ def test_report_timing(self):
                 m.y = Var(Any, dense=False)
                 xfrm.apply_to(m)
             result = out.getvalue().strip()
-            self.maxDiff = None
             for l, r in zip(result.splitlines(), ref.splitlines()):
                 self.assertRegex(str(l.strip()), str(r.strip()))
         finally:
@@ -122,7 +121,6 @@ def test_report_timing(self):
             m.y = Var(Any, dense=False)
             xfrm.apply_to(m)
             result = os.getvalue().strip()
-            self.maxDiff = None
             for l, r in zip(result.splitlines(), ref.splitlines()):
                 self.assertRegex(str(l.strip()), str(r.strip()))
         finally:
@@ -135,7 +133,6 @@ def test_report_timing(self):
             m.y = Var(Any, dense=False)
             xfrm.apply_to(m)
             result = os.getvalue().strip()
-            self.maxDiff = None
             for l, r in zip(result.splitlines(), ref.splitlines()):
                 self.assertRegex(str(l.strip()), str(r.strip()))
             self.assertEqual(buf.getvalue().strip(), "")
@@ -172,7 +169,6 @@ def test_report_timing_context_manager(self):
                     xfrm.apply_to(m)
                 self.assertEqual(OUT.getvalue(), "")
                 result = OS.getvalue().strip()
-                self.maxDiff = None
                 for l, r in zip_longest(result.splitlines(), ref.splitlines()):
                     self.assertRegex(str(l.strip()), str(r.strip()))
             # Active reporting is False: the previous log should not have changed
diff --git a/pyomo/common/unittest.py b/pyomo/common/unittest.py
index 9a21b35faa8..84d962eb784 100644
--- a/pyomo/common/unittest.py
+++ b/pyomo/common/unittest.py
@@ -498,6 +498,10 @@ class TestCase(_unittest.TestCase):
 
     __doc__ += _unittest.TestCase.__doc__
 
+    # By default, we always want to spend the time to create the full
+    # diff of the test reault and the baseline
+    maxDiff = None
+
     def assertStructuredAlmostEqual(
         self,
         first,
@@ -631,7 +635,7 @@ def initialize_dependencies(self):
         cls.package_modules = {}
         packages_used = set(sum(list(cls.package_dependencies.values()), []))
         for package_ in packages_used:
-            pack, pack_avail = attempt_import(package_, defer_check=False)
+            pack, pack_avail = attempt_import(package_, defer_import=False)
             cls.package_available[package_] = pack_avail
             cls.package_modules[package_] = pack
 
diff --git a/pyomo/contrib/pynumero/algorithms/solvers/cyipopt_solver.py b/pyomo/contrib/pynumero/algorithms/solvers/cyipopt_solver.py
index 9d24c0dd562..cdea542295b 100644
--- a/pyomo/contrib/pynumero/algorithms/solvers/cyipopt_solver.py
+++ b/pyomo/contrib/pynumero/algorithms/solvers/cyipopt_solver.py
@@ -24,6 +24,8 @@
 from pyomo.common.deprecation import relocated_module_attribute
 from pyomo.common.dependencies import attempt_import, numpy as np, numpy_available
 from pyomo.common.tee import redirect_fd, TeeStream
+from pyomo.common.modeling import unique_component_name
+from pyomo.core.base.objective import Objective
 
 # Because pynumero.interfaces requires numpy, we will leverage deferred
 # imports here so that the solver can be registered even when numpy is
@@ -332,11 +334,22 @@ def solve(self, model, **kwds):
         grey_box_blocks = list(
             model.component_data_objects(egb.ExternalGreyBoxBlock, active=True)
         )
-        if grey_box_blocks:
-            # nlp = pyomo_nlp.PyomoGreyBoxNLP(model)
-            nlp = pyomo_grey_box.PyomoNLPWithGreyBoxBlocks(model)
-        else:
-            nlp = pyomo_nlp.PyomoNLP(model)
+        # if there is no objective, add one temporarily so we can construct an NLP
+        objectives = list(model.component_data_objects(Objective, active=True))
+        if not objectives:
+            objname = unique_component_name(model, "_obj")
+            objective = model.add_component(objname, Objective(expr=0.0))
+        try:
+            if grey_box_blocks:
+                # nlp = pyomo_nlp.PyomoGreyBoxNLP(model)
+                nlp = pyomo_grey_box.PyomoNLPWithGreyBoxBlocks(model)
+            else:
+                nlp = pyomo_nlp.PyomoNLP(model)
+        finally:
+            # We only need the objective to construct the NLP, so we delete
+            # it from the model ASAP
+            if not objectives:
+                model.del_component(objective)
 
         problem = cyipopt_interface.CyIpoptNLP(
             nlp,
diff --git a/pyomo/contrib/pynumero/algorithms/solvers/tests/test_cyipopt_solver.py b/pyomo/contrib/pynumero/algorithms/solvers/tests/test_cyipopt_solver.py
index e9da31097a0..0af5a772c98 100644
--- a/pyomo/contrib/pynumero/algorithms/solvers/tests/test_cyipopt_solver.py
+++ b/pyomo/contrib/pynumero/algorithms/solvers/tests/test_cyipopt_solver.py
@@ -316,3 +316,13 @@ def test_hs071_evalerror_old_cyipopt(self):
         msg = "Error in AMPL evaluation"
         with self.assertRaisesRegex(PyNumeroEvaluationError, msg):
             res = solver.solve(m, tee=True)
+
+    def test_solve_without_objective(self):
+        m = create_model1()
+        m.o.deactivate()
+        m.x[2].fix(0.0)
+        m.x[3].fix(4.0)
+        solver = pyo.SolverFactory("cyipopt")
+        res = solver.solve(m, tee=True)
+        pyo.assert_optimal_termination(res)
+        self.assertAlmostEqual(m.x[1].value, 9.0)
diff --git a/pyomo/contrib/pynumero/dependencies.py b/pyomo/contrib/pynumero/dependencies.py
index 9e2088ffa0a..d323bd43e84 100644
--- a/pyomo/contrib/pynumero/dependencies.py
+++ b/pyomo/contrib/pynumero/dependencies.py
@@ -17,7 +17,7 @@
     'numpy',
     'Pynumero requires the optional Pyomo dependency "numpy"',
     minimum_version='1.13.0',
-    defer_check=False,
+    defer_import=False,
 )
 
 if not numpy_available:
diff --git a/pyomo/contrib/pynumero/examples/tests/test_cyipopt_examples.py b/pyomo/contrib/pynumero/examples/tests/test_cyipopt_examples.py
index 1f45f26d43b..408a0197382 100644
--- a/pyomo/contrib/pynumero/examples/tests/test_cyipopt_examples.py
+++ b/pyomo/contrib/pynumero/examples/tests/test_cyipopt_examples.py
@@ -35,7 +35,7 @@
     'One of the tests below requires a recent version of pandas for'
     ' comparing with a tolerance.',
     minimum_version='1.1.0',
-    defer_check=False,
+    defer_import=False,
 )
 
 from pyomo.contrib.pynumero.asl import AmplInterface
diff --git a/pyomo/contrib/pynumero/intrinsic.py b/pyomo/contrib/pynumero/intrinsic.py
index 84675cc4c02..34054e7ffa2 100644
--- a/pyomo/contrib/pynumero/intrinsic.py
+++ b/pyomo/contrib/pynumero/intrinsic.py
@@ -11,9 +11,7 @@
 
 from pyomo.common.dependencies import numpy as np, attempt_import
 
-block_vector = attempt_import(
-    'pyomo.contrib.pynumero.sparse.block_vector', defer_check=True
-)[0]
+block_vector = attempt_import('pyomo.contrib.pynumero.sparse.block_vector')[0]
 
 
 def norm(x, ord=None):
diff --git a/pyomo/contrib/pyros/config.py b/pyomo/contrib/pyros/config.py
index a7ca41d095f..bc2bfd591e6 100644
--- a/pyomo/contrib/pyros/config.py
+++ b/pyomo/contrib/pyros/config.py
@@ -26,71 +26,34 @@
 default_pyros_solver_logger = setup_pyros_logger()
 
 
-class LoggerType:
+def logger_domain(obj):
     """
-    Domain validator for objects castable to logging.Logger.
-    """
-
-    def __call__(self, obj):
-        """
-        Cast object to logger.
+    Domain validator for logger-type arguments.
 
-        Parameters
-        ----------
-        obj : object
-            Object to be cast.
+    This admits any object of type ``logging.Logger``,
+    or which can be cast to ``logging.Logger``.
+    """
+    if isinstance(obj, logging.Logger):
+        return obj
+    else:
+        return logging.getLogger(obj)
 
-        Returns
-        -------
-        logging.Logger
-            If `str_or_logger` is of type `logging.Logger`,then
-            `str_or_logger` is returned.
-            Otherwise, ``logging.getLogger(str_or_logger)``
-            is returned.
-        """
-        if isinstance(obj, logging.Logger):
-            return obj
-        else:
-            return logging.getLogger(obj)
 
-    def domain_name(self):
-        """Return str briefly describing domain encompassed by self."""
-        return "None, str or logging.Logger"
+logger_domain.domain_name = "None, str or logging.Logger"
 
 
-class PositiveIntOrMinusOne:
+def positive_int_or_minus_one(obj):
     """
-    Domain validator for objects castable to a
-    strictly positive int or -1.
+    Domain validator for objects castable to a strictly
+    positive int or -1.
     """
+    ans = int(obj)
+    if ans != float(obj) or (ans <= 0 and ans != -1):
+        raise ValueError(f"Expected positive int or -1, but received value {obj!r}")
+    return ans
 
-    def __call__(self, obj):
-        """
-        Cast object to positive int or -1.
 
-        Parameters
-        ----------
-        obj : object
-            Object of interest.
-
-        Returns
-        -------
-        int
-            Positive int, or -1.
-
-        Raises
-        ------
-        ValueError
-            If object not castable to positive int, or -1.
-        """
-        ans = int(obj)
-        if ans != float(obj) or (ans <= 0 and ans != -1):
-            raise ValueError(f"Expected positive int or -1, but received value {obj!r}")
-        return ans
-
-    def domain_name(self):
-        """Return str briefly describing domain encompassed by self."""
-        return "positive int or -1"
+positive_int_or_minus_one.domain_name = "positive int or -1"
 
 
 def mutable_param_validator(param_obj):
@@ -721,7 +684,7 @@ def pyros_config():
         "max_iter",
         ConfigValue(
             default=-1,
-            domain=PositiveIntOrMinusOne(),
+            domain=positive_int_or_minus_one,
             description=(
                 """
                 Iteration limit. If -1 is provided, then no iteration
@@ -766,7 +729,7 @@ def pyros_config():
         "progress_logger",
         ConfigValue(
             default=default_pyros_solver_logger,
-            domain=LoggerType(),
+            domain=logger_domain,
             doc=(
                 """
                 Logger (or name thereof) used for reporting PyROS solver
diff --git a/pyomo/contrib/pyros/tests/test_config.py b/pyomo/contrib/pyros/tests/test_config.py
index 76b9114b9e6..0f52d04135d 100644
--- a/pyomo/contrib/pyros/tests/test_config.py
+++ b/pyomo/contrib/pyros/tests/test_config.py
@@ -12,9 +12,9 @@
 from pyomo.contrib.pyros.config import (
     InputDataStandardizer,
     mutable_param_validator,
-    LoggerType,
+    logger_domain,
     SolverNotResolvable,
-    PositiveIntOrMinusOne,
+    positive_int_or_minus_one,
     pyros_config,
     SolverIterable,
     SolverResolvable,
@@ -557,16 +557,29 @@ def test_positive_int_or_minus_one(self):
         """
         Test positive int or -1 validator works as expected.
         """
-        standardizer_func = PositiveIntOrMinusOne()
-        self.assertIs(
-            standardizer_func(1.0),
+        standardizer_func = positive_int_or_minus_one
+        ans = standardizer_func(1.0)
+        self.assertEqual(
+            ans,
             1,
-            msg=(f"{PositiveIntOrMinusOne.__name__} does not standardize as expected."),
+            msg=f"{positive_int_or_minus_one.__name__} output value not as expected.",
+        )
+        self.assertIs(
+            type(ans),
+            int,
+            msg=f"{positive_int_or_minus_one.__name__} output type not as expected.",
         )
+
+        ans = standardizer_func(-1.0)
         self.assertEqual(
-            standardizer_func(-1.00),
+            ans,
             -1,
-            msg=(f"{PositiveIntOrMinusOne.__name__} does not standardize as expected."),
+            msg=f"{positive_int_or_minus_one.__name__} output value not as expected.",
+        )
+        self.assertIs(
+            type(ans),
+            int,
+            msg=f"{positive_int_or_minus_one.__name__} output type not as expected.",
         )
 
         exc_str = r"Expected positive int or -1, but received value.*"
@@ -576,26 +589,26 @@ def test_positive_int_or_minus_one(self):
             standardizer_func(0)
 
 
-class TestLoggerType(unittest.TestCase):
+class TestLoggerDomain(unittest.TestCase):
     """
-    Test logger type validator.
+    Test logger type domain validator.
     """
 
     def test_logger_type(self):
         """
         Test logger type validator.
         """
-        standardizer_func = LoggerType()
+        standardizer_func = logger_domain
         mylogger = logging.getLogger("example")
         self.assertIs(
             standardizer_func(mylogger),
             mylogger,
-            msg=f"{LoggerType.__name__} output not as expected",
+            msg=f"{standardizer_func.__name__} output not as expected",
         )
         self.assertIs(
             standardizer_func(mylogger.name),
             mylogger,
-            msg=f"{LoggerType.__name__} output not as expected",
+            msg=f"{standardizer_func.__name__} output not as expected",
         )
 
         exc_str = r"A logger name must be a string"
diff --git a/pyomo/contrib/solver/ipopt.py b/pyomo/contrib/solver/ipopt.py
index 3ac1a5ac4a2..dc632adb184 100644
--- a/pyomo/contrib/solver/ipopt.py
+++ b/pyomo/contrib/solver/ipopt.py
@@ -17,7 +17,11 @@
 
 from pyomo.common import Executable
 from pyomo.common.config import ConfigValue, document_kwargs_from_configdict, ConfigDict
-from pyomo.common.errors import PyomoException, DeveloperError
+from pyomo.common.errors import (
+    PyomoException,
+    DeveloperError,
+    InfeasibleConstraintException,
+)
 from pyomo.common.tempfiles import TempfileManager
 from pyomo.common.timing import HierarchicalTimer
 from pyomo.core.base.var import _GeneralVarData
@@ -72,11 +76,7 @@ def __init__(
             ),
         )
         self.writer_config: ConfigDict = self.declare(
-            'writer_config',
-            ConfigValue(
-                default=NLWriter.CONFIG(),
-                description="Configuration that controls options in the NL writer.",
-            ),
+            'writer_config', NLWriter.CONFIG()
         )
 
 
@@ -314,15 +314,19 @@ def solve(self, model, **kwds):
             ) as row_file, open(basename + '.col', 'w') as col_file:
                 timer.start('write_nl_file')
                 self._writer.config.set_value(config.writer_config)
-                nl_info = self._writer.write(
-                    model,
-                    nl_file,
-                    row_file,
-                    col_file,
-                    symbolic_solver_labels=config.symbolic_solver_labels,
-                )
+                try:
+                    nl_info = self._writer.write(
+                        model,
+                        nl_file,
+                        row_file,
+                        col_file,
+                        symbolic_solver_labels=config.symbolic_solver_labels,
+                    )
+                    proven_infeasible = False
+                except InfeasibleConstraintException:
+                    proven_infeasible = True
                 timer.stop('write_nl_file')
-            if len(nl_info.variables) > 0:
+            if not proven_infeasible and len(nl_info.variables) > 0:
                 # Get a copy of the environment to pass to the subprocess
                 env = os.environ.copy()
                 if nl_info.external_function_libraries:
@@ -361,11 +365,17 @@ def solve(self, model, **kwds):
                     timer.stop('subprocess')
                     # This is the stuff we need to parse to get the iterations
                     # and time
-                    iters, ipopt_time_nofunc, ipopt_time_func, ipopt_total_time = (
+                    (iters, ipopt_time_nofunc, ipopt_time_func, ipopt_total_time) = (
                         self._parse_ipopt_output(ostreams[0])
                     )
 
-            if len(nl_info.variables) == 0:
+            if proven_infeasible:
+                results = Results()
+                results.termination_condition = TerminationCondition.provenInfeasible
+                results.solution_loader = SolSolutionLoader(None, None)
+                results.iteration_count = 0
+                results.timing_info.total_seconds = 0
+            elif len(nl_info.variables) == 0:
                 if len(nl_info.eliminated_vars) == 0:
                     results = Results()
                     results.termination_condition = TerminationCondition.emptyModel
@@ -457,7 +467,7 @@ def solve(self, model, **kwds):
                 )
 
         results.solver_configuration = config
-        if len(nl_info.variables) > 0:
+        if not proven_infeasible and len(nl_info.variables) > 0:
             results.solver_log = ostreams[0].getvalue()
 
         # Capture/record end-time / wall-time
diff --git a/pyomo/contrib/solver/tests/solvers/test_solvers.py b/pyomo/contrib/solver/tests/solvers/test_solvers.py
index cf5f6cf5c57..a4f4a3bc389 100644
--- a/pyomo/contrib/solver/tests/solvers/test_solvers.py
+++ b/pyomo/contrib/solver/tests/solvers/test_solvers.py
@@ -1508,6 +1508,71 @@ def test_bug_2(self, name: str, opt_class: Type[SolverBase], use_presolve: bool)
             res = opt.solve(m)
             self.assertAlmostEqual(res.incumbent_objective, -18, 5)
 
+    @parameterized.expand(input=_load_tests(nl_solvers))
+    def test_presolve_with_zero_coef(
+        self, name: str, opt_class: Type[SolverBase], use_presolve: bool
+    ):
+        opt: SolverBase = opt_class()
+        if not opt.available():
+            raise unittest.SkipTest(f'Solver {opt.name} not available.')
+        if use_presolve:
+            opt.config.writer_config.linear_presolve = True
+        else:
+            opt.config.writer_config.linear_presolve = False
+
+        """
+        when c2 gets presolved out, c1 becomes 
+        x - y + y = 0 which becomes
+        x - 0*y == 0 which is the zero we are testing for
+        """
+        m = pe.ConcreteModel()
+        m.x = pe.Var()
+        m.y = pe.Var()
+        m.z = pe.Var()
+        m.obj = pe.Objective(expr=m.x**2 + m.y**2 + m.z**2)
+        m.c1 = pe.Constraint(expr=m.x == m.y + m.z + 1.5)
+        m.c2 = pe.Constraint(expr=m.z == -m.y)
+
+        res = opt.solve(m)
+        self.assertAlmostEqual(res.incumbent_objective, 2.25)
+        self.assertAlmostEqual(m.x.value, 1.5)
+        self.assertAlmostEqual(m.y.value, 0)
+        self.assertAlmostEqual(m.z.value, 0)
+
+        m.x.setlb(2)
+        res = opt.solve(
+            m, load_solutions=False, raise_exception_on_nonoptimal_result=False
+        )
+        if use_presolve:
+            exp = TerminationCondition.provenInfeasible
+        else:
+            exp = TerminationCondition.locallyInfeasible
+        self.assertEqual(res.termination_condition, exp)
+
+        m = pe.ConcreteModel()
+        m.w = pe.Var()
+        m.x = pe.Var()
+        m.y = pe.Var()
+        m.z = pe.Var()
+        m.obj = pe.Objective(expr=m.x**2 + m.y**2 + m.z**2 + m.w**2)
+        m.c1 = pe.Constraint(expr=m.x + m.w == m.y + m.z)
+        m.c2 = pe.Constraint(expr=m.z == -m.y)
+        m.c3 = pe.Constraint(expr=m.x == -m.w)
+
+        res = opt.solve(m)
+        self.assertAlmostEqual(res.incumbent_objective, 0)
+        self.assertAlmostEqual(m.w.value, 0)
+        self.assertAlmostEqual(m.x.value, 0)
+        self.assertAlmostEqual(m.y.value, 0)
+        self.assertAlmostEqual(m.z.value, 0)
+
+        del m.c1
+        m.c1 = pe.Constraint(expr=m.x + m.w == m.y + m.z + 1.5)
+        res = opt.solve(
+            m, load_solutions=False, raise_exception_on_nonoptimal_result=False
+        )
+        self.assertEqual(res.termination_condition, exp)
+
     @parameterized.expand(input=_load_tests(all_solvers))
     def test_scaling(self, name: str, opt_class: Type[SolverBase], use_presolve: bool):
         opt: SolverBase = opt_class()
diff --git a/pyomo/core/base/units_container.py b/pyomo/core/base/units_container.py
index 1bf25ffdead..fb3d28385d0 100644
--- a/pyomo/core/base/units_container.py
+++ b/pyomo/core/base/units_container.py
@@ -127,7 +127,6 @@
 
 pint_module, pint_available = attempt_import(
     'pint',
-    defer_check=True,
     error_message=(
         'The "pint" package failed to import. '
         'This package is necessary to use Pyomo units.'
diff --git a/pyomo/core/tests/unit/test_block.py b/pyomo/core/tests/unit/test_block.py
index 88646643703..71e80d90a73 100644
--- a/pyomo/core/tests/unit/test_block.py
+++ b/pyomo/core/tests/unit/test_block.py
@@ -2667,7 +2667,6 @@ def test_pprint(self):
 
 5 Declarations: a1_IDX a3_IDX c a b
 """
-        self.maxDiff = None
         self.assertEqual(ref, buf.getvalue())
 
     @unittest.skipIf(not 'glpk' in solvers, "glpk solver is not available")
diff --git a/pyomo/core/tests/unit/test_numeric_expr.py b/pyomo/core/tests/unit/test_numeric_expr.py
index c073ee0f726..c1066c292d7 100644
--- a/pyomo/core/tests/unit/test_numeric_expr.py
+++ b/pyomo/core/tests/unit/test_numeric_expr.py
@@ -1424,7 +1424,6 @@ def test_sumOf_nestedTrivialProduct2(self):
         e1 = m.a * m.p
         e2 = m.b - m.c
         e = e2 - e1
-        self.maxDiff = None
         self.assertExpressionsEqual(
             e,
             LinearExpression(
diff --git a/pyomo/core/tests/unit/test_numvalue.py b/pyomo/core/tests/unit/test_numvalue.py
index bd784d655e8..1cccd3863ea 100644
--- a/pyomo/core/tests/unit/test_numvalue.py
+++ b/pyomo/core/tests/unit/test_numvalue.py
@@ -18,6 +18,7 @@
 
 import pyomo.common.unittest as unittest
 from pyomo.common.dependencies import numpy, numpy_available
+from pyomo.core.base.units_container import pint_available
 
 from pyomo.environ import (
     value,
@@ -50,7 +51,16 @@ def __init__(self, val=0):
 
 class MyBogusNumericType(MyBogusType):
     def __add__(self, other):
-        return MyBogusNumericType(self.val + float(other))
+        if other.__class__ in native_numeric_types:
+            return MyBogusNumericType(self.val + float(other))
+        else:
+            return NotImplemented
+
+    def __le__(self, other):
+        if other.__class__ in native_numeric_types:
+            return self.val <= float(other)
+        else:
+            return NotImplemented
 
     def __lt__(self, other):
         return self.val < float(other)
@@ -534,6 +544,8 @@ def test_unknownNumericType(self):
         try:
             val = as_numeric(ref)
             self.assertEqual(val().val, 42.0)
+            self.assertIn(MyBogusNumericType, native_numeric_types)
+            self.assertIn(MyBogusNumericType, native_types)
         finally:
             native_numeric_types.remove(MyBogusNumericType)
             native_types.remove(MyBogusNumericType)
@@ -562,10 +574,43 @@ def test_numpy_basic_bool_registration(self):
     @unittest.skipUnless(numpy_available, "This test requires NumPy")
     def test_automatic_numpy_registration(self):
         cmd = (
-            'import pyomo; from pyomo.core.base import Var, Param; '
-            'from pyomo.core.base.units_container import units; import numpy as np; '
-            'print(np.float64 in pyomo.common.numeric_types.native_numeric_types); '
-            '%s; print(np.float64 in pyomo.common.numeric_types.native_numeric_types)'
+            'from pyomo.common.numeric_types import native_numeric_types as nnt; '
+            'print("float64" in [_.__name__ for _ in nnt]); '
+            'import numpy; '
+            'print("float64" in [_.__name__ for _ in nnt])'
+        )
+
+        rc = subprocess.run(
+            [sys.executable, '-c', cmd],
+            stdout=subprocess.PIPE,
+            stderr=subprocess.STDOUT,
+            text=True,
+        )
+        self.assertEqual((rc.returncode, rc.stdout), (0, "False\nTrue\n"))
+
+        cmd = (
+            'import numpy; '
+            'from pyomo.common.numeric_types import native_numeric_types as nnt; '
+            'print("float64" in [_.__name__ for _ in nnt])'
+        )
+
+        rc = subprocess.run(
+            [sys.executable, '-c', cmd],
+            stdout=subprocess.PIPE,
+            stderr=subprocess.STDOUT,
+            text=True,
+        )
+        self.assertEqual((rc.returncode, rc.stdout), (0, "True\n"))
+
+    def test_unknownNumericType_expr_registration(self):
+        cmd = (
+            'import pyomo; '
+            'from pyomo.core.base import Var, Param; '
+            'from pyomo.core.base.units_container import units; '
+            'from pyomo.common.numeric_types import native_numeric_types as nnt; '
+            f'from {__name__} import MyBogusNumericType; '
+            'ref = MyBogusNumericType(42); '
+            'print(MyBogusNumericType in nnt); %s; print(MyBogusNumericType in nnt); '
         )
 
         def _tester(expr):
@@ -575,19 +620,32 @@ def _tester(expr):
                 stderr=subprocess.STDOUT,
                 text=True,
             )
-            self.assertEqual((rc.returncode, rc.stdout), (0, "False\nTrue\n"))
-
-        _tester('Var() <= np.float64(5)')
-        _tester('np.float64(5) <= Var()')
-        _tester('np.float64(5) + Var()')
-        _tester('Var() + np.float64(5)')
-        _tester('v = Var(); v.construct(); v.value = np.float64(5)')
-        _tester('p = Param(mutable=True); p.construct(); p.value = np.float64(5)')
-        _tester('v = Var(units=units.m); v.construct(); v.value = np.float64(5)')
-        _tester(
-            'p = Param(mutable=True, units=units.m); p.construct(); '
-            'p.value = np.float64(5)'
-        )
+            self.assertEqual(
+                (rc.returncode, rc.stdout),
+                (
+                    0,
+                    '''False
+WARNING: Dynamically registering the following numeric type:
+        pyomo.core.tests.unit.test_numvalue.MyBogusNumericType
+    Dynamic registration is supported for convenience, but there are known
+    limitations to this approach.  We recommend explicitly registering numeric
+    types using RegisterNumericType() or RegisterIntegerType().
+True
+''',
+                ),
+            )
+
+        _tester('Var() <= ref')
+        _tester('ref <= Var()')
+        _tester('ref + Var()')
+        _tester('Var() + ref')
+        _tester('v = Var(); v.construct(); v.value = ref')
+        _tester('p = Param(mutable=True); p.construct(); p.value = ref')
+        if pint_available:
+            _tester('v = Var(units=units.m); v.construct(); v.value = ref')
+            _tester(
+                'p = Param(mutable=True, units=units.m); p.construct(); p.value = ref'
+            )
 
 
 if __name__ == "__main__":
diff --git a/pyomo/core/tests/unit/test_reference.py b/pyomo/core/tests/unit/test_reference.py
index 287ff204f9e..cfd9b99f945 100644
--- a/pyomo/core/tests/unit/test_reference.py
+++ b/pyomo/core/tests/unit/test_reference.py
@@ -1280,7 +1280,6 @@ def test_contains_with_nonflattened(self):
             normalize_index.flatten = _old_flatten
 
     def test_pprint_nonfinite_sets(self):
-        self.maxDiff = None
         m = ConcreteModel()
         m.v = Var(NonNegativeIntegers, dense=False)
         m.ref = Reference(m.v)
@@ -1322,7 +1321,6 @@ def test_pprint_nonfinite_sets(self):
 
     def test_pprint_nonfinite_sets_ctypeNone(self):
         # test issue #2039
-        self.maxDiff = None
         m = ConcreteModel()
         m.v = Var(NonNegativeIntegers, dense=False)
         m.ref = Reference(m.v, ctype=None)
diff --git a/pyomo/core/tests/unit/test_set.py b/pyomo/core/tests/unit/test_set.py
index 1ad08ba025c..4bbac6ecaa0 100644
--- a/pyomo/core/tests/unit/test_set.py
+++ b/pyomo/core/tests/unit/test_set.py
@@ -6267,7 +6267,6 @@ def test_issue_835(self):
 
     @unittest.skipIf(NamedTuple is None, "typing module not available")
     def test_issue_938(self):
-        self.maxDiff = None
         NodeKey = NamedTuple('NodeKey', [('id', int)])
         ArcKey = NamedTuple('ArcKey', [('node_from', NodeKey), ('node_to', NodeKey)])
 
diff --git a/pyomo/environ/tests/test_environ.py b/pyomo/environ/tests/test_environ.py
index 9c89fd135d5..9811b412af7 100644
--- a/pyomo/environ/tests/test_environ.py
+++ b/pyomo/environ/tests/test_environ.py
@@ -140,6 +140,7 @@ def test_tpl_import_time(self):
             'cPickle',
             'csv',
             'ctypes',  # mandatory import in core/base/external.py; TODO: fix this
+            'datetime',  # imported by contrib.solver
             'decimal',
             'gc',  # Imported on MacOS, Windows; Linux in 3.10
             'glob',
diff --git a/pyomo/gdp/plugins/bigm.py b/pyomo/gdp/plugins/bigm.py
index 1f9f561b192..3f450dbbd4f 100644
--- a/pyomo/gdp/plugins/bigm.py
+++ b/pyomo/gdp/plugins/bigm.py
@@ -13,6 +13,7 @@
 
 import logging
 
+from pyomo.common.autoslots import AutoSlots
 from pyomo.common.collections import ComponentMap
 from pyomo.common.config import ConfigDict, ConfigValue
 from pyomo.common.gc_manager import PauseGC
@@ -58,6 +59,26 @@
 logger = logging.getLogger('pyomo.gdp.bigm')
 
 
+class _BigMData(AutoSlots.Mixin):
+    __slots__ = ('bigm_src',)
+
+    def __init__(self):
+        # we will keep a map of constraints (hashable, ha!) to a tuple to
+        # indicate what their M value is and where it came from, of the form:
+        # ((lower_value, lower_source, lower_key), (upper_value, upper_source,
+        # upper_key)), where the first tuple is the information for the lower M,
+        # the second tuple is the info for the upper M, source is the Suffix or
+        # argument dictionary and None if the value was calculated, and key is
+        # the key in the Suffix or argument dictionary, and None if it was
+        # calculated. (Note that it is possible the lower or upper is
+        # user-specified and the other is not, hence the need to store
+        # information for both.)
+        self.bigm_src = {}
+
+
+Block.register_private_data_initializer(_BigMData)
+
+
 @TransformationFactory.register(
     'gdp.bigm', doc="Relax disjunctive model using big-M terms."
 )
@@ -94,15 +115,8 @@ class BigM_Transformation(GDP_to_MIP_Transformation, _BigM_MixIn):
     name beginning "_pyomo_gdp_bigm_reformulation".  That Block will
     contain an indexed Block named "relaxedDisjuncts", which will hold
     the relaxed disjuncts.  This block is indexed by an integer
-    indicating the order in which the disjuncts were relaxed.
-    Each block has a dictionary "_constraintMap":
-
-        'srcConstraints': ComponentMap(<transformed constraint>:
-                                       <src constraint>)
-        'transformedConstraints': ComponentMap(<src constraint>:
-                                               <transformed constraint>)
-
-    All transformed Disjuncts will have a pointer to the block their transformed
+    indicating the order in which the disjuncts were relaxed. All
+    transformed Disjuncts will have a pointer to the block their transformed
     constraints are on, and all transformed Disjunctions will have a
     pointer to the corresponding 'Or' or 'ExactlyOne' constraint.
 
@@ -247,18 +261,6 @@ def _transform_disjunct(self, obj, bigM, transBlock):
 
         relaxationBlock = self._get_disjunct_transformation_block(obj, transBlock)
 
-        # we will keep a map of constraints (hashable, ha!) to a tuple to
-        # indicate what their M value is and where it came from, of the form:
-        # ((lower_value, lower_source, lower_key), (upper_value, upper_source,
-        # upper_key)), where the first tuple is the information for the lower M,
-        # the second tuple is the info for the upper M, source is the Suffix or
-        # argument dictionary and None if the value was calculated, and key is
-        # the key in the Suffix or argument dictionary, and None if it was
-        # calculated. (Note that it is possible the lower or upper is
-        # user-specified and the other is not, hence the need to store
-        # information for both.)
-        relaxationBlock.bigm_src = {}
-
         # This is crazy, but if the disjunction has been previously
         # relaxed, the disjunct *could* be deactivated.  This is a big
         # deal for Hull, as it uses the component_objects /
@@ -278,8 +280,8 @@ def _transform_constraint(
     ):
         # add constraint to the transformation block, we'll transform it there.
         transBlock = disjunct._transformation_block()
-        bigm_src = transBlock.bigm_src
-        constraintMap = transBlock._constraintMap
+        bigm_src = transBlock.private_data().bigm_src
+        constraint_map = transBlock.private_data('pyomo.gdp')
 
         disjunctionRelaxationBlock = transBlock.parent_block()
 
@@ -346,7 +348,7 @@ def _transform_constraint(
             bigm_src[c] = (lower, upper)
 
             self._add_constraint_expressions(
-                c, i, M, disjunct.binary_indicator_var, newConstraint, constraintMap
+                c, i, M, disjunct.binary_indicator_var, newConstraint, constraint_map
             )
 
             # deactivate because we relaxed
@@ -409,7 +411,7 @@ def _update_M_from_suffixes(self, constraint, suffix_list, lower, upper):
     def get_m_value_src(self, constraint):
         transBlock = _get_constraint_transBlock(constraint)
         ((lower_val, lower_source, lower_key), (upper_val, upper_source, upper_key)) = (
-            transBlock.bigm_src[constraint]
+            transBlock.private_data().bigm_src[constraint]
         )
 
         if (
@@ -464,7 +466,7 @@ def get_M_value_src(self, constraint):
         transBlock = _get_constraint_transBlock(constraint)
         # This is a KeyError if it fails, but it is also my fault if it
         # fails... (That is, it's a bug in the mapping.)
-        return transBlock.bigm_src[constraint]
+        return transBlock.private_data().bigm_src[constraint]
 
     def get_M_value(self, constraint):
         """Returns the M values used to transform constraint. Return is a tuple:
@@ -479,7 +481,7 @@ def get_M_value(self, constraint):
         transBlock = _get_constraint_transBlock(constraint)
         # This is a KeyError if it fails, but it is also my fault if it
         # fails... (That is, it's a bug in the mapping.)
-        lower, upper = transBlock.bigm_src[constraint]
+        lower, upper = transBlock.private_data().bigm_src[constraint]
         return (lower[0], upper[0])
 
     def get_all_M_values_by_constraint(self, model):
@@ -499,9 +501,8 @@ def get_all_M_values_by_constraint(self, model):
             # First check if it was transformed at all.
             if transBlock is not None:
                 # If it was transformed with BigM, we get the M values.
-                if hasattr(transBlock, 'bigm_src'):
-                    for cons in transBlock.bigm_src:
-                        m_values[cons] = self.get_M_value(cons)
+                for cons in transBlock.private_data().bigm_src:
+                    m_values[cons] = self.get_M_value(cons)
         return m_values
 
     def get_largest_M_value(self, model):
diff --git a/pyomo/gdp/plugins/bigm_mixin.py b/pyomo/gdp/plugins/bigm_mixin.py
index ad6e6dcad86..510b36b5102 100644
--- a/pyomo/gdp/plugins/bigm_mixin.py
+++ b/pyomo/gdp/plugins/bigm_mixin.py
@@ -232,7 +232,7 @@ def _estimate_M(self, expr, constraint):
         return tuple(M)
 
     def _add_constraint_expressions(
-        self, c, i, M, indicator_var, newConstraint, constraintMap
+        self, c, i, M, indicator_var, newConstraint, constraint_map
     ):
         # Since we are both combining components from multiple blocks and using
         # local names, we need to make sure that the first index for
@@ -253,8 +253,10 @@ def _add_constraint_expressions(
                 )
             M_expr = M[0] * (1 - indicator_var)
             newConstraint.add((name, i, 'lb'), c.lower <= c.body - M_expr)
-            constraintMap['transformedConstraints'][c] = [newConstraint[name, i, 'lb']]
-            constraintMap['srcConstraints'][newConstraint[name, i, 'lb']] = c
+            constraint_map.transformed_constraints[c].append(
+                newConstraint[name, i, 'lb']
+            )
+            constraint_map.src_constraint[newConstraint[name, i, 'lb']] = c
         if c.upper is not None:
             if M[1] is None:
                 raise GDP_Error(
@@ -263,13 +265,7 @@ def _add_constraint_expressions(
                 )
             M_expr = M[1] * (1 - indicator_var)
             newConstraint.add((name, i, 'ub'), c.body - M_expr <= c.upper)
-            transformed = constraintMap['transformedConstraints'].get(c)
-            if transformed is not None:
-                constraintMap['transformedConstraints'][c].append(
-                    newConstraint[name, i, 'ub']
-                )
-            else:
-                constraintMap['transformedConstraints'][c] = [
-                    newConstraint[name, i, 'ub']
-                ]
-            constraintMap['srcConstraints'][newConstraint[name, i, 'ub']] = c
+            constraint_map.transformed_constraints[c].append(
+                newConstraint[name, i, 'ub']
+            )
+            constraint_map.src_constraint[newConstraint[name, i, 'ub']] = c
diff --git a/pyomo/gdp/plugins/binary_multiplication.py b/pyomo/gdp/plugins/binary_multiplication.py
index d68f7efe76f..bea33580ed6 100644
--- a/pyomo/gdp/plugins/binary_multiplication.py
+++ b/pyomo/gdp/plugins/binary_multiplication.py
@@ -121,7 +121,7 @@ def _transform_disjunct(self, obj, transBlock):
     def _transform_constraint(self, obj, disjunct):
         # add constraint to the transformation block, we'll transform it there.
         transBlock = disjunct._transformation_block()
-        constraintMap = transBlock._constraintMap
+        constraint_map = transBlock.private_data('pyomo.gdp')
 
         disjunctionRelaxationBlock = transBlock.parent_block()
 
@@ -137,14 +137,14 @@ def _transform_constraint(self, obj, disjunct):
                 continue
 
             self._add_constraint_expressions(
-                c, i, disjunct.binary_indicator_var, newConstraint, constraintMap
+                c, i, disjunct.binary_indicator_var, newConstraint, constraint_map
             )
 
             # deactivate because we relaxed
             c.deactivate()
 
     def _add_constraint_expressions(
-        self, c, i, indicator_var, newConstraint, constraintMap
+        self, c, i, indicator_var, newConstraint, constraint_map
     ):
         # Since we are both combining components from multiple blocks and using
         # local names, we need to make sure that the first index for
@@ -156,21 +156,21 @@ def _add_constraint_expressions(
         # over the constraint indices, but I don't think it matters a lot.)
         unique = len(newConstraint)
         name = c.local_name + "_%s" % unique
-        transformed = constraintMap['transformedConstraints'][c] = []
+        transformed = constraint_map.transformed_constraints[c]
 
         lb, ub = c.lower, c.upper
         if (c.equality or lb is ub) and lb is not None:
             # equality
             newConstraint.add((name, i, 'eq'), (c.body - lb) * indicator_var == 0)
             transformed.append(newConstraint[name, i, 'eq'])
-            constraintMap['srcConstraints'][newConstraint[name, i, 'eq']] = c
+            constraint_map.src_constraint[newConstraint[name, i, 'eq']] = c
         else:
             # inequality
             if lb is not None:
                 newConstraint.add((name, i, 'lb'), 0 <= (c.body - lb) * indicator_var)
                 transformed.append(newConstraint[name, i, 'lb'])
-                constraintMap['srcConstraints'][newConstraint[name, i, 'lb']] = c
+                constraint_map.src_constraint[newConstraint[name, i, 'lb']] = c
             if ub is not None:
                 newConstraint.add((name, i, 'ub'), (c.body - ub) * indicator_var <= 0)
                 transformed.append(newConstraint[name, i, 'ub'])
-                constraintMap['srcConstraints'][newConstraint[name, i, 'ub']] = c
+                constraint_map.src_constraint[newConstraint[name, i, 'ub']] = c
diff --git a/pyomo/gdp/plugins/gdp_to_mip_transformation.py b/pyomo/gdp/plugins/gdp_to_mip_transformation.py
index 59cb221321a..8dcd22b292a 100644
--- a/pyomo/gdp/plugins/gdp_to_mip_transformation.py
+++ b/pyomo/gdp/plugins/gdp_to_mip_transformation.py
@@ -11,7 +11,8 @@
 
 from functools import wraps
 
-from pyomo.common.collections import ComponentMap
+from pyomo.common.autoslots import AutoSlots
+from pyomo.common.collections import ComponentMap, DefaultComponentMap
 from pyomo.common.log import is_debug_set
 from pyomo.common.modeling import unique_component_name
 
@@ -48,6 +49,17 @@
 from weakref import ref as weakref_ref
 
 
+class _GDPTransformationData(AutoSlots.Mixin):
+    __slots__ = ('src_constraint', 'transformed_constraints')
+
+    def __init__(self):
+        self.src_constraint = ComponentMap()
+        self.transformed_constraints = DefaultComponentMap(list)
+
+
+Block.register_private_data_initializer(_GDPTransformationData, scope='pyomo.gdp')
+
+
 class GDP_to_MIP_Transformation(Transformation):
     """
     Base class for transformations from GDP to MIP
@@ -243,14 +255,7 @@ def _get_disjunct_transformation_block(self, disjunct, transBlock):
         relaxationBlock = relaxedDisjuncts[len(relaxedDisjuncts)]
 
         relaxationBlock.transformedConstraints = Constraint(Any)
-
         relaxationBlock.localVarReferences = Block()
-        # add the map that will link back and forth between transformed
-        # constraints and their originals.
-        relaxationBlock._constraintMap = {
-            'srcConstraints': ComponentMap(),
-            'transformedConstraints': ComponentMap(),
-        }
 
         # add mappings to source disjunct (so we'll know we've relaxed)
         disjunct._transformation_block = weakref_ref(relaxationBlock)
diff --git a/pyomo/gdp/plugins/hull.py b/pyomo/gdp/plugins/hull.py
index d1c38bde039..5b9d2ad08a9 100644
--- a/pyomo/gdp/plugins/hull.py
+++ b/pyomo/gdp/plugins/hull.py
@@ -58,12 +58,18 @@
 
 
 class _HullTransformationData(AutoSlots.Mixin):
-    __slots__ = ('disaggregated_var_map', 'original_var_map', 'bigm_constraint_map')
+    __slots__ = (
+        'disaggregated_var_map',
+        'original_var_map',
+        'bigm_constraint_map',
+        'disaggregation_constraint_map',
+    )
 
     def __init__(self):
         self.disaggregated_var_map = DefaultComponentMap(ComponentMap)
         self.original_var_map = ComponentMap()
         self.bigm_constraint_map = DefaultComponentMap(ComponentMap)
+        self.disaggregation_constraint_map = DefaultComponentMap(ComponentMap)
 
 
 Block.register_private_data_initializer(_HullTransformationData)
@@ -95,24 +101,11 @@ class Hull_Reformulation(GDP_to_MIP_Transformation):
     The transformation will create a new Block with a unique
     name beginning "_pyomo_gdp_hull_reformulation". It will contain an
     indexed Block named "relaxedDisjuncts" that will hold the relaxed
-    disjuncts.  This block is indexed by an integer indicating the order
-    in which the disjuncts were relaxed. Each block has a dictionary
-    "_constraintMap":
-
-        'srcConstraints': ComponentMap(<transformed constraint>:
-                                       <src constraint>),
-        'transformedConstraints':
-            ComponentMap(<src constraint container> :
-                         <transformed constraint container>,
-                         <src constraintData> : [<transformed constraintDatas>])
-
-    All transformed Disjuncts will have a pointer to the block their transformed
-    constraints are on, and all transformed Disjunctions will have a
-    pointer to the corresponding OR or XOR constraint.
-
-    The _pyomo_gdp_hull_reformulation block will have a ComponentMap
-    "_disaggregationConstraintMap":
-        <src var>:ComponentMap(<srcDisjunction>: <disaggregation constraint>)
+    disjuncts. This block is indexed by an integer indicating the order
+    in which the disjuncts were relaxed. All transformed Disjuncts will
+    have a pointer to the block their transformed constraints are on,
+    and all transformed Disjunctions will have a pointer to the
+    corresponding OR or XOR constraint.
     """
 
     CONFIG = cfg.ConfigDict('gdp.hull')
@@ -294,10 +287,6 @@ def _add_transformation_block(self, to_block):
         # Disjunctions we transform onto this block here.
         transBlock.disaggregationConstraints = Constraint(NonNegativeIntegers)
 
-        # This will map from srcVar to a map of srcDisjunction to the
-        # disaggregation constraint corresponding to srcDisjunction
-        transBlock._disaggregationConstraintMap = ComponentMap()
-
         # we are going to store some of the disaggregated vars directly here
         # when we have vars that don't appear in every disjunct
         transBlock._disaggregatedVars = Var(NonNegativeIntegers, dense=False)
@@ -330,7 +319,9 @@ def _transform_disjunctionData(
         )
 
         disaggregationConstraint = transBlock.disaggregationConstraints
-        disaggregationConstraintMap = transBlock._disaggregationConstraintMap
+        disaggregationConstraintMap = (
+            transBlock.private_data().disaggregation_constraint_map
+        )
         disaggregatedVars = transBlock._disaggregatedVars
         disaggregated_var_bounds = transBlock._boundsConstraints
 
@@ -499,13 +490,7 @@ def _transform_disjunctionData(
             # and update the map so that we can find this later. We index by
             # variable and the particular disjunction because there is a
             # different one for each disjunction
-            if var in disaggregationConstraintMap:
-                disaggregationConstraintMap[var][obj] = disaggregationConstraint[
-                    cons_idx
-                ]
-            else:
-                thismap = disaggregationConstraintMap[var] = ComponentMap()
-                thismap[obj] = disaggregationConstraint[cons_idx]
+            disaggregationConstraintMap[var][obj] = disaggregationConstraint[cons_idx]
 
         # deactivate for the writers
         obj.deactivate()
@@ -675,7 +660,7 @@ def _transform_constraint(
     ):
         # we will put a new transformed constraint on the relaxation block.
         relaxationBlock = disjunct._transformation_block()
-        constraintMap = relaxationBlock._constraintMap
+        constraint_map = relaxationBlock.private_data('pyomo.gdp')
 
         # We will make indexes from ({obj.local_name} x obj.index_set() x ['lb',
         # 'ub']), but don't bother construct that set here, as taking Cartesian
@@ -757,19 +742,19 @@ def _transform_constraint(
                         # this variable, so I'm going to return
                         # it. Alternatively we could return an empty list, but I
                         # think I like this better.
-                        constraintMap['transformedConstraints'][c] = [v[0]]
+                        constraint_map.transformed_constraints[c].append(v[0])
                         # Reverse map also (this is strange)
-                        constraintMap['srcConstraints'][v[0]] = c
+                        constraint_map.src_constraint[v[0]] = c
                         continue
                     newConsExpr = expr - (1 - y) * h_0 == c.lower * y
 
                 if obj.is_indexed():
                     newConstraint.add((name, i, 'eq'), newConsExpr)
                     # map the _ConstraintDatas (we mapped the container above)
-                    constraintMap['transformedConstraints'][c] = [
+                    constraint_map.transformed_constraints[c].append(
                         newConstraint[name, i, 'eq']
-                    ]
-                    constraintMap['srcConstraints'][newConstraint[name, i, 'eq']] = c
+                    )
+                    constraint_map.src_constraint[newConstraint[name, i, 'eq']] = c
                 else:
                     newConstraint.add((name, 'eq'), newConsExpr)
                     # map to the _ConstraintData (And yes, for
@@ -779,10 +764,10 @@ def _transform_constraint(
                     # IndexedConstraints, we can map the container to the
                     # container, but more importantly, we are mapping the
                     # _ConstraintDatas to each other above)
-                    constraintMap['transformedConstraints'][c] = [
+                    constraint_map.transformed_constraints[c].append(
                         newConstraint[name, 'eq']
-                    ]
-                    constraintMap['srcConstraints'][newConstraint[name, 'eq']] = c
+                    )
+                    constraint_map.src_constraint[newConstraint[name, 'eq']] = c
 
                 continue
 
@@ -797,16 +782,16 @@ def _transform_constraint(
 
                 if obj.is_indexed():
                     newConstraint.add((name, i, 'lb'), newConsExpr)
-                    constraintMap['transformedConstraints'][c] = [
+                    constraint_map.transformed_constraints[c].append(
                         newConstraint[name, i, 'lb']
-                    ]
-                    constraintMap['srcConstraints'][newConstraint[name, i, 'lb']] = c
+                    )
+                    constraint_map.src_constraint[newConstraint[name, i, 'lb']] = c
                 else:
                     newConstraint.add((name, 'lb'), newConsExpr)
-                    constraintMap['transformedConstraints'][c] = [
+                    constraint_map.transformed_constraints[c].append(
                         newConstraint[name, 'lb']
-                    ]
-                    constraintMap['srcConstraints'][newConstraint[name, 'lb']] = c
+                    )
+                    constraint_map.src_constraint[newConstraint[name, 'lb']] = c
 
             if c.upper is not None:
                 if self._generate_debug_messages:
@@ -821,24 +806,16 @@ def _transform_constraint(
                     newConstraint.add((name, i, 'ub'), newConsExpr)
                     # map (have to account for fact we might have created list
                     # above
-                    transformed = constraintMap['transformedConstraints'].get(c)
-                    if transformed is not None:
-                        transformed.append(newConstraint[name, i, 'ub'])
-                    else:
-                        constraintMap['transformedConstraints'][c] = [
-                            newConstraint[name, i, 'ub']
-                        ]
-                    constraintMap['srcConstraints'][newConstraint[name, i, 'ub']] = c
+                    constraint_map.transformed_constraints[c].append(
+                        newConstraint[name, i, 'ub']
+                    )
+                    constraint_map.src_constraint[newConstraint[name, i, 'ub']] = c
                 else:
                     newConstraint.add((name, 'ub'), newConsExpr)
-                    transformed = constraintMap['transformedConstraints'].get(c)
-                    if transformed is not None:
-                        transformed.append(newConstraint[name, 'ub'])
-                    else:
-                        constraintMap['transformedConstraints'][c] = [
-                            newConstraint[name, 'ub']
-                        ]
-                    constraintMap['srcConstraints'][newConstraint[name, 'ub']] = c
+                    constraint_map.transformed_constraints[c].append(
+                        newConstraint[name, 'ub']
+                    )
+                    constraint_map.src_constraint[newConstraint[name, 'ub']] = c
 
         # deactivate now that we have transformed
         obj.deactivate()
@@ -931,9 +908,11 @@ def get_disaggregation_constraint(
             )
 
         try:
-            cons = transBlock.parent_block()._disaggregationConstraintMap[original_var][
-                disjunction
-            ]
+            cons = (
+                transBlock.parent_block()
+                .private_data()
+                .disaggregation_constraint_map[original_var][disjunction]
+            )
         except:
             if raise_exception:
                 logger.error(
diff --git a/pyomo/gdp/plugins/multiple_bigm.py b/pyomo/gdp/plugins/multiple_bigm.py
index 6177de3c037..4dffd4e9f9a 100644
--- a/pyomo/gdp/plugins/multiple_bigm.py
+++ b/pyomo/gdp/plugins/multiple_bigm.py
@@ -359,7 +359,7 @@ def _transform_disjunct(self, obj, transBlock, active_disjuncts, Ms):
     def _transform_constraint(self, obj, disjunct, active_disjuncts, Ms):
         # we will put a new transformed constraint on the relaxation block.
         relaxationBlock = disjunct._transformation_block()
-        constraintMap = relaxationBlock._constraintMap
+        constraint_map = relaxationBlock.private_data('pyomo.gdp')
         transBlock = relaxationBlock.parent_block()
 
         # Though rare, it is possible to get naming conflicts here
@@ -378,7 +378,7 @@ def _transform_constraint(self, obj, disjunct, active_disjuncts, Ms):
                 continue
 
             if not self._config.only_mbigm_bound_constraints:
-                transformed = []
+                transformed = constraint_map.transformed_constraints[c]
                 if c.lower is not None:
                     rhs = sum(
                         Ms[c, disj][0] * disj.indicator_var.get_associated_binary()
@@ -397,8 +397,7 @@ def _transform_constraint(self, obj, disjunct, active_disjuncts, Ms):
                     newConstraint.add((i, 'ub'), c.body - c.upper <= rhs)
                     transformed.append(newConstraint[i, 'ub'])
                 for c_new in transformed:
-                    constraintMap['srcConstraints'][c_new] = [c]
-                constraintMap['transformedConstraints'][c] = transformed
+                    constraint_map.src_constraint[c_new] = [c]
             else:
                 lower = (None, None, None)
                 upper = (None, None, None)
@@ -427,7 +426,7 @@ def _transform_constraint(self, obj, disjunct, active_disjuncts, Ms):
                     M,
                     disjunct.indicator_var.get_associated_binary(),
                     newConstraint,
-                    constraintMap,
+                    constraint_map,
                 )
 
             # deactivate now that we have transformed
@@ -496,6 +495,7 @@ def _transform_bound_constraints(self, active_disjuncts, transBlock, Ms):
                 relaxationBlock = self._get_disjunct_transformation_block(
                     disj, transBlock
                 )
+                constraint_map = relaxationBlock.private_data('pyomo.gdp')
                 if len(lower_dict) > 0:
                     M = lower_dict.get(disj, None)
                     if M is None:
@@ -527,39 +527,24 @@ def _transform_bound_constraints(self, active_disjuncts, transBlock, Ms):
             idx = i + offset
             if len(lower_dict) > 0:
                 transformed.add((idx, 'lb'), v >= lower_rhs)
-                relaxationBlock._constraintMap['srcConstraints'][
-                    transformed[idx, 'lb']
-                ] = []
+                constraint_map.src_constraint[transformed[idx, 'lb']] = []
                 for c, disj in lower_bound_constraints_by_var[v]:
-                    relaxationBlock._constraintMap['srcConstraints'][
-                        transformed[idx, 'lb']
-                    ].append(c)
-                    disj.transformation_block._constraintMap['transformedConstraints'][
-                        c
-                    ] = [transformed[idx, 'lb']]
+                    constraint_map.src_constraint[transformed[idx, 'lb']].append(c)
+                    disj.transformation_block.private_data(
+                        'pyomo.gdp'
+                    ).transformed_constraints[c].append(transformed[idx, 'lb'])
             if len(upper_dict) > 0:
                 transformed.add((idx, 'ub'), v <= upper_rhs)
-                relaxationBlock._constraintMap['srcConstraints'][
-                    transformed[idx, 'ub']
-                ] = []
+                constraint_map.src_constraint[transformed[idx, 'ub']] = []
                 for c, disj in upper_bound_constraints_by_var[v]:
-                    relaxationBlock._constraintMap['srcConstraints'][
-                        transformed[idx, 'ub']
-                    ].append(c)
+                    constraint_map.src_constraint[transformed[idx, 'ub']].append(c)
                     # might already be here if it had an upper bound
-                    if (
-                        c
-                        in disj.transformation_block._constraintMap[
-                            'transformedConstraints'
-                        ]
-                    ):
-                        disj.transformation_block._constraintMap[
-                            'transformedConstraints'
-                        ][c].append(transformed[idx, 'ub'])
-                    else:
-                        disj.transformation_block._constraintMap[
-                            'transformedConstraints'
-                        ][c] = [transformed[idx, 'ub']]
+                    disj_constraint_map = disj.transformation_block.private_data(
+                        'pyomo.gdp'
+                    )
+                    disj_constraint_map.transformed_constraints[c].append(
+                        transformed[idx, 'ub']
+                    )
 
         return transformed_constraints
 
diff --git a/pyomo/gdp/tests/test_bigm.py b/pyomo/gdp/tests/test_bigm.py
index 00efcb46485..2383d4587f5 100644
--- a/pyomo/gdp/tests/test_bigm.py
+++ b/pyomo/gdp/tests/test_bigm.py
@@ -1316,18 +1316,10 @@ def test_do_not_transform_deactivated_constraintDatas(self):
         bigm.apply_to(m)
 
         # the real test: This wasn't transformed
-        log = StringIO()
-        with LoggingIntercept(log, 'pyomo.gdp', logging.ERROR):
-            self.assertRaisesRegex(
-                KeyError,
-                r".*b.simpledisj1.c\[1\]",
-                bigm.get_transformed_constraints,
-                m.b.simpledisj1.c[1],
-            )
-        self.assertRegex(
-            log.getvalue(),
-            r".*Constraint 'b.simpledisj1.c\[1\]' has not been transformed.",
-        )
+        with self.assertRaisesRegex(
+            GDP_Error, r"Constraint 'b.simpledisj1.c\[1\]' has not been transformed."
+        ):
+            bigm.get_transformed_constraints(m.b.simpledisj1.c[1])
 
         # and the rest of the container was transformed
         cons_list = bigm.get_transformed_constraints(m.b.simpledisj1.c[2])
@@ -2272,18 +2264,12 @@ def check_all_but_evil1_b_anotherblock_constraint_transformed(self, m):
         self.assertEqual(len(evil1), 2)
         self.assertIs(evil1[0].parent_block(), disjBlock[1])
         self.assertIs(evil1[1].parent_block(), disjBlock[1])
-        out = StringIO()
-        with LoggingIntercept(out, 'pyomo.gdp', logging.ERROR):
-            self.assertRaisesRegex(
-                KeyError,
-                r".*.evil\[1\].b.anotherblock.c",
-                bigm.get_transformed_constraints,
-                m.evil[1].b.anotherblock.c,
-            )
-        self.assertRegex(
-            out.getvalue(),
-            r".*Constraint 'evil\[1\].b.anotherblock.c' has not been transformed.",
-        )
+        with self.assertRaisesRegex(
+            GDP_Error,
+            r"Constraint 'evil\[1\].b.anotherblock.c' has not been transformed.",
+        ):
+            bigm.get_transformed_constraints(m.evil[1].b.anotherblock.c)
+
         evil1 = bigm.get_transformed_constraints(m.evil[1].bb[1].c)
         self.assertEqual(len(evil1), 2)
         self.assertIs(evil1[0].parent_block(), disjBlock[1])
diff --git a/pyomo/gdp/tests/test_hull.py b/pyomo/gdp/tests/test_hull.py
index 858764759ee..55edf244731 100644
--- a/pyomo/gdp/tests/test_hull.py
+++ b/pyomo/gdp/tests/test_hull.py
@@ -897,18 +897,10 @@ def test_do_not_transform_deactivated_constraintDatas(self):
         hull = TransformationFactory('gdp.hull')
         hull.apply_to(m)
         # can't ask for simpledisj1.c[1]: it wasn't transformed
-        log = StringIO()
-        with LoggingIntercept(log, 'pyomo.gdp', logging.ERROR):
-            self.assertRaisesRegex(
-                KeyError,
-                r".*b.simpledisj1.c\[1\]",
-                hull.get_transformed_constraints,
-                m.b.simpledisj1.c[1],
-            )
-        self.assertRegex(
-            log.getvalue(),
-            r".*Constraint 'b.simpledisj1.c\[1\]' has not been transformed.",
-        )
+        with self.assertRaisesRegex(
+            GDP_Error, r"Constraint 'b.simpledisj1.c\[1\]' has not been transformed."
+        ):
+            hull.get_transformed_constraints(m.b.simpledisj1.c[1])
 
         # this fixes a[2] to 0, so we should get the disggregated var
         transformed = hull.get_transformed_constraints(m.b.simpledisj1.c[2])
@@ -2314,8 +2306,7 @@ def test_mapping_method_errors(self):
         with LoggingIntercept(log, 'pyomo.gdp.hull', logging.ERROR):
             self.assertRaisesRegex(
                 KeyError,
-                r".*_pyomo_gdp_hull_reformulation.relaxedDisjuncts\[1\]."
-                r"disaggregatedVars.w",
+                r".*disjunction",
                 hull.get_disaggregation_constraint,
                 m.d[1].transformation_block.disaggregatedVars.w,
                 m.disjunction,
diff --git a/pyomo/gdp/util.py b/pyomo/gdp/util.py
index 57eef29eded..fe11975954d 100644
--- a/pyomo/gdp/util.py
+++ b/pyomo/gdp/util.py
@@ -477,16 +477,17 @@ def get_src_constraint(transformedConstraint):
     a transformation block
     """
     transBlock = transformedConstraint.parent_block()
+    src_constraints = transBlock.private_data('pyomo.gdp').src_constraint
     # This should be our block, so if it's not, the user messed up and gave
     # us the wrong thing. If they happen to also have a _constraintMap then
     # the world is really against us.
-    if not hasattr(transBlock, "_constraintMap"):
+    if transformedConstraint not in src_constraints:
         raise GDP_Error(
             "Constraint '%s' is not a transformed constraint"
             % transformedConstraint.name
         )
     # if something goes wrong here, it's a bug in the mappings.
-    return transBlock._constraintMap['srcConstraints'][transformedConstraint]
+    return src_constraints[transformedConstraint]
 
 
 def _find_parent_disjunct(constraint):
@@ -537,11 +538,15 @@ def get_transformed_constraints(srcConstraint):
             "from any of its _ComponentDatas.)"
         )
     transBlock = _get_constraint_transBlock(srcConstraint)
-    try:
-        return transBlock._constraintMap['transformedConstraints'][srcConstraint]
-    except:
-        logger.error("Constraint '%s' has not been transformed." % srcConstraint.name)
-        raise
+    transformed_constraints = transBlock.private_data(
+        'pyomo.gdp'
+    ).transformed_constraints
+    if srcConstraint in transformed_constraints:
+        return transformed_constraints[srcConstraint]
+    else:
+        raise GDP_Error(
+            "Constraint '%s' has not been transformed." % srcConstraint.name
+        )
 
 
 def _warn_for_active_disjunct(innerdisjunct, outerdisjunct):
diff --git a/pyomo/repn/plugins/nl_writer.py b/pyomo/repn/plugins/nl_writer.py
index f3ff94ea8c9..a256cd1b900 100644
--- a/pyomo/repn/plugins/nl_writer.py
+++ b/pyomo/repn/plugins/nl_writer.py
@@ -1760,7 +1760,7 @@ def _linear_presolve(self, comp_by_linear_var, lcon_by_linear_nnz, var_bounds):
                 id2_isdiscrete = var_map[id2].domain.isdiscrete()
                 if var_map[_id].domain.isdiscrete() ^ id2_isdiscrete:
                     # if only one variable is discrete, then we need to
-                    # substiitute out the other
+                    # substitute out the other
                     if id2_isdiscrete:
                         _id, id2 = id2, _id
                         coef, coef2 = coef2, coef
@@ -1820,10 +1820,15 @@ def _linear_presolve(self, comp_by_linear_var, lcon_by_linear_nnz, var_bounds):
                 # appropriately (that expr_info is persisting in the
                 # eliminated_vars dict - and we will use that to
                 # update other linear expressions later.)
+                old_nnz = len(expr_info.linear)
                 c = expr_info.linear.pop(_id, 0)
+                nnz = old_nnz - 1
                 expr_info.const += c * b
                 if x in expr_info.linear:
                     expr_info.linear[x] += c * a
+                    if expr_info.linear[x] == 0:
+                        nnz -= 1
+                        coef = expr_info.linear.pop(x)
                 elif a:
                     expr_info.linear[x] = c * a
                     # replacing _id with x... NNZ is not changing,
@@ -1831,10 +1836,17 @@ def _linear_presolve(self, comp_by_linear_var, lcon_by_linear_nnz, var_bounds):
                     # this constraint
                     comp_by_linear_var[x].append((con_id, expr_info))
                     continue
-                # NNZ has been reduced by 1
-                nnz = len(expr_info.linear)
-                _old = lcon_by_linear_nnz[nnz + 1]
+                _old = lcon_by_linear_nnz[old_nnz]
                 if con_id in _old:
+                    if not nnz:
+                        if abs(expr_info.const) > TOL:
+                            # constraint is trivially infeasible
+                            raise InfeasibleConstraintException(
+                                "model contains a trivially infeasible constraint "
+                                f"{expr_info.const} == {coef}*{var_map[x]}"
+                            )
+                        # constraint is trivially feasible
+                        eliminated_cons.add(con_id)
                     lcon_by_linear_nnz[nnz][con_id] = _old.pop(con_id)
             # If variables were replaced by the variable that
             # we are currently eliminating, then we need to update
diff --git a/pyomo/repn/tests/ampl/test_nlv2.py b/pyomo/repn/tests/ampl/test_nlv2.py
index 8b95fc03bdb..86eb43d9a37 100644
--- a/pyomo/repn/tests/ampl/test_nlv2.py
+++ b/pyomo/repn/tests/ampl/test_nlv2.py
@@ -1096,7 +1096,6 @@ def test_log_timing(self):
         m.c1 = Constraint([1, 2], rule=lambda m, i: sum(m.x.values()) == 1)
         m.c2 = Constraint(expr=m.p * m.x[1] ** 2 + m.x[2] ** 3 <= 100)
 
-        self.maxDiff = None
         OUT = io.StringIO()
         with capture_output() as LOG:
             with report_timing(level=logging.DEBUG):
@@ -1683,6 +1682,131 @@ def test_presolve_named_expressions(self):
 G0 2	#obj
 0 0
 1 0
+""",
+                OUT.getvalue(),
+            )
+        )
+
+    def test_presolve_zero_coef(self):
+        m = ConcreteModel()
+        m.x = Var()
+        m.y = Var()
+        m.z = Var()
+        m.obj = Objective(expr=m.x**2 + m.y**2 + m.z**2)
+        m.c1 = Constraint(expr=m.x == m.y + m.z + 1.5)
+        m.c2 = Constraint(expr=m.z == -m.y)
+
+        OUT = io.StringIO()
+        with LoggingIntercept() as LOG:
+            nlinfo = nl_writer.NLWriter().write(
+                m, OUT, symbolic_solver_labels=True, linear_presolve=True
+            )
+        self.assertEqual(LOG.getvalue(), "")
+
+        self.assertEqual(nlinfo.eliminated_vars[0], (m.x, 1.5))
+        self.assertIs(nlinfo.eliminated_vars[1][0], m.y)
+        self.assertExpressionsEqual(
+            nlinfo.eliminated_vars[1][1], LinearExpression([-1.0 * m.z])
+        )
+
+        self.assertEqual(
+            *nl_diff(
+                """g3 1 1 0	# problem unknown
+ 1 0 1 0 0	#vars, constraints, objectives, ranges, eqns
+ 0 1 0 0 0 0	#nonlinear constrs, objs; ccons: lin, nonlin, nd, nzlb
+ 0 0	#network constraints: nonlinear, linear
+ 0 1 0	#nonlinear vars in constraints, objectives, both
+ 0 0 0 1	#linear network variables; functions; arith, flags
+ 0 0 0 0 0	#discrete variables: binary, integer, nonlinear (b,c,o)
+ 0 1	#nonzeros in Jacobian, obj. gradient
+ 3 1	#max name lengths: constraints, variables
+ 0 0 0 0 0	#common exprs: b,c,o,c1,o1
+O0 0	#obj
+o54	#sumlist
+3	#(n)
+o5	#^
+n1.5
+n2
+o5	#^
+o16	#-
+v0	#z
+n2
+o5	#^
+v0	#z
+n2
+x0	#initial guess
+r	#0 ranges (rhs's)
+b	#1 bounds (on variables)
+3	#z
+k0	#intermediate Jacobian column lengths
+G0 1	#obj
+0 0
+""",
+                OUT.getvalue(),
+            )
+        )
+
+        m.c3 = Constraint(expr=m.x == 2)
+        OUT = io.StringIO()
+        with LoggingIntercept() as LOG:
+            with self.assertRaisesRegex(
+                nl_writer.InfeasibleConstraintException,
+                r"model contains a trivially infeasible constraint 0.5 == 0.0\*y",
+            ):
+                nlinfo = nl_writer.NLWriter().write(
+                    m, OUT, symbolic_solver_labels=True, linear_presolve=True
+                )
+        self.assertEqual(LOG.getvalue(), "")
+
+        m.c1.set_value(m.x >= m.y + m.z + 1.5)
+        OUT = io.StringIO()
+        with LoggingIntercept() as LOG:
+            nlinfo = nl_writer.NLWriter().write(
+                m, OUT, symbolic_solver_labels=True, linear_presolve=True
+            )
+        self.assertEqual(LOG.getvalue(), "")
+
+        self.assertIs(nlinfo.eliminated_vars[0][0], m.y)
+        self.assertExpressionsEqual(
+            nlinfo.eliminated_vars[0][1], LinearExpression([-1.0 * m.z])
+        )
+        self.assertEqual(nlinfo.eliminated_vars[1], (m.x, 2))
+
+        self.assertEqual(
+            *nl_diff(
+                """g3 1 1 0	# problem unknown
+ 1 1 1 0 0	#vars, constraints, objectives, ranges, eqns
+ 0 1 0 0 0 0	#nonlinear constrs, objs; ccons: lin, nonlin, nd, nzlb
+ 0 0	#network constraints: nonlinear, linear
+ 0 1 0	#nonlinear vars in constraints, objectives, both
+ 0 0 0 1	#linear network variables; functions; arith, flags
+ 0 0 0 0 0	#discrete variables: binary, integer, nonlinear (b,c,o)
+ 0 1	#nonzeros in Jacobian, obj. gradient
+ 3 1	#max name lengths: constraints, variables
+ 0 0 0 0 0	#common exprs: b,c,o,c1,o1
+C0	#c1
+n0
+O0 0	#obj
+o54	#sumlist
+3	#(n)
+o5	#^
+n2
+n2
+o5	#^
+o16	#-
+v0	#z
+n2
+o5	#^
+v0	#z
+n2
+x0	#initial guess
+r	#1 ranges (rhs's)
+1 0.5	#c1
+b	#1 bounds (on variables)
+3	#z
+k0	#intermediate Jacobian column lengths
+G0 1	#obj
+0 0
 """,
                 OUT.getvalue(),
             )
diff --git a/pyomo/solvers/plugins/solvers/GAMS.py b/pyomo/solvers/plugins/solvers/GAMS.py
index e84cbdb441d..be3499a2f6b 100644
--- a/pyomo/solvers/plugins/solvers/GAMS.py
+++ b/pyomo/solvers/plugins/solvers/GAMS.py
@@ -41,7 +41,7 @@
 
 from pyomo.common.dependencies import attempt_import
 
-gdxcc, gdxcc_available = attempt_import('gdxcc', defer_check=True)
+gdxcc, gdxcc_available = attempt_import('gdxcc')
 
 logger = logging.getLogger('pyomo.solvers')