diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py
index 8c5f5f86ada..cbadab3de88 100644
--- a/pyomo/contrib/doe/doe.py
+++ b/pyomo/contrib/doe/doe.py
@@ -450,12 +450,13 @@ def run_doe(self, model=None, results_file=None):
         self.results["Unknown Parameters"] = self.get_unknown_parameter_values()
         self.results["Unknown Parameter Names"] = [
             str(pyo.ComponentUID(k, context=model.scenario_blocks[0]))
-            for k in model.scenario_blocks[0].unknown_parameters
+            for k in self._expanded_unknown_parameters(model.scenario_blocks[0])
         ]
         self.results["Measurement Error"] = self.get_measurement_error_values()
         self.results["Measurement Error Names"] = [
             str(pyo.ComponentUID(k, context=model.scenario_blocks[0]))
             for k in model.scenario_blocks[0].measurement_error
+            if hasattr(k, "name")
         ]
 
         self.results["Prior FIM"] = [list(row) for row in list(self.prior_FIM)]
@@ -581,6 +582,156 @@ def run_multi_doe_sequential(self, N_exp=1):
     def run_multi_doe_simultaneous(self, N_exp=1):
         raise NotImplementedError("Multiple experiment optimization not yet supported.")
 
+    def _build_meas_error_covariance_matrix(self, model):
+        """
+        Builds the full measurement-error covariance matrix
+
+        Note: The code does not automatically build the full covariance matrix
+        It only places whatever standard deviation and covariances the user
+        explicitly supplies
+        Therefore, for correlation in time, shared timepoints, or other types of
+        correlation, the user must provide all the desired covariance terms
+        They standard deviations may be constant or depend on the value (i.e.,
+        data) of the measured or input variables (e.g., be proportional to them)
+
+        Parameters
+        ----------
+        model : ConcreteModel
+            Annotated Pyomo model
+
+        Returns
+        -------
+        Sigma_y: numpy.ndarray
+            Full measurement-error covariance matrix
+        """
+        # get the output variables
+        outputs = list(model.experiment_outputs.keys())
+        outputs_name = [y_hat.name for y_hat in outputs]
+        outputs_index = {y_hat_name: i for i, y_hat_name in enumerate(outputs_name)}
+
+        # get the number of output variables
+        number_outputs = len(outputs)
+
+        # define the measurement-error covariance matrix
+        Sigma_y = np.zeros((number_outputs, number_outputs))
+
+        # check if all the values of the measurement-error standard deviation
+        # have been supplied
+        try:
+            all_known_errors = all(
+                model.measurement_error[y_hat] is not None
+                for y_hat in model.experiment_outputs
+            )
+        except KeyError:
+            raise KeyError(
+                'One or more experiment outputs are not defined in the '
+                '"measurement_error" attribute. All the variables defined '
+                'in "experiment_outputs" must be defined as keys in '
+                '"measurement_error".'
+            )
+
+        # fill the leading-diagonal elements from the standard deviation of
+        # the measurement errors
+        for y_hat in outputs:
+            # get the index of y_hat
+            i = outputs_index[y_hat.name]
+
+            if all_known_errors:
+                standard_dev = model.measurement_error[y_hat]
+                Sigma_y[i, i] = standard_dev**2
+            else:
+                raise ValueError(
+                    'One or more values are missing from "measurement_error". All '
+                    'values of the measurement errors are required to compute the '
+                    'Fisher information matrix.'
+                )
+
+        # fill the off-diagonal elements from covariance entries
+        # supplied by the user
+        for key, err_cov in model.measurement_error.items():
+            if isinstance(key, tuple) and len(key) == 2:
+                yi, yj = key
+                if yi.name not in outputs_name or yj.name not in outputs_name:
+                    raise ValueError(
+                        "Measurement-error covariance must be defined only between "
+                        "experiment output variables."
+                    )
+
+                # get the indices of yi and yj
+                i = outputs_index[yi.name]
+                j = outputs_index[yj.name]
+
+                # update the measurement-error covariance matrix which
+                # is a symmetric matrix
+                Sigma_y[i, j] = err_cov
+                Sigma_y[j, i] = err_cov
+            elif not isinstance(key, tuple) and not hasattr(key, "name"):
+                raise TypeError(
+                    "Expected a tuple of two measured variables when specifying a "
+                    "measurement-error covariance, e.g., "
+                    "measurement_error[(y1, y2)] = covariance."
+                )
+
+        return Sigma_y
+
+    def get_meas_error_covariance_matrix_inv(self, model):
+        """
+        Computes the inverse of the measurement-error covariance matrix
+
+        Parameters
+        ----------
+        model : ConcreteModel
+            Annotated Pyomo model
+
+        Returns
+        -------
+        Sigma_y_inv: numpy.ndarray
+            Inverse of the measurement-error covariance matrix
+        """
+        # get the measurement-error covariance matrix
+        Sigma_y = self._build_meas_error_covariance_matrix(model)
+
+        # compute the inverse of the measurement-error covariance matrix
+        try:
+            Sigma_y_inv = np.linalg.inv(Sigma_y)
+        except np.linalg.LinAlgError:
+            Sigma_y_inv = np.linalg.pinv(Sigma_y)
+            logger.warning(
+                "The measurement-error covariance matrix is singular. "
+                "Using pseudo-inverse instead."
+            )
+
+        return Sigma_y_inv
+
+    def _expanded_unknown_parameters(self, model):
+        """
+        Creates a list of scalar unknown parameter components
+
+        The unknown_parameters suffix may contain either scalar ComponentData
+        objects or indexed components. Indexed components are expanded to
+        their scalar data objects.
+
+        Parameters
+        ----------
+        model : ConcreteModel
+            Annotated Pyomo model
+
+        Returns
+        -------
+        params_data_object: list
+            List of scalar parameter data objects
+        """
+        params_data_object = []
+        for component in model.unknown_parameters:
+            # check if it is indexed
+            if component.is_indexed():
+                # get the parameter data objects
+                params_data_object.extend(component.values())
+            else:
+                params_data_object.append(component)
+
+        return params_data_object
+
     # Compute FIM for the DoE object
     def compute_FIM(self, model=None, method="sequential"):
         """
@@ -615,16 +766,16 @@ def compute_FIM(self, model=None, method="sequential"):
         self.check_model_labels(model=model)
 
         # Set length values for the model features
-        self.n_parameters = len(model.unknown_parameters)
-        self.n_measurement_error = len(model.measurement_error)
+        self.n_parameters = len(self._expanded_unknown_parameters(model))
+        self.n_measurement_error = len(
+            [k for k in model.measurement_error if hasattr(k, "name")]
+        )
         self.n_experiment_inputs = len(model.experiment_inputs)
         self.n_experiment_outputs = len(model.experiment_outputs)
 
         # Check FIM input, if it exists. Otherwise, set the prior_FIM attribute
         if self.prior_FIM is None:
-            self.prior_FIM = np.zeros(
-                (len(model.unknown_parameters), len(model.unknown_parameters))
-            )
+            self.prior_FIM = np.zeros((self.n_parameters, self.n_parameters))
         else:
             self.check_model_FIM(FIM=self.prior_FIM)
 
@@ -668,25 +819,28 @@ def _sequential_FIM(self, model=None):
             model.del_component(model.parameter_scenarios)
         model.parameter_scenarios = pyo.Suffix(direction=pyo.Suffix.LOCAL)
 
+        # get the parameter data objects
+        unknown_params = self._expanded_unknown_parameters(model)
+        unknown_params_val = {p.name: p.value for p in unknown_params}
+
         # Populate parameter scenarios, and scenario
         # inds based on finite difference scheme
         if self.fd_formula == FiniteDifferenceStep.central:
             model.parameter_scenarios.update(
-                (2 * ind, k) for ind, k in enumerate(model.unknown_parameters.keys())
+                (2 * ind, k) for ind, k in enumerate(unknown_params)
             )
             model.parameter_scenarios.update(
-                (2 * ind + 1, k)
-                for ind, k in enumerate(model.unknown_parameters.keys())
+                (2 * ind + 1, k) for ind, k in enumerate(unknown_params)
             )
-            model.scenarios = range(len(model.unknown_parameters) * 2)
+            model.scenarios = range(self.n_parameters * 2)
         elif self.fd_formula in [
             FiniteDifferenceStep.forward,
             FiniteDifferenceStep.backward,
         ]:
             model.parameter_scenarios.update(
-                (ind + 1, k) for ind, k in enumerate(model.unknown_parameters.keys())
+                (ind + 1, k) for ind, k in enumerate(unknown_params)
             )
-            model.scenarios = range(len(model.unknown_parameters) + 1)
+            model.scenarios = range(self.n_parameters + 1)
         else:
             raise DeveloperError(
                 "Finite difference option not recognized. Please "
@@ -721,7 +875,7 @@ def _sequential_FIM(self, model=None):
             if not skip_param_update:
                 param = model.parameter_scenarios[s]
                 # Update parameter values for the given finite difference scenario
-                param.set_value(model.unknown_parameters[param] * (1 + diff))
+                param.set_value(unknown_params_val[param.name] * (1 + diff))
             else:
                 continue
 
@@ -741,7 +895,7 @@ def _sequential_FIM(self, model=None):
 
             # Reset value of parameter to default value
             # before computing finite difference perturbation
-            param.set_value(model.unknown_parameters[param])
+            param.set_value(unknown_params_val[param.name])
 
             # Extract the measurement values for the scenario and append
             measurement_vals.append(
@@ -751,12 +905,7 @@ def _sequential_FIM(self, model=None):
         # Use the measurement outputs to make the Q matrix
         measurement_vals_np = np.array(measurement_vals).T
 
-        self.seq_jac = np.zeros(
-            (
-                len(model.experiment_outputs.items()),
-                len(model.unknown_parameters.items()),
-            )
-        )
+        self.seq_jac = np.zeros((self.n_experiment_outputs, self.n_parameters))
 
         # Counting variable for loop
         i = 0
@@ -764,7 +913,8 @@ def _sequential_FIM(self, model=None):
         # Loop over parameter values and grab correct
         # columns for finite difference calculation
 
-        for k, v in model.unknown_parameters.items():
+        for k in unknown_params:
+            v = pyo.value(k)
             curr_step = v * self.step
 
             if self.fd_formula == FiniteDifferenceStep.central:
@@ -792,18 +942,11 @@ def _sequential_FIM(self, model=None):
             # Increment the count
             i += 1
 
-        # TODO: As more complex measurement error schemes
-        #       are put in place, this needs to change
-        # Add independent (non-correlated) measurement
-        # error for FIM calculation
-        cov_y = np.zeros((len(model.measurement_error), len(model.measurement_error)))
-        count = 0
-        for k, v in model.measurement_error.items():
-            cov_y[count, count] = 1 / v**2
-            count += 1
+        # get the inverse of the measurement-error covariance matrix
+        Sigma_y_inv = self.get_meas_error_covariance_matrix_inv(model)
 
         # Compute and record FIM
-        self.seq_FIM = self.seq_jac.T @ cov_y @ self.seq_jac + self.prior_FIM
+        self.seq_FIM = self.seq_jac.T @ Sigma_y_inv @ self.seq_jac + self.prior_FIM
 
     # Use kaug to get FIM
     def _kaug_FIM(self, model=None):
@@ -834,8 +977,11 @@ def _kaug_FIM(self, model=None):
 
         self.solver.solve(model, tee=self.tee)
 
+        # get the parameter data objects
+        unknown_params = self._expanded_unknown_parameters(model)
+
         # Probe the solved model for dsdp results (sensitivities s.t. parameters)
-        params_dict = {k.name: v for k, v in model.unknown_parameters.items()}
+        params_dict = {p.name: p.value for p in unknown_params}
         params_names = list(params_dict.keys())
 
         dsdp_re, col = get_dsdp(model, params_names, params_dict, tee=self.tee)
@@ -868,12 +1014,12 @@ def _kaug_FIM(self, model=None):
         jac = [[] for k in params_names]
 
         for d in range(len(dsdp_extract)):
-            for k, v in model.unknown_parameters.items():
+            for k in unknown_params:
                 p = params_names.index(k.name)  # Index of parameter in np array
                 # if scaled by parameter value or constant value
                 sensi = dsdp_extract[d][p] * self.scale_constant_value
                 if self.scale_nominal_param_value:
-                    sensi *= v
+                    sensi *= pyo.value(k)
                 jac[p].append(sensi)
 
         # record kaug jacobian
@@ -881,23 +1027,14 @@ def _kaug_FIM(self, model=None):
 
         # Compute FIM
         if self.prior_FIM is None:
-            self.prior_FIM = np.zeros((len(params_names), len(params_names)))
+            self.prior_FIM = np.zeros(self.n_parameters, self.n_parameters)
         else:
             self.check_model_FIM(FIM=self.prior_FIM)
 
-        # Constructing the Covariance of the measurements for the FIM calculation
-        # The following assumes independent measurement error.
-        cov_y = np.zeros((len(model.measurement_error), len(model.measurement_error)))
-        count = 0
-        for k, v in model.measurement_error.items():
-            cov_y[count, count] = 1 / v**2
-            count += 1
-
-        # TODO: need to add a covariance matrix for measurements (sigma inverse)
-        # i.e., cov_y = self.cov_y or model.cov_y
-        # Still deciding where this would be best.
+        # get the inverse of the measurement-error covariance matrix
+        Sigma_y_inv = self.get_meas_error_covariance_matrix_inv(model)
 
-        self.kaug_FIM = self.kaug_jac.T @ cov_y @ self.kaug_jac + self.prior_FIM
+        self.kaug_FIM = self.kaug_jac.T @ Sigma_y_inv @ self.kaug_jac + self.prior_FIM
 
     # Create the DoE model (with ``scenarios`` from finite differencing scheme)
     def create_doe_model(self, model=None):
@@ -948,13 +1085,13 @@ def create_doe_model(self, model=None):
         scen_block_ind = min(
             [
                 k.name.split(".").index("scenario_blocks[0]")
-                for k in model.scenario_blocks[0].unknown_parameters.keys()
+                for k in self._expanded_unknown_parameters(model.scenario_blocks[0])
             ]
         )
         model.parameter_names = pyo.Set(
             initialize=[
                 ".".join(k.name.split(".")[(scen_block_ind + 1) :])
-                for k in model.scenario_blocks[0].unknown_parameters.keys()
+                for k in self._expanded_unknown_parameters(model.scenario_blocks[0])
             ]
         )
         model.output_names = pyo.Set(
@@ -1100,8 +1237,7 @@ def jacobian_rule(m, n, p):
             var_lo = cuid.find_component_on(m.scenario_blocks[s2])
 
             param = m.parameter_scenarios[max(s1, s2)]
-            param_loc = pyo.ComponentUID(param).find_component_on(m.scenario_blocks[0])
-            param_val = m.scenario_blocks[0].unknown_parameters[param_loc]
+            param_val = pyo.value(param)
             param_diff = param_val * fd_step_mult * self.step
 
             if self.scale_nominal_param_value:
@@ -1154,21 +1290,24 @@ def fim_rule(m, p, q):
                 else:
                     return m.fim[p, q] == m.fim[q, p]
             else:
-                return (
-                    m.fim[p, q]
-                    == sum(
-                        1
-                        / m.scenario_blocks[0].measurement_error[
-                            pyo.ComponentUID(n).find_component_on(m.scenario_blocks[0])
-                        ]
-                        ** 2
-                        * m.sensitivity_jacobian[n, p]
-                        * m.sensitivity_jacobian[n, q]
-                        for n in m.output_names
-                    )
-                    + m.prior_FIM[p, q]
+                # create a numpy array for the sensitivity Jacobian
+                sens_jac = np.array(
+                    [
+                        [m.sensitivity_jacobian[y, p] for p in m.parameter_names]
+                        for y in m.output_names
+                    ],
+                    dtype=object,
                 )
 
+                # get the inverse of the measurement-error covariance matrix
+                Sigma_y_inv = self.get_meas_error_covariance_matrix_inv(
+                    m.scenario_blocks[0]
+                )
+
+                fim_expr = sens_jac.T @ Sigma_y_inv @ sens_jac
+
+                return m.fim[p, q] == fim_expr[p_ind, q_ind] + m.prior_FIM[p, q]
+
         model.jacobian_constraint = pyo.Constraint(
             model.output_names, model.parameter_names, rule=jacobian_rule
         )
@@ -1216,8 +1355,10 @@ def _generate_scenario_blocks(self, model=None):
         self.check_model_labels(model=model.base_model)
 
         # Gather lengths of label structures for later use in the model build process
-        self.n_parameters = len(model.base_model.unknown_parameters)
-        self.n_measurement_error = len(model.base_model.measurement_error)
+        self.n_parameters = len(self._expanded_unknown_parameters(model.base_model))
+        self.n_measurement_error = len(
+            [k for k in model.base_model.measurement_error if hasattr(k, "name")]
+        )
         self.n_experiment_inputs = len(model.base_model.experiment_inputs)
         self.n_experiment_outputs = len(model.base_model.experiment_outputs)
 
@@ -1249,27 +1390,28 @@ def _generate_scenario_blocks(self, model=None):
         # are associated with parameters
         model.parameter_scenarios = pyo.Suffix(direction=pyo.Suffix.LOCAL)
 
+        # get the parameter data objects
+        unknown_params = self._expanded_unknown_parameters(model.base_model)
+        unknown_params_val = {p.name: p.value for p in unknown_params}
+
         # Populate parameter scenarios, and scenario
         # inds based on finite difference scheme
         if self.fd_formula == FiniteDifferenceStep.central:
             model.parameter_scenarios.update(
-                (2 * ind, k)
-                for ind, k in enumerate(model.base_model.unknown_parameters.keys())
+                (2 * ind, k) for ind, k in enumerate(unknown_params)
             )
             model.parameter_scenarios.update(
-                (2 * ind + 1, k)
-                for ind, k in enumerate(model.base_model.unknown_parameters.keys())
+                (2 * ind + 1, k) for ind, k in enumerate(unknown_params)
             )
-            model.scenarios = range(len(model.base_model.unknown_parameters) * 2)
+            model.scenarios = range(self.n_parameters * 2)
         elif self.fd_formula in [
             FiniteDifferenceStep.forward,
             FiniteDifferenceStep.backward,
         ]:
             model.parameter_scenarios.update(
-                (ind + 1, k)
-                for ind, k in enumerate(model.base_model.unknown_parameters.keys())
+                (ind + 1, k) for ind, k in enumerate(unknown_params)
             )
-            model.scenarios = range(len(model.base_model.unknown_parameters) + 1)
+            model.scenarios = range(self.n_parameters + 1)
         else:
             raise DeveloperError(
                 "Finite difference option not recognized. Please contact "
@@ -1327,7 +1469,7 @@ def build_block_scenarios(b, s):
             # Update parameter values for the given finite difference scenario
             pyo.ComponentUID(param, context=m.base_model).find_component_on(
                 b
-            ).set_value(m.base_model.unknown_parameters[param] * (1 + diff))
+            ).set_value(unknown_params_val[param.name] * (1 + diff))
 
             # Fix experiment inputs before solve (enforce square solve)
             for comp in b.experiment_inputs:
@@ -1728,7 +1870,7 @@ def check_model_labels(self, model=None):
 
         # Check that experimental inputs exist
         try:
-            outputs = [k.name for k, v in model.experiment_inputs.items()]
+            exp_inputs = [k.name for k, v in model.experiment_inputs.items()]
         except:
             raise RuntimeError(
                 "Experiment model does not have suffix " + '"experiment_inputs".'
@@ -1736,7 +1878,7 @@ def check_model_labels(self, model=None):
 
         # Check that unknown parameters exist
         try:
-            outputs = [k.name for k, v in model.unknown_parameters.items()]
+            unknown_params = self._expanded_unknown_parameters(model)
         except:
             raise RuntimeError(
                 "Experiment model does not have suffix " + '"unknown_parameters".'
@@ -1744,7 +1886,7 @@ def check_model_labels(self, model=None):
 
         # Check that measurement errors exist
         try:
-            outputs = [k.name for k, v in model.measurement_error.items()]
+            meas_error = [k.name for k in model.measurement_error if hasattr(k, "name")]
         except:
             raise RuntimeError(
                 "Experiment model does not have suffix " + '"measurement_error".'
@@ -2714,12 +2856,11 @@ def get_unknown_parameter_values(self, model=None):
                     "`get_unknown_parameter_values`"
                 )
 
-            theta_vals = [
-                pyo.value(k)
-                for k, v in model.scenario_blocks[0].unknown_parameters.items()
-            ]
+            unknown_params = self._expanded_unknown_parameters(model.scenario_blocks[0])
+            theta_vals = [pyo.value(k) for k in unknown_params]
         else:
-            theta_vals = [pyo.value(k) for k, v in model.unknown_parameters.items()]
+            unknown_params = self._expanded_unknown_parameters(model)
+            theta_vals = [pyo.value(k) for k in unknown_params]
 
         return theta_vals
 
@@ -2788,12 +2929,17 @@ def get_measurement_error_values(self, model=None):
                     "`get_measurement_error_values`"
                 )
 
-            sigma_vals = [
-                pyo.value(k)
-                for k, v in model.scenario_blocks[0].measurement_error.items()
-            ]
+            # get the measurement-error covariance matrix
+            Sigma_matrix = self._build_meas_error_covariance_matrix(
+                model.scenario_blocks[0]
+            )
+
+            # get the measurement-error standard deviation
+            sigma_vals = np.sqrt(np.diag(Sigma_matrix))
+
         else:
-            sigma_vals = [pyo.value(k) for k, v in model.measurement_error.items()]
+            Sigma_matrix = self._build_meas_error_covariance_matrix(model)
+            sigma_vals = np.sqrt(np.diag(Sigma_matrix))
 
         return sigma_vals
 
diff --git a/pyomo/contrib/doe/tests/test_doe_build.py b/pyomo/contrib/doe/tests/test_doe_build.py
index 333ac1669ef..2dd0f7ef55e 100644
--- a/pyomo/contrib/doe/tests/test_doe_build.py
+++ b/pyomo/contrib/doe/tests/test_doe_build.py
@@ -462,7 +462,7 @@ def test_get_measurement_error_without_blocks(self):
 
         count = 0
         for k, v in doe_obj.compute_FIM_model.measurement_error.items():
-            self.assertEqual(pyo.value(k), stuff[count])
+            self.assertEqual(v, stuff[count])
             count += 1
 
     def test_get_unknown_parameters_without_blocks(self):
diff --git a/pyomo/contrib/doe/tests/test_doe_solve.py b/pyomo/contrib/doe/tests/test_doe_solve.py
index 4329ddf37ab..a96ce97f620 100644
--- a/pyomo/contrib/doe/tests/test_doe_solve.py
+++ b/pyomo/contrib/doe/tests/test_doe_solve.py
@@ -70,17 +70,79 @@ def get_rooney_biegler_data():
     return data.iloc[0]
 
 
-def get_rooney_biegler_experiment():
+def get_rooney_biegler_experiment(index_vars=False):
     """Get a fresh RooneyBieglerExperiment instance for testing.
 
     Creates a new experiment instance to ensure test isolation.
     Each test gets its own instance to avoid state sharing.
