Extended Parmest Capability for weighted SSE objective #3535
Conversation
|
@slilonfe5 Here is some quick feedback
Feedback on the
|
|
@adowling2 @djlaky I also updated the calculation for the normal SSE such that we can use the user-supplied measurement error if defined; otherwise, we calculate the measurement error as usual. |
|
@slilonfe5 Once you have the tests ready, tag us for feedback. Also, I think you can skip adding this to the depreciated class. |
|
@adowling2 @djlaky I have created a separate method ( I tested these with three examples (2 steady state and 1 dynamic), and all work well. I'm yet to write the test file for these. |
Co-authored-by: Bethany Nicholson <[email protected]>
Co-authored-by: Bethany Nicholson <[email protected]>
Co-authored-by: Bethany Nicholson <[email protected]>
Co-authored-by: Bethany Nicholson <[email protected]>
There was a problem hiding this comment.
@slilonfe5 here are a few more review comments. I'm still going through parmest.py
|
@blnicho @mrmundt @jsiirola @djlaky @adowling2 I have implemented Bethany's final review. |
| (0.1, "incorrect_obj"), | ||
| ], | ||
| ) | ||
| class TestRooneyBieglerWSSE(unittest.TestCase): |
There was a problem hiding this comment.
Since this test class is parameterized to test both SSE and WSSE you might want to change the class name to be more generic.
| class TestRooneyBieglerWSSE(unittest.TestCase): | |
| class TestParmestCovEst(unittest.TestCase): |
| Covariance Matrix Estimation | ||
| ================================= | ||
| ============================ |
There was a problem hiding this comment.
It doesn't look like the functions for the "SSE" or "WSSE" objective functions are defined anywhere in the documentation. Also, what are the "weights" in the "WSSE"? You should probably add the functions here and perhaps again in the overview.rst file here:
pyomo/doc/OnlineDocs/explanation/analysis/parmest/overview.rst
Lines 43 to 55 in df215bf
| def test_cov_scipy_least_squares_comparison(self): | ||
| """ | ||
| Estimates the unknown parameters and covariance matrix from the measurement | ||
| error standard deviation using Scipy least_squares function. | ||
| """ | ||
| if self.measurement_std is None or self.objective_function == "incorrect_obj": | ||
| self.skipTest( | ||
| "This test only applies to 'SSE' and 'SSE_weighted' " | ||
| "objectives with user-supplied measurement error" | ||
| ) | ||
|
|
||
| def model(theta, t): | ||
| """ | ||
| Model to be fitted y = model(theta, t) | ||
| Arguments: | ||
| theta: vector of fitted parameters | ||
| t: independent variable [hours] | ||
|
|
||
| Returns: | ||
| y: model predictions [need to check paper for units] | ||
| """ | ||
| asymptote = theta[0] | ||
| rate_constant = theta[1] | ||
|
|
||
| return asymptote * (1 - np.exp(-rate_constant * t)) | ||
|
|
||
| def residual(theta, t, y): | ||
| """ | ||
| Calculate residuals | ||
| Arguments: | ||
| theta: vector of fitted parameters | ||
| t: independent variable [hours] | ||
| y: dependent variable [?] | ||
| """ | ||
| return y - model(theta, t) | ||
|
|
||
| # define data | ||
| t = self.data["hour"].to_numpy() | ||
| y = self.data["y"].to_numpy() | ||
|
|
||
| # define initial guess | ||
| theta_guess = np.array([15, 0.5]) | ||
|
|
||
| ## solve with optimize.least_squares | ||
| sol = scipy.optimize.least_squares( | ||
| residual, theta_guess, method="trf", args=(t, y), verbose=2 | ||
| ) | ||
| theta_hat = sol.x | ||
|
|
||
| self.assertAlmostEqual( | ||
| theta_hat[0], 19.1426, places=2 | ||
| ) # 19.1426 from the paper | ||
| self.assertAlmostEqual(theta_hat[1], 0.5311, places=2) # 0.5311 from the paper | ||
|
|
||
| # calculate the variance of the measurement error | ||
| sigre = self.measurement_std**2 | ||
|
|
||
| cov = sigre * np.linalg.inv(np.matmul(sol.jac.T, sol.jac)) | ||
|
|
||
| self.assertAlmostEqual(cov[0, 0], 0.009588, places=4) | ||
| self.assertAlmostEqual(cov[0, 1], -0.000665, places=4) | ||
| self.assertAlmostEqual(cov[1, 0], -0.000665, places=4) | ||
| self.assertAlmostEqual(cov[1, 1], 0.000063, places=4) |
There was a problem hiding this comment.
Is this test actually testing anything in parmest? I think it's an interesting example as an alternate calculation approach but it doesn't look like it's using parmest or Pyomo so I'm not sure it's appropriate to include here.
| # check if the model has all the required suffixes | ||
| _check_model_labels_helper(model, logging_level) |
There was a problem hiding this comment.
Is there any way to reach the _compute_jacobian function without previously calling either SSE or SSE_weighted? I'm wondering if this call to _check_model_labels_helper is really needed or if we can assume that the model labels have already been checked.
| try: | ||
| solver = pyo.SolverFactory(solver) | ||
| solver.solve(model, tee=tee) | ||
| except Exception as e: | ||
| raise RuntimeError( | ||
| f"Model from experiment did not solve appropriately. Make sure the " | ||
| f"model is well-posed. The original error was {e}." | ||
| ) |
There was a problem hiding this comment.
You usually want to avoid catching a generic Exception and check for the specific things you're worried about. In this case I'm guessing you're trying to catch the error that is thrown when you try to load an unconverged/failed Ipopt solution into the model. Here is an example of how you can check the solver termination condition and avoid the error:
from pyomo.opt import TerminationCondition
...
solver = pyo.SolverFactory(solver_name)
results = solver.solve(model, tee=tee, load_solutions=False)
if results.solver.termination_condition == TerminationCondition.optimal:
model.solutions.load_from(results)
else:
print ("Solution is not optimal")
# now do something about it? raise an exception? ...
| "with the `calc_cov` and `cov_n` arguments. This usage will be " | ||
| "removed in a future release. Please update to the new parmest " | ||
| "interface using `cov_est()` function for covariance calculation.", | ||
| version="6.9.3.dev0", |
There was a problem hiding this comment.
| version="6.9.3.dev0", | |
| version="6.9.4.dev0", |
| "You're using a deprecated call to the `theta_est()` function " | ||
| "with the `calc_cov` and `cov_n` arguments. This usage will be " | ||
| "removed in a future release. Please update to the new parmest " | ||
| "interface using `cov_est()` function for covariance calculation.", |
There was a problem hiding this comment.
| "interface using `cov_est()` function for covariance calculation.", | |
| "interface using the `cov_est()` function for covariance calculation.", |
| for key in self.solver_options: | ||
| solver.options[key] = self.solver_options[key] | ||
|
|
||
| solve_result = solver.solve(self.ef_instance, tee=self.tee) |
There was a problem hiding this comment.
Should this call to solve be guarded in case of solver failure?
| Argument: | ||
| method: string ``method`` object specified by the user, | ||
| e.g., 'finite_difference' | ||
| solver: string ``solver`` object specified by the user, e.g., 'ipopt' | ||
| step: float used for relative perturbation of the parameters, | ||
| e.g., step=0.02 is a 2% perturbation | ||
|
|
||
| Returns: | ||
| cov: pd.DataFrame, covariance matrix of the estimated parameters |
There was a problem hiding this comment.
Same comment about using the numpydoc format for docstrings
| try: | ||
| pyo.SolverFactory(solver).solve(model, tee=self.tee) | ||
| except Exception as e: | ||
| raise RuntimeError( | ||
| f"Model from experiment did not solve appropriately. Make sure the " | ||
| f"model is well-posed. The original error was {e}." | ||
| ) |
There was a problem hiding this comment.
Same comment about catching a generic Exception
Fixes # .
Summary/Motivation:
Currently, the Parmest SSE objective does not support measurements in different units. This work adds new capabilities to Parmest, including weighted SSE to handle measurements in different units, and more robust covariance matrix calculation methods for more accurate uncertainty quantification. This work also enables the calculation of the covariance matrix using a user-supplied measurement error standard deviation.
Changes proposed in this PR:
Legal Acknowledgement
By contributing to this software project, I have read the contribution guide and agree to the following terms and conditions for my contribution: