diff --git a/pyomo/contrib/appsi/base.py b/pyomo/contrib/appsi/base.py
index ca7255d5628..e2a694cfbd6 100644
--- a/pyomo/contrib/appsi/base.py
+++ b/pyomo/contrib/appsi/base.py
@@ -900,7 +900,7 @@ def invalidate(self):
 
 
 class PersistentBase(abc.ABC):
-    def __init__(self, only_child_vars=True):
+    def __init__(self, only_child_vars=False):
         self._model = None
         self._active_constraints = dict()  # maps constraint to (lower, body, upper)
         self._vars = dict()  # maps var id to (var, lb, ub, fixed, domain, value)
@@ -1334,9 +1334,6 @@ def update(self, timer: HierarchicalTimer = None):
         self.remove_constraints(old_cons)
         self.remove_sos_constraints(old_sos)
         timer.stop('cons')
-        timer.start('vars')
-        self.remove_variables(old_vars)
-        timer.stop('vars')
         timer.start('params')
         self.remove_params(old_params)
 
@@ -1463,6 +1460,12 @@ def update(self, timer: HierarchicalTimer = None):
             self.set_objective(pyomo_obj)
         timer.stop('objective')
 
+        # this has to be done after the objective and constraints in case the
+        # old objective/constraints use old variables
+        timer.start('vars')
+        self.remove_variables(old_vars)
+        timer.stop('vars')
+
 
 legacy_termination_condition_map = {
     TerminationCondition.unknown: LegacyTerminationCondition.unknown,
diff --git a/pyomo/contrib/appsi/solvers/cbc.py b/pyomo/contrib/appsi/solvers/cbc.py
index b31a96dbf8a..a3aae2a9213 100644
--- a/pyomo/contrib/appsi/solvers/cbc.py
+++ b/pyomo/contrib/appsi/solvers/cbc.py
@@ -64,7 +64,7 @@ def __init__(
 
 
 class Cbc(PersistentSolver):
-    def __init__(self, only_child_vars=True):
+    def __init__(self, only_child_vars=False):
         self._config = CbcConfig()
         self._solver_options = dict()
         self._writer = LPWriter(only_child_vars=only_child_vars)
diff --git a/pyomo/contrib/appsi/solvers/cplex.py b/pyomo/contrib/appsi/solvers/cplex.py
index 6c5e281ffac..f03bee6ecc5 100644
--- a/pyomo/contrib/appsi/solvers/cplex.py
+++ b/pyomo/contrib/appsi/solvers/cplex.py
@@ -67,7 +67,7 @@ def __init__(self, solver):
 class Cplex(PersistentSolver):
     _available = None
 
-    def __init__(self, only_child_vars=True):
+    def __init__(self, only_child_vars=False):
         self._config = CplexConfig()
         self._solver_options = dict()
         self._writer = LPWriter(only_child_vars=only_child_vars)
diff --git a/pyomo/contrib/appsi/solvers/gurobi.py b/pyomo/contrib/appsi/solvers/gurobi.py
index 2362612e9ee..a173c69abc6 100644
--- a/pyomo/contrib/appsi/solvers/gurobi.py
+++ b/pyomo/contrib/appsi/solvers/gurobi.py
@@ -232,7 +232,7 @@ class Gurobi(PersistentBase, PersistentSolver):
     _available = None
     _num_instances = 0
 
-    def __init__(self, only_child_vars=True):
+    def __init__(self, only_child_vars=False):
         super(Gurobi, self).__init__(only_child_vars=only_child_vars)
         self._num_instances += 1
         self._config = GurobiConfig()
diff --git a/pyomo/contrib/appsi/solvers/highs.py b/pyomo/contrib/appsi/solvers/highs.py
index 9de5accfb91..1cf60d36e10 100644
--- a/pyomo/contrib/appsi/solvers/highs.py
+++ b/pyomo/contrib/appsi/solvers/highs.py
@@ -144,7 +144,7 @@ class Highs(PersistentBase, PersistentSolver):
 
     _available = None
 
-    def __init__(self, only_child_vars=True):
+    def __init__(self, only_child_vars=False):
         super().__init__(only_child_vars=only_child_vars)
         self._config = HighsConfig()
         self._solver_options = dict()
diff --git a/pyomo/contrib/appsi/solvers/ipopt.py b/pyomo/contrib/appsi/solvers/ipopt.py
index fde4c55073d..d38a836a2ac 100644
--- a/pyomo/contrib/appsi/solvers/ipopt.py
+++ b/pyomo/contrib/appsi/solvers/ipopt.py
@@ -127,7 +127,7 @@ def __init__(
 
 
 class Ipopt(PersistentSolver):
-    def __init__(self, only_child_vars=True):
+    def __init__(self, only_child_vars=False):
         self._config = IpoptConfig()
         self._solver_options = dict()
         self._writer = NLWriter(only_child_vars=only_child_vars)
diff --git a/pyomo/contrib/appsi/solvers/tests/test_gurobi_persistent.py b/pyomo/contrib/appsi/solvers/tests/test_gurobi_persistent.py
index 6366077642d..b032f5c827e 100644
--- a/pyomo/contrib/appsi/solvers/tests/test_gurobi_persistent.py
+++ b/pyomo/contrib/appsi/solvers/tests/test_gurobi_persistent.py
@@ -373,6 +373,7 @@ def test_quadratic_constraint_attr(self):
     def test_var_attr(self):
         m = pe.ConcreteModel()
         m.x = pe.Var(within=pe.Binary)
+        m.obj = pe.Objective(expr=m.x)
 
         opt = Gurobi()
         opt.set_instance(m)
@@ -694,6 +695,8 @@ def test_update7(self):
         m.y = pe.Var()
 
         opt = self.opt
+        orig_only_child_vars = opt._only_child_vars
+        opt._only_child_vars = True
         opt.set_instance(m)
         self.assertEqual(opt._solver_model.getAttr('NumVars'), 2)
 
@@ -712,3 +715,4 @@ def test_update7(self):
         opt.remove_variables([m.x])
         opt.update()
         self.assertEqual(opt._solver_model.getAttr('NumVars'), 1)
+        opt._only_child_vars = orig_only_child_vars
diff --git a/pyomo/contrib/appsi/solvers/tests/test_persistent_solvers.py b/pyomo/contrib/appsi/solvers/tests/test_persistent_solvers.py
index bafccb3527c..1520407e294 100644
--- a/pyomo/contrib/appsi/solvers/tests/test_persistent_solvers.py
+++ b/pyomo/contrib/appsi/solvers/tests/test_persistent_solvers.py
@@ -30,6 +30,7 @@
 nlp_solvers = [('ipopt', Ipopt)]
 qcp_solvers = [('gurobi', Gurobi), ('ipopt', Ipopt), ('cplex', Cplex)]
 miqcqp_solvers = [('gurobi', Gurobi), ('cplex', Cplex)]
+only_child_vars_options = [True, False]
 
 
 """
@@ -66,12 +67,46 @@
 """
 
 
+def _load_tests(solver_list, only_child_vars_list):
+    res = list()
+    for solver_name, solver in solver_list:
+        for child_var_option in only_child_vars_list:
+            test_name = f"{solver_name}_only_child_vars_{child_var_option}"
+            res.append((test_name, solver, child_var_option))
+    return res
+
+
 @unittest.skipUnless(cmodel_available, 'appsi extensions are not available')
 @unittest.skipUnless(numpy_available, 'numpy is not available')
 class TestSolvers(unittest.TestCase):
-    @parameterized.expand(input=all_solvers)
-    def test_stale_vars(self, name: str, opt_class: Type[PersistentSolver]):
-        opt: PersistentSolver = opt_class()
+    @parameterized.expand(input=_load_tests(all_solvers, only_child_vars_options))
+    def test_remove_variable_and_objective(
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
+    ):
+        # this test is for issue #2888
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
+        if not opt.available():
+            raise unittest.SkipTest
+        m = pe.ConcreteModel()
+        m.x = pe.Var(bounds=(2, None))
+        m.obj = pe.Objective(expr=m.x)
+        res = opt.solve(m)
+        self.assertEqual(res.termination_condition, TerminationCondition.optimal)
+        self.assertAlmostEqual(m.x.value, 2)
+
+        del m.x
+        del m.obj
+        m.x = pe.Var(bounds=(2, None))
+        m.obj = pe.Objective(expr=m.x)
+        res = opt.solve(m)
+        self.assertEqual(res.termination_condition, TerminationCondition.optimal)
+        self.assertAlmostEqual(m.x.value, 2)
+
+    @parameterized.expand(input=_load_tests(all_solvers, only_child_vars_options))
+    def test_stale_vars(
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
+    ):
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
         m = pe.ConcreteModel()
@@ -112,9 +147,11 @@ def test_stale_vars(self, name: str, opt_class: Type[PersistentSolver]):
         self.assertFalse(m.y.stale)
         self.assertTrue(m.z.stale)
 
-    @parameterized.expand(input=all_solvers)
-    def test_range_constraint(self, name: str, opt_class: Type[PersistentSolver]):
-        opt: PersistentSolver = opt_class()
+    @parameterized.expand(input=_load_tests(all_solvers, only_child_vars_options))
+    def test_range_constraint(
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
+    ):
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
         m = pe.ConcreteModel()
@@ -133,9 +170,11 @@ def test_range_constraint(self, name: str, opt_class: Type[PersistentSolver]):
         duals = opt.get_duals()
         self.assertAlmostEqual(duals[m.c], 1)
 
-    @parameterized.expand(input=all_solvers)
-    def test_reduced_costs(self, name: str, opt_class: Type[PersistentSolver]):
-        opt: PersistentSolver = opt_class()
+    @parameterized.expand(input=_load_tests(all_solvers, only_child_vars_options))
+    def test_reduced_costs(
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
+    ):
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
         m = pe.ConcreteModel()
@@ -150,9 +189,11 @@ def test_reduced_costs(self, name: str, opt_class: Type[PersistentSolver]):
         self.assertAlmostEqual(rc[m.x], 3)
         self.assertAlmostEqual(rc[m.y], 4)
 
-    @parameterized.expand(input=all_solvers)
-    def test_reduced_costs2(self, name: str, opt_class: Type[PersistentSolver]):
-        opt: PersistentSolver = opt_class()
+    @parameterized.expand(input=_load_tests(all_solvers, only_child_vars_options))
+    def test_reduced_costs2(
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
+    ):
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
         m = pe.ConcreteModel()
@@ -170,9 +211,11 @@ def test_reduced_costs2(self, name: str, opt_class: Type[PersistentSolver]):
         rc = opt.get_reduced_costs()
         self.assertAlmostEqual(rc[m.x], 1)
 
-    @parameterized.expand(input=all_solvers)
-    def test_param_changes(self, name: str, opt_class: Type[PersistentSolver]):
-        opt: PersistentSolver = opt_class()
+    @parameterized.expand(input=_load_tests(all_solvers, only_child_vars_options))
+    def test_param_changes(
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
+    ):
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
         m = pe.ConcreteModel()
@@ -202,13 +245,15 @@ def test_param_changes(self, name: str, opt_class: Type[PersistentSolver]):
             self.assertAlmostEqual(duals[m.c1], (1 + a1 / (a2 - a1)))
             self.assertAlmostEqual(duals[m.c2], a1 / (a2 - a1))
 
-    @parameterized.expand(input=all_solvers)
-    def test_immutable_param(self, name: str, opt_class: Type[PersistentSolver]):
+    @parameterized.expand(input=_load_tests(all_solvers, only_child_vars_options))
+    def test_immutable_param(
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
+    ):
         """
         This test is important because component_data_objects returns immutable params as floats.
         We want to make sure we process these correctly.
         """
-        opt: PersistentSolver = opt_class()
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
         m = pe.ConcreteModel()
@@ -238,9 +283,11 @@ def test_immutable_param(self, name: str, opt_class: Type[PersistentSolver]):
             self.assertAlmostEqual(duals[m.c1], (1 + a1 / (a2 - a1)))
             self.assertAlmostEqual(duals[m.c2], a1 / (a2 - a1))
 
-    @parameterized.expand(input=all_solvers)
-    def test_equality(self, name: str, opt_class: Type[PersistentSolver]):
-        opt: PersistentSolver = opt_class()
+    @parameterized.expand(input=_load_tests(all_solvers, only_child_vars_options))
+    def test_equality(
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
+    ):
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
         m = pe.ConcreteModel()
@@ -270,9 +317,11 @@ def test_equality(self, name: str, opt_class: Type[PersistentSolver]):
             self.assertAlmostEqual(duals[m.c1], (1 + a1 / (a2 - a1)))
             self.assertAlmostEqual(duals[m.c2], -a1 / (a2 - a1))
 
-    @parameterized.expand(input=all_solvers)
-    def test_linear_expression(self, name: str, opt_class: Type[PersistentSolver]):
-        opt: PersistentSolver = opt_class()
+    @parameterized.expand(input=_load_tests(all_solvers, only_child_vars_options))
+    def test_linear_expression(
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
+    ):
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
         m = pe.ConcreteModel()
@@ -304,9 +353,11 @@ def test_linear_expression(self, name: str, opt_class: Type[PersistentSolver]):
             self.assertAlmostEqual(res.best_feasible_objective, m.y.value)
             self.assertTrue(res.best_objective_bound <= m.y.value)
 
-    @parameterized.expand(input=all_solvers)
-    def test_no_objective(self, name: str, opt_class: Type[PersistentSolver]):
-        opt: PersistentSolver = opt_class()
+    @parameterized.expand(input=_load_tests(all_solvers, only_child_vars_options))
+    def test_no_objective(
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
+    ):
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
         m = pe.ConcreteModel()
@@ -336,9 +387,11 @@ def test_no_objective(self, name: str, opt_class: Type[PersistentSolver]):
             self.assertAlmostEqual(duals[m.c1], 0)
             self.assertAlmostEqual(duals[m.c2], 0)
 
-    @parameterized.expand(input=all_solvers)
-    def test_add_remove_cons(self, name: str, opt_class: Type[PersistentSolver]):
-        opt: PersistentSolver = opt_class()
+    @parameterized.expand(input=_load_tests(all_solvers, only_child_vars_options))
+    def test_add_remove_cons(
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
+    ):
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
         m = pe.ConcreteModel()
@@ -386,9 +439,11 @@ def test_add_remove_cons(self, name: str, opt_class: Type[PersistentSolver]):
         self.assertAlmostEqual(duals[m.c1], -(1 + a1 / (a2 - a1)))
         self.assertAlmostEqual(duals[m.c2], a1 / (a2 - a1))
 
-    @parameterized.expand(input=all_solvers)
-    def test_results_infeasible(self, name: str, opt_class: Type[PersistentSolver]):
-        opt: PersistentSolver = opt_class()
+    @parameterized.expand(input=_load_tests(all_solvers, only_child_vars_options))
+    def test_results_infeasible(
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
+    ):
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
         m = pe.ConcreteModel()
@@ -430,9 +485,9 @@ def test_results_infeasible(self, name: str, opt_class: Type[PersistentSolver]):
         ):
             res.solution_loader.get_reduced_costs()
 
-    @parameterized.expand(input=all_solvers)
-    def test_duals(self, name: str, opt_class: Type[PersistentSolver]):
-        opt: PersistentSolver = opt_class()
+    @parameterized.expand(input=_load_tests(all_solvers, only_child_vars_options))
+    def test_duals(self, name: str, opt_class: Type[PersistentSolver], only_child_vars):
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
         m = pe.ConcreteModel()
@@ -453,11 +508,11 @@ def test_duals(self, name: str, opt_class: Type[PersistentSolver]):
         self.assertAlmostEqual(duals[m.c1], 0.5)
         self.assertNotIn(m.c2, duals)
 
-    @parameterized.expand(input=qcp_solvers)
+    @parameterized.expand(input=_load_tests(qcp_solvers, only_child_vars_options))
     def test_mutable_quadratic_coefficient(
-        self, name: str, opt_class: Type[PersistentSolver]
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
     ):
-        opt: PersistentSolver = opt_class()
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
         m = pe.ConcreteModel()
@@ -477,11 +532,11 @@ def test_mutable_quadratic_coefficient(
         self.assertAlmostEqual(m.x.value, 0.10256137418973625, 4)
         self.assertAlmostEqual(m.y.value, 0.0869525991355825, 4)
 
-    @parameterized.expand(input=qcp_solvers)
+    @parameterized.expand(input=_load_tests(qcp_solvers, only_child_vars_options))
     def test_mutable_quadratic_objective(
-        self, name: str, opt_class: Type[PersistentSolver]
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
     ):
-        opt: PersistentSolver = opt_class()
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
         m = pe.ConcreteModel()
@@ -504,9 +559,11 @@ def test_mutable_quadratic_objective(
         self.assertAlmostEqual(m.x.value, 0.6962249634573562, 4)
         self.assertAlmostEqual(m.y.value, 0.09227926676152151, 4)
 
-    @parameterized.expand(input=all_solvers)
-    def test_fixed_vars(self, name: str, opt_class: Type[PersistentSolver]):
-        opt: PersistentSolver = opt_class()
+    @parameterized.expand(input=_load_tests(all_solvers, only_child_vars_options))
+    def test_fixed_vars(
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
+    ):
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         for treat_fixed_vars_as_params in [True, False]:
             opt.update_config.treat_fixed_vars_as_params = treat_fixed_vars_as_params
             if not opt.available():
@@ -542,9 +599,11 @@ def test_fixed_vars(self, name: str, opt_class: Type[PersistentSolver]):
             self.assertAlmostEqual(m.x.value, 0)
             self.assertAlmostEqual(m.y.value, 2)
 
-    @parameterized.expand(input=all_solvers)
-    def test_fixed_vars_2(self, name: str, opt_class: Type[PersistentSolver]):
-        opt: PersistentSolver = opt_class()
+    @parameterized.expand(input=_load_tests(all_solvers, only_child_vars_options))
+    def test_fixed_vars_2(
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
+    ):
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         opt.update_config.treat_fixed_vars_as_params = True
         if not opt.available():
             raise unittest.SkipTest
@@ -579,9 +638,11 @@ def test_fixed_vars_2(self, name: str, opt_class: Type[PersistentSolver]):
         self.assertAlmostEqual(m.x.value, 0)
         self.assertAlmostEqual(m.y.value, 2)
 
-    @parameterized.expand(input=all_solvers)
-    def test_fixed_vars_3(self, name: str, opt_class: Type[PersistentSolver]):
-        opt: PersistentSolver = opt_class()
+    @parameterized.expand(input=_load_tests(all_solvers, only_child_vars_options))
+    def test_fixed_vars_3(
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
+    ):
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         opt.update_config.treat_fixed_vars_as_params = True
         if not opt.available():
             raise unittest.SkipTest
@@ -594,9 +655,11 @@ def test_fixed_vars_3(self, name: str, opt_class: Type[PersistentSolver]):
         res = opt.solve(m)
         self.assertAlmostEqual(m.x.value, 2)
 
-    @parameterized.expand(input=nlp_solvers)
-    def test_fixed_vars_4(self, name: str, opt_class: Type[PersistentSolver]):
-        opt: PersistentSolver = opt_class()
+    @parameterized.expand(input=_load_tests(nlp_solvers, only_child_vars_options))
+    def test_fixed_vars_4(
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
+    ):
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         opt.update_config.treat_fixed_vars_as_params = True
         if not opt.available():
             raise unittest.SkipTest
@@ -613,11 +676,11 @@ def test_fixed_vars_4(self, name: str, opt_class: Type[PersistentSolver]):
         self.assertAlmostEqual(m.x.value, 2**0.5)
         self.assertAlmostEqual(m.y.value, 2**0.5)
 
-    @parameterized.expand(input=all_solvers)
+    @parameterized.expand(input=_load_tests(all_solvers, only_child_vars_options))
     def test_mutable_param_with_range(
-        self, name: str, opt_class: Type[PersistentSolver]
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
     ):
-        opt: PersistentSolver = opt_class()
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
         try:
@@ -703,9 +766,11 @@ def test_mutable_param_with_range(
                 self.assertAlmostEqual(duals[m.con1], (1 + a1 / (a2 - a1)), 6)
                 self.assertAlmostEqual(duals[m.con2], -a1 / (a2 - a1), 6)
 
-    @parameterized.expand(input=all_solvers)
-    def test_add_and_remove_vars(self, name: str, opt_class: Type[PersistentSolver]):
-        opt = opt_class()
+    @parameterized.expand(input=_load_tests(all_solvers, only_child_vars_options))
+    def test_add_and_remove_vars(
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
+    ):
+        opt = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
         m = pe.ConcreteModel()
@@ -730,7 +795,8 @@ def test_add_and_remove_vars(self, name: str, opt_class: Type[PersistentSolver])
         b2 = 1
         m.c1 = pe.Constraint(expr=(0, m.y - a1 * m.x - b1, None))
         m.c2 = pe.Constraint(expr=(None, -m.y + a2 * m.x + b2, 0))
-        opt.add_variables([m.x])
+        if only_child_vars:
+            opt.add_variables([m.x])
         opt.add_constraints([m.c1, m.c2])
         res = opt.solve(m)
         self.assertEqual(res.termination_condition, TerminationCondition.optimal)
@@ -738,7 +804,8 @@ def test_add_and_remove_vars(self, name: str, opt_class: Type[PersistentSolver])
         self.assertAlmostEqual(m.x.value, (b2 - b1) / (a1 - a2))
         self.assertAlmostEqual(m.y.value, a1 * (b2 - b1) / (a1 - a2) + b1)
         opt.remove_constraints([m.c1, m.c2])
-        opt.remove_variables([m.x])
+        if only_child_vars:
+            opt.remove_variables([m.x])
         m.x.value = None
         res = opt.solve(m)
         self.assertEqual(res.termination_condition, TerminationCondition.optimal)
@@ -748,9 +815,9 @@ def test_add_and_remove_vars(self, name: str, opt_class: Type[PersistentSolver])
         with self.assertRaises(Exception):
             opt.load_vars([m.x])
 
-    @parameterized.expand(input=nlp_solvers)
-    def test_exp(self, name: str, opt_class: Type[PersistentSolver]):
-        opt = opt_class()
+    @parameterized.expand(input=_load_tests(nlp_solvers, only_child_vars_options))
+    def test_exp(self, name: str, opt_class: Type[PersistentSolver], only_child_vars):
+        opt = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
         m = pe.ConcreteModel()
@@ -762,9 +829,9 @@ def test_exp(self, name: str, opt_class: Type[PersistentSolver]):
         self.assertAlmostEqual(m.x.value, -0.42630274815985264)
         self.assertAlmostEqual(m.y.value, 0.6529186341994245)
 
-    @parameterized.expand(input=nlp_solvers)
-    def test_log(self, name: str, opt_class: Type[PersistentSolver]):
-        opt = opt_class()
+    @parameterized.expand(input=_load_tests(nlp_solvers, only_child_vars_options))
+    def test_log(self, name: str, opt_class: Type[PersistentSolver], only_child_vars):
+        opt = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
         m = pe.ConcreteModel()
@@ -776,9 +843,11 @@ def test_log(self, name: str, opt_class: Type[PersistentSolver]):
         self.assertAlmostEqual(m.x.value, 0.6529186341994245)
         self.assertAlmostEqual(m.y.value, -0.42630274815985264)
 
-    @parameterized.expand(input=all_solvers)
-    def test_with_numpy(self, name: str, opt_class: Type[PersistentSolver]):
-        opt: PersistentSolver = opt_class()
+    @parameterized.expand(input=_load_tests(all_solvers, only_child_vars_options))
+    def test_with_numpy(
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
+    ):
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
         m = pe.ConcreteModel()
@@ -804,9 +873,11 @@ def test_with_numpy(self, name: str, opt_class: Type[PersistentSolver]):
         self.assertAlmostEqual(m.x.value, (b2 - b1) / (a1 - a2))
         self.assertAlmostEqual(m.y.value, a1 * (b2 - b1) / (a1 - a2) + b1)
 
-    @parameterized.expand(input=all_solvers)
-    def test_bounds_with_params(self, name: str, opt_class: Type[PersistentSolver]):
-        opt: PersistentSolver = opt_class()
+    @parameterized.expand(input=_load_tests(all_solvers, only_child_vars_options))
+    def test_bounds_with_params(
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
+    ):
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
         m = pe.ConcreteModel()
@@ -836,9 +907,11 @@ def test_bounds_with_params(self, name: str, opt_class: Type[PersistentSolver]):
         res = opt.solve(m)
         self.assertAlmostEqual(m.y.value, 3)
 
-    @parameterized.expand(input=all_solvers)
-    def test_solution_loader(self, name: str, opt_class: Type[PersistentSolver]):
-        opt: PersistentSolver = opt_class()
+    @parameterized.expand(input=_load_tests(all_solvers, only_child_vars_options))
+    def test_solution_loader(
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
+    ):
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
         m = pe.ConcreteModel()
@@ -887,9 +960,11 @@ def test_solution_loader(self, name: str, opt_class: Type[PersistentSolver]):
         self.assertIn(m.c1, duals)
         self.assertAlmostEqual(duals[m.c1], 1)
 
-    @parameterized.expand(input=all_solvers)
-    def test_time_limit(self, name: str, opt_class: Type[PersistentSolver]):
-        opt: PersistentSolver = opt_class()
+    @parameterized.expand(input=_load_tests(all_solvers, only_child_vars_options))
+    def test_time_limit(
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
+    ):
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
         from sys import platform
@@ -947,9 +1022,11 @@ def test_time_limit(self, name: str, opt_class: Type[PersistentSolver]):
                 res.termination_condition, TerminationCondition.maxTimeLimit
             )
 
-    @parameterized.expand(input=all_solvers)
-    def test_objective_changes(self, name: str, opt_class: Type[PersistentSolver]):
-        opt: PersistentSolver = opt_class()
+    @parameterized.expand(input=_load_tests(all_solvers, only_child_vars_options))
+    def test_objective_changes(
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
+    ):
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
         m = pe.ConcreteModel()
@@ -1006,9 +1083,11 @@ def test_objective_changes(self, name: str, opt_class: Type[PersistentSolver]):
         res = opt.solve(m)
         self.assertAlmostEqual(res.best_feasible_objective, 4)
 
-    @parameterized.expand(input=all_solvers)
-    def test_domain(self, name: str, opt_class: Type[PersistentSolver]):
-        opt: PersistentSolver = opt_class()
+    @parameterized.expand(input=_load_tests(all_solvers, only_child_vars_options))
+    def test_domain(
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
+    ):
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
         m = pe.ConcreteModel()
@@ -1030,9 +1109,11 @@ def test_domain(self, name: str, opt_class: Type[PersistentSolver]):
         res = opt.solve(m)
         self.assertAlmostEqual(res.best_feasible_objective, 0)
 
-    @parameterized.expand(input=mip_solvers)
-    def test_domain_with_integers(self, name: str, opt_class: Type[PersistentSolver]):
-        opt: PersistentSolver = opt_class()
+    @parameterized.expand(input=_load_tests(mip_solvers, only_child_vars_options))
+    def test_domain_with_integers(
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
+    ):
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
         m = pe.ConcreteModel()
@@ -1054,9 +1135,11 @@ def test_domain_with_integers(self, name: str, opt_class: Type[PersistentSolver]
         res = opt.solve(m)
         self.assertAlmostEqual(res.best_feasible_objective, 1)
 
-    @parameterized.expand(input=all_solvers)
-    def test_fixed_binaries(self, name: str, opt_class: Type[PersistentSolver]):
-        opt: PersistentSolver = opt_class()
+    @parameterized.expand(input=_load_tests(all_solvers, only_child_vars_options))
+    def test_fixed_binaries(
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
+    ):
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
         m = pe.ConcreteModel()
@@ -1071,7 +1154,7 @@ def test_fixed_binaries(self, name: str, opt_class: Type[PersistentSolver]):
         res = opt.solve(m)
         self.assertAlmostEqual(res.best_feasible_objective, 1)
 
-        opt: PersistentSolver = opt_class()
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         opt.update_config.treat_fixed_vars_as_params = False
         m.x.fix(0)
         res = opt.solve(m)
@@ -1080,9 +1163,11 @@ def test_fixed_binaries(self, name: str, opt_class: Type[PersistentSolver]):
         res = opt.solve(m)
         self.assertAlmostEqual(res.best_feasible_objective, 1)
 
-    @parameterized.expand(input=mip_solvers)
-    def test_with_gdp(self, name: str, opt_class: Type[PersistentSolver]):
-        opt: PersistentSolver = opt_class()
+    @parameterized.expand(input=_load_tests(mip_solvers, only_child_vars_options))
+    def test_with_gdp(
+        self, name: str, opt_class: Type[PersistentSolver], only_child_vars
+    ):
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
 
@@ -1104,7 +1189,7 @@ def test_with_gdp(self, name: str, opt_class: Type[PersistentSolver]):
         self.assertAlmostEqual(m.x.value, 0)
         self.assertAlmostEqual(m.y.value, 1)
 
-        opt: PersistentSolver = opt_class()
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         opt.use_extensions = True
         res = opt.solve(m)
         self.assertAlmostEqual(res.best_feasible_objective, 1)
@@ -1173,9 +1258,9 @@ def test_variables_elsewhere2(self, name: str, opt_class: Type[PersistentSolver]
         self.assertIn(m.y, sol)
         self.assertNotIn(m.z, sol)
 
-    @parameterized.expand(input=all_solvers)
-    def test_bug_1(self, name: str, opt_class: Type[PersistentSolver]):
-        opt: PersistentSolver = opt_class()
+    @parameterized.expand(input=_load_tests(all_solvers, only_child_vars_options))
+    def test_bug_1(self, name: str, opt_class: Type[PersistentSolver], only_child_vars):
+        opt: PersistentSolver = opt_class(only_child_vars=only_child_vars)
         if not opt.available():
             raise unittest.SkipTest
 
diff --git a/pyomo/contrib/appsi/writers/lp_writer.py b/pyomo/contrib/appsi/writers/lp_writer.py
index fe440261452..8a76fa5f9eb 100644
--- a/pyomo/contrib/appsi/writers/lp_writer.py
+++ b/pyomo/contrib/appsi/writers/lp_writer.py
@@ -16,7 +16,7 @@
 
 
 class LPWriter(PersistentBase):
-    def __init__(self, only_child_vars=True):
+    def __init__(self, only_child_vars=False):
         super(LPWriter, self).__init__(only_child_vars=only_child_vars)
         self._config = WriterConfig()
         self._writer = None
diff --git a/pyomo/contrib/appsi/writers/nl_writer.py b/pyomo/contrib/appsi/writers/nl_writer.py
index 0800f7bb2ad..9c739fd6ebb 100644
--- a/pyomo/contrib/appsi/writers/nl_writer.py
+++ b/pyomo/contrib/appsi/writers/nl_writer.py
@@ -19,7 +19,7 @@
 
 
 class NLWriter(PersistentBase):
-    def __init__(self, only_child_vars=True):
+    def __init__(self, only_child_vars=False):
         super(NLWriter, self).__init__(only_child_vars=only_child_vars)
         self._config = WriterConfig()
         self._writer = None