+
+    Parameters
+    ----------
+    index_vars : Boolean
+        Specifies if the Rooney-Biegler model with indexed
+        parameters is used
     """
-    return RooneyBieglerExperiment(
-        data=get_rooney_biegler_data(),
-        theta={'asymptote': 15, 'rate_constant': 0.5},
-        measure_error=0.1,
-    )
+
+    if index_vars:
+
+        class RooneyBieglerExperimentIndexed(RooneyBieglerExperiment):
+            def create_model(self):
+                data = self.data.to_frame().transpose()
+
+                model = pyo.ConcreteModel()
+
+                model.var_names = pyo.Set(initialize=["asymptote", "rate_constant"])
+                model.theta = pyo.Var(model.var_names, initialize=self.theta)
+                model.theta["asymptote"].fix()
+                model.theta["rate_constant"].fix()
+
+                # Add the experiment inputs
+                model.hour = pyo.Var(initialize=data["hour"].iloc[0], bounds=(0, 10))
+
+                # Fix the experiment inputs
+                model.hour.fix()
+
+                # Add experiment outputs
+                model.y = pyo.Var(
+                    initialize=data['y'].iloc[0], within=pyo.PositiveReals
+                )
+
+                # Define the model equations
+                def response_rule(m):
+                    return m.y == m.theta["asymptote"] * (
+                        1 - pyo.exp(-m.theta["rate_constant"] * m.hour)
+                    )
+
+                model.response_con = pyo.Constraint(rule=response_rule)
+
+                self.model = model
+
+            def label_model(self):
+                m = self.model
+
+                m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+                m.experiment_outputs.update([(m.y, self.data["y"])])
+
+                m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+                m.unknown_parameters.update((k, pyo.ComponentUID(k)) for k in [m.theta])
+
+                m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+                m.measurement_error.update([(m.y, self.measure_error)])
+
+                m.experiment_inputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+                m.experiment_inputs.update([(m.hour, self.data['hour'])])
+
+        return RooneyBieglerExperimentIndexed(
+            data=get_rooney_biegler_data(),
+            theta={'asymptote': 15, 'rate_constant': 0.5},
+            measure_error=0.1,
+        )
+    else:
+        return RooneyBieglerExperiment(
+            data=get_rooney_biegler_data(),
+            theta={'asymptote': 15, 'rate_constant': 0.5},
+            measure_error=0.1,
+        )
 
 
 def get_FIM_Q_L(doe_obj=None):
@@ -111,9 +173,8 @@ def get_FIM_Q_L(doe_obj=None):
         for i in model.output_names
         for j in model.parameter_names
     ]
-    sigma_inv = [
-        1 / v**2 for k, v in model.scenario_blocks[0].measurement_error.items()
-    ]
+    sigma_inv = doe_obj.get_meas_error_covariance_matrix_inv(model.scenario_blocks[0])
+
     FIM_vals_np = np.array(FIM_vals).reshape((n_param, n_param))
 
     for i in range(n_param):
@@ -124,12 +185,7 @@ def get_FIM_Q_L(doe_obj=None):
     L_vals_np = np.array(L_vals).reshape((n_param, n_param))
     Q_vals_np = np.array(Q_vals).reshape((n_y, n_param))
 
-    sigma_inv_np = np.zeros((n_y, n_y))
-
-    for ind, v in enumerate(sigma_inv):
-        sigma_inv_np[ind, ind] = v
-
-    return FIM_vals_np, Q_vals_np, L_vals_np, sigma_inv_np
+    return FIM_vals_np, Q_vals_np, L_vals_np, sigma_inv
 
 
 def get_standard_args(experiment, fd_method, obj_used):
@@ -168,23 +224,32 @@ def test_rooney_biegler_fd_central_solve(self):
 
         # Use RooneyBiegler for algorithm validation (faster)
         experiment = get_rooney_biegler_experiment()
+        ind_param_experiment = get_rooney_biegler_experiment(index_vars=True)
 
         DoE_args = get_standard_args(experiment, fd_method, obj_used)
+        ind_param_DoE_args = get_standard_args(
+            ind_param_experiment, fd_method, obj_used
+        )
 
         doe_obj = DesignOfExperiments(**DoE_args)
+        ind_param_doe_obj = DesignOfExperiments(**ind_param_DoE_args)
 
         doe_obj.run_doe()
+        ind_param_doe_obj.run_doe()
 
         # assert model solves
         self.assertEqual(doe_obj.results["Solver Status"], "ok")
+        self.assertEqual(ind_param_doe_obj.results["Solver Status"], "ok")
 
         # assert that Q, F, and L are the same.
         FIM, Q, L, sigma_inv = get_FIM_Q_L(doe_obj=doe_obj)
+        FIM_ind, Q_ind, L_ind, sigma_inv_ind = get_FIM_Q_L(doe_obj=ind_param_doe_obj)
 
         # Since Trace is used, no comparison for FIM and L.T @ L
 
         # Make sure FIM and Q.T @ sigma_inv @ Q are close (alternate definition of FIM)
         self.assertTrue(np.all(np.isclose(FIM, Q.T @ sigma_inv @ Q)))
+        self.assertTrue(np.all(np.isclose(FIM_ind, Q_ind.T @ sigma_inv_ind @ Q_ind)))
 
     @unittest.skipIf(not pandas_available, "pandas is not available")
     def test_rooney_biegler_fd_forward_solve(self):
@@ -230,22 +295,31 @@ def test_rooney_biegler_fd_backward_solve(self):
 
         # Use RooneyBiegler for algorithm validation (faster)
         experiment = get_rooney_biegler_experiment()
+        ind_param_experiment = get_rooney_biegler_experiment(index_vars=True)
 
         DoE_args = get_standard_args(experiment, fd_method, obj_used)
+        ind_param_DoE_args = get_standard_args(
+            ind_param_experiment, fd_method, obj_used
+        )
 
         doe_obj = DesignOfExperiments(**DoE_args)
+        ind_param_doe_obj = DesignOfExperiments(**ind_param_DoE_args)
 
         doe_obj.run_doe()
+        ind_param_doe_obj.run_doe()
 
         self.assertEqual(doe_obj.results["Solver Status"], "ok")
+        self.assertEqual(ind_param_doe_obj.results["Solver Status"], "ok")
 
         # assert that Q, F, and L are the same.
         FIM, Q, L, sigma_inv = get_FIM_Q_L(doe_obj=doe_obj)
+        FIM_ind, Q_ind, L_ind, sigma_inv_ind = get_FIM_Q_L(doe_obj=ind_param_doe_obj)
 
         # Since Trace is used, no comparison for FIM and L.T @ L
 
         # Make sure FIM and Q.T @ sigma_inv @ Q are close (alternate definition of FIM)
         self.assertTrue(np.all(np.isclose(FIM, Q.T @ sigma_inv @ Q)))
+        self.assertTrue(np.all(np.isclose(FIM_ind, Q_ind.T @ sigma_inv_ind @ Q_ind)))
 
     @unittest.skipIf(not pandas_available, "pandas is not available")
     def test_rooney_biegler_obj_det_solve(self):
@@ -254,27 +328,43 @@ def test_rooney_biegler_obj_det_solve(self):
 
         # Use RooneyBiegler for algorithm validation (faster)
         experiment = get_rooney_biegler_experiment()
+        ind_param_experiment = get_rooney_biegler_experiment(index_vars=True)
 
         DoE_args = get_standard_args(experiment, fd_method, obj_used)
+        ind_param_DoE_args = get_standard_args(
+            ind_param_experiment, fd_method, obj_used
+        )
+
         DoE_args["scale_nominal_param_value"] = (
             False  # Vanilla determinant solve needs this
         )
         DoE_args["_Cholesky_option"] = False
         DoE_args["_only_compute_fim_lower"] = False
 
+        ind_param_DoE_args["scale_nominal_param_value"] = False
+        ind_param_DoE_args["_Cholesky_option"] = False
+        ind_param_DoE_args["_only_compute_fim_lower"] = False
+
         doe_obj = DesignOfExperiments(**DoE_args)
+        ind_param_doe_obj = DesignOfExperiments(**ind_param_DoE_args)
 
         # Increase numerical performance by adding a prior
         prior_FIM = doe_obj.compute_FIM()
+        prior_FIM_ind = ind_param_doe_obj.compute_FIM()
         doe_obj.prior_FIM = prior_FIM
+        ind_param_doe_obj.prior_FIM = prior_FIM_ind
 
         doe_obj.run_doe()
+        ind_param_doe_obj.run_doe()
 
         self.assertEqual(doe_obj.results["Solver Status"], "ok")
+        self.assertEqual(ind_param_doe_obj.results["Solver Status"], "ok")
 
         expected_design = 9.999213890476453
         actual_design = doe_obj.results["Experiment Design"][0]
+        actual_design_ind = ind_param_doe_obj.results["Experiment Design"][0]
         self.assertAlmostEqual(actual_design, expected_design, places=3)
+        self.assertAlmostEqual(actual_design_ind, expected_design, places=3)
 
     @unittest.skipIf(not pandas_available, "pandas is not available")
     def test_rooney_biegler_obj_cholesky_solve(self):
@@ -283,29 +373,44 @@ def test_rooney_biegler_obj_cholesky_solve(self):
 
         # Use RooneyBiegler for algorithm validation (faster)
         experiment = get_rooney_biegler_experiment()
+        ind_param_experiment = get_rooney_biegler_experiment(index_vars=True)
 
         DoE_args = get_standard_args(experiment, fd_method, obj_used)
+        ind_param_DoE_args = get_standard_args(
+            ind_param_experiment, fd_method, obj_used
+        )
 
         # Add prior FIM for better numerical conditioning
         # This follows the pattern in rooney_biegler_doe_example.py
         doe_obj_prior = DesignOfExperiments(**DoE_args)
+        ind_doe_obj_prior = DesignOfExperiments(**ind_param_DoE_args)
         prior_FIM = doe_obj_prior.compute_FIM()
+        prior_FIM_ind = ind_doe_obj_prior.compute_FIM()
         DoE_args['prior_FIM'] = prior_FIM
+        ind_param_DoE_args['prior_FIM'] = prior_FIM_ind
 
         doe_obj = DesignOfExperiments(**DoE_args)
+        ind_param_doe_obj = DesignOfExperiments(**ind_param_DoE_args)
 
         doe_obj.run_doe()
+        ind_param_doe_obj.run_doe()
 
         self.assertEqual(doe_obj.results["Solver Status"], "ok")
+        self.assertEqual(ind_param_doe_obj.results["Solver Status"], "ok")
 
         # assert that Q, F, and L are the same.
         FIM, Q, L, sigma_inv = get_FIM_Q_L(doe_obj=doe_obj)
+        FIM_ind, Q_ind, L_ind, sigma_inv_ind = get_FIM_Q_L(doe_obj=ind_param_doe_obj)
 
         # Since Cholesky is used, there is comparison for FIM and L.T @ L
         self.assertTrue(np.all(np.isclose(FIM, L @ L.T)))
+        self.assertTrue(np.all(np.isclose(FIM_ind, L_ind @ L_ind.T)))
 
         # Note: When using prior_FIM, the relationship FIM = Q.T @ sigma_inv @ Q + prior_FIM
         self.assertTrue(np.all(np.isclose(FIM, Q.T @ sigma_inv @ Q + prior_FIM)))
+        self.assertTrue(
+            np.all(np.isclose(FIM_ind, Q_ind.T @ sigma_inv_ind @ Q_ind + prior_FIM_ind))
+        )
 
     def DISABLE_test_reactor_obj_cholesky_solve_bad_prior(self):
         # [10/2025] This test has been disabled because it frequently
@@ -349,10 +454,15 @@ def test_compute_FIM_seq_centr(self):
 
         # Use RooneyBiegler for algorithm validation (faster)
         experiment = get_rooney_biegler_experiment()
+        ind_param_experiment = get_rooney_biegler_experiment(index_vars=True)
 
         DoE_args = get_standard_args(experiment, fd_method, obj_used)
+        ind_param_DoE_args = get_standard_args(
+            ind_param_experiment, fd_method, obj_used
+        )
 
         doe_obj = DesignOfExperiments(**DoE_args)
+        ind_param_doe_obj = DesignOfExperiments(**ind_param_DoE_args)
 
         expected_FIM = np.array(
             [[18957.7788694, 4238.27606876], [4238.27606876, 947.52577076]]
@@ -361,6 +471,13 @@ def test_compute_FIM_seq_centr(self):
         self.assertTrue(
             np.all(np.isclose(doe_obj.compute_FIM(method="sequential"), expected_FIM))
         )
+        self.assertTrue(
+            np.all(
+                np.isclose(
+                    ind_param_doe_obj.compute_FIM(method="sequential"), expected_FIM
+                )
+            )
+        )
 
     # This test ensure that compute FIM runs without error using the
     # `sequential` option with forward finite differences
@@ -371,12 +488,18 @@ def test_compute_FIM_seq_forward(self):
 
         # Use RooneyBiegler for algorithm validation (faster)
         experiment = get_rooney_biegler_experiment()
+        ind_param_experiment = get_rooney_biegler_experiment(index_vars=True)
 
         DoE_args = get_standard_args(experiment, fd_method, obj_used)
+        ind_param_DoE_args = get_standard_args(
+            ind_param_experiment, fd_method, obj_used
+        )
 
         doe_obj = DesignOfExperiments(**DoE_args)
+        ind_param_doe_obj = DesignOfExperiments(**ind_param_DoE_args)
 
         doe_obj.compute_FIM(method="sequential")
+        ind_param_doe_obj.compute_FIM(method="sequential")
 
     # This test ensure that compute FIM runs without error using the
     # `kaug` option. kaug computes the FIM directly so no finite difference
@@ -390,10 +513,15 @@ def test_compute_FIM_kaug(self):
         obj_used = "determinant"
 
         experiment = get_rooney_biegler_experiment()
+        ind_param_experiment = get_rooney_biegler_experiment(index_vars=True)
 
         DoE_args = get_standard_args(experiment, fd_method, obj_used)
+        ind_param_DoE_args = get_standard_args(
+            ind_param_experiment, fd_method, obj_used
+        )
 
         doe_obj = DesignOfExperiments(**DoE_args)
+        ind_param_doe_obj = DesignOfExperiments(**ind_param_DoE_args)
 
         expected_FIM = np.array(
             [[18957.7788694, 4238.27606876], [4238.27606876, 947.52577076]]
@@ -402,6 +530,11 @@ def test_compute_FIM_kaug(self):
         self.assertTrue(
             np.all(np.isclose(doe_obj.compute_FIM(method="kaug"), expected_FIM))
         )
+        self.assertTrue(
+            np.all(
+                np.isclose(ind_param_doe_obj.compute_FIM(method="kaug"), expected_FIM)
+            )
+        )
 
     # This test ensure that compute FIM runs without error using the
     # `sequential` option with backward finite differences
@@ -412,12 +545,18 @@ def test_compute_FIM_seq_backward(self):
 
         # Use RooneyBiegler for algorithm validation (faster)
         experiment = get_rooney_biegler_experiment()
+        ind_param_experiment = get_rooney_biegler_experiment(index_vars=True)
 
         DoE_args = get_standard_args(experiment, fd_method, obj_used)
+        ind_param_DoE_args = get_standard_args(
+            ind_param_experiment, fd_method, obj_used
+        )
 
         doe_obj = DesignOfExperiments(**DoE_args)
+        ind_param_doe_obj = DesignOfExperiments(**ind_param_DoE_args)
 
         doe_obj.compute_FIM(method="sequential")
+        ind_param_doe_obj.compute_FIM(method="sequential")
 
     @unittest.skipIf(not pandas_available, "pandas is not available")
     def test_reactor_grid_search(self):
diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py
index 7f05b19b9ba..4b710fe154a 100644
--- a/pyomo/contrib/parmest/parmest.py
+++ b/pyomo/contrib/parmest/parmest.py
@@ -103,8 +103,12 @@ def SSE(model):
 def SSE_weighted(model):
     """
     Returns an expression that is used to compute the 'SSE_weighted' objective,
