diff --git a/pyomo/contrib/parmest/examples/reaction_kinetics/simple_reaction_parmest_example.py b/pyomo/contrib/parmest/examples/reaction_kinetics/simple_reaction_parmest_example.py index 72eeb606fc2..b130df6ed34 100644 --- a/pyomo/contrib/parmest/examples/reaction_kinetics/simple_reaction_parmest_example.py +++ b/pyomo/contrib/parmest/examples/reaction_kinetics/simple_reaction_parmest_example.py @@ -52,23 +52,6 @@ def simple_reaction_model(data): # fix all of the regressed parameters model.k.fix() - # =================================================================== - # Stage-specific cost computations - def ComputeFirstStageCost_rule(model): - return 0 - - model.FirstStageCost = Expression(rule=ComputeFirstStageCost_rule) - - def AllMeasurements(m): - return (float(data['y']) - m.y) ** 2 - - model.SecondStageCost = Expression(rule=AllMeasurements) - - def total_cost_rule(m): - return m.FirstStageCost + m.SecondStageCost - - model.Total_Cost_Objective = Objective(rule=total_cost_rule, sense=minimize) - return model @@ -90,6 +73,8 @@ def label_model(self): m.experiment_outputs.update( [(m.x1, self.data['x1']), (m.x2, self.data['x2']), (m.y, self.data['y'])] ) + m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.measurement_error.update([(m.y, None), (m.x1, None), (m.x2, None)]) return m @@ -161,7 +146,7 @@ def main(): # Only estimate the parameter k[1]. The parameter k[2] will remain fixed # at its initial value - pest = parmest.Estimator(exp_list) + pest = parmest.Estimator(exp_list, obj_function="SSE") obj, theta = pest.theta_est() print(obj) print(theta) @@ -174,7 +159,7 @@ def main(): # ======================================================================= # Estimate both k1 and k2 and compute the covariance matrix - pest = parmest.Estimator(exp_list) + pest = parmest.Estimator(exp_list, obj_function="SSE") n = 15 # total number of data points used in the objective (y in 15 scenarios) obj, theta, cov = pest.theta_est(calc_cov=True, cov_n=n) print(obj) diff --git a/pyomo/contrib/parmest/examples/reactor_design/parameter_estimation_example.py b/pyomo/contrib/parmest/examples/reactor_design/parameter_estimation_example.py index 630802ddba9..451207f3af0 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/parameter_estimation_example.py +++ b/pyomo/contrib/parmest/examples/reactor_design/parameter_estimation_example.py @@ -34,9 +34,12 @@ def main(): pest = parmest.Estimator(exp_list, obj_function='SSE') # Parameter estimation with covariance - obj, theta, cov = pest.theta_est(calc_cov=True, cov_n=17) - print(obj) + obj, theta, cov = pest.theta_est(calc_cov=True, cov_n=19) + print("Least squares objective value:", obj) + print("Estimated parameters (theta):\n") print(theta) + print("Covariance matrix:\n") + print(cov) if __name__ == "__main__": diff --git a/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py b/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py index 1b3360a2d93..6172096e1ad 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py +++ b/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py @@ -22,15 +22,11 @@ def reactor_design_model(): # Create the concrete model model = pyo.ConcreteModel() - # Rate constants - model.k1 = pyo.Param( - initialize=5.0 / 6.0, within=pyo.PositiveReals, mutable=True - ) # min^-1 - model.k2 = pyo.Param( - initialize=5.0 / 3.0, within=pyo.PositiveReals, mutable=True - ) # min^-1 - model.k3 = pyo.Param( - initialize=1.0 / 6000.0, within=pyo.PositiveReals, mutable=True + # Rate constants, make unknown parameters variables + model.k1 = pyo.Var(initialize=5.0 / 6.0, within=pyo.PositiveReals) # min^-1 + model.k2 = pyo.Var(initialize=5.0 / 3.0, within=pyo.PositiveReals) # min^-1 + model.k3 = pyo.Var( + initialize=1.0 / 6000.0, within=pyo.PositiveReals ) # m^3/(gmol min) # Inlet concentration of A, gmol/m^3 @@ -119,6 +115,11 @@ def label_model(self): (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, None), (m.cb, None), (m.cc, None), (m.cd, None)] + ) + return m def get_labeled_model(self): diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py index 21138d472da..c1292a8cfe1 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -11,6 +11,7 @@ #### Wrapping mpi-sppy functionality and local option Jan 2021, Feb 2021 #### Redesign with Experiment class Dec 2023 +# Options for using mpi-sppy or local EF only used in the deprecatedEstimator class # TODO: move use_mpisppy to a Pyomo configuration option # False implies always use the EF that is local to parmest use_mpisppy = True # Use it if we can but use local if not. @@ -80,6 +81,7 @@ logger = logging.getLogger(__name__) +# Only used in the deprecatedEstimator class def ef_nonants(ef): # Wrapper to call someone's ef_nonants # (the function being called is very short, but it might be changed) @@ -89,6 +91,7 @@ def ef_nonants(ef): return local_ef.ef_nonants(ef) +# Only used in the deprecatedEstimator class def _experiment_instance_creation_callback( scenario_name, node_names=None, cb_data=None ): @@ -354,14 +357,28 @@ def _count_total_experiments(experiment_list): Returns ------- - total_number_data : int + total_data_points : int The total number of data points in the list of experiments """ - total_number_data = 0 + total_data_points = 0 + for experiment in experiment_list: - total_number_data += len(experiment.get_labeled_model().experiment_outputs) + output_vars = experiment.get_labeled_model().experiment_outputs + + # 1. Identify the first parent component + # (e.g., the 'ca' Var container itself) + first_var_key = list(output_vars.keys())[0] + first_parent = first_var_key.parent_component() + + # 2. Count only the keys that belong to this specific parent + # This filters out 'cb', 'cc', etc. + first_param_indices = [ + v for v in output_vars.keys() if v.parent_component() is first_parent + ] - return total_number_data + total_data_points += len(first_param_indices) + + return total_data_points class CovarianceMethod(Enum): @@ -969,197 +986,308 @@ def _instance_creation_callback(self, experiment_number=None, cb_data=None): model = self._create_parmest_model(experiment_number) return model - def _Q_opt( - self, - ThetaVals=None, - solver="ef_ipopt", - return_values=[], - bootlist=None, - calc_cov=NOTSET, - cov_n=NOTSET, - ): + def _create_scenario_blocks(self, bootlist=None, theta_vals=None, fix_theta=False): + # Create scenario block structure """ - Set up all thetas as first stage Vars, return resulting theta - values as well as the objective function value. + Create scenario blocks for parameter estimation + Parameters + ---------- + bootlist : list, optional + List of bootstrap experiment numbers to use. If None, use all experiments in exp_list. + Default is None. + theta_vals : dict, optional + Dictionary of theta values to set in the model. If None, use default values from experiment class. + Default is None. + fix_theta : bool, optional + If True, fix the theta values in the model. If False, leave them free. + Default is False. + Returns + ------- + model : ConcreteModel + Pyomo model with scenario blocks for parameter estimation. Contains indexed block for + each experiment in exp_list or bootlist. """ - if solver == "k_aug": - raise RuntimeError("k_aug no longer supported.") + # Utility function for updated _Q_opt + # Make an indexed block of model scenarios, one for each experiment in exp_list + + # Create a parent model to hold scenario blocks + model = self.ef_instance = self._create_parmest_model(0) + expanded_theta_names = self._expand_indexed_unknowns(model) + if fix_theta: + for name in expanded_theta_names: + theta_var = model.find_component(name) + theta_var.fix() + + # Set the number of experiments to use, either from bootlist or all experiments + self.obj_probability_constant = ( + len(bootlist) if bootlist is not None else len(self.exp_list) + ) - # (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))] + # Create indexed block for holding scenario models + model.exp_scenarios = pyo.Block(range(self.obj_probability_constant)) - # 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) + # Otherwise, use all experiments in exp_list + for i in range(self.obj_probability_constant): + # If bootlist is provided, use it to create scenario blocks for specified experiments + if bootlist is not None: + # Create parmest model for experiment i + parmest_model = self._create_parmest_model(bootlist[i]) - # 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 + # Assign parmest model to block + model.exp_scenarios[i].transfer_attributes_from(parmest_model) - 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, + # Otherwise, use all experiments in exp_list + else: + # Create parmest model for experiment i + parmest_model = self._create_parmest_model(i) + if theta_vals is not None: + # Set theta values in the block model + for name in expanded_theta_names: + # Check the name is in the parmest model + if name in theta_vals: + theta_var = parmest_model.find_component(name) + theta_var.set_value(theta_vals[name]) + if fix_theta: + theta_var.fix() + else: + theta_var.unfix() + + # parmest_model.pprint() + # Assign parmest model to block + model.exp_scenarios[i].transfer_attributes_from(parmest_model) + # model.exp_scenarios[i].pprint() + + # Add linking constraints for theta variables between blocks and parent model + for name in expanded_theta_names: + # Constrain the variable in the first block to equal the parent variable + # If fixing theta, do not add linking constraints + parent_theta_var = model.find_component(name) + if not fix_theta: + for i in range(self.obj_probability_constant): + child_theta_var = model.exp_scenarios[i].find_component(name) + model.add_component( + f"Link_{name}_Block{i}_Parent", + pyo.Constraint(expr=child_theta_var == parent_theta_var), ) - ) - - 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 + # Deactivate existing objectives in the parent model and indexed scenarios + for obj in model.component_objects(pyo.Objective): + obj.deactivate() - 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 + # 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()) + / self.obj_probability_constant + ) - # Extract number of fitted parameters - l = len(theta_vals) + model.Obj = pyo.Objective(rule=total_obj, sense=pyo.minimize) - # Assumption: Objective value is sum of squared errors - sse = obj_val + return model - '''Calculate covariance assuming experimental observation errors - are independent and follow a Gaussian distribution - with constant variance. + # Redesigned _Q_opt method using scenario blocks, and combined with + # _Q_at_theta structure. + def _Q_opt( + self, + return_values=None, + bootlist=None, + solver="ef_ipopt", + theta_vals=None, + calc_cov=NOTSET, + cov_n=NOTSET, + fix_theta=False, + ): + ''' + Making new version of _Q_opt that uses scenario blocks, similar to DoE. - The formula used in parmest was verified against equations - (7-5-15) and (7-5-16) in "Nonlinear Parameter Estimation", - Y. Bard, 1974. + 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 - 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() - ) + Parameters + ---------- + return_values : list, optional + List of variable names to return values for. Default is None. + bootlist : list, optional + List of bootstrap experiment numbers to use. If None, use all experiments in exp_list. + Default is None. + theta_vals : dict, optional + Dictionary of theta values to set in the model. If None, use default values from experiment class. + Default is None. + solver : str, optional + Solver to use for optimization. Default is "ef_ipopt". + calc_cov : bool, optional + If True, calculate covariance matrix of estimated parameters. Default is NOTSET. + cov_n : int, optional + Number of data points to use for covariance calculation. Required if calc_cov is True. Default is NOTSET. + fix_theta : bool, optional + If True, fix the theta values in the model. If False, leave them free. + Default is False. + Returns + ------- + If fix_theta is False: + obj_value : float + Objective value at optimal parameter estimates. + theta_estimates : pd.Series + Series of estimated parameter values. + If fix_theta is True: + obj_value : float + Objective value at fixed parameter values. + theta_estimates : dict + Dictionary of fixed parameter values. + WorstStatus : TerminationCondition + Solver termination condition. + + ''' + # Create extended form model with scenario blocks + model = self._create_scenario_blocks( + bootlist=bootlist, theta_vals=theta_vals, fix_theta=fix_theta + ) + expanded_theta_names = self._expand_indexed_unknowns(model) - theta_vals = pd.Series(theta_vals) + # Print model if in diagnostic mode + if self.diagnostic_mode: + print("Parmest _Q_opt model with scenario blocks:") + model.pprint() - 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 + # Check solver and set options + if solver == "k_aug": + raise RuntimeError("k_aug no longer supported.") + if solver == "ef_ipopt": + sol = SolverFactory('ipopt') + # Currently, parmest is only tested with ipopt via ef_ipopt + # No other pyomo solvers have been verified to work with parmest from current release + # to my knowledge. + + # Seeing if other solvers work here. + # else: + # raise RuntimeError("Unknown solver in Q_Opt=" + solver) + + if self.solver_options is not None: + for key in self.solver_options: + sol.options[key] = self.solver_options[key] + + # Solve model + solve_result = sol.solve(model, tee=self.tee) + + # Separate handling of termination conditions for _Q_at_theta vs _Q_opt + # If not fixing theta, ensure optimal termination of the solve to return result + if not fix_theta: + # Ensure optimal termination + assert_optimal_termination(solve_result) + # If fixing theta, capture termination condition if not optimal unless infeasible + else: + # Initialize worst_status to optimal, update if not optimal + worst_status = pyo.TerminationCondition.optimal + # Get termination condition from solve result + status = solve_result.solver.termination_condition + + # In case of fixing theta, just log a warning if not optimal + if status != pyo.TerminationCondition.optimal: + # logger.warning( + # "Solver did not terminate optimally when thetas were fixed. " + # "Termination condition: %s", + # str(status), + # ) + # Unless infeasible, update worst_status + if worst_status != pyo.TerminationCondition.infeasible: + worst_status = status + + # Extract objective value + obj_value = pyo.value(model.Obj) + theta_estimates = {} + # Extract theta estimates from parent model + for name in expanded_theta_names: + # Value returns value in suffix, which does not change after estimation + # Need to use pyo.value to get variable value + theta_estimates[name] = pyo.value(model.find_component(name)) + + self.obj_value = obj_value + self.estimated_theta = theta_estimates + + # If fixing theta, return objective value, theta estimates, and worst status + if fix_theta: + return obj_value, theta_estimates, worst_status + + # Return theta estimates as a pandas Series + theta_estimates = pd.Series(theta_estimates) + + # Extract return values if requested + # Assumes the model components are named the same in each block, and are pyo.Vars. + if return_values is not None and len(return_values) > 0: + var_values = [] + # In the scenario blocks structure, exp_scenarios is an IndexedBlock + exp_blocks = self.ef_instance.exp_scenarios.values() + # Loop over each experiment block and extract requested variable values + for exp_i in exp_blocks: + # In each block, extract requested variables + vals = {} + for var in return_values: + # Find the variable in the block + exp_i_var = exp_i.find_component(str(var)) + # Check if variable exists in the block + if exp_i_var is None: + continue + # Extract value(s) from variable + 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 + # Only append if vals is not empty + if len(vals) > 0: + var_values.append(vals) + # Convert to DataFrame + var_values = pd.DataFrame(var_values) + + # Calculate covariance if requested using cov_est() + if calc_cov is not NOTSET and calc_cov: + + # Check cov_n argument is set correctly + # Needs to be provided + assert cov_n is not NOTSET, ( + "The number of data points 'cov_n' must be provided to calculate " + "the covariance matrix." + ) + # Needs to be an integer + assert isinstance(cov_n, int), ( + f"Expected an integer for the 'cov_n' argument. " f"Got {type(cov_n)}." + ) + # Needs to equal total number of data points across all experiments + # In progress: Adjusting number_exp to be more robust. + # Can be removed in future when cov_n is no longer an input. + # assert cov_n == self.number_exp, ( + # "The number of data points 'cov_n' must equal the total number " + # "of data points across all experiments." + # ) + + # Needs to be greater than number of parameters + n = cov_n # number of data points + l = len(self.estimated_theta) # number of fitted parameters + assert n > l, ( + "The number of data points 'cov_n' must be greater than " + "the number of fitted parameters." + ) - 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 + cov = self.cov_est(method='reduced_hessian') + if return_values is not None and len(return_values) > 0: + return obj_value, theta_estimates, var_values, cov + else: + return obj_value, theta_estimates, cov + if return_values is not None and len(return_values) > 0: + return obj_value, theta_estimates, var_values else: - raise RuntimeError("Unknown solver in Q_Opt=" + solver) + return obj_value, theta_estimates + + # Removed old _Q_opt function def _cov_at_theta(self, method, solver, step): """ @@ -1184,11 +1312,13 @@ def _cov_at_theta(self, method, solver, step): if method == CovarianceMethod.reduced_hessian.value: # compute the inverse reduced hessian to be used # in the "reduced_hessian" method - # parmest makes the fitted parameters stage 1 variables + # retrieve the independent variables (i.e., estimated parameters) ind_vars = [] - for nd_name, Var, sol_val in ef_nonants(self.ef_instance): - ind_vars.append(Var) - # calculate the reduced hessian + for key, _ in self.ef_instance.unknown_parameters.items(): + name = key.name + var = self.ef_instance.find_component(name) + ind_vars.append(var) + solve_result, inv_red_hes = ( inverse_reduced_hessian.inv_reduced_hessian_barrier( self.ef_instance, @@ -1199,6 +1329,43 @@ def _cov_at_theta(self, method, solver, step): ) self.inv_red_hes = inv_red_hes + else: + # if not using the 'reduced_hessian' method, calculate the sum of squared errors + # using 'finite_difference' method or 'automatic_differentiation_kaug' + sse_vals = [] + for experiment in self.exp_list: + model = _get_labeled_model(experiment) + + # fix the value of the unknown parameters to the estimated values + for param in model.unknown_parameters: + param.fix(self.estimated_theta[param.name]) + + # re-solve the model with the estimated parameters + results = pyo.SolverFactory(solver).solve(model, tee=self.tee) + assert_optimal_termination(results) + + # choose and evaluate the sum of squared errors expression + if self.obj_function == ObjectiveType.SSE: + sse_expr = SSE(model) + elif self.obj_function == ObjectiveType.SSE_weighted: + sse_expr = SSE_weighted(model) + else: + raise ValueError( + f"Invalid objective function for covariance calculation. " + f"The covariance matrix can only be calculated using the built-in " + f"objective functions: {[e.value for e in ObjectiveType]}. Supply " + f"the Estimator object one of these built-in objectives and " + f"re-run the code." + ) + + # evaluate the numerical SSE and store it + sse_val = pyo.value(sse_expr) + sse_vals.append(sse_val) + + sse = sum(sse_vals) + logger.info( + f"The sum of squared errors at the estimated parameter(s) is: {sse}" + ) # Number of data points considered n = self.number_exp @@ -1206,42 +1373,6 @@ def _cov_at_theta(self, method, solver, step): # Extract the number of fitted parameters l = len(self.estimated_theta) - # calculate the sum of squared errors at the estimated parameter values - sse_vals = [] - for experiment in self.exp_list: - model = _get_labeled_model(experiment) - - # fix the value of the unknown parameters to the estimated values - for param in model.unknown_parameters: - param.fix(self.estimated_theta[param.name]) - - # re-solve the model with the estimated parameters - results = pyo.SolverFactory(solver).solve(model, tee=self.tee) - assert_optimal_termination(results) - - # choose and evaluate the sum of squared errors expression - if self.obj_function == ObjectiveType.SSE: - sse_expr = SSE(model) - elif self.obj_function == ObjectiveType.SSE_weighted: - sse_expr = SSE_weighted(model) - else: - raise ValueError( - f"Invalid objective function for covariance calculation. " - f"The covariance matrix can only be calculated using the built-in " - f"objective functions: {[e.value for e in ObjectiveType]}. Supply " - f"the Estimator object one of these built-in objectives and " - f"re-run the code." - ) - - # evaluate the numerical SSE and store it - sse_val = pyo.value(sse_expr) - sse_vals.append(sse_val) - - sse = sum(sse_vals) - logger.info( - f"The sum of squared errors at the estimated parameter(s) is: {sse}" - ) - """Calculate covariance assuming experimental observation errors are independent and follow a Gaussian distribution with constant variance. @@ -1264,11 +1395,11 @@ def _cov_at_theta(self, method, solver, step): # check if the user specified 'SSE' or 'SSE_weighted' as the objective function if self.obj_function == ObjectiveType.SSE: # check if the user defined the 'measurement_error' attribute - if hasattr(model, "measurement_error"): + if hasattr(self.ef_instance, "measurement_error"): # get the measurement errors meas_error = [ - model.measurement_error[y_hat] - for y_hat, y in model.experiment_outputs.items() + self.ef_instance.measurement_error[y_hat] + for y_hat, y in self.ef_instance.experiment_outputs.items() ] # check if the user supplied the values of the measurement errors @@ -1340,10 +1471,10 @@ def _cov_at_theta(self, method, solver, step): ) elif self.obj_function == ObjectiveType.SSE_weighted: # check if the user defined the 'measurement_error' attribute - if hasattr(model, "measurement_error"): + if hasattr(self.ef_instance, "measurement_error"): meas_error = [ - model.measurement_error[y_hat] - for y_hat, y in model.experiment_outputs.items() + self.ef_instance.measurement_error[y_hat] + for y_hat, y in self.ef_instance.experiment_outputs.items() ] # check if the user supplied the values for the measurement errors @@ -1380,200 +1511,16 @@ def _cov_at_theta(self, method, solver, step): raise AttributeError( 'Experiment model does not have suffix "measurement_error".' ) - - return cov - - def _Q_at_theta(self, thetavals, initialize_parmest_model=False): - """ - Return the objective function value with fixed theta values. - - Parameters - ---------- - thetavals: dict - A dictionary of theta values. - - initialize_parmest_model: boolean - If True: Solve square problem instance, build extensive form of the model for - parameter estimation, and set flag model_initialized to True. Default is False. - - Returns - ------- - objectiveval: float - The objective function value. - thetavals: dict - A dictionary of all values for theta that were input. - solvertermination: Pyomo TerminationCondition - Tries to return the "worst" solver status across the scenarios. - pyo.TerminationCondition.optimal is the best and - pyo.TerminationCondition.infeasible is the worst. - """ - - optimizer = pyo.SolverFactory('ipopt') - - if len(thetavals) > 0: - dummy_cb = { - "callback": self._instance_creation_callback, - "ThetaVals": thetavals, - "theta_names": self._return_theta_names(), - "cb_data": None, - } else: - dummy_cb = { - "callback": self._instance_creation_callback, - "theta_names": self._return_theta_names(), - "cb_data": None, - } - - if self.diagnostic_mode: - if len(thetavals) > 0: - print(' Compute objective at theta = ', str(thetavals)) - else: - print(' Compute objective at initial theta') - - # start block of code to deal with models with no constraints - # (ipopt will crash or complain on such problems without special care) - instance = _experiment_instance_creation_callback("FOO0", None, dummy_cb) - try: # deal with special problems so Ipopt will not crash - first = next(instance.component_objects(pyo.Constraint, active=True)) - active_constraints = True - except: - active_constraints = False - # end block of code to deal with models with no constraints - - WorstStatus = pyo.TerminationCondition.optimal - totobj = 0 - scenario_numbers = list(range(len(self.exp_list))) - if initialize_parmest_model: - # create dictionary to store pyomo model instances (scenarios) - scen_dict = dict() - - for snum in scenario_numbers: - sname = "scenario_NODE" + str(snum) - instance = _experiment_instance_creation_callback(sname, None, dummy_cb) - model_theta_names = self._expand_indexed_unknowns(instance) - - if initialize_parmest_model: - # list to store fitted parameter names that will be unfixed - # after initialization - theta_init_vals = [] - # use appropriate theta_names member - theta_ref = model_theta_names - - for i, theta in enumerate(theta_ref): - # Use parser in ComponentUID to locate the component - var_cuid = ComponentUID(theta) - var_validate = var_cuid.find_component_on(instance) - if var_validate is None: - logger.warning( - "theta_name %s was not found on the model", (theta) - ) - else: - try: - if len(thetavals) == 0: - var_validate.fix() - else: - var_validate.fix(thetavals[theta]) - theta_init_vals.append(var_validate) - except: - logger.warning( - 'Unable to fix model parameter value for %s (not a Pyomo model Var)', - (theta), - ) - - if active_constraints: - 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 - ) - ) - print( - " status_obj, solved, iters, time, regularization_stat = ", - str(status_obj), - str(solved), - str(iters), - str(time), - str(regu), - ) - - results = optimizer.solve(instance) - if self.diagnostic_mode: - print( - 'standard solve solver termination condition=', - str(results.solver.termination_condition), - ) - - if ( - results.solver.termination_condition - != pyo.TerminationCondition.optimal - ): - # DLW: Aug2018: not distinguishing "middlish" conditions - if WorstStatus != pyo.TerminationCondition.infeasible: - WorstStatus = results.solver.termination_condition - if initialize_parmest_model: - if self.diagnostic_mode: - print( - "Scenario {:d} infeasible with initialized parameter values".format( - snum - ) - ) - else: - if initialize_parmest_model: - if self.diagnostic_mode: - print( - "Scenario {:d} initialization successful with initial parameter values".format( - snum - ) - ) - if initialize_parmest_model: - # unfix parameters after initialization - for theta in theta_init_vals: - theta.unfix() - scen_dict[sname] = instance - else: - if initialize_parmest_model: - # unfix parameters after initialization - for theta in theta_init_vals: - theta.unfix() - scen_dict[sname] = instance - - objobject = getattr(instance, self._second_stage_cost_exp) - objval = pyo.value(objobject) - totobj += objval - - retval = totobj / len(scenario_numbers) # -1?? - if initialize_parmest_model and not hasattr(self, 'ef_instance'): - # create extensive form of the model using scenario dictionary - if len(scen_dict) > 0: - for scen in scen_dict.values(): - scen._mpisppy_probability = 1 / len(scen_dict) - - if use_mpisppy: - EF_instance = sputils._create_EF_from_scen_dict( - scen_dict, - EF_name="_Q_at_theta", - # suppress_warnings=True - ) - else: - EF_instance = local_ef._create_EF_from_scen_dict( - scen_dict, EF_name="_Q_at_theta", nonant_for_fixed_vars=True - ) - - self.ef_instance = EF_instance - # set self.model_initialized flag to True to skip extensive form model - # creation using theta_est() - self.model_initialized = True - - # return initialized theta values - if len(thetavals) == 0: - # use appropriate theta_names member - theta_ref = self._return_theta_names() - for i, theta in enumerate(theta_ref): - thetavals[theta] = theta_init_vals[i]() + raise ValueError( + f"Invalid objective function for covariance calculation. " + f"The covariance matrix can only be calculated using the built-in " + f"objective functions: {[e.value for e in ObjectiveType]}. Supply " + f"the Estimator object one of these built-in objectives and " + f"re-run the code." + ) - return retval, thetavals, WorstStatus + return cov def _get_sample_list(self, samplesize, num_samples, replacement=True): samplelist = list() @@ -1944,9 +1891,11 @@ def leaveNout_bootstrap_test( return results + # Updated version that uses _Q_opt def objective_at_theta(self, theta_values=None, initialize_parmest_model=False): """ - Objective value for each theta + Objective value for each theta, solving extensive form problem with + fixed theta values. Parameters ---------- @@ -1958,7 +1907,6 @@ def objective_at_theta(self, theta_values=None, initialize_parmest_model=False): of the model for parameter estimation, and set flag model_initialized to True. Default is False. - Returns ------- obj_at_theta: pd.DataFrame @@ -1973,65 +1921,66 @@ def objective_at_theta(self, theta_values=None, initialize_parmest_model=False): initialize_parmest_model=initialize_parmest_model, ) - if len(self.estimator_theta_names) == 0: - pass # skip assertion if model has no fitted parameters - else: - # create a local instance of the pyomo model to access model variables and parameters - model_temp = self._create_parmest_model(0) - model_theta_list = self._expand_indexed_unknowns(model_temp) - - # if self.estimator_theta_names is not the same as temp model_theta_list, - # create self.theta_names_updated - if set(self.estimator_theta_names) == set(model_theta_list) and len( - self.estimator_theta_names - ) == len(set(model_theta_list)): - pass - else: - self.theta_names_updated = model_theta_list + if initialize_parmest_model: + # Print deprecation warning, that this option will be removed in + # future releases. + deprecation_warning( + "The `initialize_parmest_model` option in `objective_at_theta()` is " + "deprecated and will be removed in future releases. Please ensure the" + "model is initialized within the Experiment class definition.", + version="6.9.5", + ) if theta_values is None: all_thetas = {} # dictionary to store fitted variables # use appropriate theta names member - theta_names = model_theta_list + # Get theta names from fresh parmest model, assuming this can be called + # directly after creating Estimator. + theta_names = self._expand_indexed_unknowns(self._create_parmest_model(0)) else: assert isinstance(theta_values, pd.DataFrame) # for parallel code we need to use lists and dicts in the loop theta_names = theta_values.columns # # check if theta_names are in model - for theta in list(theta_names): - theta_temp = theta.replace("'", "") # cleaning quotes from theta_names - assert theta_temp in [ - t.replace("'", "") for t in model_theta_list - ], "Theta name {} in 'theta_values' not in 'theta_names' {}".format( - theta_temp, model_theta_list + # Clean names, ignore quotes, and compare sets + clean_provided = [t.replace("'", "") for t in theta_names] + clean_expected = [ + t.replace("'", "") + for t in self._expand_indexed_unknowns(self._create_parmest_model(0)) + ] + # If they do not match, raise error + if set(clean_provided) != set(clean_expected): + raise ValueError( + f"Provided theta values {clean_provided} do not match expected theta names {clean_expected}." ) + # Rename columns using cleaned names + if set(clean_provided) != set(theta_names): + theta_values.columns = clean_provided - assert len(list(theta_names)) == len(model_theta_list) - + # Convert to list of dicts for parallel processing all_thetas = theta_values.to_dict('records') + # Initialize task manager + num_tasks = len(all_thetas) if all_thetas else 1 + task_mgr = utils.ParallelTaskManager(num_tasks) + + # Use local theta values for each task if all_thetas is provided, else empty list if all_thetas: - task_mgr = utils.ParallelTaskManager(len(all_thetas)) local_thetas = task_mgr.global_to_local_data(all_thetas) - else: - if initialize_parmest_model: - task_mgr = utils.ParallelTaskManager( - 1 - ) # initialization performed using just 1 set of theta values + elif initialize_parmest_model: + local_thetas = [] + # walk over the mesh, return objective function all_obj = list() if len(all_thetas) > 0: for Theta in local_thetas: - obj, thetvals, worststatus = self._Q_at_theta( - Theta, initialize_parmest_model=initialize_parmest_model + obj, thetvals, worststatus = self._Q_opt( + theta_vals=Theta, fix_theta=True ) if worststatus != pyo.TerminationCondition.infeasible: all_obj.append(list(Theta.values()) + [obj]) - # DLW, Aug2018: should we also store the worst solver status? else: - obj, thetvals, worststatus = self._Q_at_theta( - thetavals={}, initialize_parmest_model=initialize_parmest_model - ) + obj, thetvals, worststatus = self._Q_opt(theta_vals=None, fix_theta=True) if worststatus != pyo.TerminationCondition.infeasible: all_obj.append(list(thetvals.values()) + [obj]) diff --git a/pyomo/contrib/parmest/tests/test_parmest.py b/pyomo/contrib/parmest/tests/test_parmest.py index dfeab769e71..d514de01eae 100644 --- a/pyomo/contrib/parmest/tests/test_parmest.py +++ b/pyomo/contrib/parmest/tests/test_parmest.py @@ -15,6 +15,8 @@ from pyomo.common.unittest import pytest from parameterized import parameterized, parameterized_class import pyomo.common.unittest as unittest +from pyomo.contrib.mpc import data +from pyomo.contrib.mpc.examples.cstr import model import pyomo.contrib.parmest.parmest as parmest import pyomo.contrib.parmest.graphics as graphics import pyomo.contrib.parmest as parmestbase @@ -30,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 @@ -508,40 +511,6 @@ def test_parallel_parmest(self): retcode = subprocess.call(rlist) self.assertEqual(retcode, 0) - @unittest.skipIf(not pynumero_ASL_available, "pynumero_ASL is not available") - def test_theta_est_cov(self): - objval, thetavals, cov = self.pest.theta_est(calc_cov=True, cov_n=6) - - self.assertAlmostEqual(objval, 4.3317112, places=2) - self.assertAlmostEqual( - thetavals["asymptote"], 19.1426, places=2 - ) # 19.1426 from the paper - self.assertAlmostEqual( - thetavals["rate_constant"], 0.5311, places=2 - ) # 0.5311 from the paper - - # Covariance matrix - self.assertAlmostEqual( - cov["asymptote"]["asymptote"], 6.155892, places=2 - ) # 6.22864 from paper - self.assertAlmostEqual( - cov["asymptote"]["rate_constant"], -0.425232, places=2 - ) # -0.4322 from paper - self.assertAlmostEqual( - cov["rate_constant"]["asymptote"], -0.425232, places=2 - ) # -0.4322 from paper - self.assertAlmostEqual( - cov["rate_constant"]["rate_constant"], 0.040571, places=2 - ) # 0.04124 from paper - - """ Why does the covariance matrix from parmest not match the paper? Parmest is - calculating the exact reduced Hessian. The paper (Rooney and Bielger, 2001) likely - employed the first order approximation common for nonlinear regression. The paper - values were verified with Scipy, which uses the same first order approximation. - The formula used in parmest was verified against equations (7-5-15) and (7-5-16) in - "Nonlinear Parameter Estimation", Y. Bard, 1974. - """ - def test_cov_scipy_least_squares_comparison(self): """ Scipy results differ in the 3rd decimal place from the paper. It is possible @@ -655,20 +624,27 @@ def setUp(self): columns=["hour", "y"], ) + # Updated models to use Vars for experiment output, and Constraints def rooney_biegler_params(data): model = pyo.ConcreteModel() model.asymptote = pyo.Param(initialize=15, mutable=True) model.rate_constant = pyo.Param(initialize=0.5, mutable=True) - model.hour = pyo.Param(within=pyo.PositiveReals, mutable=True) - model.y = pyo.Param(within=pyo.PositiveReals, mutable=True) + # Add the experiment inputs + model.h = pyo.Var(initialize=data["hour"].iloc[0], bounds=(0, 10)) - def response_rule(m, h): - expr = m.asymptote * (1 - pyo.exp(-m.rate_constant * h)) - return expr + # Fix the experiment inputs + model.h.fix() - model.response_function = pyo.Expression(data.hour, rule=response_rule) + # 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.asymptote * (1 - pyo.exp(-m.rate_constant * m.h)) + + model.response_con = pyo.Constraint(rule=response_rule) return model @@ -683,14 +659,14 @@ def label_model(self): m = self.model m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL) - m.experiment_outputs.update( - [(m.hour, self.data["hour"]), (m.y, self.data["y"])] - ) + 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.asymptote, m.rate_constant] ) + m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.measurement_error.update([(m.y, None)]) rooney_biegler_params_exp_list = [] for i in range(self.data.shape[0]): @@ -701,24 +677,30 @@ def label_model(self): def rooney_biegler_indexed_params(data): model = pyo.ConcreteModel() + # Define the indexed parameters model.param_names = pyo.Set(initialize=["asymptote", "rate_constant"]) model.theta = pyo.Param( model.param_names, initialize={"asymptote": 15, "rate_constant": 0.5}, mutable=True, ) + # Add the experiment inputs + model.h = pyo.Var(initialize=data["hour"].iloc[0], bounds=(0, 10)) - model.hour = pyo.Param(within=pyo.PositiveReals, mutable=True) - model.y = pyo.Param(within=pyo.PositiveReals, mutable=True) + # Fix the experiment inputs + model.h.fix() - def response_rule(m, h): - expr = m.theta["asymptote"] * ( - 1 - pyo.exp(-m.theta["rate_constant"] * h) - ) - return expr + # Add experiment outputs + model.y = pyo.Var(initialize=data['y'].iloc[0], within=pyo.PositiveReals) - model.response_function = pyo.Expression(data.hour, rule=response_rule) + # Define the model equations + def response_rule(m): + return m.y == m.theta["asymptote"] * ( + 1 - pyo.exp(-m.theta["rate_constant"] * m.h) + ) + # Add the model equations to the model + model.response_con = pyo.Constraint(rule=response_rule) return model class RooneyBieglerExperimentIndexedParams(RooneyBieglerExperiment): @@ -732,13 +714,14 @@ def label_model(self): m = self.model m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL) - m.experiment_outputs.update( - [(m.hour, self.data["hour"]), (m.y, self.data["y"])] - ) + 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, None)]) + rooney_biegler_indexed_params_exp_list = [] for i in range(self.data.shape[0]): rooney_biegler_indexed_params_exp_list.append( @@ -753,14 +736,20 @@ def rooney_biegler_vars(data): model.asymptote.fixed = True # parmest will unfix theta variables model.rate_constant.fixed = True - model.hour = pyo.Param(within=pyo.PositiveReals, mutable=True) - model.y = pyo.Param(within=pyo.PositiveReals, mutable=True) + # Add the experiment inputs + model.h = pyo.Var(initialize=data["hour"].iloc[0], bounds=(0, 10)) - def response_rule(m, h): - expr = m.asymptote * (1 - pyo.exp(-m.rate_constant * h)) - return expr + # Fix the experiment inputs + model.h.fix() - model.response_function = pyo.Expression(data.hour, rule=response_rule) + # 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.asymptote * (1 - pyo.exp(-m.rate_constant * m.h)) + + model.response_con = pyo.Constraint(rule=response_rule) return model @@ -775,14 +764,14 @@ def label_model(self): m = self.model m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL) - m.experiment_outputs.update( - [(m.hour, self.data["hour"]), (m.y, self.data["y"])] - ) + 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.asymptote, m.rate_constant] ) + m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.measurement_error.update([(m.y, None)]) rooney_biegler_vars_exp_list = [] for i in range(self.data.shape[0]): @@ -802,16 +791,22 @@ def rooney_biegler_indexed_vars(data): ) model.theta["rate_constant"].fixed = True - model.hour = pyo.Param(within=pyo.PositiveReals, mutable=True) - model.y = pyo.Param(within=pyo.PositiveReals, mutable=True) + # Add the experiment inputs + model.h = pyo.Var(initialize=data["hour"].iloc[0], bounds=(0, 10)) - def response_rule(m, h): - expr = m.theta["asymptote"] * ( - 1 - pyo.exp(-m.theta["rate_constant"] * h) + # Fix the experiment inputs + model.h.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.h) ) - return expr - model.response_function = pyo.Expression(data.hour, rule=response_rule) + model.response_con = pyo.Constraint(rule=response_rule) return model @@ -826,28 +821,21 @@ def label_model(self): m = self.model m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL) - m.experiment_outputs.update( - [(m.hour, self.data["hour"]), (m.y, self.data["y"])] - ) + 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, None)]) + rooney_biegler_indexed_vars_exp_list = [] for i in range(self.data.shape[0]): rooney_biegler_indexed_vars_exp_list.append( RooneyBieglerExperimentIndexedVars(self.data.loc[i, :]) ) - # Sum of squared error function - def SSE(model): - expr = ( - model.experiment_outputs[model.y] - - model.response_function[model.experiment_outputs[model.hour]] - ) ** 2 - return expr - - self.objective_function = SSE + self.objective_function = 'SSE' theta_vals = pd.DataFrame([20, 1], index=["asymptote", "rate_constant"]).T theta_vals_index = pd.DataFrame( @@ -899,16 +887,16 @@ def check_rooney_biegler_results(self, objval, cov): self.assertAlmostEqual(objval, 4.3317112, places=2) self.assertAlmostEqual( - cov.iloc[asymptote_index, asymptote_index], 6.30579403, places=2 + cov.iloc[asymptote_index, asymptote_index], 6.155892, places=2 ) # 6.22864 from paper self.assertAlmostEqual( - cov.iloc[asymptote_index, rate_constant_index], -0.4395341, places=2 + 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.4395341, places=2 + 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.04193591, places=2 + 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') @@ -1107,7 +1095,8 @@ def _dccrate(m, t): def ComputeFirstStageCost_rule(m): return 0 - m.FirstStageCost = pyo.Expression(rule=ComputeFirstStageCost_rule) + # Model used in + m.FirstStage = pyo.Expression(rule=ComputeFirstStageCost_rule) def ComputeSecondStageCost_rule(m): return sum( @@ -1117,14 +1106,12 @@ def ComputeSecondStageCost_rule(m): for t in meas_t ) - m.SecondStageCost = pyo.Expression(rule=ComputeSecondStageCost_rule) + m.SecondStage = pyo.Expression(rule=ComputeSecondStageCost_rule) def total_cost_rule(model): - return model.FirstStageCost + model.SecondStageCost + return model.FirstStage + model.SecondStage - m.Total_Cost_Objective = pyo.Objective( - rule=total_cost_rule, sense=pyo.minimize - ) + m.Total_Cost = pyo.Objective(rule=total_cost_rule, sense=pyo.minimize) disc = pyo.TransformationFactory("dae.collocation") disc.apply_to(m, nfe=20, ncp=2) @@ -1165,6 +1152,10 @@ def label_model(self): 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], None) for t in meas_time_points) + m.measurement_error.update((m.cb[t], None) for t in meas_time_points) + m.measurement_error.update((m.cc[t], None) for t in meas_time_points) def get_labeled_model(self): self.create_model() @@ -1210,8 +1201,8 @@ def get_labeled_model(self): exp_list_df = [ReactorDesignExperimentDAE(data_df)] exp_list_dict = [ReactorDesignExperimentDAE(data_dict)] - self.pest_df = parmest.Estimator(exp_list_df) - self.pest_dict = parmest.Estimator(exp_list_dict) + self.pest_df = parmest.Estimator(exp_list_df, obj_function="SSE") + self.pest_dict = parmest.Estimator(exp_list_dict, obj_function="SSE") # Estimator object with multiple scenarios exp_list_df_multiple = [ @@ -1223,8 +1214,12 @@ def get_labeled_model(self): ReactorDesignExperimentDAE(data_dict), ] - self.pest_df_multiple = parmest.Estimator(exp_list_df_multiple) - self.pest_dict_multiple = parmest.Estimator(exp_list_dict_multiple) + self.pest_df_multiple = parmest.Estimator( + exp_list_df_multiple, obj_function="SSE" + ) + self.pest_dict_multiple = parmest.Estimator( + exp_list_dict_multiple, obj_function="SSE" + ) # Create an instance of the model self.m_df = ABC_model(data_df) @@ -1305,6 +1300,7 @@ def test_return_continuous_set_multiple_datasets(self): self.assertAlmostEqual(return_vals1["time"].loc[1][18], 2.368, places=3) self.assertAlmostEqual(return_vals2["time"].loc[1][18], 2.368, places=3) + # Currently failing, _count_total_experiments problem @unittest.skipUnless(pynumero_ASL_available, 'pynumero_ASL is not available') def test_covariance(self): from pyomo.contrib.interior_point.inverse_reduced_hessian import ( @@ -1315,7 +1311,7 @@ def test_covariance(self): # 3 data components (ca, cb, cc), 20 timesteps, 1 scenario = 60 # In this example, this is the number of data points in data_df, but that's # only because the data is indexed by time and contains no additional information. - n = 60 + n = 20 # Compute covariance using parmest obj, theta, cov = self.pest_df.theta_est(calc_cov=True, cov_n=n) diff --git a/pyomo/contrib/parmest/utils/create_ef.py b/pyomo/contrib/parmest/utils/create_ef.py index 56048fced10..17336114668 100644 --- a/pyomo/contrib/parmest/utils/create_ef.py +++ b/pyomo/contrib/parmest/utils/create_ef.py @@ -21,6 +21,7 @@ from pyomo.core import Objective +# File no longer used in parmest; retained for possible future use. def get_objs(scenario_instance): """return the list of objective functions for scenario_instance""" scenario_objs = scenario_instance.component_data_objects( diff --git a/pyomo/contrib/parmest/utils/mpi_utils.py b/pyomo/contrib/parmest/utils/mpi_utils.py index 82376062c45..30560a6a77d 100644 --- a/pyomo/contrib/parmest/utils/mpi_utils.py +++ b/pyomo/contrib/parmest/utils/mpi_utils.py @@ -10,6 +10,8 @@ from collections import OrderedDict import importlib +# ParallelTaskManager is used, MPI Interface is not. + """ This module is a collection of classes that provide a friendlier interface to MPI (through mpi4py). They help