diff --git a/pyomo/repn/plugins/nl_writer.py b/pyomo/repn/plugins/nl_writer.py
index 59a83e61730..f79b5009122 100644
--- a/pyomo/repn/plugins/nl_writer.py
+++ b/pyomo/repn/plugins/nl_writer.py
@@ -1691,6 +1691,15 @@ def _linear_presolve(self, comp_by_linear_var, lcon_by_linear_nnz, var_bounds):
         if not self.config.linear_presolve:
             return eliminated_cons, eliminated_vars
 
+        # We need to record all named expressions with linear components
+        # so that any eliminated variables are removed from them.
+        for expr, info, _ in self.subexpression_cache.values():
+            if not info.linear:
+                continue
+            expr_id = id(expr)
+            for _id in info.linear:
+                comp_by_linear_var[_id].append((expr_id, info))
+
         fixed_vars = [
             _id for _id, (lb, ub) in var_bounds.items() if lb == ub and lb is not None
         ]
diff --git a/pyomo/repn/tests/ampl/test_nlv2.py b/pyomo/repn/tests/ampl/test_nlv2.py
index f0bc5262183..fe5f422d323 100644
--- a/pyomo/repn/tests/ampl/test_nlv2.py
+++ b/pyomo/repn/tests/ampl/test_nlv2.py
@@ -1481,7 +1481,6 @@ def test_presolve_almost_lower_triangular_nonlinear(self):
         # Note: bounds on x[1] are:
         #   min(22/3, 82/17, 23/4, -39/-6) == 4.823529411764706
         #   max(2/3, 62/17, 3/4, -19/-6) == 3.6470588235294117
-        print(OUT.getvalue())
         self.assertEqual(
             *nl_diff(
                 """g3 1 1 0	# problem unknown
@@ -1558,6 +1557,83 @@ def test_presolve_lower_triangular_out_of_bounds(self):
                 nlinfo = nl_writer.NLWriter().write(m, OUT, linear_presolve=True)
         self.assertEqual(LOG.getvalue(), "")
 
+    def test_presolve_named_expressions(self):
+        # Test from #3055
+        m = pyo.ConcreteModel()
+        m.x = pyo.Var([1, 2, 3], initialize=1, bounds=(0, 10))
+        m.subexpr = pyo.Expression(pyo.Integers)
+        m.subexpr[1] = m.x[1] + m.x[2]
+        m.eq = pyo.Constraint(pyo.Integers)
+        m.eq[1] = m.x[1] == 7
+        m.eq[2] = m.x[3] == 0.1 * m.subexpr[1] * m.x[2]
+        m.obj = pyo.Objective(expr=m.x[1] ** 2 + m.x[2] ** 2 + m.x[3] ** 3)
+
+        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, [(m.x[1], nl_writer.AMPLRepn(7, {}, None))]
+        )
+
+        self.assertEqual(
+            *nl_diff(
+                """g3 1 1 0	# problem unknown
+ 2 1 1 0 1 	# vars, constraints, objectives, ranges, eqns
+ 1 1 0 0 0 0	# nonlinear constrs, objs; ccons: lin, nonlin, nd, nzlb
+ 0 0	# network constraints: nonlinear, linear
+ 1 2 1 	# 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)
+ 2 2 	# nonzeros in Jacobian, obj. gradient
+ 5 4	# max name lengths: constraints, variables
+ 0 0 0 1 0	# common exprs: b,c,o,c1,o1
+V2 1 1	#subexpr[1]
+0 1
+n7.0
+C0	#eq[2]
+o16	#-
+o2	#*
+o2	#*
+n0.1
+v2	#subexpr[1]
+v0	#x[2]
+O0 0	#obj
+o54	# sumlist
+3	# (n)
+o5	#^
+n7.0
+n2
+o5	#^
+v0	#x[2]
+n2
+o5	#^
+v1	#x[3]
+n3
+x2	# initial guess
+0 1	#x[2]
+1 1	#x[3]
+r	#1 ranges (rhs's)
+4 0	#eq[2]
+b	#2 bounds (on variables)
+0 0 10	#x[2]
+0 0 10	#x[3]
+k1	#intermediate Jacobian column lengths
+1
+J0 2	#eq[2]
+0 0
+1 1
+G0 2	#obj
+0 0
+1 0
+""",
+                OUT.getvalue(),
+            )
+        )
+
     def test_scaling(self):
         m = pyo.ConcreteModel()
         m.x = pyo.Var(initialize=0)
@@ -1665,7 +1741,7 @@ def test_scaling(self):
         self.assertEqual(LOG.getvalue(), "")
 
         nl2 = OUT.getvalue()
-        print(nl2)
+
         self.assertEqual(
             *nl_diff(
                 """g3 1 1 0	# problem unknown
@@ -1759,7 +1835,7 @@ def test_named_expressions(self):
 
         OUT = io.StringIO()
         nl_writer.NLWriter().write(m, OUT, symbolic_solver_labels=True)
-        print(OUT.getvalue())
+
         self.assertEqual(
             *nl_diff(
                 """g3 1 1 0	# problem unknown