-    assuming Gaussian i.i.d. errors, with measurement error standard deviation
-    defined in the annotated Pyomo model
+    assuming Gaussian correlated or i.i.d. errors, with the standard deviation
+    or covariance of the measurement errors defined in the annotated Pyomo model
+
+    This objective function is applicable to both homoscedastic
+    (constant-variance) and heteroskedastic (non-constant-variance, e.g.,
+    proportional-error) measurement-error models
 
     Parameters
     ----------
@@ -122,26 +126,34 @@ def SSE_weighted(model):
             'objective.'
         )
 
-    # check if all the values of the measurement error standard deviation
+    # check if all the values of the measurement-error standard deviation
     # have been supplied
-    all_known_errors = all(
-        model.measurement_error[y_hat] is not None for y_hat in model.experiment_outputs
-    )
+    try:
+        all_known_errors = all(
+            model.measurement_error[y_hat] is not None
+            for y_hat in model.experiment_outputs
+        )
+    except KeyError:
+        raise KeyError(
+            'One or more experiment outputs are not defined in the '
+            '"measurement_error" attribute. All the variables defined '
+            'in "experiment_outputs" must be defined as keys in '
+            '"measurement_error".'
+        )
 
     if all_known_errors:
+        # calculate the residuals between the model predictions and data
+        prediction_resid = np.array(
+            [y - y_hat for y_hat, y in model.experiment_outputs.items()]
+        ).reshape(1, -1)
+
+        # get the inverse of the measurement-error covariance matrix
+        Sigma_y_inv = get_meas_error_covariance_matrix_inv(model)
+
         # calculate the weighted SSE between the prediction
         # and observation of the measured variables
