diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py
index a13e6785bf4..8c5f5f86ada 100644
--- a/pyomo/contrib/doe/doe.py
+++ b/pyomo/contrib/doe/doe.py
@@ -52,9 +52,10 @@
 
 import pyomo.environ as pyo
 from pyomo.contrib.doe.utils import (
+    _SMALL_TOLERANCE_DEFINITENESS,
     check_FIM,
     compute_FIM_metrics,
-    _SMALL_TOLERANCE_DEFINITENESS,
+    regularize_fim_for_cholesky,
 )
 
 from pyomo.opt import SolverStatus
@@ -281,7 +282,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 +376,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 +394,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 +487,92 @@ def run_doe(self, model=None, results_file=None):
             with open(results_file, "w") as file:
                 json.dump(self.results, file)
 
+    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"):
+            # The model doesn't have the Cholesky variables, so we can't initialize them.
+            # This happens if the function is called with a model using GreyBox.
+            return
+
+        fim_np = self._get_fim_numpy(model)
+        fim_pd, _ = regularize_fim_for_cholesky(fim_np)
+
+        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 +934,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 "
@@ -1372,29 +1409,14 @@ def create_objective_function(self, model=None):
         model.obj_cons = pyo.Block()
 
         # Assemble the FIM matrix. This is helpful for initialization!
-        fim_vals = [
-            model.fim[bu, un].value
-            for i, bu in enumerate(model.parameter_names)
-            for j, un in enumerate(model.parameter_names)
-        ]
-        fim = np.array(fim_vals).reshape(
-            len(model.parameter_names), len(model.parameter_names)
-        )
+        fim = self._get_fim_numpy(model)
 
         ### Initialize the Cholesky decomposition matrix
         if self.Cholesky_option and self.objective_option in (
             ObjectiveLib.determinant,
             ObjectiveLib.trace,
         ):
-            # Calculate the eigenvalues of the FIM matrix
-            eig = np.linalg.eigvals(fim)
-
-            # If the smallest eigenvalue is (practically) negative,
-            # add a diagonal matrix to make it positive definite
-            if min(eig) < small_number:
-                fim = fim + np.eye(len(model.parameter_names)) * (
-                    small_number - min(eig)
-                )
+            fim, _ = regularize_fim_for_cholesky(fim)
 
             # Compute the Cholesky decomposition of the FIM matrix
             L = np.linalg.cholesky(fim)
@@ -1663,6 +1685,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..84eb6187a71 100644
--- a/pyomo/contrib/doe/grey_box_utilities.py
+++ b/pyomo/contrib/doe/grey_box_utilities.py
@@ -108,9 +108,11 @@ 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 equation S5.
+        # https://arxiv.org/abs/2604.03354v1
         upt_FIM = self._input_values
 
         # Create FIM in the correct way
@@ -164,6 +166,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 +200,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 +237,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 +276,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_build.py b/pyomo/contrib/doe/tests/test_doe_build.py
index 7ba299dc4d0..333ac1669ef 100644
--- a/pyomo/contrib/doe/tests/test_doe_build.py
+++ b/pyomo/contrib/doe/tests/test_doe_build.py
@@ -25,9 +25,6 @@
 
 if scipy_available:
     from pyomo.contrib.doe import DesignOfExperiments
-    from pyomo.contrib.doe.examples.reactor_example import (
-        ReactorExperiment as FullReactorExperiment,
-    )
     from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import (
         RooneyBieglerExperiment,
     )
@@ -98,10 +95,6 @@ def get_FIM_FIMPrior_Q_L(doe_obj=None):
         for j in model.parameter_names
     ]
     sigma_inv = [1 / v for k, v in model.scenario_blocks[0].measurement_error.items()]
-    param_vals = np.array(
-        [[v for k, v in model.scenario_blocks[0].unknown_parameters.items()]]
-    )
-
     FIM_vals_np = np.array(FIM_vals).reshape((n_param, n_param))
     FIM_prior_vals_np = np.array(FIM_prior_vals).reshape((n_param, n_param))
 
diff --git a/pyomo/contrib/doe/tests/test_doe_errors.py b/pyomo/contrib/doe/tests/test_doe_errors.py
index 10165aa6267..b5587760148 100644
--- a/pyomo/contrib/doe/tests/test_doe_errors.py
+++ b/pyomo/contrib/doe/tests/test_doe_errors.py
@@ -27,9 +27,6 @@
         BadExperiment,
         RooneyBieglerExperimentFlag,
     )
