diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py
index ea9dfc00640..0350cebd1ed 100644
--- a/pyomo/contrib/parmest/parmest.py
+++ b/pyomo/contrib/parmest/parmest.py
@@ -234,6 +234,37 @@ def SSE(model):
     expr = sum((y - y_hat) ** 2 for y, y_hat in model.experiment_outputs.items())
     return expr
 
+def regularize_term(model, prior_FIM, theta_ref):
+    """
+    Regularization term for the objective function, which is used to penalize deviation from a 
+    reference theta
+    (theta - theta_ref).transpose() * prior_FIM * (theta - theta_ref)
+
+    theta_ref: Reference parameter value, element of matrix
+    prior_FIM: Fisher Information Matrix from prior experimental design, matrix
+    theta: Parameter value, matrix
+
+    Added to SSE objective function
+    """
+    # Check if prior_FIM is a square matrix
+    if prior_FIM.shape[0] != prior_FIM.shape[1]:
+        raise ValueError("prior_FIM must be a square matrix")
+    
+    # Check if theta_ref is a vector of the same size as prior_FIM
+    if len(theta_ref) != prior_FIM.shape[0]:
+        raise ValueError("theta_ref must be a vector of the same size as prior_FIM")
+    
+    # (theta - theta_ref).transpose() * prior_FIM * (theta - theta_ref)
+    expr = np.zeros(len(theta_ref))
+
+    for i in range(len(theta_ref)):
+        if theta_ref[i] is None:
+            raise ValueError("theta_ref must not contain None values")
+        expr[i] = (model.unknown_parameters[i] - theta_ref[i]).transpose() * prior_FIM[i] * (model.unknown_parameters[i] - theta_ref[i])
+    return sum(expr)**2
+    
+    return expr
+
 
 class Estimator(object):
     """
@@ -270,6 +301,8 @@ def __init__(
         self,
         experiment_list,
         obj_function=None,
+        prior_FIM=None,
+        theta_ref=None,
         tee=False,
         diagnostic_mode=False,
         solver_options=None,
@@ -296,6 +329,8 @@ def __init__(
 
         # populate keyword argument options
         self.obj_function = obj_function
+        self.prior_FIM = prior_FIM
+        self.theta_ref = theta_ref
         self.tee = tee
         self.diagnostic_mode = diagnostic_mode
         self.solver_options = solver_options
@@ -427,7 +462,18 @@ def _create_parmest_model(self, experiment_number):
             # TODO, this needs to be turned into an enum class of options that still support
             # custom functions
             if self.obj_function == 'SSE':
-                second_stage_rule = SSE
+                
+                if self.prior_FIM and self.theta_ref is not None:
+                    # Regularize the objective function
+                    second_stage_rule = SSE + regularize_term(model = self.model_initialized, prior_FIM = self.prior_FIM, theta_ref = self.theta_ref)
+                elif self.prior_FIM:
+                    theta_ref = model.unknown_parameters.values()
+                    second_stage_rule = SSE + regularize_term(prior_FIM = self.prior_FIM, theta_ref = theta_ref)
+
+                else:
+                    # Sum of squared errors
+                    second_stage_rule = SSE
+            
             else:
                 # A custom function uses model.experiment_outputs as data
                 second_stage_rule = self.obj_function