-        try:
-            expr = (1 / 2) * sum(
-                ((y - y_hat) / model.measurement_error[y_hat]) ** 2
-                for y_hat, y in model.experiment_outputs.items()
-            )
-            return expr
-        except ZeroDivisionError:
-            raise ValueError(
-                'Division by zero encountered in the "SSE_weighted" objective. '
-                'One or more values of the measurement error are zero.'
-            )
+        expr = (1 / 2) * prediction_resid @ Sigma_y_inv @ prediction_resid.T
+        return expr[0, 0]
     else:
         raise ValueError(
             'One or more values are missing from "measurement_error". All values of '
@@ -149,6 +161,127 @@ def SSE_weighted(model):
         )
 
 
+def _build_meas_error_covariance_matrix(model, estimated_var=None):
+    """
+    Builds the full measurement-error covariance matrix
+
+    Note: The code does not automatically build the full covariance matrix
+    It only places whatever covariances the user explicitly supplies
+    Therefore, for correlation in time, shared timepoints, or other types of
+    correlation, the user must provide all the desired covariance terms
+    The diagonal elements can be constructed automatically or the standard
+    deviations can be specified by the user. They standard deviations may be
+    constant or depend on the value (i.e., data) of the measured or input
+    variables (e.g., be proportional to them)
+
+    Parameters
+    ----------
+    model : ConcreteModel
+        Annotated Pyomo model
+    estimated_var: float or int, optional
+        Value of the estimated variance of the measurement error
+        in cases where the user does not supply the
+        measurement-error standard deviation
+
+    Returns
+    -------
+    Sigma_y: numpy.ndarray
+        Full measurement-error covariance matrix
+    """
+    # get the output variables
+    outputs = list(model.experiment_outputs.keys())
+    outputs_name = [y_hat.name for y_hat in outputs]
+    outputs_index = {y_hat_name: i for i, y_hat_name in enumerate(outputs_name)}
+
+    # get the number of output variables
+    number_outputs = len(outputs)
+
+    # define the measurement-error covariance matrix
+    Sigma_y = np.zeros((number_outputs, number_outputs))
+
+    if hasattr(model, "measurement_error"):
+        # check if all the measurement-error standard deviations
+        # have been supplied
+        all_known_errors = all(
+            model.measurement_error[y_hat] is not None
+            for y_hat in model.experiment_outputs
+        )
+
+        # fill the leading-diagonal elements from the standard deviation of
+        # the measurement errors
+        for y_hat in outputs:
+            # get the index of y_hat
+            i = outputs_index[y_hat.name]
+
+            if all_known_errors:
+                standard_dev = model.measurement_error[y_hat]
+                Sigma_y[i, i] = standard_dev**2
+            else:
+                Sigma_y[i, i] = estimated_var
+
+        # fill the off-diagonal elements from covariance entries
+        # supplied by the user
+        for key, err_cov in model.measurement_error.items():
+            if isinstance(key, tuple) and len(key) == 2:
+                yi, yj = key
+                if yi.name not in outputs_name or yj.name not in outputs_name:
+                    raise ValueError(
+                        "Measurement-error covariance must be defined only between "
+                        "experiment output variables."
+                    )
+
+                # get the indices of yi and yj
+                i = outputs_index[yi.name]
+                j = outputs_index[yj.name]
+
+                # update the measurement-error covariance matrix which
+                # is a symmetric matrix
+                Sigma_y[i, j] = err_cov
+                Sigma_y[j, i] = err_cov
+            elif not isinstance(key, tuple) and not hasattr(key, "name"):
+                raise TypeError(
+                    "Expected a tuple of two measured variables when specifying a "
+                    "measurement-error covariance, e.g., "
+                    "measurement_error[(y1, y2)] = covariance."
+                )
+
+    return Sigma_y
+
+
+def get_meas_error_covariance_matrix_inv(model, estimated_var=None):
+    """
+    Computes the inverse of the measurement-error covariance matrix
+
+    Parameters
+    ----------
+    model : ConcreteModel
+        Annotated Pyomo model
+    estimated_var: float or int, optional
+        Value of the estimated variance of the measurement error
+        in cases where the user does not supply the
+        measurement-error standard deviation
+
+    Returns
+    -------
+    Sigma_y_inv: numpy.ndarray
+        Inverse of the measurement-error covariance matrix
+    """
+    # get the measurement-error covariance matrix
+    Sigma_y = _build_meas_error_covariance_matrix(model, estimated_var)
+
+    # compute the inverse of the measurement-error covariance matrix
+    try:
+        Sigma_y_inv = np.linalg.inv(Sigma_y)
+    except np.linalg.LinAlgError:
+        Sigma_y_inv = np.linalg.pinv(Sigma_y)
+        logger.warning(
+            "The measurement-error covariance matrix is singular. "
+            "Using pseudo-inverse instead."
+        )
+
+    return Sigma_y_inv
+
+
 def _validate_prior_FIM(prior_FIM, require_psd=True):
     """
     Validate user-supplied prior Fisher Information Matrix.
@@ -570,8 +703,11 @@ def _compute_jacobian(experiment, theta_vals, step, solver, tee, solver_options)
     # grab the model
     model = _get_labeled_model(experiment)
 
+    # get the parameter data objects
+    param_data_objects, _, _ = _expanded_unknown_parameter_info(model)
+
     # fix the value of the unknown parameters to the estimated values
-    for param in model.unknown_parameters:
+    for param in param_data_objects:
         param.fix(theta_vals[param.name])
 
     # re-solve the model with the estimated parameters
@@ -585,7 +721,7 @@ def _compute_jacobian(experiment, theta_vals, step, solver, tee, solver_options)
     assert_optimal_termination(results)
 
     # get the estimated parameter values
-    param_values = [p.value for p in model.unknown_parameters]
+    param_values = [p.value for p in param_data_objects]
 
     # get the number of parameters and measured variables
     n_params = len(param_values)
@@ -594,7 +730,7 @@ def _compute_jacobian(experiment, theta_vals, step, solver, tee, solver_options)
     # compute the sensitivity of the measured variables w.r.t the parameters
     J = np.zeros((n_outputs, n_params))
 
-    for i, param in enumerate(model.unknown_parameters):
+    for i, param in enumerate(param_data_objects):
         # store original value of the parameter
         orig_value = param_values[i]
 
@@ -813,36 +949,12 @@ def _finite_difference_FIM(
     # grab the model
     model = _get_labeled_model(experiment)
 
-    # extract the measured variables and measurement errors
-    y_hat_list = [y_hat for y_hat, y in model.experiment_outputs.items()]
-
-    # check if the model has a 'measurement_error' attribute and
-    # the measurement error standard deviation has been supplied
-    all_known_errors = all(
-        model.measurement_error[y_hat] is not None for y_hat in model.experiment_outputs
-    )
-
-    if hasattr(model, "measurement_error") and all_known_errors:
-        error_list = [
-            model.measurement_error[y_hat] for y_hat in model.experiment_outputs
-        ]
-
-        # check if the dimension of error_list is the same with that of y_hat_list
-        if len(error_list) != len(y_hat_list):
-            raise ValueError(
-                "Experiment outputs and measurement errors are not the same length."
-            )
-
-        # compute the matrix of the inverse of the measurement error variance
-        # the following assumes independent and identically distributed
-        # measurement errors
-        W = np.diag([1 / (err**2) for err in error_list])
+    # get the inverse of the measurement-error covariance matrix
+    Sigma_y_inv = get_meas_error_covariance_matrix_inv(model, estimated_var)
 
-        # calculate the FIM using the formula in our future paper
-        # Lilonfe et al. (2025)
-        FIM = J.T @ W @ J
-    else:
-        FIM = (1 / estimated_var) * (J.T @ J)
+    # calculate the FIM using the formula in our future paper
+    # Lilonfe and Dowling. (2026)
+    FIM = J.T @ Sigma_y_inv @ J
 
     return FIM
 
@@ -900,8 +1012,11 @@ def _kaug_FIM(
     # add the built-in objective function selected by the user
     model.objective = pyo.Objective(expr=obj_function, sense=pyo.minimize)
 
+    # get the parameter data objects
+    param_data_objects, _, _ = _expanded_unknown_parameter_info(model)
+
     # fix the parameter values to the estimated values
-    for param in model.unknown_parameters:
+    for param in param_data_objects:
         param.fix(theta_vals[param.name])
 
     solver = pyo.SolverFactory(solver)
@@ -912,7 +1027,7 @@ def _kaug_FIM(
     assert_optimal_termination(results)
 
     # Probe the solved model for dsdp results (sensitivities s.t. parameters)
-    params_dict = {k.name: v for k, v in model.unknown_parameters.items()}
+    params_dict = {k.name: theta_vals[k.name] for k in param_data_objects}
     params_names = list(params_dict.keys())
 
     dsdp_re, col = get_dsdp(model, params_names, params_dict, tee=tee)
@@ -946,7 +1061,7 @@ def _kaug_FIM(
     jac = [[] for _ in params_names]
 
     for d in range(len(dsdp_extract)):
-        for k, v in model.unknown_parameters.items():
+        for k in param_data_objects:
             p = params_names.index(k.name)  # Index of parameter in np array
             sensi = dsdp_extract[d][p]
             jac[p].append(sensi)
@@ -954,24 +1069,11 @@ def _kaug_FIM(
     # record kaug jacobian
     kaug_jac = np.array(jac).T
 
-    # compute FIM
-    # compute the matrix of the inverse of the measurement error variance
-    # the following assumes independent and identically distributed
-    # measurement errors
-    W = np.zeros((len(model.measurement_error), len(model.measurement_error)))
-    all_known_errors = all(
-        model.measurement_error[y_hat] is not None for y_hat in model.experiment_outputs
-    )
-
-    count = 0
-    for k, v in model.measurement_error.items():
-        if all_known_errors:
-            W[count, count] = 1 / (v**2)
-        else:
-            W[count, count] = 1 / estimated_var
-        count += 1
+    # get the inverse of the measurement-error covariance matrix
+    Sigma_y_inv = get_meas_error_covariance_matrix_inv(model, estimated_var)
 
-    FIM = kaug_jac.T @ W @ kaug_jac
+    # compute the FIM
+    FIM = kaug_jac.T @ Sigma_y_inv @ kaug_jac
 
     return FIM
 
@@ -1698,13 +1800,12 @@ def _cov_at_theta(self, method, solver, step):
             for experiment in self.exp_list:
                 model = _get_labeled_model(experiment)
 
+                # get the parameter data objects
+                param_data_objects, _, _ = _expanded_unknown_parameter_info(model)
+
                 # fix the value of the unknown parameters to the estimated values
-                for param in model.unknown_parameters:
-                    if param.is_indexed():
-                        for idx in param:
-                            param[idx].fix(self.estimated_theta[param[idx].name])
-                    else:
-                        param.fix(self.estimated_theta[param.name])
+                for param in param_data_objects:
+                    param.fix(self.estimated_theta[param.name])
 
                 # re-solve the model with the estimated parameters
                 results = pyo.SolverFactory(solver).solve(model, tee=self.tee)
@@ -1763,10 +1864,18 @@ def _cov_at_theta(self, method, solver, step):
             # check if the user defined the 'measurement_error' attribute
             if hasattr(ref_model, "measurement_error"):
                 # get the measurement errors
-                meas_error = [
-                    ref_model.measurement_error[y_hat]
-                    for y_hat, y in ref_model.experiment_outputs.items()
-                ]
+                try:
+                    meas_error = [
+                        ref_model.measurement_error[y_hat]
+                        for y_hat, y in ref_model.experiment_outputs.items()
+                    ]
+                except KeyError:
+                    raise KeyError(
+                        'One or more experiment outputs are not defined in the '
+                        '"measurement_error" attribute. All the variables defined '
+                        'in "experiment_outputs" must be defined as keys in '
+                        '"measurement_error".'
+                    )
 
                 # check if the user supplied the values of the measurement errors
                 if all(item is None for item in meas_error):
diff --git a/pyomo/contrib/parmest/tests/test_parmest.py b/pyomo/contrib/parmest/tests/test_parmest.py
index af8fc9ca757..7faa668d3b3 100644
--- a/pyomo/contrib/parmest/tests/test_parmest.py
+++ b/pyomo/contrib/parmest/tests/test_parmest.py
@@ -146,10 +146,7 @@ def check_rooney_biegler_covariance(
         ][0]
 
         if measurement_error is None and obj_function == "SSE":
-            if (
-                cov_method == "finite_difference"
-                or cov_method == "automatic_differentiation_kaug"
-            ):
+            if cov_method in ("finite_difference", "automatic_differentiation_kaug"):
                 self.assertAlmostEqual(
                     cov.iloc[asymptote_index, asymptote_index], 6.229612, places=2
                 )  # 6.22864 from paper
@@ -180,10 +177,7 @@ def check_rooney_biegler_covariance(
                     places=2,
                 )  # 0.04124 from paper
         elif measurement_error is not None and obj_function in ("SSE", "SSE_weighted"):
-            if (
-                cov_method == "finite_difference"
-                or cov_method == "automatic_differentiation_kaug"
-            ):
+            if cov_method in ("finite_difference", "automatic_differentiation_kaug"):
                 self.assertAlmostEqual(
                     cov.iloc[asymptote_index, asymptote_index], 0.009588, places=4
                 )
@@ -845,9 +839,9 @@ def label_model(self):
                 m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL)
                 m.measurement_error.update([(m.y, None)])
 
-        rooney_biegler_indexed_vars_exp_list = []
+        self.rooney_biegler_indexed_vars_exp_list = []
         for i in range(self.data.shape[0]):
-            rooney_biegler_indexed_vars_exp_list.append(
+            self.rooney_biegler_indexed_vars_exp_list.append(
                 RooneyBieglerExperimentIndexedVars(self.data.loc[i, :])
             )
 
@@ -875,24 +869,24 @@ def label_model(self):
                 "theta_vals": theta_vals,
             },
             "vars_index": {
-                "exp_list": rooney_biegler_indexed_vars_exp_list,
+                "exp_list": self.rooney_biegler_indexed_vars_exp_list,
                 "theta_names": ["theta"],
                 "theta_vals": theta_vals_index,
             },
             "vars_quoted_index": {
-                "exp_list": rooney_biegler_indexed_vars_exp_list,
+                "exp_list": self.rooney_biegler_indexed_vars_exp_list,
                 "theta_names": ["theta['asymptote']", "theta['rate_constant']"],
                 "theta_vals": theta_vals_index,
             },
             "vars_str_index": {
-                "exp_list": rooney_biegler_indexed_vars_exp_list,
+                "exp_list": self.rooney_biegler_indexed_vars_exp_list,
                 "theta_names": ["theta[asymptote]", "theta[rate_constant]"],
                 "theta_vals": theta_vals_index,
             },
         }
 
     @unittest.skipIf(not pynumero_ASL_available, "pynumero_ASL is not available")
-    def check_rooney_biegler_results(self, objval, cov):
+    def check_rooney_biegler_results(self, objval, cov, cov_method="reduced_hessian"):
 
         # get indices in covariance matrix
         cov_cols = cov.columns.to_list()
@@ -902,18 +896,33 @@ def check_rooney_biegler_results(self, objval, cov):
         ][0]
 
         self.assertAlmostEqual(objval, 4.3317112, places=2)
-        self.assertAlmostEqual(
-            cov.iloc[asymptote_index, asymptote_index], 6.155892, places=2
-        )  # 6.22864 from paper
-        self.assertAlmostEqual(
-            cov.iloc[asymptote_index, rate_constant_index], -0.425232, places=2
-        )  # -0.4322 from paper
-        self.assertAlmostEqual(
-            cov.iloc[rate_constant_index, asymptote_index], -0.425232, places=2
-        )  # -0.4322 from paper
-        self.assertAlmostEqual(
-            cov.iloc[rate_constant_index, rate_constant_index], 0.040571, places=2
-        )  # 0.04124 from paper
+
+        if cov_method in ("finite_difference", "automatic_differentiation_kaug"):
+            self.assertAlmostEqual(
+                cov.iloc[asymptote_index, asymptote_index], 6.229612, places=2
+            )  # 6.22864 from paper
+            self.assertAlmostEqual(
+                cov.iloc[asymptote_index, rate_constant_index], -0.432265, places=2
+            )  # -0.4322 from paper
+            self.assertAlmostEqual(
+                cov.iloc[rate_constant_index, asymptote_index], -0.432265, places=2
+            )  # -0.4322 from paper
+            self.assertAlmostEqual(
+                cov.iloc[rate_constant_index, rate_constant_index], 0.041242, places=2
+            )  # 0.04124 from paper
+        else:
+            self.assertAlmostEqual(
+                cov.iloc[asymptote_index, asymptote_index], 6.155892, places=2
+            )  # 6.22864 from paper
+            self.assertAlmostEqual(
+                cov.iloc[asymptote_index, rate_constant_index], -0.425232, places=2
+            )  # -0.4322 from paper
+            self.assertAlmostEqual(
+                cov.iloc[rate_constant_index, asymptote_index], -0.425232, places=2
+            )  # -0.4322 from paper
+            self.assertAlmostEqual(
+                cov.iloc[rate_constant_index, rate_constant_index], 0.040571, places=2
+            )  # 0.04124 from paper
 
     @unittest.skipUnless(pynumero_ASL_available, 'pynumero_ASL is not available')
     def test_parmest_basics(self):
@@ -985,6 +994,32 @@ def test_parmest_basics_with_square_problem_solve_no_theta_vals(self):
             cov = pest.cov_est(method="reduced_hessian")
             self.check_rooney_biegler_results(objval, cov)
 
+    @unittest.skipUnless(pynumero_ASL_available, 'pynumero_ASL is not available')
+    def test_parmest_indexed_vars_finite_difference_cov(self):
+
+        pest = parmest.Estimator(
+            self.rooney_biegler_indexed_vars_exp_list,
+            obj_function=self.objective_function,
+        )
+
+        objval, thetavals = pest.theta_est()
+        cov_method = "finite_difference"
+        cov = pest.cov_est(method=cov_method)
+        self.check_rooney_biegler_results(objval, cov, cov_method)
+
+    @unittest.skipUnless(pynumero_ASL_available, 'pynumero_ASL is not available')
+    def test_parmest_indexed_vars_auto_differentiation_cov(self):
+
+        pest = parmest.Estimator(
+            self.rooney_biegler_indexed_vars_exp_list,
+            obj_function=self.objective_function,
+        )
+
+        objval, thetavals = pest.theta_est()
+        cov_method = "automatic_differentiation_kaug"
+        cov = pest.cov_est(method=cov_method)
+        self.check_rooney_biegler_results(objval, cov, cov_method)
+
 
 @unittest.skipIf(
     not parmest.parmest_available,
@@ -1034,6 +1069,54 @@ def setUp(self):
             exp_list, obj_function="SSE", solver_options=solver_options
         )
 
+        # create an inherited class to test the results when a
+        # full measurement-error covariance matrix is used
+        class ReactorFullErrorCov(ReactorDesignExperiment):
+            def label_model(self):
+                m = self.model
+
+                m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+                m.experiment_outputs.update(
+                    [
+                        (m.ca, self.data_i['ca']),
+                        (m.cb, self.data_i['cb']),
+                        (m.cc, self.data_i['cc']),
+                        (m.cd, self.data_i['cd']),
+                    ]
+                )
+
+                m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+                m.unknown_parameters.update(
+                    (k, pyo.ComponentUID(k)) for k in [m.k1, m.k2, m.k3]
+                )
+
+                m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+                m.measurement_error.update(
+                    [(m.ca, 1), (m.cb, 0.1), (m.cc, 0.1), (m.cd, 0.1)]
+                )
+                m.measurement_error.update(
+                    [
+                        ((m.ca, m.cb), 0),
+                        ((m.ca, m.cc), 0),
+                        ((m.ca, m.cd), 0),
+                        ((m.cb, m.cc), 0),
+                        ((m.cb, m.cd), 0),
+                        ((m.cc, m.cd), 0),
+                    ]
+                )
+
+                return m
+
+        exp_list_full_error = []
+        for i in range(data.shape[0]):
+            exp_list_full_error.append(ReactorFullErrorCov(data, i))
+
+        self.pest_full_error = parmest.Estimator(
+            exp_list_full_error,
+            obj_function="SSE_weighted",
+            solver_options=solver_options,
+        )
+
     def test_theta_est(self):
         # used in data reconciliation
         objval, thetavals = self.pest.theta_est()
@@ -1042,6 +1125,12 @@ def test_theta_est(self):
         self.assertAlmostEqual(thetavals["k2"], 5.0 / 3.0, places=4)
         self.assertAlmostEqual(thetavals["k3"], 1.0 / 6000.0, places=7)
 
+        objval2, thetavals2 = self.pest_full_error.theta_est()
+
+        self.assertAlmostEqual(thetavals2["k1"], 5.0 / 6.0, places=4)
+        self.assertAlmostEqual(thetavals2["k2"], 5.0 / 3.0, places=4)
+        self.assertAlmostEqual(thetavals2["k3"], 1.0 / 6000.0, places=7)
+
     def test_return_values(self):
         objval, thetavals, data_rec = self.pest.theta_est(
             return_values=["ca", "cb", "cc", "cd", "caf"]
@@ -1251,7 +1340,7 @@ def get_labeled_model(self):
 
         # create an instance of the ReactorDesignExperimentDAE class
         # without the "unknown_parameters" attribute
-        class ReactorDesignExperimentException(ReactorDesignExperimentDAE):
+        class ReactorDesignParameterException(ReactorDesignExperimentDAE):
             def label_model(self):
 
                 m = self.model
@@ -1273,11 +1362,118 @@ def label_model(self):
                 )
 
         # create an experiment list without the "unknown_parameters" attribute
-        exp_list_df_no_params = [ReactorDesignExperimentException(data_df)]
-        exp_list_dict_no_params = [ReactorDesignExperimentException(data_dict)]
+        self.exp_list_df_no_params = [ReactorDesignParameterException(data_df)]
+        self.exp_list_dict_no_params = [ReactorDesignParameterException(data_dict)]
+
+        # create instances of the ReactorDesignExperimentDAE class
+        # with incorrect definition of the measurement-error covariance
+        class ReactorErrorCovarianceException1(ReactorDesignExperimentDAE):
+            def label_model(self):
+
+                m = self.model
+
+                if isinstance(self.data, pd.DataFrame):
+                    meas_time_points = self.data.index
+                else:
+                    meas_time_points = list(self.data["ca"].keys())
+
+                m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+                m.experiment_outputs.update(
+                    (m.ca[t], self.data["ca"][t]) for t in meas_time_points
+                )
+                m.experiment_outputs.update(
+                    (m.cb[t], self.data["cb"][t]) for t in meas_time_points
+                )
+                m.experiment_outputs.update(
+                    (m.cc[t], self.data["cc"][t]) for t in meas_time_points
+                )
+
+                m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+                m.unknown_parameters.update(
+                    (k, pyo.ComponentUID(k)) for k in [m.k1, m.k2]
+                )
+
+                m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+                m.measurement_error.update((m.ca[t], 0.01) for t in meas_time_points)
+                m.measurement_error.update((m.cb[t], 0.01) for t in meas_time_points)
+                m.measurement_error.update((m.cc[t], 0.01) for t in meas_time_points)
+                m.measurement_error.update(
+                    ([m.ca[t], m.cb[t]], 0.1) for t in meas_time_points
+                )
 
-        self.exp_list_df_no_params = exp_list_df_no_params
-        self.exp_list_dict_no_params = exp_list_dict_no_params
+        class ReactorErrorCovarianceException2(ReactorDesignExperimentDAE):
+            def label_model(self):
+
+                m = self.model
+
+                if isinstance(self.data, pd.DataFrame):
+                    meas_time_points = self.data.index
+                else:
+                    meas_time_points = list(self.data["ca"].keys())
+
+                m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+                m.experiment_outputs.update(
+                    (m.ca[t], self.data["ca"][t]) for t in meas_time_points
+                )
+                m.experiment_outputs.update(
+                    (m.cb[t], self.data["cb"][t]) for t in meas_time_points
+                )
+                m.experiment_outputs.update(
+                    (m.cc[t], self.data["cc"][t]) for t in meas_time_points
+                )
+
+                m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+                m.unknown_parameters.update(
+                    (k, pyo.ComponentUID(k)) for k in [m.k1, m.k2]
+                )
+
+                m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+                m.measurement_error.update((m.ca[t], 0.01) for t in meas_time_points)
+                m.measurement_error.update((m.cb[t], 0.01) for t in meas_time_points)
+                m.measurement_error.update((m.cc[t], 0.01) for t in meas_time_points)
+                m.measurement_error.update(
+                    ((m.ca[t], m.k1), 0.1) for t in meas_time_points
+                )
+
+        class ReactorIncompleteErrorException(ReactorDesignExperimentDAE):
+            def label_model(self):
+
+                m = self.model
+
+                if isinstance(self.data, pd.DataFrame):
+                    meas_time_points = self.data.index
+                else:
+                    meas_time_points = list(self.data["ca"].keys())
+
+                m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+                m.experiment_outputs.update(
+                    (m.ca[t], self.data["ca"][t]) for t in meas_time_points
+                )
+                m.experiment_outputs.update(
+                    (m.cb[t], self.data["cb"][t]) for t in meas_time_points
+                )
+                m.experiment_outputs.update(
+                    (m.cc[t], self.data["cc"][t]) for t in meas_time_points
+                )
+
+                m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+                m.unknown_parameters.update(
+                    (k, pyo.ComponentUID(k)) for k in [m.k1, m.k2]
+                )
+
+                m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL)
+                m.measurement_error.update((m.ca[t], 0.01) for t in meas_time_points)
+                m.measurement_error.update((m.cb[t], 0.01) for t in meas_time_points)
+
+        # create an experiment list with the incorrect definition of the
+        # measurement-error covariance
+        self.exp_list_df_incorrect_err_cov1 = [
+            ReactorErrorCovarianceException1(data_df)
+        ]
+        self.exp_list_df_incorrect_err_cov2 = [
+            ReactorErrorCovarianceException2(data_df)
+        ]
+        self.exp_list_df_incomplete_err = [ReactorIncompleteErrorException(data_df)]
 
     def test_unknown_parameters_exception(self):
         """