-    from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import (
-        RooneyBieglerExperiment,
-    )
 
 from pyomo.contrib.doe.examples.rooney_biegler_doe_example import run_rooney_biegler_doe
 from pyomo.opt import SolverFactory
@@ -95,7 +92,7 @@ def test_experiment_none_error(self):
             # Experiment provided as None
             DoE_args = get_standard_args(None, fd_method, obj_used, flag_val)
 
-            doe_obj = DesignOfExperiments(**DoE_args)
+            DesignOfExperiments(**DoE_args)
 
     def test_reactor_check_no_get_labeled_model(self):
         fd_method = "central"
@@ -110,7 +107,7 @@ def test_reactor_check_no_get_labeled_model(self):
         ):
             DoE_args = get_standard_args(experiment, fd_method, obj_used, flag_val)
 
-            doe_obj = DesignOfExperiments(**DoE_args)
+            DesignOfExperiments(**DoE_args)
 
     def test_reactor_check_no_experiment_outputs(self):
         fd_method = "central"
@@ -781,6 +778,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..82390383cb4
--- /dev/null
+++ b/pyomo/contrib/doe/tests/test_doe_initialization.py
@@ -0,0 +1,163 @@
+# ____________________________________________________________________________________
+#
+# 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.parmest.experiment import Experiment
+from pyomo.contrib.doe.utils import (
+    _SMALL_TOLERANCE_DEFINITENESS,
+    regularize_fim_for_cholesky,
+)
+from pyomo.opt import SolverFactory
+
+ipopt_available = SolverFactory("ipopt").available()
+
+
+class _StructurallyUnidentifiableExperiment(Experiment):
+    """
+    Minimal real experiment with a rank-deficient FIM.
+
+    The single measured response depends only on ``theta1 + theta2``:
+    ``y = (theta1 + theta2) * x``. With one output and two parameters, both
+    sensitivities are identical, so the analytical FIM is
+    ``[[1, 1], [1, 1]]`` when ``x = 1`` and the measurement error is 1.
+    """
+
+    def get_labeled_model(self):
+        m = pyo.ConcreteModel()
+        m.x = pyo.Var(initialize=1.0, bounds=(1.0, 1.0))
+        m.theta1 = pyo.Var(initialize=2.0)
+        m.theta2 = pyo.Var(initialize=3.0)
+        m.theta1.fix()
+        m.theta2.fix()
+        m.y = pyo.Var(initialize=5.0)
+        m.response = pyo.Constraint(expr=m.y == (m.theta1 + m.theta2) * m.x)
+
+        m.experiment_inputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+        m.experiment_inputs[m.x] = 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.y] = 5.0
+        m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+        m.measurement_error[m.y] = 1.0
+        return m
+
+
+def _make_unidentifiable_doe_object(objective_option="trace"):
+    """
+    Build a DoE object for the structurally unidentifiable two-parameter model.
+    """
+    return DesignOfExperiments(
+        experiment=_StructurallyUnidentifiableExperiment(),
+        fd_formula="central",
+        step=1e-3,
+        objective_option=objective_option,
+        solver=SolverFactory("ipopt"),
+        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_regularize_fim_for_cholesky_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.
+        """
+        fim = np.array([[2.0, -1.0], [-1.0, 0.0]])
+        fim_pd, jitter = regularize_fim_for_cholesky(fim)
+        self.assertGreater(jitter, 0.0)
+        min_eig = np.min(np.linalg.eigvalsh(fim_pd))
+        self.assertAlmostEqual(
+            min_eig,
+            _SMALL_TOLERANCE_DEFINITENESS,
+            delta=_SMALL_TOLERANCE_DEFINITENESS / 10,
+        )
+
+    def test_regularize_fim_for_cholesky_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.
+        """
+        fim = np.array([[1.0e-2]])
+        fim_pd, jitter = regularize_fim_for_cholesky(fim)
+        self.assertEqual(jitter, 0.0)
+        self.assertAlmostEqual(fim_pd[0, 0], 1.0e-2)
+
+    @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available")
+    def test_trace_initialization_regularizes_structurally_unidentifiable_fim(self):
+        """
+        Verify trace-mode initialization regularizes a singular real-model FIM.
+
+        This experiment is structurally unidentifiable because the single
+        response depends only on ``theta1 + theta2``. Its analytical FIM is
+        rank-deficient, so Cholesky-based initialization must add diagonal
+        regularization before computing the Cholesky factors.
+
+        We intentionally do not call ``run_doe()`` here. The behavior under
+        test is the Cholesky/FIM synchronization step after a singular FIM is
+        available. ``run_doe()`` continues on to additional solver-driven
+        optimization steps and does not expose a clean pause point between
+        "the FIM values are now known" and "initialize the Cholesky-related
+        variables from those values." Instead, this test uses the real DoE
+        model-building machinery, writes the analytically known singular FIM
+        for this experiment onto ``model.fim``, and then directly exercises
+        ``_initialize_cholesky_from_fim()``.
+        """
+        doe_obj = _make_unidentifiable_doe_object(objective_option="trace")
+        doe_obj.create_doe_model()
+        doe_obj.create_objective_function()
+
+        model = doe_obj.model
+        params = list(model.parameter_names)
+        expected_fim = np.array([[1.0, 1.0], [1.0, 1.0]])
+        # Pyomo.DoE stores only the lower triangle when
+        # ``only_compute_fim_lower`` is enabled. Populate ``model.fim`` using
+        # the analytical singular FIM for this experiment so the initialization
+        # helper sees the exact pathological case we want to regularize.
+        for i, p in enumerate(params):
+            for j, q in enumerate(params):
+                if doe_obj.only_compute_fim_lower and i < j:
+                    model.fim[p, q].set_value(0.0)
+                else:
+                    model.fim[p, q].set_value(expected_fim[i, j])
+
+        doe_obj._initialize_cholesky_from_fim()
+
+        L = np.array([[pyo.value(model.L[i, j]) for j in params] for i in params])
+
+        expected_fim_pd, jitter = regularize_fim_for_cholesky(expected_fim)
+        expected_L = np.linalg.cholesky(expected_fim_pd)
+        reconstructed_fim = L @ L.T
+
+        self.assertGreater(jitter, 0.0)
+        self.assertTrue(np.allclose(L, expected_L))
+        self.assertTrue(np.allclose(reconstructed_fim, expected_fim_pd))
+        self.assertAlmostEqual(
+            np.min(np.linalg.eigvalsh(reconstructed_fim)),
+            _SMALL_TOLERANCE_DEFINITENESS,
+            delta=_SMALL_TOLERANCE_DEFINITENESS / 10,
+        )
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/pyomo/contrib/doe/tests/test_doe_solve.py b/pyomo/contrib/doe/tests/test_doe_solve.py
index 9e31accf894..4329ddf37ab 100644
--- a/pyomo/contrib/doe/tests/test_doe_solve.py
+++ b/pyomo/contrib/doe/tests/test_doe_solve.py
@@ -33,7 +33,6 @@
 
 if scipy_available:
     from pyomo.contrib.doe import DesignOfExperiments
-    from pyomo.contrib.doe.examples.reactor_experiment import ReactorExperiment
     from pyomo.contrib.doe.examples.reactor_example import (
         ReactorExperiment as FullReactorExperiment,
         run_reactor_doe,
@@ -115,10 +114,6 @@ def get_FIM_Q_L(doe_obj=None):
     sigma_inv = [
         1 / v**2 for k, v in model.scenario_blocks[0].measurement_error.items()
     ]
-    param_vals = np.array(
-        [[v for k, v in model.scenario_blocks[0].unknown_parameters.items()]]
-    )
-
     FIM_vals_np = np.array(FIM_vals).reshape((n_param, n_param))
 
     for i in range(n_param):
diff --git a/pyomo/contrib/doe/tests/test_greybox.py b/pyomo/contrib/doe/tests/test_greybox.py
index 2fd007e6d52..f1bcbbe9f64 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
 
@@ -376,7 +426,7 @@ def make_greybox_and_doe_objects_rooney_biegler(objective_option):
         cyipopt_call_working = not (
             bad_message in doe_object.results["Termination Message"]
         )
-    except:
+    except Exception:
         cyipopt_call_working = False
 
 
@@ -389,12 +439,13 @@ class TestFIMExternalGreyBox(unittest.TestCase):
     # set the inputs for the
     # Grey Box object
     def test_set_inputs(self):
+        """Verify grey-box packed inputs reconstruct the full symmetric FIM."""
         objective_option = "trace"
         doe_obj, grey_box_object = make_greybox_and_doe_objects(
             objective_option=objective_option
         )
 
-        # Set input values to the random testing matrix
+        # Set input values to the random testing matrix (upper triangular part)
         grey_box_object.set_input_values(testing_matrix[masking_matrix > 0])
 
         # Grab the values from get_FIM
@@ -402,9 +453,8 @@ def test_set_inputs(self):
 
         self.assertTrue(np.all(np.isclose(grey_box_FIM, testing_matrix)))
 
-    # Testing that getting the
-    # input names works properly
     def test_input_names(self):
+        """Verify grey-box input names follow the expected packed triangular order."""
         objective_option = "trace"
         doe_obj, grey_box_object = make_greybox_and_doe_objects(
             objective_option=objective_option
@@ -432,6 +482,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 +490,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 +517,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 +531,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 +546,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 +579,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 +596,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 +614,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 +638,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 +687,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 +712,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 +736,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 +780,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 +827,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 +850,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 +873,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 +887,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 +901,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 +916,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,42 +934,27 @@ 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)
+        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()
 
-        # Check to see if each component exists
-        all_exist = True
-
         # Check output and value
         # FIM Initial will be the prior FIM
         # added with the identity matrix.
         A_opt_val = np.trace(np.linalg.inv(testing_matrix + np.eye(4)))
-
-        try:
-            A_opt_val_gb = doe_obj.model.obj_cons.egb_fim_block.outputs["A-opt"].value
-        except:
-            A_opt_val_gb = -10.0  # Trace should never be negative
-            all_exist = False
-
-        # Intermediate check for output existence
-        self.assertTrue(all_exist)
+        A_opt_val_gb = doe_obj.model.obj_cons.egb_fim_block.outputs["A-opt"].value
         self.assertAlmostEqual(A_opt_val, A_opt_val_gb)
 
         # Check inputs and values
-        try:
-            input_values = []
-            for i in _.input_names():
-                input_values.append(doe_obj.model.obj_cons.egb_fim_block.inputs[i]())
-        except:
-            input_values = np.zeros_like(testing_matrix)
-            all_exist = False
-
-        # Final check on existence of inputs
-        self.assertTrue(all_exist)
+        input_values = []
+        for i in grey_box_object.input_names():
+            input_values.append(doe_obj.model.obj_cons.egb_fim_block.inputs[i]())
         # Rebuild the current FIM from the input
         # values taken from the egb_fim_block
         current_FIM = np.zeros_like(testing_matrix)
@@ -823,44 +964,33 @@ 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)
+        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()
 
-        # Check to see if each component exists
-        all_exist = True
+        self.assertTrue(doe_obj.only_compute_fim_lower)
+        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)
 
         # Check output and value
         # FIM Initial will be the prior FIM
         # added with the identity matrix.
         D_opt_val = np.log(np.linalg.det(testing_matrix + np.eye(4)))
