diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py index 12599de7c60..b288670d2d0 100644 --- a/pyomo/contrib/doe/doe.py +++ b/pyomo/contrib/doe/doe.py @@ -49,6 +49,19 @@ from pyomo.opt import SolverStatus +# This small and positive tolerance is used when checking +# if the prior is negative definite or approximately +# indefinite. It is defined as a tolerance here to ensure +# consistency between the code below and the tests. The +# user should not need to adjust it. +_SMALL_TOLERANCE_DEFINITENESS = 1e-6 + +# This small and positive tolerance is used to check +# the FIM is approximately symmetric. It is defined as +# a tolerance here to ensure consistency between the code +# below and the tests. The user should not need to adjust it. +_SMALL_TOLERANCE_SYMMETRY = 1e-6 + class ObjectiveLib(Enum): determinant = "determinant" @@ -309,7 +322,31 @@ def run_doe(self, model=None, results_file=None): (len(model.parameter_names), len(model.parameter_names)) ) - L_vals_sq = np.linalg.cholesky(fim_np) + # Need to compute the full FIM before initializing the Cholesky factorization + if self.only_compute_fim_lower: + fim_np = fim_np + fim_np.T - np.diag(np.diag(fim_np)) + + # Check if the FIM is positive definite + # If not, add jitter to the diagonal + # to ensure positive definiteness + min_eig = np.min(np.linalg.eigvals(fim_np)) + + if min_eig < _SMALL_TOLERANCE_DEFINITENESS: + # Raise the minimum eigenvalue to at least _SMALL_TOLERANCE_DEFINITENESS + jitter = np.min( + [ + -min_eig + _SMALL_TOLERANCE_DEFINITENESS, + _SMALL_TOLERANCE_DEFINITENESS, + ] + ) + else: + # No jitter needed + jitter = 0 + + # Add jitter to the diagonal to ensure positive definiteness + L_vals_sq = np.linalg.cholesky( + fim_np + jitter * np.eye(len(model.parameter_names)) + ) for i, c in enumerate(model.parameter_names): for j, d in enumerate(model.parameter_names): model.L[c, d].value = L_vals_sq[i, j] @@ -1338,7 +1375,28 @@ def check_model_FIM(self, model=None, FIM=None): ) ) - self.logger.info("FIM provided matches expected dimensions from model.") + # Compute the eigenvalues of the FIM + evals = np.linalg.eigvals(FIM) + + # Check if the FIM is positive definite + if np.min(evals) < -_SMALL_TOLERANCE_DEFINITENESS: + raise ValueError( + "FIM provided is not positive definite. It has one or more negative eigenvalue(s) less than -{:.1e}".format( + _SMALL_TOLERANCE_DEFINITENESS + ) + ) + + # Check if the FIM is symmetric + if not np.allclose(FIM, FIM.T, atol=_SMALL_TOLERANCE_SYMMETRY): + raise ValueError( + "FIM provided is not symmetric using absolute tolerance {}".format( + _SMALL_TOLERANCE_SYMMETRY + ) + ) + + self.logger.info( + "FIM provided matches expected dimensions from model and is approximately positive (semi) definite." + ) # Check the jacobian shape against what is expected from the model. def check_model_jac(self, jac=None): diff --git a/pyomo/contrib/doe/tests/test_doe_errors.py b/pyomo/contrib/doe/tests/test_doe_errors.py index 611cc7c31e9..4c9823a251d 100644 --- a/pyomo/contrib/doe/tests/test_doe_errors.py +++ b/pyomo/contrib/doe/tests/test_doe_errors.py @@ -165,6 +165,55 @@ def test_reactor_check_bad_prior_size(self): ): doe_obj.create_doe_model() + def test_reactor_check_bad_prior_negative_eigenvalue(self): + from pyomo.contrib.doe.doe import _SMALL_TOLERANCE_DEFINITENESS + + fd_method = "central" + obj_used = "trace" + flag_val = 0 # Value for faulty model build mode - 0: full model + + prior_FIM = -np.ones((4, 4)) + + experiment = FullReactorExperiment(data_ex, 10, 3) + + DoE_args = get_standard_args(experiment, fd_method, obj_used, flag_val) + DoE_args['prior_FIM'] = prior_FIM + + doe_obj = DesignOfExperiments(**DoE_args) + + with self.assertRaisesRegex( + ValueError, + r"FIM provided is not positive definite. It has one or more negative eigenvalue\(s\) less than -{:.1e}".format( + _SMALL_TOLERANCE_DEFINITENESS + ), + ): + doe_obj.create_doe_model() + + def test_reactor_check_bad_prior_not_symmetric(self): + from pyomo.contrib.doe.doe import _SMALL_TOLERANCE_SYMMETRY + + fd_method = "central" + obj_used = "trace" + flag_val = 0 # Value for faulty model build mode - 0: full model + + prior_FIM = np.zeros((4, 4)) + prior_FIM[0, 1] = 1e-3 + + experiment = FullReactorExperiment(data_ex, 10, 3) + + DoE_args = get_standard_args(experiment, fd_method, obj_used, flag_val) + DoE_args['prior_FIM'] = prior_FIM + + doe_obj = DesignOfExperiments(**DoE_args) + + with self.assertRaisesRegex( + ValueError, + "FIM provided is not symmetric using absolute tolerance {}".format( + _SMALL_TOLERANCE_SYMMETRY + ), + ): + doe_obj.create_doe_model() + def test_reactor_check_bad_jacobian_init_size(self): fd_method = "central" obj_used = "trace" diff --git a/pyomo/contrib/doe/tests/test_doe_solve.py b/pyomo/contrib/doe/tests/test_doe_solve.py index 1c4019bb3cd..83884be0b2c 100644 --- a/pyomo/contrib/doe/tests/test_doe_solve.py +++ b/pyomo/contrib/doe/tests/test_doe_solve.py @@ -230,6 +230,37 @@ def test_reactor_obj_cholesky_solve(self): # Make sure FIM and Q.T @ sigma_inv @ Q are close (alternate definition of FIM) self.assertTrue(np.all(np.isclose(FIM, Q.T @ sigma_inv @ Q))) + def test_reactor_obj_cholesky_solve_bad_prior(self): + + from pyomo.contrib.doe.doe import _SMALL_TOLERANCE_DEFINITENESS + + fd_method = "central" + obj_used = "determinant" + + experiment = FullReactorExperiment(data_ex, 10, 3) + + DoE_args = get_standard_args(experiment, fd_method, obj_used) + + # Specify a prior that is slightly negative definite + # Because it is less than the tolerance, it should be adjusted to be positive definite + # No error should be thrown + DoE_args['prior_FIM'] = -(_SMALL_TOLERANCE_DEFINITENESS / 100) * np.eye(4) + + doe_obj = DesignOfExperiments(**DoE_args) + + doe_obj.run_doe() + + self.assertEqual(doe_obj.results["Solver Status"], "ok") + + # assert that Q, F, and L are the same. + FIM, Q, L, sigma_inv = get_FIM_Q_L(doe_obj=doe_obj) + + # Since Cholesky is used, there is comparison for FIM and L.T @ L + self.assertTrue(np.all(np.isclose(FIM, L @ L.T))) + + # Make sure FIM and Q.T @ sigma_inv @ Q are close (alternate definition of FIM) + self.assertTrue(np.all(np.isclose(FIM, Q.T @ sigma_inv @ Q))) + # This test ensure that compute FIM runs without error using the # `sequential` option with central finite differences def test_compute_FIM_seq_centr(self):