diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py
index a13e6785bf4..4dc6f42718b 100644
--- a/pyomo/contrib/doe/doe.py
+++ b/pyomo/contrib/doe/doe.py
@@ -281,7 +281,6 @@ def run_doe(self, model=None, results_file=None):
         results_file: string name of the file path to save the results
                       to in the form of a .json file
                       default: None --> don't save
-
         """
         # Check results file name
         if results_file is not None:
@@ -376,6 +375,11 @@ def run_doe(self, model=None, results_file=None):
             if self.objective_option == ObjectiveLib.trace:
                 trace_val = np.trace(np.linalg.pinv(self.get_FIM()))
                 model.obj_cons.egb_fim_block.outputs["A-opt"].set_value(trace_val)
+            elif self.objective_option == ObjectiveLib.pseudo_trace:
+                pseudo_trace_val = np.trace(np.array(self.get_FIM()))
+                model.obj_cons.egb_fim_block.outputs["pseudo-A-opt"].set_value(
+                    pseudo_trace_val
+                )
             elif self.objective_option == ObjectiveLib.determinant:
                 det_val = np.linalg.det(np.array(self.get_FIM()))
                 model.obj_cons.egb_fim_block.outputs["log-D-opt"].set_value(
@@ -389,63 +393,8 @@ def run_doe(self, model=None, results_file=None):
                 cond_number = np.log(np.abs(np.max(eig) / np.min(eig)))
                 model.obj_cons.egb_fim_block.outputs["ME-opt"].set_value(cond_number)
 
-        # If the model has L, initialize it with the solved FIM
-        if hasattr(model, "L"):
-            # Get the FIM values
-            fim_vals = [
-                pyo.value(model.fim[i, j])
-                for i in model.parameter_names
-                for j in model.parameter_names
-            ]
-            fim_np = np.array(fim_vals).reshape(
-                (len(model.parameter_names), len(model.parameter_names))
-            )
-
-            # Need to compute the full FIM before
-            # initializing the Cholesky factorization
-            if self.only_compute_fim_lower:
-                fim_np = fim_np + fim_np.T - np.diag(np.diag(fim_np))
-
-            # Check if the FIM is positive definite
-            # If not, add jitter to the diagonal
-            # to ensure positive definiteness
-            min_eig = np.min(np.linalg.eigvals(fim_np))
-
-            if min_eig < _SMALL_TOLERANCE_DEFINITENESS:
-                # Raise the minimum eigenvalue to at
-                # least _SMALL_TOLERANCE_DEFINITENESS
-                jitter = np.min(
-                    [
-                        -min_eig + _SMALL_TOLERANCE_DEFINITENESS,
-                        _SMALL_TOLERANCE_DEFINITENESS,
-                    ]
-                )
-            else:
-                # No jitter needed
-                jitter = 0
-
-            # Add jitter to the diagonal to ensure positive definiteness
-            L_vals_sq = np.linalg.cholesky(
-                fim_np + jitter * np.eye(len(model.parameter_names))
-            )
-            for i, c in enumerate(model.parameter_names):
-                for j, d in enumerate(model.parameter_names):
-                    model.L[c, d].value = L_vals_sq[i, j]
-
-            # Initialize the inverse of L if it exists
-            if hasattr(model, "L_inv"):
-                L_inv_vals = np.linalg.inv(L_vals_sq)
-
-                for i, c in enumerate(model.parameter_names):
-                    for j, d in enumerate(model.parameter_names):
-                        if i >= j:
-                            model.L_inv[c, d].value = L_inv_vals[i, j]
-                        else:
-                            model.L_inv[c, d].value = 0.0
-                # Initialize the cov_trace if it exists
-                if hasattr(model, "cov_trace"):
-                    initial_cov_trace = np.sum(L_inv_vals**2)
-                    model.cov_trace.value = initial_cov_trace
+        # Keep Cholesky-related variables synchronized with current FIM values
+        self._initialize_cholesky_from_fim(model=model)
 
         if hasattr(model, "determinant"):
             model.determinant.value = np.linalg.det(np.array(self.get_FIM()))
@@ -537,6 +486,109 @@ def run_doe(self, model=None, results_file=None):
             with open(results_file, "w") as file:
                 json.dump(self.results, file)
 
+    def _compute_cholesky_jitter(self, min_eig):
+        """
+        Compute diagonal regularization for Cholesky initialization.
+
+        Parameters
+        ----------
+        min_eig: float
+            Minimum eigenvalue of the current FIM estimate.
+
+        Returns
+        -------
+        float
+            Nonnegative diagonal shift needed so the minimum eigenvalue
+            is at least ``_SMALL_TOLERANCE_DEFINITENESS``.
+        """
+        return max(0.0, _SMALL_TOLERANCE_DEFINITENESS - float(min_eig))
+
+    def _get_fim_numpy(self, model):
+        """
+        Assemble the current FIM variable values into a NumPy array.
+
+        Parameters
+        ----------
+        model: ConcreteModel
+            DoE model containing variable ``fim``.
+
+        Returns
+        -------
+        ndarray
+            Dense FIM array. If ``only_compute_fim_lower`` is True, the
+            returned array is symmetrized from the lower triangle.
+        """
+        fim_vals = [
+            pyo.value(model.fim[i, j])
+            for i in model.parameter_names
+            for j in model.parameter_names
+        ]
+        fim_np = np.array(fim_vals, dtype=float).reshape(
+            (len(model.parameter_names), len(model.parameter_names))
+        )
+        if self.only_compute_fim_lower:
+            fim_np = fim_np + fim_np.T - np.diag(np.diag(fim_np))
+        return fim_np
+
+    def _initialize_cholesky_from_fim(self, model=None):
+        """
+        Synchronize Cholesky-related variables using the current FIM.
+
+        Parameters
+        ----------
+        model: ConcreteModel, optional
+            DoE model to update. Defaults to ``self.model``.
+
+        Returns
+        -------
+        None
+            Updates model values in place for available variables:
+            ``L``, ``L_inv``, ``fim_inv``, and ``cov_trace``.
+        """
+        if model is None:
+            model = self.model
+        if not hasattr(model, "L"):
+            return
+
+        fim_np = self._get_fim_numpy(model)
+        min_eig = float(np.min(np.linalg.eigvalsh(fim_np)))
+        jitter = self._compute_cholesky_jitter(min_eig)
+        fim_pd = fim_np + jitter * np.eye(len(model.parameter_names))
+
+        L_vals = np.linalg.cholesky(fim_pd)
+        for i, c in enumerate(model.parameter_names):
+            for j, d in enumerate(model.parameter_names):
+                if i >= j:
+                    model.L[c, d].value = L_vals[i, j]
+                else:
+                    model.L[c, d].value = 0.0
+
+        if hasattr(model, "L_inv"):
+            L_inv_vals = np.linalg.inv(L_vals)
+            for i, c in enumerate(model.parameter_names):
+                for j, d in enumerate(model.parameter_names):
+                    if i >= j:
+                        model.L_inv[c, d].value = L_inv_vals[i, j]
+                    else:
+                        model.L_inv[c, d].value = 0.0
+
+        if hasattr(model, "fim_inv"):
+            # Use the pseudo-inverse here rather than the strict inverse.
+            # The jittered matrix should be positive definite, but ``pinv``
+            # is safer for borderline ill-conditioned cases and matches the
+            # defensive approach already used when initializing ``fim_inv``
+            # from user-provided starting values.
+            fim_inv_vals = np.linalg.pinv(fim_pd)
+            for i, c in enumerate(model.parameter_names):
+                for j, d in enumerate(model.parameter_names):
+                    if self.only_compute_fim_lower and i < j:
+                        model.fim_inv[c, d].value = 0.0
+                    else:
+                        model.fim_inv[c, d].value = fim_inv_vals[i, j]
+
+        if hasattr(model, "cov_trace"):
+            model.cov_trace.value = np.trace(fim_inv_vals)
+
     # Perform multi-experiment doe (sequential, or ``greedy`` approach)
     def run_multi_doe_sequential(self, N_exp=1):
         raise NotImplementedError("Multiple experiment optimization not yet supported.")
@@ -898,6 +950,7 @@ def create_doe_model(self, model=None):
             self.only_compute_fim_lower
             and self.objective_option == ObjectiveLib.determinant
             and not self.Cholesky_option
+            and not self.use_grey_box
         ):
             raise ValueError(
                 "Cannot compute determinant with explicit formula "
@@ -1663,6 +1716,11 @@ def FIM_egb_cons(m, p1, p2):
             model.objective = pyo.Objective(
                 expr=model.obj_cons.egb_fim_block.outputs["A-opt"], sense=pyo.minimize
             )
+        elif self.objective_option == ObjectiveLib.pseudo_trace:
+            model.objective = pyo.Objective(
+                expr=model.obj_cons.egb_fim_block.outputs["pseudo-A-opt"],
+                sense=pyo.maximize,
+            )
         elif self.objective_option == ObjectiveLib.determinant:
             model.objective = pyo.Objective(
                 expr=model.obj_cons.egb_fim_block.outputs["log-D-opt"],
diff --git a/pyomo/contrib/doe/grey_box_utilities.py b/pyomo/contrib/doe/grey_box_utilities.py
index 41e5cf2d78c..645cd46f04c 100644
--- a/pyomo/contrib/doe/grey_box_utilities.py
+++ b/pyomo/contrib/doe/grey_box_utilities.py
@@ -108,9 +108,10 @@ def __init__(self, doe_object, objective_option="determinant", logger_level=None
     def _get_FIM(self):
         # Grabs the current FIM subject
         # to the input values.
-        # This function currently assumes
-        # that we use a lower triangular
-        # FIM.
+        # Inputs store one triangular half
+        # of a symmetric FIM. Reconstruct
+        # the full symmetric matrix here,
+        # consistent with manuscript S5.
         upt_FIM = self._input_values
 
         # Create FIM in the correct way
@@ -164,6 +165,8 @@ def output_names(self):
 
         if self.objective_option == ObjectiveLib.trace:
             obj_name = "A-opt"
+        elif self.objective_option == ObjectiveLib.pseudo_trace:
+            obj_name = "pseudo-A-opt"
         elif self.objective_option == ObjectiveLib.determinant:
             obj_name = "log-D-opt"
         elif self.objective_option == ObjectiveLib.minimum_eigenvalue:
@@ -196,6 +199,8 @@ def evaluate_outputs(self):
 
         if self.objective_option == ObjectiveLib.trace:
             obj_value = np.trace(np.linalg.pinv(M))
+        elif self.objective_option == ObjectiveLib.pseudo_trace:
+            obj_value = np.trace(M)
         elif self.objective_option == ObjectiveLib.determinant:
             sign, logdet = np.linalg.slogdet(M)
             obj_value = logdet
@@ -231,6 +236,8 @@ def finalize_block_construction(self, pyomo_block):
         # objective function.
         if self.objective_option == ObjectiveLib.trace:
             pyomo_block.outputs["A-opt"] = output_value
+        elif self.objective_option == ObjectiveLib.pseudo_trace:
+            pyomo_block.outputs["pseudo-A-opt"] = output_value
         elif self.objective_option == ObjectiveLib.determinant:
             pyomo_block.outputs["log-D-opt"] = output_value
         elif self.objective_option == ObjectiveLib.minimum_eigenvalue:
@@ -268,6 +275,8 @@ def evaluate_jacobian_outputs(self):
             # is -inv(FIM) @ inv(FIM). Add reference to
             # pyomo.DoE 2.0 manuscript S.I.
             jac_M = -Minv @ Minv
+        elif self.objective_option == ObjectiveLib.pseudo_trace:
+            jac_M = np.eye(self._n_params, dtype=np.float64)
         elif self.objective_option == ObjectiveLib.determinant:
             Minv = np.linalg.pinv(M)
             # Derivative formula derived using tensor
diff --git a/pyomo/contrib/doe/tests/test_doe_errors.py b/pyomo/contrib/doe/tests/test_doe_errors.py
index 10165aa6267..3f8b538d51e 100644
--- a/pyomo/contrib/doe/tests/test_doe_errors.py
+++ b/pyomo/contrib/doe/tests/test_doe_errors.py
@@ -781,6 +781,25 @@ def test_invalid_trace_without_cholesky(self):
         ):
             doe_obj.create_objective_function()
 
+    def test_invalid_determinant_without_cholesky(self):
+        fd_method = "central"
+        obj_used = "determinant"
+
+        experiment = get_rooney_biegler_experiment_flag()
+
+        DoE_args = get_standard_args(experiment, fd_method, obj_used, flag=None)
+        DoE_args["_Cholesky_option"] = False
+
+        doe_obj = DesignOfExperiments(**DoE_args)
+
+        # The explicit determinant formulation needs the full FIM, so we keep
+        # this regression test focused on the early validation guard.
+        with self.assertRaisesRegex(
+            ValueError,
+            "Cannot compute determinant with explicit formula if only_compute_fim_lower is True.",
+        ):
+            doe_obj.create_doe_model()
+
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/pyomo/contrib/doe/tests/test_doe_initialization.py b/pyomo/contrib/doe/tests/test_doe_initialization.py
new file mode 100644
index 00000000000..c0a3fa35f6d
--- /dev/null
+++ b/pyomo/contrib/doe/tests/test_doe_initialization.py
@@ -0,0 +1,194 @@
+# ____________________________________________________________________________________
+#
+# Pyomo: Python Optimization Modeling Objects
+# Copyright (c) 2008-2026 National Technology and Engineering Solutions of Sandia, LLC
+# Under the terms of Contract DE-NA0003525 with National Technology and Engineering
+# Solutions of Sandia, LLC, the U.S. Government retains certain rights in this
+# software.  This software is distributed under the 3-clause BSD License.
+# ____________________________________________________________________________________
+
+from pyomo.common.dependencies import numpy as np, numpy_available
+import pyomo.common.unittest as unittest
+import pyomo.environ as pyo
+
+from pyomo.contrib.doe import DesignOfExperiments
+from pyomo.contrib.doe.doe import _SMALL_TOLERANCE_DEFINITENESS
+
+
+class _FakeResult:
+    """Minimal fake solver result object matching attributes used by DoE."""
+
+    class _SolverData:
+        status = "ok"
+        termination_condition = "optimal"
+        message = "fake solve"
+
+    solver = _SolverData()
+
+
+class _MutatingRecordingSolver:
+    """
+    Fake solver that perturbs FIM during dummy-objective initialization solve.
+
+    The key behavior we want to exercise is not optimization success, but the
+    fact that the square initialization solve can change the FIM after the
+    object has already created Cholesky-related variables. That is exactly the
+    bug this regression test protects.
+    """
+
+    def __init__(self):
+        self.options = {}
+
+    def solve(self, model, tee=False):
+        """Mutate 2x2 FIM values during init-stage solve."""
+        if (
+            hasattr(model, "dummy_obj")
+            and model.dummy_obj.active
+            and hasattr(model, "fim")
+            and len(list(model.parameter_names)) == 2
+        ):
+            p1, p2 = list(model.parameter_names)
+            model.fim[p1, p1].set_value(16.0)
+            model.fim[p2, p1].set_value(4.0)
+            if (p1, p2) in model.fim:
+                model.fim[p1, p2].set_value(0.0)
+            model.fim[p2, p2].set_value(9.0)
+        return _FakeResult()
+
+
+class _TwoParamExperiment:
+    """
+    Minimal two-parameter experiment fixture for trace/Cholesky tests.
+
+    Two parameters is the smallest useful case here because it creates a real
+    matrix-valued FIM while still keeping the expected values easy to reason
+    about by hand. That makes it a good target for the initialization
+    synchronization regression.
+    """
+
+    def get_labeled_model(self):
+        m = pyo.ConcreteModel()
+        m.x1 = pyo.Var(initialize=1.0)
+        m.x2 = pyo.Var(initialize=1.0)
+        m.theta1 = pyo.Var(initialize=2.0)
+        m.theta2 = pyo.Var(initialize=3.0)
+        m.y1 = pyo.Var(initialize=2.0)
+        m.y2 = pyo.Var(initialize=3.0)
+        m.eq1 = pyo.Constraint(expr=m.y1 == m.theta1 * m.x1)
+        m.eq2 = pyo.Constraint(expr=m.y2 == m.theta2 * m.x2)
+
+        m.experiment_inputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+        m.experiment_inputs[m.x1] = 1.0
+        m.experiment_inputs[m.x2] = 1.0
+        m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+        m.unknown_parameters[m.theta1] = 2.0
+        m.unknown_parameters[m.theta2] = 3.0
+        m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+        m.experiment_outputs[m.y1] = 2.0
+        m.experiment_outputs[m.y2] = 3.0
+        m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+        m.measurement_error[m.y1] = 1.0
+        m.measurement_error[m.y2] = 1.0
+        return m
+
+
+def _make_trace_doe_object():
+    """
+    Build a trace-objective DoE object with a FIM-mutating fake solver.
+
+    The trace objective is the path that introduces ``L``, ``L_inv``,
+    ``fim_inv``, and ``cov_trace``, so it is the right place to check that the
+    helper re-synchronizes all of those values after the init solve.
+    """
+    return DesignOfExperiments(
+        experiment=_TwoParamExperiment(),
+        fd_formula="central",
+        step=1e-3,
+        objective_option="trace",
+        solver=_MutatingRecordingSolver(),
+        tee=False,
+    )
+
+
+@unittest.skipIf(not numpy_available, "Pyomo.DoE needs numpy to run tests")
+class TestCholeskyInitialization(unittest.TestCase):
+    """Regression tests for DoE initialization and Cholesky sync."""
+
+    def test_compute_cholesky_jitter_raises_negative_eigenvalue(self):
+        """
+        Negative minimum eigenvalue should produce positive corrective jitter.
+
+        This is a direct unit test for the helper used by the initialization
+        sync path, so we can catch arithmetic regressions independently of the
+        solver/modeled example below.
+        """
+        doe_obj = _make_trace_doe_object()
+        min_eig = -1.0e-3
+        jitter = doe_obj._compute_cholesky_jitter(min_eig)
+        self.assertAlmostEqual(
+            jitter, _SMALL_TOLERANCE_DEFINITENESS - min_eig, places=14
+        )
+
+    def test_compute_cholesky_jitter_zero_when_not_needed(self):
+        """
+        Positive minimum eigenvalue above tolerance should yield zero jitter.
+
+        We keep this separate from the negative-eigenvalue case because it
+        verifies the helper does not add unnecessary regularization when the
+        FIM is already well-conditioned.
+        """
+        doe_obj = _make_trace_doe_object()
+        min_eig = 1.0e-2
+        jitter = doe_obj._compute_cholesky_jitter(min_eig)
+        self.assertEqual(jitter, 0.0)
+
+    def test_trace_initialization_resynchronizes_fim_inverse_variables(self):
+        """
+        Verify trace-mode initialization re-synchronizes inverse-related vars.
+
+        The fake solver mutates the FIM during the square solve. After that
+        solve, the DoE code should rebuild the dependent quantities from the
+        updated FIM so that:
+        - ``L`` matches the Cholesky factor of the current FIM,
+        - ``L_inv`` is the inverse of that Cholesky factor,
+        - ``fim_inv`` matches the inverse FIM, and
+        - ``cov_trace`` matches the trace of the inverse FIM.
+        """
+        doe_obj = _make_trace_doe_object()
+        doe_obj.run_doe()
+
+        model = doe_obj.model
+        params = list(model.parameter_names)
+        fim = np.array(
+            [[pyo.value(model.fim[i, j]) for j in params] for i in params], dtype=float
+        )
+        L = np.array([[pyo.value(model.L[i, j]) for j in params] for i in params])
+        L_inv = np.array(
+            [[pyo.value(model.L_inv[i, j]) for j in params] for i in params]
+        )
+        fim_inv = np.array(
+            [[pyo.value(model.fim_inv[i, j]) for j in params] for i in params]
+        )
+
+        expected_fim = np.array([[16.0, 4.0], [4.0, 9.0]])
+        expected_fim_inv = np.linalg.inv(expected_fim)
+        expected_L = np.linalg.cholesky(expected_fim)
+        expected_L_inv = np.linalg.inv(expected_L)
+        fim_sym = fim + fim.T - np.diag(np.diag(fim))
+        fim_inv_sym = fim_inv + fim_inv.T - np.diag(np.diag(fim_inv))
+
+        # The model stores only the lower triangle when only_compute_fim_lower=True.
+        # Reconstruct the symmetric FIM before checking the initialization math,
+        # because Cholesky and inverse quantities are defined on the full matrix.
+        self.assertTrue(np.allclose(fim_sym, expected_fim))
+        self.assertTrue(np.allclose(L, expected_L))
+        self.assertTrue(np.allclose(L @ L_inv, np.eye(len(params)), atol=1e-8))
+        self.assertTrue(np.allclose(fim_inv_sym, expected_fim_inv))
+        self.assertTrue(np.allclose(L_inv, expected_L_inv))
+        self.assertAlmostEqual(
+            pyo.value(model.cov_trace), float(np.trace(expected_fim_inv)), places=10
+        )
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/pyomo/contrib/doe/tests/test_greybox.py b/pyomo/contrib/doe/tests/test_greybox.py
index 2fd007e6d52..a0064493fa7 100644
--- a/pyomo/contrib/doe/tests/test_greybox.py
+++ b/pyomo/contrib/doe/tests/test_greybox.py
@@ -69,6 +69,15 @@
 
 
 def get_numerical_derivative(grey_box_object=None):
+    """
+    Compute finite-difference first derivatives of the selected DoE objective.
+
+    Parameters
+    ----------
+    grey_box_object : FIMExternalGreyBox
+        Grey-box objective wrapper that provides the current FIM and objective
+        type (trace, determinant, minimum-eigenvalue, or condition-number).
+    """
     # Internal import to avoid circular imports
     from pyomo.contrib.doe import ObjectiveLib
 
@@ -120,6 +129,19 @@ def get_numerical_derivative(grey_box_object=None):
 
 
 def get_numerical_second_derivative(grey_box_object=None, return_reduced=True):
+    """
+    Compute finite-difference second derivatives for Hessian validation.
+
+    Parameters
+    ----------
+    grey_box_object : FIMExternalGreyBox
+        Grey-box objective wrapper that provides the current FIM and objective
+        type.
+    return_reduced : bool, optional
+        If True, return the 10x10 reduced upper-triangular Hessian form used by
+        the external grey-box interface for a 4-parameter case; otherwise
+        return the full 4-D derivative tensor.
+    """
     # Internal import to avoid circular imports
     from pyomo.contrib.doe import ObjectiveLib
 
@@ -243,6 +265,18 @@ def get_numerical_second_derivative(grey_box_object=None, return_reduced=True):
 
 
 def get_standard_args(experiment, fd_method, obj_used):
+    """
+    Build common DesignOfExperiments keyword arguments used in this test module.
+
+    Parameters
+    ----------
+    experiment : object
+        Experiment fixture implementing ``get_labeled_model``.
+    fd_method : str
+        Finite-difference mode passed to DoE (for example, ``"central"``).
+    obj_used : str
+        Objective option name (for example, ``"trace"`` or ``"determinant"``).
+    """
     args = {}
     args['experiment'] = experiment
     args['fd_formula'] = fd_method
@@ -279,6 +313,14 @@ def get_standard_args(experiment, fd_method, obj_used):
 
 
 def make_greybox_and_doe_objects(objective_option):
+    """
+    Create reactor-based DoE and grey-box objects for derivative/build tests.
+
+    Parameters
+    ----------
+    objective_option : str
+        Objective option to validate with the grey-box wrapper.
+    """
     fd_method = "central"
     obj_used = objective_option
 
@@ -299,6 +341,14 @@ def make_greybox_and_doe_objects(objective_option):
 
 
 def make_greybox_and_doe_objects_rooney_biegler(objective_option):
+    """
+    Create Rooney-Biegler DoE/grey-box objects with a data-driven prior FIM.
+
+    Parameters
+    ----------
+    objective_option : str
+        Objective option to validate and/or solve in grey-box mode.
+    """
     fd_method = "central"
     obj_used = objective_option
 
@@ -389,6 +439,7 @@ class TestFIMExternalGreyBox(unittest.TestCase):
     # set the inputs for the
     # Grey Box object
     def test_set_inputs(self):
+        """Verify packed upper-triangular inputs reconstruct the full FIM."""
         objective_option = "trace"
         doe_obj, grey_box_object = make_greybox_and_doe_objects(
             objective_option=objective_option
@@ -405,6 +456,7 @@ def test_set_inputs(self):
     # Testing that getting the
     # input names works properly
     def test_input_names(self):
+        """Verify the grey-box input ordering matches the expected FIM entries."""
         objective_option = "trace"
         doe_obj, grey_box_object = make_greybox_and_doe_objects(
             objective_option=objective_option
@@ -432,6 +484,7 @@ def test_input_names(self):
     # Testing that getting the
     # output names works properly
     def test_output_names(self):
+        """Verify output labels for all supported grey-box objective variants."""
         # Need to test for each objective type
         # A-opt
         objective_option = "trace"
@@ -439,6 +492,12 @@ def test_output_names(self):
             objective_option=objective_option
         )
 
+        # Pseudo A-opt
+        objective_option = "pseudo_trace"
+        doe_obj_P, grey_box_object_P = make_greybox_and_doe_objects(
+            objective_option=objective_option
+        )
+
         # D-opt
         objective_option = "determinant"
         doe_obj_D, grey_box_object_D = make_greybox_and_doe_objects(
@@ -460,11 +519,12 @@ def test_output_names(self):
         # Hard-coded names of the outputs
         # There is one element per
         # objective type
-        output_names = ['A-opt', 'log-D-opt', 'E-opt', 'ME-opt']
+        output_names = ['A-opt', 'pseudo-A-opt', 'log-D-opt', 'E-opt', 'ME-opt']
 
         # Grabbing input names from grey box object
         output_names_gb = []
         output_names_gb.extend(grey_box_object_A.output_names())
+        output_names_gb.extend(grey_box_object_P.output_names())
         output_names_gb.extend(grey_box_object_D.output_names())
         output_names_gb.extend(grey_box_object_E.output_names())
         output_names_gb.extend(grey_box_object_ME.output_names())
@@ -473,6 +533,7 @@ def test_output_names(self):
 
     # Testing output computation
     def test_outputs_A_opt(self):
+        """Check A-opt output equals trace of inverse FIM."""
         objective_option = "trace"
         doe_obj, grey_box_object = make_greybox_and_doe_objects(
             objective_option=objective_option
@@ -487,7 +548,24 @@ def test_outputs_A_opt(self):
 
         self.assertTrue(np.isclose(grey_box_A_opt, A_opt))
 
+    def test_outputs_pseudo_A_opt(self):
+        """Check pseudo-A-opt output equals trace of the FIM itself."""
+        objective_option = "pseudo_trace"
+        doe_obj, grey_box_object = make_greybox_and_doe_objects(
+            objective_option=objective_option
+        )
+
+        # Set input values to the random testing matrix
+        grey_box_object.set_input_values(testing_matrix[masking_matrix > 0])
+
+        grey_box_pseudo_A_opt = grey_box_object.evaluate_outputs()
+
+        pseudo_A_opt = np.trace(testing_matrix)
+
+        self.assertTrue(np.isclose(grey_box_pseudo_A_opt, pseudo_A_opt))
+
     def test_outputs_D_opt(self):
+        """Check D-opt output equals log-determinant of FIM."""
         objective_option = "determinant"
         doe_obj, grey_box_object = make_greybox_and_doe_objects(
             objective_option=objective_option
@@ -503,6 +581,7 @@ def test_outputs_D_opt(self):
         self.assertTrue(np.isclose(grey_box_D_opt, D_opt))
 
     def test_outputs_E_opt(self):
+        """Check E-opt output equals minimum eigenvalue of FIM."""
         objective_option = "minimum_eigenvalue"
         doe_obj, grey_box_object = make_greybox_and_doe_objects(
             objective_option=objective_option
@@ -519,6 +598,7 @@ def test_outputs_E_opt(self):
         self.assertTrue(np.isclose(grey_box_E_opt, E_opt))
 
     def test_outputs_ME_opt(self):
+        """Check ME-opt output equals log condition number of FIM."""
         objective_option = "condition_number"
         doe_obj, grey_box_object = make_greybox_and_doe_objects(
             objective_option=objective_option
@@ -536,6 +616,7 @@ def test_outputs_ME_opt(self):
 
     # Testing Jacobian computation
     def test_jacobian_A_opt(self):
+        """Compare A-opt Jacobian against finite-difference derivatives."""
         objective_option = "trace"
         doe_obj, grey_box_object = make_greybox_and_doe_objects(
             objective_option=objective_option
@@ -559,7 +640,31 @@ def test_jacobian_A_opt(self):
         # assert that each component is close
         self.assertTrue(np.all(np.isclose(jac, jac_FD, rtol=1e-4, atol=1e-4)))
 
+    def test_jacobian_pseudo_A_opt(self):
+        """Pseudo-A-opt is linear, so its Jacobian should be the identity."""
+        objective_option = "pseudo_trace"
+        doe_obj, grey_box_object = make_greybox_and_doe_objects(
+            objective_option=objective_option
+        )
+
+        # Set input values to the random testing matrix
+        grey_box_object.set_input_values(testing_matrix[masking_matrix > 0])
+
+        # Grab the Jacobian values
+        utri_vals_jac = grey_box_object.evaluate_jacobian_outputs().toarray()
+
+        # Recover the Jacobian in Matrix Form
+        jac = np.zeros_like(grey_box_object._get_FIM())
+        jac[np.triu_indices_from(jac)] = utri_vals_jac
+        jac += jac.transpose()
+        jac = jac / 2
+
+        # For pseudo-trace, the objective is trace(FIM), so the full Jacobian
+        # is exactly the identity matrix.
+        self.assertTrue(np.allclose(jac, np.eye(jac.shape[0])))
+
     def test_jacobian_D_opt(self):
+        """Compare D-opt Jacobian against finite-difference derivatives."""
         objective_option = "determinant"
         doe_obj, grey_box_object = make_greybox_and_doe_objects(
             objective_option=objective_option
@@ -584,6 +689,7 @@ def test_jacobian_D_opt(self):
         self.assertTrue(np.all(np.isclose(jac, jac_FD, rtol=1e-4, atol=1e-4)))
 
     def test_jacobian_E_opt(self):
+        """Compare E-opt Jacobian against finite-difference derivatives."""
         objective_option = "minimum_eigenvalue"
         doe_obj, grey_box_object = make_greybox_and_doe_objects(
             objective_option=objective_option
@@ -608,6 +714,7 @@ def test_jacobian_E_opt(self):
         self.assertTrue(np.all(np.isclose(jac, jac_FD, rtol=1e-4, atol=1e-4)))
 
     def test_jacobian_ME_opt(self):
+        """Compare ME-opt Jacobian against finite-difference derivatives."""
         objective_option = "condition_number"
         doe_obj, grey_box_object = make_greybox_and_doe_objects(
             objective_option=objective_option
@@ -631,8 +738,29 @@ def test_jacobian_ME_opt(self):
         # assert that each component is close
         self.assertTrue(np.all(np.isclose(jac, jac_FD, rtol=1e-4, atol=1e-4)))
 
+    def test_hessian_pseudo_A_opt(self):
+        """Pseudo-A-opt is linear, so its Hessian should be identically zero."""
+        objective_option = "pseudo_trace"
+        doe_obj, grey_box_object = make_greybox_and_doe_objects(
+            objective_option=objective_option
+        )
+
+        # Set input values to the random testing matrix
+        grey_box_object.set_input_values(testing_matrix[masking_matrix > 0])
+
+        # Grab the Hessian values
+        hess_vals_from_gb = grey_box_object.evaluate_hessian_outputs().toarray()
+
+        # Recover the Hessian in Matrix Form
+        hess_gb = hess_vals_from_gb
+        hess_gb += hess_gb.transpose() - np.diag(np.diag(hess_gb))
+
+        # Pseudo-trace is linear in the FIM, so the second derivative is zero.
+        self.assertTrue(np.allclose(hess_gb, np.zeros_like(hess_gb)))
+
     # Testing Hessian Computation
     def test_hessian_A_opt(self):
+        """Compare A-opt Hessian against finite-difference second derivatives."""
         objective_option = "trace"
         doe_obj, grey_box_object = make_greybox_and_doe_objects(
             objective_option=objective_option
@@ -654,7 +782,31 @@ def test_hessian_A_opt(self):
         # assert that each component is close
         self.assertTrue(np.all(np.isclose(hess_gb, hess_FD, rtol=1e-4, atol=1e-4)))
 
+    def test_pseudo_A_opt_greybox_build(self):
+        """Validate pseudo-A-opt grey-box wiring and initialized values."""
+        objective_option = "pseudo_trace"
+        doe_obj, grey_box_object = make_greybox_and_doe_objects(
+            objective_option=objective_option
+        )
+
+        # Build the greybox objective block on the DoE object
+        doe_obj.create_grey_box_objective_function()
+
+        # Pseudo-trace uses the trace of the FIM itself, so the initial value
+        # should match the current symmetric FIM rather than its inverse.
+        pseudo_A_opt_val = np.trace(testing_matrix + np.eye(4))
+
+        self.assertTrue(hasattr(doe_obj.model.obj_cons, "egb_fim_block"))
+        self.assertIn("pseudo-A-opt", grey_box_object.output_names())
+        self.assertIsInstance(doe_obj.model.objective, pyo.Objective)
+        self.assertEqual(doe_obj.model.objective.sense, pyo.maximize)
+        self.assertAlmostEqual(
+            doe_obj.model.obj_cons.egb_fim_block.outputs["pseudo-A-opt"].value,
+            pseudo_A_opt_val,
+        )
+
     def test_hessian_D_opt(self):
+        """Compare D-opt Hessian against finite-difference second derivatives."""
         objective_option = "determinant"
         doe_obj, grey_box_object = make_greybox_and_doe_objects(
             objective_option=objective_option
@@ -677,6 +829,7 @@ def test_hessian_D_opt(self):
         self.assertTrue(np.all(np.isclose(hess_gb, hess_FD, rtol=1e-4, atol=1e-4)))
 
     def test_hessian_E_opt(self):
+        """Compare E-opt Hessian against finite-difference second derivatives."""
         objective_option = "minimum_eigenvalue"
         doe_obj, grey_box_object = make_greybox_and_doe_objects(
             objective_option=objective_option
@@ -699,6 +852,7 @@ def test_hessian_E_opt(self):
         self.assertTrue(np.all(np.isclose(hess_gb, hess_FD, rtol=1e-4, atol=1e-4)))
 
     def test_hessian_ME_opt(self):
+        """Compare ME-opt Hessian against finite-difference second derivatives."""
         objective_option = "condition_number"
         doe_obj, grey_box_object = make_greybox_and_doe_objects(
             objective_option=objective_option
@@ -721,6 +875,7 @@ def test_hessian_ME_opt(self):
         self.assertTrue(np.all(np.isclose(hess_gb, hess_FD, rtol=1e-4, atol=1e-4)))
 
     def test_equality_constraint_names(self):
+        """Confirm the FIM grey-box objective exposes no equality constraints."""
         objective_option = "condition_number"
         doe_obj, grey_box_object = make_greybox_and_doe_objects(
             objective_option=objective_option
@@ -734,6 +889,7 @@ def test_equality_constraint_names(self):
         self.assertListEqual(eq_con_names_gb, [])
 
     def test_evaluate_equality_constraints(self):
+        """Confirm equality-constraint values are ``None`` for this interface."""
         objective_option = "condition_number"
         doe_obj, grey_box_object = make_greybox_and_doe_objects(
             objective_option=objective_option
@@ -747,6 +903,7 @@ def test_evaluate_equality_constraints(self):
         self.assertIsNone(eq_con_vals_gb)
 
     def test_evaluate_jacobian_equality_constraints(self):
+        """Confirm equality-constraint Jacobian is ``None`` for this interface."""
         objective_option = "condition_number"
         doe_obj, grey_box_object = make_greybox_and_doe_objects(
             objective_option=objective_option
@@ -761,6 +918,7 @@ def test_evaluate_jacobian_equality_constraints(self):
         self.assertIsNone(jac_eq_con_vals_gb)
 
     def test_evaluate_hessian_equality_constraints(self):
+        """Confirm equality-constraint Hessian is ``None`` for this interface."""
         objective_option = "condition_number"
         doe_obj, grey_box_object = make_greybox_and_doe_objects(
             objective_option=objective_option
@@ -778,6 +936,7 @@ def test_evaluate_hessian_equality_constraints(self):
     # the DoE problem with grey box is built
     # properly.
     def test_A_opt_greybox_build(self):
+        """Validate A-opt grey-box block wiring and initialized values on DoE model."""
         objective_option = "trace"
         doe_obj, _ = make_greybox_and_doe_objects(objective_option=objective_option)
 
@@ -823,6 +982,7 @@ def test_A_opt_greybox_build(self):
         self.assertTrue(np.all(np.isclose(current_FIM, testing_matrix + np.eye(4))))
 
     def test_D_opt_greybox_build(self):
+        """Validate D-opt grey-box block wiring and initialized values on DoE model."""
         objective_option = "determinant"
         doe_obj, _ = make_greybox_and_doe_objects(objective_option=objective_option)
 
@@ -869,7 +1029,28 @@ def test_D_opt_greybox_build(self):
 
         self.assertTrue(np.all(np.isclose(current_FIM, testing_matrix + np.eye(4))))
 
+    def test_D_opt_greybox_build_with_triangular_fim(self):
+        """Determinant grey-box build accepts triangular FIM storage."""
+        objective_option = "determinant"
+        doe_obj, grey_box_object = make_greybox_and_doe_objects(
+            objective_option=objective_option
+        )
+
+        self.assertTrue(doe_obj.only_compute_fim_lower)
+        # Grey-box objective construction should not depend on the Cholesky
+        # build path. This fixture may still carry the shared default flag, so
+        # we do not treat it as a grey-box invariant here.
+        self.assertTrue(doe_obj.use_grey_box)
+
+        doe_obj.create_grey_box_objective_function()
+
+        self.assertTrue(hasattr(doe_obj.model.obj_cons, "egb_fim_block"))
+        self.assertIn("log-D-opt", list(grey_box_object.output_names()))
+        self.assertIsInstance(doe_obj.model.objective, pyo.Objective)
+        self.assertEqual(doe_obj.model.objective.sense, pyo.maximize)
+
     def test_E_opt_greybox_build(self):
+        """Validate E-opt grey-box block wiring and initialized values on DoE model."""
         objective_option = "minimum_eigenvalue"
         doe_obj, _ = make_greybox_and_doe_objects(objective_option=objective_option)
 
@@ -916,6 +1097,7 @@ def test_E_opt_greybox_build(self):
         self.assertTrue(np.all(np.isclose(current_FIM, testing_matrix + np.eye(4))))
 
     def test_ME_opt_greybox_build(self):
+        """Validate ME-opt grey-box block wiring and initialized values on DoE model."""
         objective_option = "condition_number"
         doe_obj, _ = make_greybox_and_doe_objects(objective_option=objective_option)
 
@@ -963,6 +1145,7 @@ def test_ME_opt_greybox_build(self):
 
     # Testing all the error messages
     def test_constructor_doe_object_error(self):
+        """Ensure constructor rejects missing DoE object dependencies."""
         with self.assertRaisesRegex(
             ValueError,
             "DoE Object must be provided to build external grey box of the FIM.",
@@ -970,6 +1153,7 @@ def test_constructor_doe_object_error(self):
             grey_box_object = FIMExternalGreyBox(doe_object=None)
 
     def test_constructor_objective_lib_error(self):
+        """Ensure constructor rejects unsupported objective option strings."""
         objective_option = "trace"
         doe_object, grey_box_object = make_greybox_and_doe_objects(
             objective_option=objective_option
@@ -982,6 +1166,7 @@ def test_constructor_objective_lib_error(self):
             )
 
     def test_output_names_obj_lib_error(self):
+        """Ensure output-name query raises on invalid objective option."""
         objective_option = "trace"
         doe_object, grey_box_object = make_greybox_and_doe_objects(
             objective_option=objective_option
@@ -995,6 +1180,7 @@ def test_output_names_obj_lib_error(self):
             grey_box_object.output_names()
 
     def test_evaluate_outputs_obj_lib_error(self):
+        """Ensure output evaluation raises on invalid objective option."""
         objective_option = "trace"
         doe_object, grey_box_object = make_greybox_and_doe_objects(
             objective_option=objective_option
@@ -1008,6 +1194,7 @@ def test_evaluate_outputs_obj_lib_error(self):
             grey_box_object.evaluate_outputs()
 
     def test_evaluate_jacobian_outputs_obj_lib_error(self):
+        """Ensure Jacobian evaluation raises on invalid objective option."""
         objective_option = "trace"
         doe_object, grey_box_object = make_greybox_and_doe_objects(
             objective_option=objective_option
@@ -1021,6 +1208,7 @@ def test_evaluate_jacobian_outputs_obj_lib_error(self):
             grey_box_object.evaluate_jacobian_outputs()
 
     def test_evaluate_hessian_outputs_obj_lib_error(self):
+        """Ensure Hessian evaluation raises on invalid objective option."""
         objective_option = "trace"
         doe_object, grey_box_object = make_greybox_and_doe_objects(
             objective_option=objective_option
@@ -1039,6 +1227,12 @@ def test_evaluate_hessian_outputs_obj_lib_error(self):
         not cyipopt_call_working, "cyipopt is not properly accessing linear solvers"
     )
     def test_solve_D_optimality_log_determinant(self):
+        """
+        Solve Rooney-Biegler D-opt design and verify one known local optimum.
+
+        The assertion accepts either known local design in ``(time, log10(det))``
+        space because the NLP can converge to either basin.
+        """
         # Two locally optimal design points exist
         # (time, optimal objective value)
         # Here, the objective value is
@@ -1080,6 +1274,12 @@ def test_solve_D_optimality_log_determinant(self):
         not cyipopt_call_working, "cyipopt is not properly accessing linear solvers"
     )
     def test_solve_A_optimality_trace_of_inverse(self):
+        """
+        Solve Rooney-Biegler A-opt design and verify one known local optimum.
+
+        The assertion accepts either known local design in ``(time, trace(FIM^-1))``
+        space because the NLP can converge to either basin.
+        """
         # Two locally optimal design points exist
         # (time, optimal objective value)
         # Here, the objective value is
@@ -1117,11 +1317,58 @@ def test_solve_A_optimality_trace_of_inverse(self):
             )
         )
 
+    @unittest.skipIf(
+        not cyipopt_call_working, "cyipopt is not properly accessing linear solvers"
+    )
+    def test_solve_pseudo_A_optimality_trace(self):
+        """
+        Solve Rooney-Biegler pseudo-A-opt design and verify the expected basin.
+
+        GreyBox exposes pseudo-A-opt as the raw trace(FIM), so this regression
+        ensures the full GreyBox solve path still converges to the same optimum.
+        """
+        objective_option = "pseudo_trace"
+        doe_object, grey_box_object = make_greybox_and_doe_objects_rooney_biegler(
+            objective_option=objective_option
+        )
+
+        # Set to use the grey box objective
+        doe_object.use_grey_box = True
+
+        # Solve the model
+        doe_object.run_doe()
+
+        self.assertEqual(doe_object.results["Solver Status"], "ok")
+
+        optimal_time_val = doe_object.results["Experiment Design"][0]
+        optimal_obj_val = doe_object.model.objective()
+
+        optimal_design_np_array = np.array([optimal_time_val, optimal_obj_val])
+        optimal_experimental_design = np.array([10.00, 812.9024487361])
+
+        self.assertTrue(
+            np.all(
+                np.isclose(optimal_design_np_array, optimal_experimental_design, 1e-1)
+            )
+        )
+        self.assertAlmostEqual(
+            doe_object.results["log10 pseudo A-opt"],
+            np.log10(optimal_obj_val),
+            places=6,
+        )
+
     @unittest.skipIf(
         not cyipopt_call_working, "cyipopt is not properly accessing linear solvers"
     )
     @unittest.skipIf(not pandas_available, "pandas is not available")
     def test_solve_E_optimality_minimum_eigenvalue(self):
+        """
+        Solve Rooney-Biegler E-opt design and verify one known local optimum.
+
+        The assertion accepts either known local design in
+        ``(time, min_eigenvalue(FIM))`` space because the NLP can converge to
+        either basin.
+        """
         # Two locally optimal design points exist
         # (time, optimal objective value)
         # Here, the objective value is
@@ -1164,6 +1411,13 @@ def test_solve_E_optimality_minimum_eigenvalue(self):
     )
     @unittest.skipIf(not pandas_available, "pandas is not available")
     def test_solve_ME_optimality_condition_number(self):
+        """
+        Solve Rooney-Biegler ME-opt design and verify one known local optimum.
+
+        The assertion accepts either known local design in
+        ``(time, log(condition_number(FIM)))`` space because the NLP can
+        converge to either basin.
+        """
         # Two locally optimal design points exist
         # (time, optimal objective value)
         # Here, the objective value is
diff --git a/pyomo/contrib/doe/tests/test_utils.py b/pyomo/contrib/doe/tests/test_utils.py
index a6f6cb38a23..6b935932890 100644
--- a/pyomo/contrib/doe/tests/test_utils.py
+++ b/pyomo/contrib/doe/tests/test_utils.py
@@ -13,6 +13,7 @@
     check_FIM,
     compute_FIM_metrics,
     get_FIM_metrics,
+    rescale_FIM,
     _SMALL_TOLERANCE_DEFINITENESS,
     _SMALL_TOLERANCE_SYMMETRY,
     _SMALL_TOLERANCE_IMG,
@@ -60,7 +61,29 @@ def test_check_FIM_non_symmetric(self):
         ):
             check_FIM(FIM)
 
-    """Test the compute_FIM_metrics() from utils.py."""
+    def test_rescale_FIM_bad_type(self):
+        """Reject scalar input because `param_vals` must be a list or array."""
+        FIM = np.array([[4, 1], [1, 3]])
+
+        # Keep this focused on the public input contract, not the numerical result.
+        with self.assertRaisesRegex(
+            ValueError,
+            "param_vals should be a list or numpy array of dimensions: 1 by `n_params`",
+        ):
+            rescale_FIM(FIM, 1.0)
+
+    def test_rescale_FIM_bad_shape(self):
+        """Reject multi-row arrays because scaling needs a single parameter row."""
+        FIM = np.array([[4, 1], [1, 3]])
+        param_vals = np.ones((2, 2))
+
+        # A 2x2 input would silently broadcast the wrong way if we did not check it.
+        with self.assertRaisesRegex(
+            ValueError,
+            r"param_vals should be a vector of dimensions: 1 by `n_params`\. "
+            r"The shape you provided is \(2, 2\)\.",
+        ):
+            rescale_FIM(FIM, param_vals)
 
     ### Helper methods for test cases
     # Sample FIM for testing
@@ -85,6 +108,7 @@ def _get_expected_fim_results(self):
         }
 
     def test_compute_FIM_metrics(self):
+        """Verify `compute_FIM_metrics` returns expected metrics for a known FIM."""
         # Create a sample Fisher Information Matrix (FIM)
         FIM = self._get_test_fim()
         # expected results
@@ -116,6 +140,7 @@ def test_compute_FIM_metrics(self):
         self.assertAlmostEqual(ME_opt, expected['ME_opt'])
 
     def test_FIM_eigenvalue_warning(self):
+        """Verify warning is logged when eigenvalues have non-negligible imaginary part."""
         # Create a matrix with an imaginary component large enough
         # to trigger the warning
         FIM = np.array([[6, 5j], [5j, 7]])
@@ -128,9 +153,8 @@ def test_FIM_eigenvalue_warning(self):
             )
             self.assertIn(expected_warning, cm.output[0])
 
-    """Test the get_FIM_metrics() from utils.py."""
-
     def test_get_FIM_metrics(self):
+        """Verify `get_FIM_metrics` dictionary entries match expected values."""
         # Create a sample Fisher Information Matrix (FIM)
         FIM = self._get_test_fim()
         # expected results