From b3486ca7ccf142b561a3d4fc425ea6faa712bd95 Mon Sep 17 00:00:00 2001 From: Stephen Cini <114932899+sscini@users.noreply.github.com> Date: Thu, 20 Mar 2025 16:28:58 -0400 Subject: [PATCH 01/11] Add regularize term function, and modified SSE to add option --- pyomo/contrib/parmest/parmest.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py index ea9dfc00640..bfa6703a3cf 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -234,6 +234,21 @@ def SSE(model): expr = sum((y - y_hat) ** 2 for y, y_hat in model.experiment_outputs.items()) return expr +def regularize_term(model, FIM, theta_ref): + """ + Regularization term for the objective function, which is used to penalize deviation from a + reference theta + (theta - theta_ref).transpose() * FIM * (theta - theta_ref) + + theta_ref: Reference parameter value, element of matrix + FIM: Fisher Information Matrix, matrix + theta: Parameter value, matrix + + Added to SSE objective function + """ + expr = ((theta - theta_ref).transpose() * FIM * (model - theta_ref) for theta in model.unknown_parameters.items()) + return expr + class Estimator(object): """ @@ -270,6 +285,8 @@ def __init__( self, experiment_list, obj_function=None, + FIM=None, + theta_ref=None, tee=False, diagnostic_mode=False, solver_options=None, @@ -428,6 +445,9 @@ def _create_parmest_model(self, experiment_number): # custom functions if self.obj_function == 'SSE': second_stage_rule = SSE + if self.FIM and self.theta_ref is not None: + + else: # A custom function uses model.experiment_outputs as data second_stage_rule = self.obj_function From 76ffb2272947c0f94388f07290220757e1111f57 Mon Sep 17 00:00:00 2001 From: Stephen Cini <114932899+sscini@users.noreply.github.com> Date: Thu, 20 Mar 2025 16:33:44 -0400 Subject: [PATCH 02/11] Fixed term in function --- pyomo/contrib/parmest/parmest.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py index bfa6703a3cf..28b69406900 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -246,7 +246,7 @@ def regularize_term(model, FIM, theta_ref): Added to SSE objective function """ - expr = ((theta - theta_ref).transpose() * FIM * (model - theta_ref) for theta in model.unknown_parameters.items()) + expr = ((theta - theta_ref).transpose() * FIM * (theta - theta_ref) for theta in model.unknown_parameters.items()) return expr @@ -446,6 +446,8 @@ def _create_parmest_model(self, experiment_number): if self.obj_function == 'SSE': second_stage_rule = SSE if self.FIM and self.theta_ref is not None: + # Regularize the objective function + second_stage_rule = SSE + regularize_term else: From 055f14cbdde28cc4c45fa3ec42a9640053f3ddd1 Mon Sep 17 00:00:00 2001 From: Stephen Cini <114932899+sscini@users.noreply.github.com> Date: Wed, 2 Apr 2025 00:16:11 -0400 Subject: [PATCH 03/11] Update parmest.py --- pyomo/contrib/parmest/parmest.py | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py index 28b69406900..626980dd7df 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -234,11 +234,11 @@ def SSE(model): expr = sum((y - y_hat) ** 2 for y, y_hat in model.experiment_outputs.items()) return expr -def regularize_term(model, FIM, theta_ref): +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() * FIM * (theta - theta_ref) + (theta - theta_ref).transpose() * prior_FIM * (theta - theta_ref) theta_ref: Reference parameter value, element of matrix FIM: Fisher Information Matrix, matrix @@ -246,7 +246,7 @@ def regularize_term(model, FIM, theta_ref): Added to SSE objective function """ - expr = ((theta - theta_ref).transpose() * FIM * (theta - theta_ref) for theta in model.unknown_parameters.items()) + expr = ((theta - theta_ref).transpose() * prior_FIM * (theta - theta_ref) for theta in model.unknown_parameters.items()) return expr @@ -285,7 +285,7 @@ def __init__( self, experiment_list, obj_function=None, - FIM=None, + prior_FIM=None, theta_ref=None, tee=False, diagnostic_mode=False, @@ -444,10 +444,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.FIM and self.theta_ref is not None: # Regularize the objective function - second_stage_rule = SSE + regularize_term + second_stage_rule = SSE + regularize_term(prior_FIM = self.prior_FIM, theta_ref = self.theta_ref) + elif self.FIM: + theta_ref = model.unknown_parameters.values() + second_stage_rule = SSE + regularize_term(prior_FIM = self.prior_FIM, theta_ref = self.theta_ref) + + else: + # Sum of squared errors + second_stage_rule = SSE + else: From 626c1b9ba76a2f46e1304226130ebd7a50aad51b Mon Sep 17 00:00:00 2001 From: Stephen Cini <114932899+sscini@users.noreply.github.com> Date: Wed, 2 Apr 2025 00:44:44 -0400 Subject: [PATCH 04/11] Revert "Update parmest.py" This reverts commit 055f14cbdde28cc4c45fa3ec42a9640053f3ddd1. --- pyomo/contrib/parmest/parmest.py | 20 ++++++-------------- 1 file changed, 6 insertions(+), 14 deletions(-) diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py index 626980dd7df..28b69406900 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -234,11 +234,11 @@ 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): +def regularize_term(model, 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 - theta_ref).transpose() * FIM * (theta - theta_ref) theta_ref: Reference parameter value, element of matrix FIM: Fisher Information Matrix, matrix @@ -246,7 +246,7 @@ def regularize_term(model, prior_FIM, theta_ref): Added to SSE objective function """ - expr = ((theta - theta_ref).transpose() * prior_FIM * (theta - theta_ref) for theta in model.unknown_parameters.items()) + expr = ((theta - theta_ref).transpose() * FIM * (theta - theta_ref) for theta in model.unknown_parameters.items()) return expr @@ -285,7 +285,7 @@ def __init__( self, experiment_list, obj_function=None, - prior_FIM=None, + FIM=None, theta_ref=None, tee=False, diagnostic_mode=False, @@ -444,18 +444,10 @@ 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.FIM and self.theta_ref is not None: # Regularize the objective function - second_stage_rule = SSE + regularize_term(prior_FIM = self.prior_FIM, theta_ref = self.theta_ref) - elif self.FIM: - theta_ref = model.unknown_parameters.values() - second_stage_rule = SSE + regularize_term(prior_FIM = self.prior_FIM, theta_ref = self.theta_ref) - - else: - # Sum of squared errors - second_stage_rule = SSE - + second_stage_rule = SSE + regularize_term else: From e853cc2eaa4e74cda9cd5e5379166a3e7d2619f5 Mon Sep 17 00:00:00 2001 From: Stephen Cini <114932899+sscini@users.noreply.github.com> Date: Wed, 2 Apr 2025 00:45:11 -0400 Subject: [PATCH 05/11] Reapply "Update parmest.py" This reverts commit 626c1b9ba76a2f46e1304226130ebd7a50aad51b. --- pyomo/contrib/parmest/parmest.py | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py index 28b69406900..626980dd7df 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -234,11 +234,11 @@ def SSE(model): expr = sum((y - y_hat) ** 2 for y, y_hat in model.experiment_outputs.items()) return expr -def regularize_term(model, FIM, theta_ref): +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() * FIM * (theta - theta_ref) + (theta - theta_ref).transpose() * prior_FIM * (theta - theta_ref) theta_ref: Reference parameter value, element of matrix FIM: Fisher Information Matrix, matrix @@ -246,7 +246,7 @@ def regularize_term(model, FIM, theta_ref): Added to SSE objective function """ - expr = ((theta - theta_ref).transpose() * FIM * (theta - theta_ref) for theta in model.unknown_parameters.items()) + expr = ((theta - theta_ref).transpose() * prior_FIM * (theta - theta_ref) for theta in model.unknown_parameters.items()) return expr @@ -285,7 +285,7 @@ def __init__( self, experiment_list, obj_function=None, - FIM=None, + prior_FIM=None, theta_ref=None, tee=False, diagnostic_mode=False, @@ -444,10 +444,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.FIM and self.theta_ref is not None: # Regularize the objective function - second_stage_rule = SSE + regularize_term + second_stage_rule = SSE + regularize_term(prior_FIM = self.prior_FIM, theta_ref = self.theta_ref) + elif self.FIM: + theta_ref = model.unknown_parameters.values() + second_stage_rule = SSE + regularize_term(prior_FIM = self.prior_FIM, theta_ref = self.theta_ref) + + else: + # Sum of squared errors + second_stage_rule = SSE + else: From 77c28328f9608f75a005f5f3edd8293757422f97 Mon Sep 17 00:00:00 2001 From: Stephen Cini <114932899+sscini@users.noreply.github.com> Date: Wed, 2 Apr 2025 00:51:16 -0400 Subject: [PATCH 06/11] Add theta_ref and prior_FIM to keyword args --- pyomo/contrib/parmest/parmest.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py index 626980dd7df..2ba8e066f2b 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -313,6 +313,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 From 7d0c9562c9984ba4e074ce5fa2a4b73168e97b44 Mon Sep 17 00:00:00 2001 From: Stephen Cini <114932899+sscini@users.noreply.github.com> Date: Wed, 2 Apr 2025 01:14:52 -0400 Subject: [PATCH 07/11] Adjust naming convention to ensure consistency. --- pyomo/contrib/parmest/parmest.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py index 2ba8e066f2b..5fd2429623c 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -241,7 +241,7 @@ def regularize_term(model, prior_FIM, theta_ref): (theta - theta_ref).transpose() * prior_FIM * (theta - theta_ref) theta_ref: Reference parameter value, element of matrix - FIM: Fisher Information Matrix, matrix + prior_FIM: Fisher Information Matrix from prior experimental design, matrix theta: Parameter value, matrix Added to SSE objective function @@ -447,18 +447,16 @@ def _create_parmest_model(self, experiment_number): # custom functions if self.obj_function == 'SSE': - if self.FIM and self.theta_ref is not None: + if self.prior_FIM and self.theta_ref is not None: # Regularize the objective function second_stage_rule = SSE + regularize_term(prior_FIM = self.prior_FIM, theta_ref = self.theta_ref) - elif self.FIM: + elif self.prior_FIM: theta_ref = model.unknown_parameters.values() second_stage_rule = SSE + regularize_term(prior_FIM = self.prior_FIM, theta_ref = self.theta_ref) else: # Sum of squared errors second_stage_rule = SSE - - else: # A custom function uses model.experiment_outputs as data From fae8500a69779ebb30dffd4998a9f998994a0322 Mon Sep 17 00:00:00 2001 From: Stephen Cini <114932899+sscini@users.noreply.github.com> Date: Wed, 2 Apr 2025 11:01:15 -0400 Subject: [PATCH 08/11] Update parmest.py, DoE meeting work --- pyomo/contrib/parmest/parmest.py | 22 +++++++++++++++++++--- 1 file changed, 19 insertions(+), 3 deletions(-) diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py index 5fd2429623c..0350cebd1ed 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -246,7 +246,23 @@ def regularize_term(model, prior_FIM, theta_ref): Added to SSE objective function """ - expr = ((theta - theta_ref).transpose() * prior_FIM * (theta - theta_ref) for theta in model.unknown_parameters.items()) + # 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 @@ -449,10 +465,10 @@ def _create_parmest_model(self, experiment_number): if self.prior_FIM and self.theta_ref is not None: # Regularize the objective function - second_stage_rule = SSE + regularize_term(prior_FIM = self.prior_FIM, theta_ref = self.theta_ref) + 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 = self.theta_ref) + second_stage_rule = SSE + regularize_term(prior_FIM = self.prior_FIM, theta_ref = theta_ref) else: # Sum of squared errors From 2d2ec2d1c36f125bb3bc53b638cb256aca00c316 Mon Sep 17 00:00:00 2001 From: Stephen Cini <114932899+sscini@users.noreply.github.com> Date: Tue, 15 Jul 2025 12:43:07 -0400 Subject: [PATCH 09/11] Updated for next round of review --- pyomo/contrib/parmest/parmest.py | 78 +++++++++++++++++++++++--------- 1 file changed, 57 insertions(+), 21 deletions(-) diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py index de83e252b63..5e5ed638de9 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -234,36 +234,61 @@ 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): + +# TODO: Waiting to merge with PR #3535 to follow Enum structure +def SSE_with_L2_regularization( + model, prior_FIM, theta_ref=None, regularization_weight=None +): """ - Regularization term for the objective function, which is used to penalize deviation from a - reference theta + SSE with an added L2 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 + Parameters + ---------- + model: Pyomo model + theta_ref: Reference parameter values, matrix, optional prior_FIM: Fisher Information Matrix from prior experimental design, matrix - theta: Parameter value, matrix + regularization_weight: Multiplier for regularization term, float, optional - Added to SSE objective function """ + + # Calculate sum of squared errors + SSE_expr = SSE(model) + + # Construct L2 regularization term # 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 + expr[i] = ( + (model.unknown_parameters[i] - theta_ref[i]).transpose() + * prior_FIM[i] + * (model.unknown_parameters[i] - theta_ref[i]) + ) + + # Combine SSE and regularization terms + expr_reg_L2 = sum(expr) ** 2 + + # If no regularization weight is not provided, + # scale the regularization term to be equivalent to the SSE term + if regularization_weight is None: + regularization_weight = SSE_expr / expr_reg_L2 + + expr_reg_L2 *= regularization_weight + expr_SSE_reg_L2 = SSE_expr + expr_reg_L2 + + return expr_SSE_reg_L2 class Estimator(object): @@ -458,21 +483,32 @@ def _create_parmest_model(self, experiment_number): for obj in model.component_objects(pyo.Objective): obj.deactivate() + # Completed in PR #3535, this is a temporary solution # TODO, this needs to be turned into an enum class of options that still support # custom functions if self.obj_function == 'SSE': - + # Sum of squared errors + second_stage_rule = SSE + + # Added L2 regularization option + elif self.obj_function == 'SSE_with_L2_regularization': + + # Prior FIM is required for L2 regularization + # If prior_FIM and theta_ref are provided, use them 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) + second_stage_rule = SSE_with_L2_regularization( + model=self.model_initialized, + prior_FIM=self.prior_FIM, + theta_ref=self.theta_ref, + ) + # If prior_FIM is provided but theta_ref is not, use + # unknown_parameters values as reference 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 - + theta_ref_none_provided = model.unknown_parameters.values() + second_stage_rule = SSE_with_L2_regularization( + prior_FIM=self.prior_FIM, theta_ref=theta_ref_none_provided + ) else: # A custom function uses model.experiment_outputs as data second_stage_rule = self.obj_function From d1ab692a43075d8eb165005e7d0406f9febb7241 Mon Sep 17 00:00:00 2001 From: Stephen Cini <114932899+sscini@users.noreply.github.com> Date: Tue, 15 Jul 2025 12:48:58 -0400 Subject: [PATCH 10/11] Added regularization weight to exp class initialize --- pyomo/contrib/parmest/parmest.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py index 5e5ed638de9..741c49e2d37 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -328,6 +328,7 @@ def __init__( obj_function=None, prior_FIM=None, theta_ref=None, + regularization_weight=None, tee=False, diagnostic_mode=False, solver_options=None, @@ -354,12 +355,15 @@ 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 + # Added keyword arguments for L2 regularization + self.prior_FIM = prior_FIM + self.theta_ref = theta_ref + self.regularization_weight = regularization_weight + # TODO: delete this when the deprecated interface is removed self.pest_deprecated = None From bd212f458e5c9c8f5f4ce42631080bf39167a8f1 Mon Sep 17 00:00:00 2001 From: Stephen Cini <114932899+sscini@users.noreply.github.com> Date: Tue, 15 Jul 2025 13:31:45 -0400 Subject: [PATCH 11/11] Still debugging, progress as of 7/15/2025 --- pyomo/contrib/parmest/parmest.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py index 741c49e2d37..b579864ec10 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -499,7 +499,7 @@ def _create_parmest_model(self, experiment_number): # Prior FIM is required for L2 regularization # If prior_FIM and theta_ref are provided, use them - if self.prior_FIM and self.theta_ref is not None: + if self.theta_ref.any() is not None: # Regularize the objective function second_stage_rule = SSE_with_L2_regularization( model=self.model_initialized,