-
-        try:
-            D_opt_val_gb = doe_obj.model.obj_cons.egb_fim_block.outputs[
-                "log-D-opt"
-            ].value
-        except:
-            D_opt_val_gb = -100.0  # Determinant should never be negative beyond -64
-            all_exist = False
-
-        # Intermediate check for output existence
-        self.assertTrue(all_exist)
+        D_opt_val_gb = doe_obj.model.obj_cons.egb_fim_block.outputs["log-D-opt"].value
         self.assertAlmostEqual(D_opt_val, D_opt_val_gb)
 
         # Check inputs and values
-        try:
-            input_values = []
-            for i in _.input_names():
-                input_values.append(doe_obj.model.obj_cons.egb_fim_block.inputs[i]())
-        except:
-            input_values = np.zeros_like(testing_matrix)
-            all_exist = False
-
-        # Final check on existence of inputs
-        self.assertTrue(all_exist)
+        input_values = []
+        for i in grey_box_object.input_names():
+            input_values.append(doe_obj.model.obj_cons.egb_fim_block.inputs[i]())
         # Rebuild the current FIM from the input
         # values taken from the egb_fim_block
         current_FIM = np.zeros_like(testing_matrix)
@@ -870,43 +1000,28 @@ def test_D_opt_greybox_build(self):
         self.assertTrue(np.all(np.isclose(current_FIM, testing_matrix + np.eye(4))))
 
     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)
+        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()
 
-        # Check to see if each component exists
-        all_exist = True
-
         # Check output and value
         # FIM Initial will be the prior FIM
         # added with the identity matrix.
         vals, vecs = np.linalg.eig(testing_matrix + np.eye(4))
         E_opt_val = np.min(vals)
-
-        try:
-            E_opt_val_gb = doe_obj.model.obj_cons.egb_fim_block.outputs["E-opt"].value
-        except:
-            E_opt_val_gb = -10.0  # Determinant should never be negative
-            all_exist = False
-
-        # Intermediate check for output existence
-        self.assertTrue(all_exist)
+        E_opt_val_gb = doe_obj.model.obj_cons.egb_fim_block.outputs["E-opt"].value
         self.assertAlmostEqual(E_opt_val, E_opt_val_gb)
 
         # Check inputs and values
-        try:
-            input_values = []
-            for i in _.input_names():
-                input_values.append(doe_obj.model.obj_cons.egb_fim_block.inputs[i]())
-        except:
-            input_values = np.zeros_like(testing_matrix)
-            all_exist = False
-
-        # Final check on existence of inputs
-        self.assertTrue(all_exist)
+        input_values = []
+        for i in grey_box_object.input_names():
+            input_values.append(doe_obj.model.obj_cons.egb_fim_block.inputs[i]())
         # Rebuild the current FIM from the input
         # values taken from the egb_fim_block
         current_FIM = np.zeros_like(testing_matrix)
@@ -916,43 +1031,28 @@ 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)
+        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()
 
-        # Check to see if each component exists
-        all_exist = True
-
         # Check output and value
         # FIM Initial will be the prior FIM
         # added with the identity matrix.
         vals, vecs = np.linalg.eig(testing_matrix + np.eye(4))
         ME_opt_val = np.log(np.abs(np.max(vals) / np.min(vals)))
-
-        try:
-            ME_opt_val_gb = doe_obj.model.obj_cons.egb_fim_block.outputs["ME-opt"].value
-        except:
-            ME_opt_val_gb = -10.0  # Condition number should not be negative
-            all_exist = False
-
-        # Intermediate check for output existence
-        self.assertTrue(all_exist)
+        ME_opt_val_gb = doe_obj.model.obj_cons.egb_fim_block.outputs["ME-opt"].value
         self.assertAlmostEqual(ME_opt_val, ME_opt_val_gb)
 
         # Check inputs and values
-        try:
-            input_values = []
-            for i in _.input_names():
-                input_values.append(doe_obj.model.obj_cons.egb_fim_block.inputs[i]())
-        except:
-            input_values = np.zeros_like(testing_matrix)
-            all_exist = False
-
-        # Final check on existence of inputs
-        self.assertTrue(all_exist)
+        input_values = []
+        for i in grey_box_object.input_names():
+            input_values.append(doe_obj.model.obj_cons.egb_fim_block.inputs[i]())
         # Rebuild the current FIM from the input
         # values taken from the egb_fim_block
         current_FIM = np.zeros_like(testing_matrix)
@@ -963,13 +1063,15 @@ 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.",
         ):
-            grey_box_object = FIMExternalGreyBox(doe_object=None)
+            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
@@ -977,11 +1079,12 @@ def test_constructor_objective_lib_error(self):
         with self.assertRaisesRegex(
             ValueError, "'Bad Objective Option' is not a valid ObjectiveLib"
         ):
-            bad_grey_box_object = FIMExternalGreyBox(
+            FIMExternalGreyBox(
                 doe_object=doe_object, objective_option="Bad Objective Option"
             )
 
     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 +1098,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 +1112,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 +1126,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 +1145,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 +1192,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 +1235,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 +1329,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
diff --git a/pyomo/contrib/doe/utils.py b/pyomo/contrib/doe/utils.py
index c85fbbe4c0b..24d2ecebd56 100644
--- a/pyomo/contrib/doe/utils.py
+++ b/pyomo/contrib/doe/utils.py
@@ -91,6 +91,28 @@ def rescale_FIM(FIM, param_vals):
     return scaled_FIM
 
 
+def regularize_fim_for_cholesky(FIM):
+    """
+    Add the smallest diagonal shift needed for a stable Cholesky factorization.
+
+    Parameters
+    ----------
+    FIM: ndarray
+        Fisher information matrix to regularize before Cholesky factorization.
+
+    Returns
+    -------
+    tuple
+        Pair ``(fim_pd, jitter)`` where ``fim_pd`` is the regularized matrix and
+        ``jitter`` is the nonnegative diagonal shift used to make the smallest
+        symmetric eigenvalue at least ``_SMALL_TOLERANCE_DEFINITENESS``.
+    """
+    min_eig = float(np.min(np.linalg.eigvalsh(FIM)))
+    jitter = max(0.0, _SMALL_TOLERANCE_DEFINITENESS - min_eig)
+    fim_pd = FIM + jitter * np.eye(FIM.shape[0])
+    return fim_pd, jitter
+
+
 def check_FIM(FIM):
     """
     Checks that the FIM is square, positive definite, and symmetric.