From 10f6b37696180f624eb7a0bcfde7aa0386f949c3 Mon Sep 17 00:00:00 2001 From: Stephen Cini <114932899+sscini@users.noreply.github.com> Date: Mon, 17 Nov 2025 21:15:05 -0500 Subject: [PATCH 01/15] First pass at code redesign, still need to figure out more --- pyomo/contrib/parmest/parmest.py | 154 +++++++++++++++++++++++++++++++ 1 file changed, 154 insertions(+) diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py index a9dee248a85..82e8ca9698c 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -971,6 +971,91 @@ def _instance_creation_callback(self, experiment_number=None, cb_data=None): model = self._create_parmest_model(experiment_number) return model + + def _create_scenario_blocks(self): + # Create scenario block structure + # Utility function for _Q_opt_simple + # Make a block of model scenarios, one for each experiment in exp_list + + # Create a parent model to hold scenario blocks + model = pyo.ConcreteModel() + model.Blocks = pyo.Block(range(len(self.exp_list))) + for i in range(len(self.exp_list)): + # Create parmest model for experiment i + parmest_model = self._create_parmest_model(i) + # Assign parmest model to block + model.Blocks[i].model = parmest_model + + # Define objective for the block + def block_obj_rule(b): + return b.model.Total_Cost_Objective + + model.Blocks[i].obj = pyo.Objective(rule=block_obj_rule, sense=pyo.minimize) + + # Make an objective that sums over all scenario blocks + def total_obj(m): + return sum(block.obj for block in m.Blocks.values()) + + model.Obj = pyo.Objective(rule=total_obj, sense=pyo.minimize) + + # Make sure all the parameters are linked across blocks + # for name in self.estimator_theta_names: + # first_block_param = getattr(model.Blocks[0].model, name) + # for i in range(1, len(self.exp_list)): + # block_param = getattr(model.Blocks[i].model, name) + # model.Blocks[i].model.add_constraint( + # pyo.Constraint(expr=block_param == first_block_param) + # ) + + return model + + + + # Redesigning simpler version of _Q_opt + def _Q_opt_simple( + self, + return_values=None, + bootlist=None, + ThetaVals=None, + solver="ipopt", + calc_cov=NOTSET, + cov_n=NOTSET, + ): + ''' + Making new version of _Q_opt that uses scenario blocks, similar to DoE. + + Steps: + 1. Load model - parmest model should be labeled + 2. Create scenario blocks (biggest redesign) - clone model to have one per experiment + 3. Define objective and constraints for the block + 4. Solve the block as a single problem + 5. Analyze results and extract parameter estimates + + ''' + + # Create scenario blocks using utility function + model = self._create_scenario_blocks() + + solver_instance = pyo.SolverFactory(solver) + for k, v in self.solver_options.items(): + solver_instance.options[k] = v + + solver_instance.solve(model, tee=self.tee) + + assert_optimal_termination(solver_instance) + + # Extract objective value + obj_value = pyo.value(model.Obj) + theta_estimates = {} + # Extract theta estimates from first block + first_block = model.Blocks[0].model + for name in self.estimator_theta_names: + theta_var = getattr(first_block, name) + theta_estimates[name] = pyo.value(theta_var) + + return obj_value, theta_estimates + + def _Q_opt( self, ThetaVals=None, @@ -1683,6 +1768,75 @@ def theta_est( cov_n=cov_n, ) + def theta_est_simple( + self, solver="ipopt", return_values=[], calc_cov=NOTSET, cov_n=NOTSET + ): + """ + Parameter estimation using all scenarios in the data + + Parameters + ---------- + solver: str, optional + Currently only "ef_ipopt" is supported. Default is "ef_ipopt". + return_values: list, optional + List of Variable names, used to return values from the model + for data reconciliation + calc_cov: boolean, optional + DEPRECATED. + + If True, calculate and return the covariance matrix + (only for "ef_ipopt" solver). Default is NOTSET + cov_n: int, optional + DEPRECATED. + + If calc_cov=True, then the user needs to supply the number of datapoints + that are used in the objective function. Default is NOTSET + + Returns + ------- + obj_val: float + The objective function value + theta_vals: pd.Series + Estimated values for theta + var_values: pd.DataFrame + Variable values for each variable name in + return_values (only for solver='ipopt') + """ + assert isinstance(solver, str) + assert isinstance(return_values, list) + assert (calc_cov is NOTSET) or isinstance(calc_cov, bool) + + if calc_cov is not NOTSET: + deprecation_warning( + "theta_est(): `calc_cov` and `cov_n` are deprecated options and " + "will be removed in the future. Please use the `cov_est()` function " + "for covariance calculation.", + version="6.9.5", + ) + else: + calc_cov = False + + # check if we are using deprecated parmest + if self.pest_deprecated is not None and calc_cov: + return self.pest_deprecated.theta_est( + solver=solver, + return_values=return_values, + calc_cov=calc_cov, + cov_n=cov_n, + ) + elif self.pest_deprecated is not None and not calc_cov: + return self.pest_deprecated.theta_est( + solver=solver, return_values=return_values + ) + + return self._Q_opt_simple( + solver=solver, + return_values=return_values, + bootlist=None, + calc_cov=calc_cov, + cov_n=cov_n, + ) + def cov_est(self, method="finite_difference", solver="ipopt", step=1e-3): """ Covariance matrix calculation using all scenarios in the data From 3e95e91718b7b853d5b965b16ea2f2d38e511d18 Mon Sep 17 00:00:00 2001 From: Stephen Cini <114932899+sscini@users.noreply.github.com> Date: Thu, 20 Nov 2025 11:52:20 -0500 Subject: [PATCH 02/15] Added comments where I have question --- pyomo/contrib/parmest/parmest.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py index 82e8ca9698c..7b1285458fb 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -974,6 +974,7 @@ def _instance_creation_callback(self, experiment_number=None, cb_data=None): def _create_scenario_blocks(self): # Create scenario block structure + # Code is still heavily hypothetical and needs to be thought over and debugged. # Utility function for _Q_opt_simple # Make a block of model scenarios, one for each experiment in exp_list @@ -1012,6 +1013,7 @@ def total_obj(m): # Redesigning simpler version of _Q_opt + # Still work in progress def _Q_opt_simple( self, return_values=None, @@ -1768,6 +1770,8 @@ def theta_est( cov_n=cov_n, ) + # Replicate of theta_est for testing simplified _Q_opt + # Still work in progress def theta_est_simple( self, solver="ipopt", return_values=[], calc_cov=NOTSET, cov_n=NOTSET ): From 3982e1b4019e4ab6a39d2921d39c580732e33880 Mon Sep 17 00:00:00 2001 From: Stephen Cini <114932899+sscini@users.noreply.github.com> Date: Fri, 21 Nov 2025 01:35:19 -0500 Subject: [PATCH 03/15] Got preliminary _Q_opt simple working with example! --- pyomo/contrib/parmest/parmest.py | 51 ++++++++++++++++---------------- 1 file changed, 26 insertions(+), 25 deletions(-) diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py index 7b1285458fb..7d576792a75 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -980,33 +980,36 @@ def _create_scenario_blocks(self): # Create a parent model to hold scenario blocks model = pyo.ConcreteModel() - model.Blocks = pyo.Block(range(len(self.exp_list))) + model.exp_scenarios = pyo.Block(range(len(self.exp_list))) for i in range(len(self.exp_list)): # Create parmest model for experiment i parmest_model = self._create_parmest_model(i) # Assign parmest model to block - model.Blocks[i].model = parmest_model - - # Define objective for the block - def block_obj_rule(b): - return b.model.Total_Cost_Objective - - model.Blocks[i].obj = pyo.Objective(rule=block_obj_rule, sense=pyo.minimize) + model.exp_scenarios[i].transfer_attributes_from(parmest_model) # Make an objective that sums over all scenario blocks def total_obj(m): - return sum(block.obj for block in m.Blocks.values()) + return sum(block.Total_Cost_Objective for block in m.exp_scenarios.values())/len(self.exp_list) model.Obj = pyo.Objective(rule=total_obj, sense=pyo.minimize) # Make sure all the parameters are linked across blocks - # for name in self.estimator_theta_names: - # first_block_param = getattr(model.Blocks[0].model, name) - # for i in range(1, len(self.exp_list)): - # block_param = getattr(model.Blocks[i].model, name) - # model.Blocks[i].model.add_constraint( - # pyo.Constraint(expr=block_param == first_block_param) - # ) + for name in self.estimator_theta_names: + # Get the variable from the first block + ref_var = getattr(model.exp_scenarios[0], name) + for i in range(1, len(self.exp_list)): + curr_var = getattr(model.exp_scenarios[i], name) + # Constrain current variable to equal reference variable + model.add_component( + f"Link_{name}_Block0_Block{i}", + pyo.Constraint(expr=curr_var == ref_var) + ) + + # Deactivate the objective in each block to avoid double counting + for i in range(len(self.exp_list)): + model.exp_scenarios[i].Total_Cost_Objective.deactivate() + + model.pprint() return model @@ -1038,22 +1041,20 @@ def _Q_opt_simple( # Create scenario blocks using utility function model = self._create_scenario_blocks() - solver_instance = pyo.SolverFactory(solver) - for k, v in self.solver_options.items(): - solver_instance.options[k] = v - - solver_instance.solve(model, tee=self.tee) + solver = SolverFactory('ipopt') + if self.solver_options is not None: + for key in self.solver_options: + solver.options[key] = self.solver_options[key] - assert_optimal_termination(solver_instance) + solve_result = solver.solve(model, tee=self.tee) + assert_optimal_termination(solve_result) # Extract objective value obj_value = pyo.value(model.Obj) theta_estimates = {} # Extract theta estimates from first block - first_block = model.Blocks[0].model for name in self.estimator_theta_names: - theta_var = getattr(first_block, name) - theta_estimates[name] = pyo.value(theta_var) + theta_estimates[name] = pyo.value(getattr(model.exp_scenarios[0], name)) return obj_value, theta_estimates From e829344e29df84194e749ac3a536f121bab84d9e Mon Sep 17 00:00:00 2001 From: Stephen Cini <114932899+sscini@users.noreply.github.com> Date: Fri, 21 Nov 2025 01:36:35 -0500 Subject: [PATCH 04/15] Ran black --- pyomo/contrib/parmest/parmest.py | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py index 7d576792a75..e3be94e5092 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -971,12 +971,11 @@ def _instance_creation_callback(self, experiment_number=None, cb_data=None): model = self._create_parmest_model(experiment_number) return model - def _create_scenario_blocks(self): # Create scenario block structure # Code is still heavily hypothetical and needs to be thought over and debugged. # Utility function for _Q_opt_simple - # Make a block of model scenarios, one for each experiment in exp_list + # Make a block of model scenarios, one for each experiment in exp_list # Create a parent model to hold scenario blocks model = pyo.ConcreteModel() @@ -989,8 +988,10 @@ def _create_scenario_blocks(self): # Make an objective that sums over all scenario blocks def total_obj(m): - return sum(block.Total_Cost_Objective for block in m.exp_scenarios.values())/len(self.exp_list) - + return sum( + block.Total_Cost_Objective for block in m.exp_scenarios.values() + ) / len(self.exp_list) + model.Obj = pyo.Objective(rule=total_obj, sense=pyo.minimize) # Make sure all the parameters are linked across blocks @@ -1002,7 +1003,7 @@ def total_obj(m): # Constrain current variable to equal reference variable model.add_component( f"Link_{name}_Block0_Block{i}", - pyo.Constraint(expr=curr_var == ref_var) + pyo.Constraint(expr=curr_var == ref_var), ) # Deactivate the objective in each block to avoid double counting @@ -1013,8 +1014,6 @@ def total_obj(m): return model - - # Redesigning simpler version of _Q_opt # Still work in progress def _Q_opt_simple( @@ -1025,7 +1024,7 @@ def _Q_opt_simple( solver="ipopt", calc_cov=NOTSET, cov_n=NOTSET, - ): + ): ''' Making new version of _Q_opt that uses scenario blocks, similar to DoE. @@ -1037,8 +1036,8 @@ def _Q_opt_simple( 5. Analyze results and extract parameter estimates ''' - - # Create scenario blocks using utility function + + # Create scenario blocks using utility function model = self._create_scenario_blocks() solver = SolverFactory('ipopt') @@ -1058,7 +1057,6 @@ def _Q_opt_simple( return obj_value, theta_estimates - def _Q_opt( self, ThetaVals=None, @@ -1841,7 +1839,7 @@ def theta_est_simple( calc_cov=calc_cov, cov_n=cov_n, ) - + def cov_est(self, method="finite_difference", solver="ipopt", step=1e-3): """ Covariance matrix calculation using all scenarios in the data From e46409797c9461c11edc61f79649bccc507bc670 Mon Sep 17 00:00:00 2001 From: Stephen Cini <114932899+sscini@users.noreply.github.com> Date: Fri, 21 Nov 2025 09:36:03 -0500 Subject: [PATCH 05/15] Changed name to _Q_opt_blocks --- pyomo/contrib/parmest/parmest.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py index e3be94e5092..c0fc5f1213f 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -986,7 +986,7 @@ def _create_scenario_blocks(self): # Assign parmest model to block model.exp_scenarios[i].transfer_attributes_from(parmest_model) - # Make an objective that sums over all scenario blocks + # Make an objective that sums over all scenario blocks and divides by number of experiments def total_obj(m): return sum( block.Total_Cost_Objective for block in m.exp_scenarios.values() @@ -1016,7 +1016,7 @@ def total_obj(m): # Redesigning simpler version of _Q_opt # Still work in progress - def _Q_opt_simple( + def _Q_opt_blocks( self, return_values=None, bootlist=None, @@ -1771,7 +1771,7 @@ def theta_est( # Replicate of theta_est for testing simplified _Q_opt # Still work in progress - def theta_est_simple( + def theta_est_blocks( self, solver="ipopt", return_values=[], calc_cov=NOTSET, cov_n=NOTSET ): """ From dc5ee767eac4aba5e41a81e1aecbea80ca702918 Mon Sep 17 00:00:00 2001 From: Stephen Cini <114932899+sscini@users.noreply.github.com> Date: Fri, 21 Nov 2025 10:48:13 -0500 Subject: [PATCH 06/15] Update parmest.py --- pyomo/contrib/parmest/parmest.py | 62 +++++++++++++++++++++++--------- 1 file changed, 46 insertions(+), 16 deletions(-) diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py index c0fc5f1213f..877dccaebe5 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -971,21 +971,43 @@ def _instance_creation_callback(self, experiment_number=None, cb_data=None): model = self._create_parmest_model(experiment_number) return model - def _create_scenario_blocks(self): + def _create_scenario_blocks(self, bootlist=None,): # Create scenario block structure - # Code is still heavily hypothetical and needs to be thought over and debugged. - # Utility function for _Q_opt_simple + # Utility function for _Q_opt_blocks # Make a block of model scenarios, one for each experiment in exp_list # Create a parent model to hold scenario blocks model = pyo.ConcreteModel() - model.exp_scenarios = pyo.Block(range(len(self.exp_list))) + + if bootlist is not None: + model.exp_scenarios = pyo.Block(range(len(bootlist))) + else: + model.exp_scenarios = pyo.Block(range(len(self.exp_list))) + for i in range(len(self.exp_list)): # Create parmest model for experiment i parmest_model = self._create_parmest_model(i) # Assign parmest model to block model.exp_scenarios[i].transfer_attributes_from(parmest_model) + # Transfer all the unknown parameters to the parent model + for name in self.estimator_theta_names: + # Get the variable from the first block + ref_var = getattr(model.exp_scenarios[0], name) + # Create a variable in the parent model with same bounds and initialization + parent_var = pyo.Var( + bounds=ref_var.bounds, + initialize=pyo.value(ref_var), + ) + setattr(model, name, parent_var) + # Constrain the variable in the first block to equal the parent variable + model.add_component( + f"Link_{name}_Block0_Parent", + pyo.Constraint( + expr=getattr(model.exp_scenarios[0], name) == parent_var + ), + ) + # Make an objective that sums over all scenario blocks and divides by number of experiments def total_obj(m): return sum( @@ -996,14 +1018,13 @@ def total_obj(m): # Make sure all the parameters are linked across blocks for name in self.estimator_theta_names: - # Get the variable from the first block - ref_var = getattr(model.exp_scenarios[0], name) for i in range(1, len(self.exp_list)): - curr_var = getattr(model.exp_scenarios[i], name) - # Constrain current variable to equal reference variable model.add_component( - f"Link_{name}_Block0_Block{i}", - pyo.Constraint(expr=curr_var == ref_var), + f"Link_{name}_Block{i}_Parent", + pyo.Constraint( + expr=getattr(model.exp_scenarios[i], name) + == getattr(model, name) + ), ) # Deactivate the objective in each block to avoid double counting @@ -1014,8 +1035,8 @@ def total_obj(m): return model - # Redesigning simpler version of _Q_opt - # Still work in progress + # Redesigning version of _Q_opt that uses scenario blocks + # Works, but still adding features from old _Q_opt def _Q_opt_blocks( self, return_values=None, @@ -1038,14 +1059,15 @@ def _Q_opt_blocks( ''' # Create scenario blocks using utility function - model = self._create_scenario_blocks() + model = self._create_scenario_blocks(bootlist=bootlist) - solver = SolverFactory('ipopt') + if solver == "ipopt": + sol = SolverFactory('ipopt') if self.solver_options is not None: for key in self.solver_options: solver.options[key] = self.solver_options[key] - solve_result = solver.solve(model, tee=self.tee) + solve_result = sol.solve(model, tee=self.tee) assert_optimal_termination(solve_result) # Extract objective value @@ -1055,6 +1077,14 @@ def _Q_opt_blocks( for name in self.estimator_theta_names: theta_estimates[name] = pyo.value(getattr(model.exp_scenarios[0], name)) + # Check they are equal to the second block + for name in self.estimator_theta_names: + val_block1 = pyo.value(getattr(model.exp_scenarios[1], name)) + assert theta_estimates[name] == val_block1, ( + f"Parameter {name} estimate differs between blocks: " + f"{theta_estimates[name]} vs {val_block1}" + ) + return obj_value, theta_estimates def _Q_opt( @@ -1832,7 +1862,7 @@ def theta_est_blocks( solver=solver, return_values=return_values ) - return self._Q_opt_simple( + return self._Q_opt_blocks( solver=solver, return_values=return_values, bootlist=None, From 63558185c5147df6d42baea0ade6530dcfe88a12 Mon Sep 17 00:00:00 2001 From: Stephen Cini <114932899+sscini@users.noreply.github.com> Date: Fri, 21 Nov 2025 10:48:30 -0500 Subject: [PATCH 07/15] Ran black --- pyomo/contrib/parmest/parmest.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py index 877dccaebe5..f4878be6660 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -971,7 +971,7 @@ def _instance_creation_callback(self, experiment_number=None, cb_data=None): model = self._create_parmest_model(experiment_number) return model - def _create_scenario_blocks(self, bootlist=None,): + def _create_scenario_blocks(self, bootlist=None): # Create scenario block structure # Utility function for _Q_opt_blocks # Make a block of model scenarios, one for each experiment in exp_list @@ -983,7 +983,7 @@ def _create_scenario_blocks(self, bootlist=None,): model.exp_scenarios = pyo.Block(range(len(bootlist))) else: model.exp_scenarios = pyo.Block(range(len(self.exp_list))) - + for i in range(len(self.exp_list)): # Create parmest model for experiment i parmest_model = self._create_parmest_model(i) @@ -995,10 +995,7 @@ def _create_scenario_blocks(self, bootlist=None,): # Get the variable from the first block ref_var = getattr(model.exp_scenarios[0], name) # Create a variable in the parent model with same bounds and initialization - parent_var = pyo.Var( - bounds=ref_var.bounds, - initialize=pyo.value(ref_var), - ) + parent_var = pyo.Var(bounds=ref_var.bounds, initialize=pyo.value(ref_var)) setattr(model, name, parent_var) # Constrain the variable in the first block to equal the parent variable model.add_component( From 099f541626269c50a44f68c43d60fd4666ea58e7 Mon Sep 17 00:00:00 2001 From: Stephen Cini <114932899+sscini@users.noreply.github.com> Date: Tue, 25 Nov 2025 12:07:05 -0500 Subject: [PATCH 08/15] Added in case for bootlist, works with example --- pyomo/contrib/parmest/parmest.py | 193 ++++++++++++++++++++++++++++--- 1 file changed, 176 insertions(+), 17 deletions(-) diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py index f4878be6660..14c42dd6f89 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -922,6 +922,7 @@ def _create_parmest_model(self, experiment_number): model.parmest_dummy_var = pyo.Var(initialize=1.0) # Add objective function (optional) + # @Reviewers What is the purpose of the reserved_names? Can we discuss this in a meeting? if self.obj_function: # Check for component naming conflicts reserved_names = [ @@ -981,14 +982,23 @@ def _create_scenario_blocks(self, bootlist=None): if bootlist is not None: model.exp_scenarios = pyo.Block(range(len(bootlist))) + + for i in range(len(bootlist)): + # Create parmest model for experiment i + parmest_model = self._create_parmest_model(bootlist[i]) + # Assign parmest model to block + model.exp_scenarios[i].transfer_attributes_from(parmest_model) + else: model.exp_scenarios = pyo.Block(range(len(self.exp_list))) - for i in range(len(self.exp_list)): - # Create parmest model for experiment i - parmest_model = self._create_parmest_model(i) - # Assign parmest model to block - model.exp_scenarios[i].transfer_attributes_from(parmest_model) + for i in range(len(self.exp_list)): + # Create parmest model for experiment i + parmest_model = self._create_parmest_model(i) + # parmest_model.pprint() + # Assign parmest model to block + model.exp_scenarios[i].transfer_attributes_from(parmest_model) + # model.exp_scenarios[i].pprint() # Transfer all the unknown parameters to the parent model for name in self.estimator_theta_names: @@ -1015,20 +1025,33 @@ def total_obj(m): # Make sure all the parameters are linked across blocks for name in self.estimator_theta_names: - for i in range(1, len(self.exp_list)): - model.add_component( - f"Link_{name}_Block{i}_Parent", - pyo.Constraint( - expr=getattr(model.exp_scenarios[i], name) - == getattr(model, name) - ), - ) + if bootlist is not None: + for i in range(1, len(bootlist)): + model.add_component( + f"Link_{name}_Block{i}_Parent", + pyo.Constraint( + expr=getattr(model.exp_scenarios[i], name) + == getattr(model, name) + ), + ) + # Deactivate the objective in each block to avoid double counting + for i in range(len(bootlist)): + model.exp_scenarios[i].Total_Cost_Objective.deactivate() + else: + for i in range(1, len(self.exp_list)): + model.add_component( + f"Link_{name}_Block{i}_Parent", + pyo.Constraint( + expr=getattr(model.exp_scenarios[i], name) + == getattr(model, name) + ), + ) - # Deactivate the objective in each block to avoid double counting - for i in range(len(self.exp_list)): - model.exp_scenarios[i].Total_Cost_Objective.deactivate() + # Deactivate the objective in each block to avoid double counting + for i in range(len(self.exp_list)): + model.exp_scenarios[i].Total_Cost_Objective.deactivate() - model.pprint() + # model.pprint() return model @@ -1989,6 +2012,81 @@ def theta_est_bootstrap( del bootstrap_theta['samples'] return bootstrap_theta + + # Add theta_est_bootstrap_blocks + def theta_est_bootstrap_blocks( + self, + bootstrap_samples, + samplesize=None, + replacement=True, + seed=None, + return_samples=False, + ): + """ + Parameter estimation using bootstrap resampling of the data + + Parameters + ---------- + bootstrap_samples: int + Number of bootstrap samples to draw from the data + samplesize: int or None, optional + Size of each bootstrap sample. If samplesize=None, samplesize will be + set to the number of samples in the data + replacement: bool, optional + Sample with or without replacement. Default is True. + seed: int or None, optional + Random seed + return_samples: bool, optional + Return a list of sample numbers used in each bootstrap estimation. + Default is False. + + Returns + ------- + bootstrap_theta: pd.DataFrame + Theta values for each sample and (if return_samples = True) + the sample numbers used in each estimation + """ + + # check if we are using deprecated parmest + if self.pest_deprecated is not None: + return self.pest_deprecated.theta_est_bootstrap( + bootstrap_samples, + samplesize=samplesize, + replacement=replacement, + seed=seed, + return_samples=return_samples, + ) + + assert isinstance(bootstrap_samples, int) + assert isinstance(samplesize, (type(None), int)) + assert isinstance(replacement, bool) + assert isinstance(seed, (type(None), int)) + assert isinstance(return_samples, bool) + + if samplesize is None: + samplesize = len(self.exp_list) + + if seed is not None: + np.random.seed(seed) + + global_list = self._get_sample_list(samplesize, bootstrap_samples, replacement) + + task_mgr = utils.ParallelTaskManager(bootstrap_samples) + local_list = task_mgr.global_to_local_data(global_list) + + bootstrap_theta = list() + for idx, sample in local_list: + objval, thetavals = self._Q_opt_blocks(bootlist=list(sample)) + thetavals['samples'] = sample + bootstrap_theta.append(thetavals) + + global_bootstrap_theta = task_mgr.allgather_global_data(bootstrap_theta) + bootstrap_theta = pd.DataFrame(global_bootstrap_theta) + + if not return_samples: + del bootstrap_theta['samples'] + + return bootstrap_theta def theta_est_leaveNout( self, lNo, lNo_samples=None, seed=None, return_samples=False @@ -2051,6 +2149,67 @@ def theta_est_leaveNout( return lNo_theta + def theta_est_leaveNout_blocks( + self, lNo, lNo_samples=None, seed=None, return_samples=False + ): + """ + Parameter estimation where N data points are left out of each sample + + Parameters + ---------- + lNo: int + Number of data points to leave out for parameter estimation + lNo_samples: int + Number of leave-N-out samples. If lNo_samples=None, the maximum + number of combinations will be used + seed: int or None, optional + Random seed + return_samples: bool, optional + Return a list of sample numbers that were left out. Default is False. + + Returns + ------- + lNo_theta: pd.DataFrame + Theta values for each sample and (if return_samples = True) + the sample numbers left out of each estimation + """ + + # check if we are using deprecated parmest + if self.pest_deprecated is not None: + return self.pest_deprecated.theta_est_leaveNout( + lNo, lNo_samples=lNo_samples, seed=seed, return_samples=return_samples + ) + + assert isinstance(lNo, int) + assert isinstance(lNo_samples, (type(None), int)) + assert isinstance(seed, (type(None), int)) + assert isinstance(return_samples, bool) + + samplesize = len(self.exp_list) - lNo + + if seed is not None: + np.random.seed(seed) + + global_list = self._get_sample_list(samplesize, lNo_samples, replacement=False) + + task_mgr = utils.ParallelTaskManager(len(global_list)) + local_list = task_mgr.global_to_local_data(global_list) + + lNo_theta = list() + for idx, sample in local_list: + objval, thetavals = self._Q_opt_blocks(bootlist=list(sample)) + lNo_s = list(set(range(len(self.exp_list))) - set(sample)) + thetavals['lNo'] = np.sort(lNo_s) + lNo_theta.append(thetavals) + + global_bootstrap_theta = task_mgr.allgather_global_data(lNo_theta) + lNo_theta = pd.DataFrame(global_bootstrap_theta) + + if not return_samples: + del lNo_theta['lNo'] + + return lNo_theta + def leaveNout_bootstrap_test( self, lNo, lNo_samples, bootstrap_samples, distribution, alphas, seed=None ): From 76ee05ec0b9e203c078fda98a4323a3ef0a610cf Mon Sep 17 00:00:00 2001 From: Stephen Cini <114932899+sscini@users.noreply.github.com> Date: Tue, 25 Nov 2025 12:08:16 -0500 Subject: [PATCH 09/15] Ran black --- 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 14c42dd6f89..d8dcc7839c9 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -2012,7 +2012,7 @@ def theta_est_bootstrap( del bootstrap_theta['samples'] return bootstrap_theta - + # Add theta_est_bootstrap_blocks def theta_est_bootstrap_blocks( self, From d91ce3f8ba3b139edf4ad4a9906e384d624b1ba7 Mon Sep 17 00:00:00 2001 From: Stephen Cini <114932899+sscini@users.noreply.github.com> Date: Fri, 28 Nov 2025 19:10:59 -0500 Subject: [PATCH 10/15] Simplified structure, ran black --- pyomo/contrib/parmest/parmest.py | 37 +++++++++++--------------------- 1 file changed, 13 insertions(+), 24 deletions(-) diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py index d8dcc7839c9..4fe1e12b5e9 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -981,6 +981,7 @@ def _create_scenario_blocks(self, bootlist=None): model = pyo.ConcreteModel() if bootlist is not None: + n_scenarios = len(bootlist) model.exp_scenarios = pyo.Block(range(len(bootlist))) for i in range(len(bootlist)): @@ -990,6 +991,7 @@ def _create_scenario_blocks(self, bootlist=None): model.exp_scenarios[i].transfer_attributes_from(parmest_model) else: + n_scenarios = len(self.exp_list) model.exp_scenarios = pyo.Block(range(len(self.exp_list))) for i in range(len(self.exp_list)): @@ -1025,31 +1027,18 @@ def total_obj(m): # Make sure all the parameters are linked across blocks for name in self.estimator_theta_names: - if bootlist is not None: - for i in range(1, len(bootlist)): - model.add_component( - f"Link_{name}_Block{i}_Parent", - pyo.Constraint( - expr=getattr(model.exp_scenarios[i], name) - == getattr(model, name) - ), - ) - # Deactivate the objective in each block to avoid double counting - for i in range(len(bootlist)): - model.exp_scenarios[i].Total_Cost_Objective.deactivate() - else: - for i in range(1, len(self.exp_list)): - model.add_component( - f"Link_{name}_Block{i}_Parent", - pyo.Constraint( - expr=getattr(model.exp_scenarios[i], name) - == getattr(model, name) - ), - ) + for i in range(1, n_scenarios): + model.add_component( + f"Link_{name}_Block{i}_Parent", + pyo.Constraint( + expr=getattr(model.exp_scenarios[i], name) + == getattr(model, name) + ), + ) - # Deactivate the objective in each block to avoid double counting - for i in range(len(self.exp_list)): - model.exp_scenarios[i].Total_Cost_Objective.deactivate() + # Deactivate the objective in each block to avoid double counting + for i in range(n_scenarios): + model.exp_scenarios[i].Total_Cost_Objective.deactivate() # model.pprint() From 1aea99f4f2ea82a46d2b62459a67cd30693e78a0 Mon Sep 17 00:00:00 2001 From: Stephen Cini <114932899+sscini@users.noreply.github.com> Date: Sat, 13 Dec 2025 15:09:33 -0500 Subject: [PATCH 11/15] Removed _Q_opt, and replicate functions, only using _Q_opt_blocks --- pyomo/contrib/parmest/parmest.py | 461 +++---------------------------- 1 file changed, 33 insertions(+), 428 deletions(-) diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py index 4fe1e12b5e9..98523eb219c 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -786,7 +786,6 @@ def __init__( diagnostic_mode=False, solver_options=None, ): - # check that we have a (non-empty) list of experiments assert isinstance(experiment_list, list) self.exp_list = experiment_list @@ -850,7 +849,6 @@ def _deprecated_init( diagnostic_mode=False, solver_options=None, ): - deprecation_warning( "You're using the deprecated parmest interface (model_function, " "data, theta_names). This interface will be removed in a future release, " @@ -873,26 +871,22 @@ def _return_theta_names(self): """ # check for deprecated inputs if self.pest_deprecated: - # if fitted model parameter names differ from theta_names # created when Estimator object is created if hasattr(self, 'theta_names_updated'): return self.pest_deprecated.theta_names_updated else: - # default theta_names, created when Estimator object is created return self.pest_deprecated.theta_names else: - # if fitted model parameter names differ from theta_names # created when Estimator object is created if hasattr(self, 'theta_names_updated'): return self.theta_names_updated else: - # default theta_names, created when Estimator object is created return self.estimator_theta_names @@ -1046,6 +1040,7 @@ def total_obj(m): # Redesigning version of _Q_opt that uses scenario blocks # Works, but still adding features from old _Q_opt + # @Reviewers: Trying to find best way to integrate the ability to fix thetas def _Q_opt_blocks( self, return_values=None, @@ -1096,198 +1091,6 @@ def _Q_opt_blocks( return obj_value, theta_estimates - def _Q_opt( - self, - ThetaVals=None, - solver="ef_ipopt", - return_values=[], - bootlist=None, - calc_cov=NOTSET, - cov_n=NOTSET, - ): - """ - Set up all thetas as first stage Vars, return resulting theta - values as well as the objective function value. - - """ - if solver == "k_aug": - raise RuntimeError("k_aug no longer supported.") - - # (Bootstrap scenarios will use indirection through the bootlist) - if bootlist is None: - scenario_numbers = list(range(len(self.exp_list))) - scen_names = ["Scenario{}".format(i) for i in scenario_numbers] - else: - scen_names = ["Scenario{}".format(i) for i in range(len(bootlist))] - - # get the probability constant that is applied to the objective function - # parmest solves the estimation problem by applying equal probabilities to - # the objective function of all the scenarios from the experiment list - self.obj_probability_constant = len(scen_names) - - # tree_model.CallbackModule = None - outer_cb_data = dict() - outer_cb_data["callback"] = self._instance_creation_callback - if ThetaVals is not None: - outer_cb_data["ThetaVals"] = ThetaVals - if bootlist is not None: - outer_cb_data["BootList"] = bootlist - outer_cb_data["cb_data"] = None # None is OK - outer_cb_data["theta_names"] = self.estimator_theta_names - - options = {"solver": "ipopt"} - scenario_creator_options = {"cb_data": outer_cb_data} - if use_mpisppy: - ef = sputils.create_EF( - scen_names, - _experiment_instance_creation_callback, - EF_name="_Q_opt", - suppress_warnings=True, - scenario_creator_kwargs=scenario_creator_options, - ) - else: - ef = local_ef.create_EF( - scen_names, - _experiment_instance_creation_callback, - EF_name="_Q_opt", - suppress_warnings=True, - scenario_creator_kwargs=scenario_creator_options, - ) - self.ef_instance = ef - - # Solve the extensive form with ipopt - if solver == "ef_ipopt": - if calc_cov is NOTSET or not calc_cov: - # Do not calculate the reduced hessian - - solver = SolverFactory('ipopt') - if self.solver_options is not None: - for key in self.solver_options: - solver.options[key] = self.solver_options[key] - - solve_result = solver.solve(self.ef_instance, tee=self.tee) - assert_optimal_termination(solve_result) - elif calc_cov is not NOTSET and calc_cov: - # parmest makes the fitted parameters stage 1 variables - ind_vars = [] - for nd_name, Var, sol_val in ef_nonants(ef): - ind_vars.append(Var) - # calculate the reduced hessian - (solve_result, inv_red_hes) = ( - inverse_reduced_hessian.inv_reduced_hessian_barrier( - self.ef_instance, - independent_variables=ind_vars, - solver_options=self.solver_options, - tee=self.tee, - ) - ) - - if self.diagnostic_mode: - print( - ' Solver termination condition = ', - str(solve_result.solver.termination_condition), - ) - - # assume all first stage are thetas... - theta_vals = {} - for nd_name, Var, sol_val in ef_nonants(ef): - # process the name - # the scenarios are blocks, so strip the scenario name - var_name = Var.name[Var.name.find(".") + 1 :] - theta_vals[var_name] = sol_val - - obj_val = pyo.value(ef.EF_Obj) - self.obj_value = obj_val - self.estimated_theta = theta_vals - - if calc_cov is not NOTSET and calc_cov: - # Calculate the covariance matrix - - if not isinstance(cov_n, int): - raise TypeError( - f"Expected an integer for the 'cov_n' argument. " - f"Got {type(cov_n)}." - ) - num_unknowns = max( - [ - len(experiment.get_labeled_model().unknown_parameters) - for experiment in self.exp_list - ] - ) - assert cov_n > num_unknowns, ( - "The number of datapoints must be greater than the " - "number of parameters to estimate." - ) - - # Number of data points considered - n = cov_n - - # Extract number of fitted parameters - l = len(theta_vals) - - # Assumption: Objective value is sum of squared errors - sse = obj_val - - '''Calculate covariance assuming experimental observation errors - are independent and follow a Gaussian distribution - with constant variance. - - The formula used in parmest was verified against equations - (7-5-15) and (7-5-16) in "Nonlinear Parameter Estimation", - Y. Bard, 1974. - - This formula is also applicable if the objective is scaled by a - constant; the constant cancels out. - (was scaled by 1/n because it computes an expected value.) - ''' - cov = 2 * sse / (n - l) * inv_red_hes - cov = pd.DataFrame( - cov, index=theta_vals.keys(), columns=theta_vals.keys() - ) - - theta_vals = pd.Series(theta_vals) - - if len(return_values) > 0: - var_values = [] - if len(scen_names) > 1: # multiple scenarios - block_objects = self.ef_instance.component_objects( - Block, descend_into=False - ) - else: # single scenario - block_objects = [self.ef_instance] - for exp_i in block_objects: - vals = {} - for var in return_values: - exp_i_var = exp_i.find_component(str(var)) - if ( - exp_i_var is None - ): # we might have a block such as _mpisppy_data - continue - # if value to return is ContinuousSet - if type(exp_i_var) == ContinuousSet: - temp = list(exp_i_var) - else: - temp = [pyo.value(_) for _ in exp_i_var.values()] - if len(temp) == 1: - vals[var] = temp[0] - else: - vals[var] = temp - if len(vals) > 0: - var_values.append(vals) - var_values = pd.DataFrame(var_values) - if calc_cov is not NOTSET and calc_cov: - return obj_val, theta_vals, var_values, cov - elif calc_cov is NOTSET or not calc_cov: - return obj_val, theta_vals, var_values - - if calc_cov is not NOTSET and calc_cov: - return obj_val, theta_vals, cov - elif calc_cov is NOTSET or not calc_cov: - return obj_val, theta_vals - - else: - raise RuntimeError("Unknown solver in Q_Opt=" + solver) - def _cov_at_theta(self, method, solver, step): """ Covariance matrix calculation using all scenarios in the data @@ -1316,13 +1119,14 @@ def _cov_at_theta(self, method, solver, step): for nd_name, Var, sol_val in ef_nonants(self.ef_instance): ind_vars.append(Var) # calculate the reduced hessian - (solve_result, inv_red_hes) = ( - inverse_reduced_hessian.inv_reduced_hessian_barrier( - self.ef_instance, - independent_variables=ind_vars, - solver_options=self.solver_options, - tee=self.tee, - ) + ( + solve_result, + inv_red_hes, + ) = inverse_reduced_hessian.inv_reduced_hessian_barrier( + self.ef_instance, + independent_variables=ind_vars, + solver_options=self.solver_options, + tee=self.tee, ) self.inv_red_hes = inv_red_hes @@ -1611,10 +1415,14 @@ def _Q_at_theta(self, thetavals, initialize_parmest_model=False): if self.diagnostic_mode: print(' Experiment = ', snum) print(' First solve with special diagnostics wrapper') - (status_obj, solved, iters, time, regu) = ( - utils.ipopt_solve_with_stats( - instance, optimizer, max_iter=500, max_cpu_time=120 - ) + ( + status_obj, + solved, + iters, + time, + regu, + ) = utils.ipopt_solve_with_stats( + instance, optimizer, max_iter=500, max_cpu_time=120 ) print( " status_obj, solved, iters, time, regularization_stat = ", @@ -1777,77 +1585,6 @@ def theta_est( assert isinstance(return_values, list) assert (calc_cov is NOTSET) or isinstance(calc_cov, bool) - if calc_cov is not NOTSET: - deprecation_warning( - "theta_est(): `calc_cov` and `cov_n` are deprecated options and " - "will be removed in the future. Please use the `cov_est()` function " - "for covariance calculation.", - version="6.9.5", - ) - else: - calc_cov = False - - # check if we are using deprecated parmest - if self.pest_deprecated is not None and calc_cov: - return self.pest_deprecated.theta_est( - solver=solver, - return_values=return_values, - calc_cov=calc_cov, - cov_n=cov_n, - ) - elif self.pest_deprecated is not None and not calc_cov: - return self.pest_deprecated.theta_est( - solver=solver, return_values=return_values - ) - - return self._Q_opt( - solver=solver, - return_values=return_values, - bootlist=None, - calc_cov=calc_cov, - cov_n=cov_n, - ) - - # Replicate of theta_est for testing simplified _Q_opt - # Still work in progress - def theta_est_blocks( - self, solver="ipopt", return_values=[], calc_cov=NOTSET, cov_n=NOTSET - ): - """ - Parameter estimation using all scenarios in the data - - Parameters - ---------- - solver: str, optional - Currently only "ef_ipopt" is supported. Default is "ef_ipopt". - return_values: list, optional - List of Variable names, used to return values from the model - for data reconciliation - calc_cov: boolean, optional - DEPRECATED. - - If True, calculate and return the covariance matrix - (only for "ef_ipopt" solver). Default is NOTSET - cov_n: int, optional - DEPRECATED. - - If calc_cov=True, then the user needs to supply the number of datapoints - that are used in the objective function. Default is NOTSET - - Returns - ------- - obj_val: float - The objective function value - theta_vals: pd.Series - Estimated values for theta - var_values: pd.DataFrame - Variable values for each variable name in - return_values (only for solver='ipopt') - """ - assert isinstance(solver, str) - assert isinstance(return_values, list) - assert (calc_cov is NOTSET) or isinstance(calc_cov, bool) - if calc_cov is not NOTSET: deprecation_warning( "theta_est(): `calc_cov` and `cov_n` are deprecated options and " @@ -1988,81 +1725,6 @@ def theta_est_bootstrap( task_mgr = utils.ParallelTaskManager(bootstrap_samples) local_list = task_mgr.global_to_local_data(global_list) - bootstrap_theta = list() - for idx, sample in local_list: - objval, thetavals = self._Q_opt(bootlist=list(sample)) - thetavals['samples'] = sample - bootstrap_theta.append(thetavals) - - global_bootstrap_theta = task_mgr.allgather_global_data(bootstrap_theta) - bootstrap_theta = pd.DataFrame(global_bootstrap_theta) - - if not return_samples: - del bootstrap_theta['samples'] - - return bootstrap_theta - - # Add theta_est_bootstrap_blocks - def theta_est_bootstrap_blocks( - self, - bootstrap_samples, - samplesize=None, - replacement=True, - seed=None, - return_samples=False, - ): - """ - Parameter estimation using bootstrap resampling of the data - - Parameters - ---------- - bootstrap_samples: int - Number of bootstrap samples to draw from the data - samplesize: int or None, optional - Size of each bootstrap sample. If samplesize=None, samplesize will be - set to the number of samples in the data - replacement: bool, optional - Sample with or without replacement. Default is True. - seed: int or None, optional - Random seed - return_samples: bool, optional - Return a list of sample numbers used in each bootstrap estimation. - Default is False. - - Returns - ------- - bootstrap_theta: pd.DataFrame - Theta values for each sample and (if return_samples = True) - the sample numbers used in each estimation - """ - - # check if we are using deprecated parmest - if self.pest_deprecated is not None: - return self.pest_deprecated.theta_est_bootstrap( - bootstrap_samples, - samplesize=samplesize, - replacement=replacement, - seed=seed, - return_samples=return_samples, - ) - - assert isinstance(bootstrap_samples, int) - assert isinstance(samplesize, (type(None), int)) - assert isinstance(replacement, bool) - assert isinstance(seed, (type(None), int)) - assert isinstance(return_samples, bool) - - if samplesize is None: - samplesize = len(self.exp_list) - - if seed is not None: - np.random.seed(seed) - - global_list = self._get_sample_list(samplesize, bootstrap_samples, replacement) - - task_mgr = utils.ParallelTaskManager(bootstrap_samples) - local_list = task_mgr.global_to_local_data(global_list) - bootstrap_theta = list() for idx, sample in local_list: objval, thetavals = self._Q_opt_blocks(bootlist=list(sample)) @@ -2123,67 +1785,6 @@ def theta_est_leaveNout( task_mgr = utils.ParallelTaskManager(len(global_list)) local_list = task_mgr.global_to_local_data(global_list) - lNo_theta = list() - for idx, sample in local_list: - objval, thetavals = self._Q_opt(bootlist=list(sample)) - lNo_s = list(set(range(len(self.exp_list))) - set(sample)) - thetavals['lNo'] = np.sort(lNo_s) - lNo_theta.append(thetavals) - - global_bootstrap_theta = task_mgr.allgather_global_data(lNo_theta) - lNo_theta = pd.DataFrame(global_bootstrap_theta) - - if not return_samples: - del lNo_theta['lNo'] - - return lNo_theta - - def theta_est_leaveNout_blocks( - self, lNo, lNo_samples=None, seed=None, return_samples=False - ): - """ - Parameter estimation where N data points are left out of each sample - - Parameters - ---------- - lNo: int - Number of data points to leave out for parameter estimation - lNo_samples: int - Number of leave-N-out samples. If lNo_samples=None, the maximum - number of combinations will be used - seed: int or None, optional - Random seed - return_samples: bool, optional - Return a list of sample numbers that were left out. Default is False. - - Returns - ------- - lNo_theta: pd.DataFrame - Theta values for each sample and (if return_samples = True) - the sample numbers left out of each estimation - """ - - # check if we are using deprecated parmest - if self.pest_deprecated is not None: - return self.pest_deprecated.theta_est_leaveNout( - lNo, lNo_samples=lNo_samples, seed=seed, return_samples=return_samples - ) - - assert isinstance(lNo, int) - assert isinstance(lNo_samples, (type(None), int)) - assert isinstance(seed, (type(None), int)) - assert isinstance(return_samples, bool) - - samplesize = len(self.exp_list) - lNo - - if seed is not None: - np.random.seed(seed) - - global_list = self._get_sample_list(samplesize, lNo_samples, replacement=False) - - task_mgr = utils.ParallelTaskManager(len(global_list)) - local_list = task_mgr.global_to_local_data(global_list) - lNo_theta = list() for idx, sample in local_list: objval, thetavals = self._Q_opt_blocks(bootlist=list(sample)) @@ -2263,7 +1864,6 @@ def leaveNout_bootstrap_test( results = [] for idx, sample in global_list: - obj, theta = self.theta_est() bootstrap_theta = self.theta_est_bootstrap(bootstrap_samples, seed=seed) @@ -2825,13 +2425,14 @@ def _Q_opt( for ndname, Var, solval in ef_nonants(ef): ind_vars.append(Var) # calculate the reduced hessian - (solve_result, inv_red_hes) = ( - inverse_reduced_hessian.inv_reduced_hessian_barrier( - self.ef_instance, - independent_variables=ind_vars, - solver_options=self.solver_options, - tee=self.tee, - ) + ( + solve_result, + inv_red_hes, + ) = inverse_reduced_hessian.inv_reduced_hessian_barrier( + self.ef_instance, + independent_variables=ind_vars, + solver_options=self.solver_options, + tee=self.tee, ) if self.diagnostic_mode: @@ -3021,10 +2622,14 @@ def _Q_at_theta(self, thetavals, initialize_parmest_model=False): if self.diagnostic_mode: print(' Experiment = ', snum) print(' First solve with special diagnostics wrapper') - (status_obj, solved, iters, time, regu) = ( - utils.ipopt_solve_with_stats( - instance, optimizer, max_iter=500, max_cpu_time=120 - ) + ( + status_obj, + solved, + iters, + time, + regu, + ) = utils.ipopt_solve_with_stats( + instance, optimizer, max_iter=500, max_cpu_time=120 ) print( " status_obj, solved, iters, time, regularization_stat = ", From 7d93cc0c05ef6e99ed56ccc03c5576f48b2b7f1a Mon Sep 17 00:00:00 2001 From: Stephen Cini <114932899+sscini@users.noreply.github.com> Date: Mon, 15 Dec 2025 09:45:01 -0500 Subject: [PATCH 12/15] Ran black on mac --- pyomo/contrib/parmest/parmest.py | 54 +++++++++++++------------------- 1 file changed, 22 insertions(+), 32 deletions(-) diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py index 98523eb219c..8122d93d28f 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -1119,14 +1119,13 @@ def _cov_at_theta(self, method, solver, step): for nd_name, Var, sol_val in ef_nonants(self.ef_instance): ind_vars.append(Var) # calculate the reduced hessian - ( - solve_result, - inv_red_hes, - ) = inverse_reduced_hessian.inv_reduced_hessian_barrier( - self.ef_instance, - independent_variables=ind_vars, - solver_options=self.solver_options, - tee=self.tee, + (solve_result, inv_red_hes) = ( + inverse_reduced_hessian.inv_reduced_hessian_barrier( + self.ef_instance, + independent_variables=ind_vars, + solver_options=self.solver_options, + tee=self.tee, + ) ) self.inv_red_hes = inv_red_hes @@ -1415,14 +1414,10 @@ def _Q_at_theta(self, thetavals, initialize_parmest_model=False): if self.diagnostic_mode: print(' Experiment = ', snum) print(' First solve with special diagnostics wrapper') - ( - status_obj, - solved, - iters, - time, - regu, - ) = utils.ipopt_solve_with_stats( - instance, optimizer, max_iter=500, max_cpu_time=120 + (status_obj, solved, iters, time, regu) = ( + utils.ipopt_solve_with_stats( + instance, optimizer, max_iter=500, max_cpu_time=120 + ) ) print( " status_obj, solved, iters, time, regularization_stat = ", @@ -2425,14 +2420,13 @@ def _Q_opt( for ndname, Var, solval in ef_nonants(ef): ind_vars.append(Var) # calculate the reduced hessian - ( - solve_result, - inv_red_hes, - ) = inverse_reduced_hessian.inv_reduced_hessian_barrier( - self.ef_instance, - independent_variables=ind_vars, - solver_options=self.solver_options, - tee=self.tee, + (solve_result, inv_red_hes) = ( + inverse_reduced_hessian.inv_reduced_hessian_barrier( + self.ef_instance, + independent_variables=ind_vars, + solver_options=self.solver_options, + tee=self.tee, + ) ) if self.diagnostic_mode: @@ -2622,14 +2616,10 @@ def _Q_at_theta(self, thetavals, initialize_parmest_model=False): if self.diagnostic_mode: print(' Experiment = ', snum) print(' First solve with special diagnostics wrapper') - ( - status_obj, - solved, - iters, - time, - regu, - ) = utils.ipopt_solve_with_stats( - instance, optimizer, max_iter=500, max_cpu_time=120 + (status_obj, solved, iters, time, regu) = ( + utils.ipopt_solve_with_stats( + instance, optimizer, max_iter=500, max_cpu_time=120 + ) ) print( " status_obj, solved, iters, time, regularization_stat = ", From d7d22143f83e232f15888750fe25bebfcb953d37 Mon Sep 17 00:00:00 2001 From: Stephen Cini <114932899+sscini@users.noreply.github.com> Date: Mon, 15 Dec 2025 10:09:46 -0500 Subject: [PATCH 13/15] Revert "Removed _Q_opt, and replicate functions, only using _Q_opt_blocks" This reverts commit 1aea99f4f2ea82a46d2b62459a67cd30693e78a0. --- pyomo/contrib/parmest/parmest.py | 407 ++++++++++++++++++++++++++++++- 1 file changed, 406 insertions(+), 1 deletion(-) diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py index 8122d93d28f..4fe1e12b5e9 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -786,6 +786,7 @@ def __init__( diagnostic_mode=False, solver_options=None, ): + # check that we have a (non-empty) list of experiments assert isinstance(experiment_list, list) self.exp_list = experiment_list @@ -849,6 +850,7 @@ def _deprecated_init( diagnostic_mode=False, solver_options=None, ): + deprecation_warning( "You're using the deprecated parmest interface (model_function, " "data, theta_names). This interface will be removed in a future release, " @@ -871,22 +873,26 @@ def _return_theta_names(self): """ # check for deprecated inputs if self.pest_deprecated: + # if fitted model parameter names differ from theta_names # created when Estimator object is created if hasattr(self, 'theta_names_updated'): return self.pest_deprecated.theta_names_updated else: + # default theta_names, created when Estimator object is created return self.pest_deprecated.theta_names else: + # if fitted model parameter names differ from theta_names # created when Estimator object is created if hasattr(self, 'theta_names_updated'): return self.theta_names_updated else: + # default theta_names, created when Estimator object is created return self.estimator_theta_names @@ -1040,7 +1046,6 @@ def total_obj(m): # Redesigning version of _Q_opt that uses scenario blocks # Works, but still adding features from old _Q_opt - # @Reviewers: Trying to find best way to integrate the ability to fix thetas def _Q_opt_blocks( self, return_values=None, @@ -1091,6 +1096,198 @@ def _Q_opt_blocks( return obj_value, theta_estimates + def _Q_opt( + self, + ThetaVals=None, + solver="ef_ipopt", + return_values=[], + bootlist=None, + calc_cov=NOTSET, + cov_n=NOTSET, + ): + """ + Set up all thetas as first stage Vars, return resulting theta + values as well as the objective function value. + + """ + if solver == "k_aug": + raise RuntimeError("k_aug no longer supported.") + + # (Bootstrap scenarios will use indirection through the bootlist) + if bootlist is None: + scenario_numbers = list(range(len(self.exp_list))) + scen_names = ["Scenario{}".format(i) for i in scenario_numbers] + else: + scen_names = ["Scenario{}".format(i) for i in range(len(bootlist))] + + # get the probability constant that is applied to the objective function + # parmest solves the estimation problem by applying equal probabilities to + # the objective function of all the scenarios from the experiment list + self.obj_probability_constant = len(scen_names) + + # tree_model.CallbackModule = None + outer_cb_data = dict() + outer_cb_data["callback"] = self._instance_creation_callback + if ThetaVals is not None: + outer_cb_data["ThetaVals"] = ThetaVals + if bootlist is not None: + outer_cb_data["BootList"] = bootlist + outer_cb_data["cb_data"] = None # None is OK + outer_cb_data["theta_names"] = self.estimator_theta_names + + options = {"solver": "ipopt"} + scenario_creator_options = {"cb_data": outer_cb_data} + if use_mpisppy: + ef = sputils.create_EF( + scen_names, + _experiment_instance_creation_callback, + EF_name="_Q_opt", + suppress_warnings=True, + scenario_creator_kwargs=scenario_creator_options, + ) + else: + ef = local_ef.create_EF( + scen_names, + _experiment_instance_creation_callback, + EF_name="_Q_opt", + suppress_warnings=True, + scenario_creator_kwargs=scenario_creator_options, + ) + self.ef_instance = ef + + # Solve the extensive form with ipopt + if solver == "ef_ipopt": + if calc_cov is NOTSET or not calc_cov: + # Do not calculate the reduced hessian + + solver = SolverFactory('ipopt') + if self.solver_options is not None: + for key in self.solver_options: + solver.options[key] = self.solver_options[key] + + solve_result = solver.solve(self.ef_instance, tee=self.tee) + assert_optimal_termination(solve_result) + elif calc_cov is not NOTSET and calc_cov: + # parmest makes the fitted parameters stage 1 variables + ind_vars = [] + for nd_name, Var, sol_val in ef_nonants(ef): + ind_vars.append(Var) + # calculate the reduced hessian + (solve_result, inv_red_hes) = ( + inverse_reduced_hessian.inv_reduced_hessian_barrier( + self.ef_instance, + independent_variables=ind_vars, + solver_options=self.solver_options, + tee=self.tee, + ) + ) + + if self.diagnostic_mode: + print( + ' Solver termination condition = ', + str(solve_result.solver.termination_condition), + ) + + # assume all first stage are thetas... + theta_vals = {} + for nd_name, Var, sol_val in ef_nonants(ef): + # process the name + # the scenarios are blocks, so strip the scenario name + var_name = Var.name[Var.name.find(".") + 1 :] + theta_vals[var_name] = sol_val + + obj_val = pyo.value(ef.EF_Obj) + self.obj_value = obj_val + self.estimated_theta = theta_vals + + if calc_cov is not NOTSET and calc_cov: + # Calculate the covariance matrix + + if not isinstance(cov_n, int): + raise TypeError( + f"Expected an integer for the 'cov_n' argument. " + f"Got {type(cov_n)}." + ) + num_unknowns = max( + [ + len(experiment.get_labeled_model().unknown_parameters) + for experiment in self.exp_list + ] + ) + assert cov_n > num_unknowns, ( + "The number of datapoints must be greater than the " + "number of parameters to estimate." + ) + + # Number of data points considered + n = cov_n + + # Extract number of fitted parameters + l = len(theta_vals) + + # Assumption: Objective value is sum of squared errors + sse = obj_val + + '''Calculate covariance assuming experimental observation errors + are independent and follow a Gaussian distribution + with constant variance. + + The formula used in parmest was verified against equations + (7-5-15) and (7-5-16) in "Nonlinear Parameter Estimation", + Y. Bard, 1974. + + This formula is also applicable if the objective is scaled by a + constant; the constant cancels out. + (was scaled by 1/n because it computes an expected value.) + ''' + cov = 2 * sse / (n - l) * inv_red_hes + cov = pd.DataFrame( + cov, index=theta_vals.keys(), columns=theta_vals.keys() + ) + + theta_vals = pd.Series(theta_vals) + + if len(return_values) > 0: + var_values = [] + if len(scen_names) > 1: # multiple scenarios + block_objects = self.ef_instance.component_objects( + Block, descend_into=False + ) + else: # single scenario + block_objects = [self.ef_instance] + for exp_i in block_objects: + vals = {} + for var in return_values: + exp_i_var = exp_i.find_component(str(var)) + if ( + exp_i_var is None + ): # we might have a block such as _mpisppy_data + continue + # if value to return is ContinuousSet + if type(exp_i_var) == ContinuousSet: + temp = list(exp_i_var) + else: + temp = [pyo.value(_) for _ in exp_i_var.values()] + if len(temp) == 1: + vals[var] = temp[0] + else: + vals[var] = temp + if len(vals) > 0: + var_values.append(vals) + var_values = pd.DataFrame(var_values) + if calc_cov is not NOTSET and calc_cov: + return obj_val, theta_vals, var_values, cov + elif calc_cov is NOTSET or not calc_cov: + return obj_val, theta_vals, var_values + + if calc_cov is not NOTSET and calc_cov: + return obj_val, theta_vals, cov + elif calc_cov is NOTSET or not calc_cov: + return obj_val, theta_vals + + else: + raise RuntimeError("Unknown solver in Q_Opt=" + solver) + def _cov_at_theta(self, method, solver, step): """ Covariance matrix calculation using all scenarios in the data @@ -1580,6 +1777,77 @@ def theta_est( assert isinstance(return_values, list) assert (calc_cov is NOTSET) or isinstance(calc_cov, bool) + if calc_cov is not NOTSET: + deprecation_warning( + "theta_est(): `calc_cov` and `cov_n` are deprecated options and " + "will be removed in the future. Please use the `cov_est()` function " + "for covariance calculation.", + version="6.9.5", + ) + else: + calc_cov = False + + # check if we are using deprecated parmest + if self.pest_deprecated is not None and calc_cov: + return self.pest_deprecated.theta_est( + solver=solver, + return_values=return_values, + calc_cov=calc_cov, + cov_n=cov_n, + ) + elif self.pest_deprecated is not None and not calc_cov: + return self.pest_deprecated.theta_est( + solver=solver, return_values=return_values + ) + + return self._Q_opt( + solver=solver, + return_values=return_values, + bootlist=None, + calc_cov=calc_cov, + cov_n=cov_n, + ) + + # Replicate of theta_est for testing simplified _Q_opt + # Still work in progress + def theta_est_blocks( + self, solver="ipopt", return_values=[], calc_cov=NOTSET, cov_n=NOTSET + ): + """ + Parameter estimation using all scenarios in the data + + Parameters + ---------- + solver: str, optional + Currently only "ef_ipopt" is supported. Default is "ef_ipopt". + return_values: list, optional + List of Variable names, used to return values from the model + for data reconciliation + calc_cov: boolean, optional + DEPRECATED. + + If True, calculate and return the covariance matrix + (only for "ef_ipopt" solver). Default is NOTSET + cov_n: int, optional + DEPRECATED. + + If calc_cov=True, then the user needs to supply the number of datapoints + that are used in the objective function. Default is NOTSET + + Returns + ------- + obj_val: float + The objective function value + theta_vals: pd.Series + Estimated values for theta + var_values: pd.DataFrame + Variable values for each variable name in + return_values (only for solver='ipopt') + """ + assert isinstance(solver, str) + assert isinstance(return_values, list) + assert (calc_cov is NOTSET) or isinstance(calc_cov, bool) + if calc_cov is not NOTSET: deprecation_warning( "theta_est(): `calc_cov` and `cov_n` are deprecated options and " @@ -1720,6 +1988,81 @@ def theta_est_bootstrap( task_mgr = utils.ParallelTaskManager(bootstrap_samples) local_list = task_mgr.global_to_local_data(global_list) + bootstrap_theta = list() + for idx, sample in local_list: + objval, thetavals = self._Q_opt(bootlist=list(sample)) + thetavals['samples'] = sample + bootstrap_theta.append(thetavals) + + global_bootstrap_theta = task_mgr.allgather_global_data(bootstrap_theta) + bootstrap_theta = pd.DataFrame(global_bootstrap_theta) + + if not return_samples: + del bootstrap_theta['samples'] + + return bootstrap_theta + + # Add theta_est_bootstrap_blocks + def theta_est_bootstrap_blocks( + self, + bootstrap_samples, + samplesize=None, + replacement=True, + seed=None, + return_samples=False, + ): + """ + Parameter estimation using bootstrap resampling of the data + + Parameters + ---------- + bootstrap_samples: int + Number of bootstrap samples to draw from the data + samplesize: int or None, optional + Size of each bootstrap sample. If samplesize=None, samplesize will be + set to the number of samples in the data + replacement: bool, optional + Sample with or without replacement. Default is True. + seed: int or None, optional + Random seed + return_samples: bool, optional + Return a list of sample numbers used in each bootstrap estimation. + Default is False. + + Returns + ------- + bootstrap_theta: pd.DataFrame + Theta values for each sample and (if return_samples = True) + the sample numbers used in each estimation + """ + + # check if we are using deprecated parmest + if self.pest_deprecated is not None: + return self.pest_deprecated.theta_est_bootstrap( + bootstrap_samples, + samplesize=samplesize, + replacement=replacement, + seed=seed, + return_samples=return_samples, + ) + + assert isinstance(bootstrap_samples, int) + assert isinstance(samplesize, (type(None), int)) + assert isinstance(replacement, bool) + assert isinstance(seed, (type(None), int)) + assert isinstance(return_samples, bool) + + if samplesize is None: + samplesize = len(self.exp_list) + + if seed is not None: + np.random.seed(seed) + + global_list = self._get_sample_list(samplesize, bootstrap_samples, replacement) + + task_mgr = utils.ParallelTaskManager(bootstrap_samples) + local_list = task_mgr.global_to_local_data(global_list) + bootstrap_theta = list() for idx, sample in local_list: objval, thetavals = self._Q_opt_blocks(bootlist=list(sample)) @@ -1780,6 +2123,67 @@ def theta_est_leaveNout( task_mgr = utils.ParallelTaskManager(len(global_list)) local_list = task_mgr.global_to_local_data(global_list) + lNo_theta = list() + for idx, sample in local_list: + objval, thetavals = self._Q_opt(bootlist=list(sample)) + lNo_s = list(set(range(len(self.exp_list))) - set(sample)) + thetavals['lNo'] = np.sort(lNo_s) + lNo_theta.append(thetavals) + + global_bootstrap_theta = task_mgr.allgather_global_data(lNo_theta) + lNo_theta = pd.DataFrame(global_bootstrap_theta) + + if not return_samples: + del lNo_theta['lNo'] + + return lNo_theta + + def theta_est_leaveNout_blocks( + self, lNo, lNo_samples=None, seed=None, return_samples=False + ): + """ + Parameter estimation where N data points are left out of each sample + + Parameters + ---------- + lNo: int + Number of data points to leave out for parameter estimation + lNo_samples: int + Number of leave-N-out samples. If lNo_samples=None, the maximum + number of combinations will be used + seed: int or None, optional + Random seed + return_samples: bool, optional + Return a list of sample numbers that were left out. Default is False. + + Returns + ------- + lNo_theta: pd.DataFrame + Theta values for each sample and (if return_samples = True) + the sample numbers left out of each estimation + """ + + # check if we are using deprecated parmest + if self.pest_deprecated is not None: + return self.pest_deprecated.theta_est_leaveNout( + lNo, lNo_samples=lNo_samples, seed=seed, return_samples=return_samples + ) + + assert isinstance(lNo, int) + assert isinstance(lNo_samples, (type(None), int)) + assert isinstance(seed, (type(None), int)) + assert isinstance(return_samples, bool) + + samplesize = len(self.exp_list) - lNo + + if seed is not None: + np.random.seed(seed) + + global_list = self._get_sample_list(samplesize, lNo_samples, replacement=False) + + task_mgr = utils.ParallelTaskManager(len(global_list)) + local_list = task_mgr.global_to_local_data(global_list) + lNo_theta = list() for idx, sample in local_list: objval, thetavals = self._Q_opt_blocks(bootlist=list(sample)) @@ -1859,6 +2263,7 @@ def leaveNout_bootstrap_test( results = [] for idx, sample in global_list: + obj, theta = self.theta_est() bootstrap_theta = self.theta_est_bootstrap(bootstrap_samples, seed=seed) From 7f21344e16e525a1141683a6e1895c6701e94d16 Mon Sep 17 00:00:00 2001 From: Stephen Cini <114932899+sscini@users.noreply.github.com> Date: Wed, 17 Dec 2025 11:26:19 -0500 Subject: [PATCH 14/15] Added testing statement --- pyomo/contrib/parmest/parmest.py | 2 +- pyomo/contrib/parmest/tests/test_parmest.py | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py index 4fe1e12b5e9..ffb7873b574 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -2266,7 +2266,7 @@ def leaveNout_bootstrap_test( obj, theta = self.theta_est() - bootstrap_theta = self.theta_est_bootstrap(bootstrap_samples, seed=seed) + bootstrap_theta = self.theta_est_bootstrap_blocks(bootstrap_samples, seed=seed) training, test = self.confidence_region_test( bootstrap_theta, diff --git a/pyomo/contrib/parmest/tests/test_parmest.py b/pyomo/contrib/parmest/tests/test_parmest.py index db71d280f7c..0baf481e035 100644 --- a/pyomo/contrib/parmest/tests/test_parmest.py +++ b/pyomo/contrib/parmest/tests/test_parmest.py @@ -32,6 +32,7 @@ pynumero_ASL_available = AmplInterface.available() testdir = this_file_dir() +# TESTS HERE WILL BE MODIFIED FOR _Q_OPT_BLOCKS LATER # Set the global seed for random number generation in tests _RANDOM_SEED_FOR_TESTING = 524 From 32d8d414dd7e552b4b3f0eb68c4458354993c35e Mon Sep 17 00:00:00 2001 From: Stephen Cini <114932899+sscini@users.noreply.github.com> Date: Wed, 17 Dec 2025 13:10:57 -0500 Subject: [PATCH 15/15] Ran black --- 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 ffb7873b574..4fe1e12b5e9 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -2266,7 +2266,7 @@ def leaveNout_bootstrap_test( obj, theta = self.theta_est() - bootstrap_theta = self.theta_est_bootstrap_blocks(bootstrap_samples, seed=seed) + bootstrap_theta = self.theta_est_bootstrap(bootstrap_samples, seed=seed) training, test = self.confidence_region_test( bootstrap_theta,