@@ -1294,6 +1490,61 @@ def test_unknown_parameters_exception(self):
 
         self.assertIn("unknown_parameters", str(context.exception))
 
+    def test_incorrect_error_covariance_exception(self):
+        """
+        Test the exception raised by parmest when the measurement-error
+        covariance is defined incorrectly
+        """
+        pest1 = parmest.Estimator(
+            self.exp_list_df_incorrect_err_cov1, obj_function="SSE"
+        )
+
+        obj1, theta1 = pest1.theta_est()
+        with pytest.raises(
+            TypeError,
+            match=r"Expected a tuple of two measured variables when specifying a "
+            r"measurement-error covariance, e\.g\., "
+            r"measurement_error\[\(y1, y2\)\] = covariance\.",
+        ):
+            pest1.cov_est()
+
+        pest2 = parmest.Estimator(
+            self.exp_list_df_incorrect_err_cov2, obj_function="SSE"
+        )
+
+        obj2, theta2 = pest2.theta_est()
+        with pytest.raises(
+            ValueError,
+            match=r"Measurement-error covariance must be defined only between "
+            r"experiment output variables\.",
+        ):
+            pest2.cov_est()
+
+        pest3 = parmest.Estimator(self.exp_list_df_incomplete_err, obj_function="SSE")
+
+        obj3, theta3 = pest3.theta_est()
+        with pytest.raises(
+            KeyError,
+            match='One or more experiment outputs are not defined in the '
+            '"measurement_error" attribute. All the variables defined '
+            'in "experiment_outputs" must be defined as keys in '
+            '"measurement_error".',
+        ):
+            pest3.cov_est()
+
+        pest4 = parmest.Estimator(
+            self.exp_list_df_incomplete_err, obj_function="SSE_weighted"
+        )
+
+        with pytest.raises(
+            KeyError,
+            match='One or more experiment outputs are not defined in the '
+            '"measurement_error" attribute. All the variables defined '
+            'in "experiment_outputs" must be defined as keys in '
+            '"measurement_error".',
+        ):
+            obj4, theta4 = pest4.theta_est()
+
     def test_dataformats(self):
         obj1, theta1 = self.pest_df.theta_est()
         obj2, theta2 = self.pest_dict.theta_est()