From 39f4b095e1a65658855340d849e11f834ef79469 Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Wed, 1 Oct 2025 09:42:04 -0500 Subject: [PATCH 01/29] Add absolute value to condition number --- pyomo/contrib/doe/grey_box_utilities.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyomo/contrib/doe/grey_box_utilities.py b/pyomo/contrib/doe/grey_box_utilities.py index 8341cbf89df..299c884e138 100644 --- a/pyomo/contrib/doe/grey_box_utilities.py +++ b/pyomo/contrib/doe/grey_box_utilities.py @@ -185,7 +185,7 @@ def evaluate_outputs(self): obj_value = np.min(eig) elif self.objective_option == ObjectiveLib.condition_number: eig, _ = np.linalg.eig(M) - obj_value = np.max(eig) / np.min(eig) + obj_value = np.abs(np.max(eig) / np.min(eig)) else: ObjectiveLib(self.objective_option) From c42e093c70763fca80159c9bbcdd6953796ea943 Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Wed, 1 Oct 2025 09:45:38 -0500 Subject: [PATCH 02/29] Add log to condition number, begin hessian comp --- pyomo/contrib/doe/grey_box_utilities.py | 69 ++++++++++++++++++++++++- 1 file changed, 67 insertions(+), 2 deletions(-) diff --git a/pyomo/contrib/doe/grey_box_utilities.py b/pyomo/contrib/doe/grey_box_utilities.py index 299c884e138..afb684b8d18 100644 --- a/pyomo/contrib/doe/grey_box_utilities.py +++ b/pyomo/contrib/doe/grey_box_utilities.py @@ -185,7 +185,7 @@ def evaluate_outputs(self): obj_value = np.min(eig) elif self.objective_option == ObjectiveLib.condition_number: eig, _ = np.linalg.eig(M) - obj_value = np.abs(np.max(eig) / np.min(eig)) + obj_value = np.log(np.abs(np.max(eig) / np.min(eig))) else: ObjectiveLib(self.objective_option) @@ -519,7 +519,72 @@ def evaluate_hessian_outputs(self, FIM=None): multiplier = ((i != j) + (k != l)) + ((i != j) and (k != l)) * (i != k) hess_vals[-1] += multiplier * curr_hess_val elif self.objective_option == ObjectiveLib.condition_number: - pass + # Hessian for log condition number has 4 terms + # + # Grab eigenvalues and eigenvectors + # Also need the min location + all_eig_vals, all_eig_vecs = np.linalg.eig(M) + min_eig_loc = np.argmin(all_eig_vals) + + # Grabbing min eigenvalue and corresponding + # eigenvector + min_eig = all_eig_vals[min_eig_loc] + min_eig_vec = np.array([all_eig_vecs[:, min_eig_loc]]) + + for current_differential in input_differentials_2D: + # Row, Col and i, j, k, l values are + # obtained identically as in the trace + # for loop above. + d1, d2 = current_differential + row = self.input_names().index(d2) + col = self.input_names().index(d1) + + i = self._param_names.index(d1[0]) + j = self._param_names.index(d1[1]) + k = self._param_names.index(d2[0]) + l = self._param_names.index(d2[1]) + + # For lop to iterate over all + # eigenvalues/vectors + curr_hess_val = 0 + for curr_eig in range(len(all_eig_vals)): + # Skip if we are at the minimum + # eigenvalue. Denominator is + # zero. + if curr_eig == min_eig_loc: + continue + + # Formula derived in Pyomo.DoE Paper + curr_hess_val += ( + 1 + * ( + min_eig_vec[0, i] + * all_eig_vecs[j, curr_eig] + * min_eig_vec[0, l] + * all_eig_vecs[k, curr_eig] + ) + / (min_eig - all_eig_vals[curr_eig]) + ) + curr_hess_val += ( + 1 + * ( + min_eig_vec[0, k] + * all_eig_vecs[i, curr_eig] + * min_eig_vec[0, j] + * all_eig_vecs[l, curr_eig] + ) + / (min_eig - all_eig_vals[curr_eig]) + ) + hess_vals.append(curr_hess_val) + hess_rows.append(row) + hess_cols.append(col) + + # See explanation above in trace hessian code + # for more brief information. Also, you may + # consult the documentation for a detailed + # description. + multiplier = ((i != j) + (k != l)) + ((i != j) and (k != l)) * (i != k) + hess_vals[-1] += multiplier * curr_hess_val else: ObjectiveLib(self.objective_option) From 42e32dea0a0a9ee3ace870fdacf07bf947212b9a Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Wed, 1 Oct 2025 10:07:34 -0500 Subject: [PATCH 03/29] Added log-condition number hessian --- pyomo/contrib/doe/grey_box_utilities.py | 96 ++++++++++++++++++++++--- 1 file changed, 88 insertions(+), 8 deletions(-) diff --git a/pyomo/contrib/doe/grey_box_utilities.py b/pyomo/contrib/doe/grey_box_utilities.py index afb684b8d18..af14451c25e 100644 --- a/pyomo/contrib/doe/grey_box_utilities.py +++ b/pyomo/contrib/doe/grey_box_utilities.py @@ -477,7 +477,7 @@ def evaluate_hessian_outputs(self, FIM=None): k = self._param_names.index(d2[0]) l = self._param_names.index(d2[1]) - # For lop to iterate over all + # For loop to iterate over all # eigenvalues/vectors curr_hess_val = 0 for curr_eig in range(len(all_eig_vals)): @@ -519,18 +519,30 @@ def evaluate_hessian_outputs(self, FIM=None): multiplier = ((i != j) + (k != l)) + ((i != j) and (k != l)) * (i != k) hess_vals[-1] += multiplier * curr_hess_val elif self.objective_option == ObjectiveLib.condition_number: - # Hessian for log condition number has 4 terms + # Hessian for log condition number has 4 + # terms. Two are multiples of the second + # derivative of the maximum and minimum + # eigenvalues. The other two are tensor products + # of the first derivative of the maximum + # eigenvalue with itself, and minimum + # eigenvalue with itself. # # Grab eigenvalues and eigenvectors - # Also need the min location + # Also need the max and min locations all_eig_vals, all_eig_vecs = np.linalg.eig(M) min_eig_loc = np.argmin(all_eig_vals) + max_eig_loc = np.argmax(all_eig_vals) # Grabbing min eigenvalue and corresponding # eigenvector min_eig = all_eig_vals[min_eig_loc] min_eig_vec = np.array([all_eig_vecs[:, min_eig_loc]]) + # Grabbing max eigenvalue and corresponding + # eigenvector + max_eig = all_eig_vals[max_eig_loc] + max_eig_vec = np.array([all_eig_vecs[:, max_eig_loc]]) + for current_differential in input_differentials_2D: # Row, Col and i, j, k, l values are # obtained identically as in the trace @@ -544,9 +556,45 @@ def evaluate_hessian_outputs(self, FIM=None): k = self._param_names.index(d2[0]) l = self._param_names.index(d2[1]) - # For lop to iterate over all - # eigenvalues/vectors - curr_hess_val = 0 + # For loop to iterate over all + # eigenvalues/vectors for first + # term (second derivative of + # maximum eigenvalue) + log_cond_term_1 = 0 + for curr_eig in range(len(all_eig_vals)): + # Skip if we are at the minimum + # eigenvalue. Denominator is + # zero. + if curr_eig == max_eig_loc: + continue + + # Formula derived in Pyomo.DoE Paper + log_cond_term_1 += ( + 1 + * ( + max_eig_vec[0, i] + * all_eig_vecs[j, curr_eig] + * max_eig_vec[0, l] + * all_eig_vecs[k, curr_eig] + ) + / (max_eig - all_eig_vals[curr_eig]) + ) + log_cond_term_1 += ( + 1 + * ( + max_eig_vec[0, k] + * all_eig_vecs[i, curr_eig] + * max_eig_vec[0, j] + * all_eig_vecs[l, curr_eig] + ) + / (max_eig - all_eig_vals[curr_eig]) + ) + + # For loop to iterate over all + # eigenvalues/vectors for third + # term (second derivative of + # minimum eigenvalue) + log_cond_term_3 = 0 for curr_eig in range(len(all_eig_vals)): # Skip if we are at the minimum # eigenvalue. Denominator is @@ -555,7 +603,7 @@ def evaluate_hessian_outputs(self, FIM=None): continue # Formula derived in Pyomo.DoE Paper - curr_hess_val += ( + log_cond_term_3 += ( 1 * ( min_eig_vec[0, i] @@ -565,7 +613,7 @@ def evaluate_hessian_outputs(self, FIM=None): ) / (min_eig - all_eig_vals[curr_eig]) ) - curr_hess_val += ( + log_cond_term_3 += ( 1 * ( min_eig_vec[0, k] @@ -575,6 +623,38 @@ def evaluate_hessian_outputs(self, FIM=None): ) / (min_eig - all_eig_vals[curr_eig]) ) + + # Computing each term of the hessian formula + # Second derivative of max eigenvalue term + log_cond_term_1 = 1 / max_eig * log_cond_term_1 + + # First derivative of max eigenvalue term + log_cond_term_2 = ( + 1 + / (max_eig**2) + * (max_eig_vec[0, l] * max_eig_vec[0, k]) + * (max_eig_vec[0, j] * max_eig_vec[0, i]) + ) + + # Second derivative of min eigenvalue term + log_cond_term_3 = 1 / min_eig * log_cond_term_3 + + # First derivative of min eigenvalue term + log_cond_term_4 = ( + 1 + / (min_eig**2) + * (min_eig_vec[0, l] * min_eig_vec[0, k]) + * (min_eig_vec[0, j] * min_eig_vec[0, i]) + ) + + # Combining all the components + curr_hess_val = ( + log_cond_term_1 + - log_cond_term_2 + - log_cond_term_3 + + log_cond_term_4 + ) + hess_vals.append(curr_hess_val) hess_rows.append(row) hess_cols.append(col) From cbd7324368b472fa48da5aa40d235dc012b91c26 Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Wed, 1 Oct 2025 10:08:20 -0500 Subject: [PATCH 04/29] ran black --- pyomo/contrib/doe/grey_box_utilities.py | 64 ++++++++++++------------- 1 file changed, 32 insertions(+), 32 deletions(-) diff --git a/pyomo/contrib/doe/grey_box_utilities.py b/pyomo/contrib/doe/grey_box_utilities.py index af14451c25e..e10eb39f449 100644 --- a/pyomo/contrib/doe/grey_box_utilities.py +++ b/pyomo/contrib/doe/grey_box_utilities.py @@ -570,24 +570,24 @@ def evaluate_hessian_outputs(self, FIM=None): # Formula derived in Pyomo.DoE Paper log_cond_term_1 += ( - 1 - * ( - max_eig_vec[0, i] - * all_eig_vecs[j, curr_eig] - * max_eig_vec[0, l] - * all_eig_vecs[k, curr_eig] - ) - / (max_eig - all_eig_vals[curr_eig]) + 1 + * ( + max_eig_vec[0, i] + * all_eig_vecs[j, curr_eig] + * max_eig_vec[0, l] + * all_eig_vecs[k, curr_eig] + ) + / (max_eig - all_eig_vals[curr_eig]) ) log_cond_term_1 += ( - 1 - * ( - max_eig_vec[0, k] - * all_eig_vecs[i, curr_eig] - * max_eig_vec[0, j] - * all_eig_vecs[l, curr_eig] - ) - / (max_eig - all_eig_vals[curr_eig]) + 1 + * ( + max_eig_vec[0, k] + * all_eig_vecs[i, curr_eig] + * max_eig_vec[0, j] + * all_eig_vecs[l, curr_eig] + ) + / (max_eig - all_eig_vals[curr_eig]) ) # For loop to iterate over all @@ -604,24 +604,24 @@ def evaluate_hessian_outputs(self, FIM=None): # Formula derived in Pyomo.DoE Paper log_cond_term_3 += ( - 1 - * ( - min_eig_vec[0, i] - * all_eig_vecs[j, curr_eig] - * min_eig_vec[0, l] - * all_eig_vecs[k, curr_eig] - ) - / (min_eig - all_eig_vals[curr_eig]) + 1 + * ( + min_eig_vec[0, i] + * all_eig_vecs[j, curr_eig] + * min_eig_vec[0, l] + * all_eig_vecs[k, curr_eig] + ) + / (min_eig - all_eig_vals[curr_eig]) ) log_cond_term_3 += ( - 1 - * ( - min_eig_vec[0, k] - * all_eig_vecs[i, curr_eig] - * min_eig_vec[0, j] - * all_eig_vecs[l, curr_eig] - ) - / (min_eig - all_eig_vals[curr_eig]) + 1 + * ( + min_eig_vec[0, k] + * all_eig_vecs[i, curr_eig] + * min_eig_vec[0, j] + * all_eig_vecs[l, curr_eig] + ) + / (min_eig - all_eig_vals[curr_eig]) ) # Computing each term of the hessian formula From 9f00f0406c7cdf24fb37378fd25910eaa209bb1a Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Wed, 1 Oct 2025 10:14:40 -0500 Subject: [PATCH 05/29] Compute log condition number, and use the eigenvalue method --- pyomo/contrib/doe/tests/test_greybox.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyomo/contrib/doe/tests/test_greybox.py b/pyomo/contrib/doe/tests/test_greybox.py index 7a0d13fe465..de8fa80701a 100644 --- a/pyomo/contrib/doe/tests/test_greybox.py +++ b/pyomo/contrib/doe/tests/test_greybox.py @@ -506,7 +506,8 @@ def test_outputs_ME_opt(self): grey_box_ME_opt = grey_box_object.evaluate_outputs() - ME_opt = np.linalg.cond(testing_matrix) + vals, vecs = np.linalg.eig(testing_matrix) + ME_opt = np.log(np.abs(np.max(vals) / np.min(vals))) self.assertTrue(np.isclose(grey_box_ME_opt, ME_opt)) From 9f04e500a08f8e3327838a0cb600db8aaff651ee Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Wed, 1 Oct 2025 10:19:05 -0500 Subject: [PATCH 06/29] Update another test with new log-cond formalism --- pyomo/contrib/doe/tests/test_greybox.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyomo/contrib/doe/tests/test_greybox.py b/pyomo/contrib/doe/tests/test_greybox.py index de8fa80701a..9a340e1d6c0 100644 --- a/pyomo/contrib/doe/tests/test_greybox.py +++ b/pyomo/contrib/doe/tests/test_greybox.py @@ -884,7 +884,8 @@ def test_ME_opt_greybox_build(self): # Check output and value # FIM Initial will be the prior FIM # added with the identity matrix. - ME_opt_val = np.linalg.cond(testing_matrix + np.eye(4)) + vals, vecs = np.linalg.eig(testing_matrix + np.eye(4)) + ME_opt_val = np.log(np.abs(np.max(vals) / np.min(vals))) try: ME_opt_val_gb = doe_obj.model.obj_cons.egb_fim_block.outputs["ME-opt"].value From 621740cabbeccd5e46ef10d87cba6647b2277b18 Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Wed, 1 Oct 2025 10:20:38 -0500 Subject: [PATCH 07/29] Update last test with log-cond formalism --- pyomo/contrib/doe/tests/test_greybox.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyomo/contrib/doe/tests/test_greybox.py b/pyomo/contrib/doe/tests/test_greybox.py index 9a340e1d6c0..d91217d7190 100644 --- a/pyomo/contrib/doe/tests/test_greybox.py +++ b/pyomo/contrib/doe/tests/test_greybox.py @@ -1122,8 +1122,8 @@ def test_solve_ME_optimality_condition_number(self): # Here, the objective value is # condition number of the FIM optimal_experimental_designs = [ - np.array([0.943, 13.524]), - np.array([10.00, 27.675]), + np.array([0.943, np.log(13.524)]), + np.array([10.00, np.log(27.675)]), ] objective_option = "condition_number" doe_object, grey_box_object = make_greybox_and_doe_objects_rooney_biegler( From 06246fe3ea343329c9f1ae01798685fdcb6e9de4 Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Wed, 1 Oct 2025 10:43:02 -0500 Subject: [PATCH 08/29] Fixing initialization for condition number to be log --- pyomo/contrib/doe/doe.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py index 5d16e327561..c060afbc582 100644 --- a/pyomo/contrib/doe/doe.py +++ b/pyomo/contrib/doe/doe.py @@ -370,7 +370,8 @@ def run_doe(self, model=None, results_file=None): eig, _ = np.linalg.eig(np.array(self.get_FIM())) model.obj_cons.egb_fim_block.outputs["E-opt"].set_value(np.min(eig)) elif self.objective_option == ObjectiveLib.condition_number: - cond_number = np.linalg.cond(np.array(self.get_FIM())) + eig, _ = np.linalg.eig(np.array(self.get_FIM())) + cond_number = np.log(np.abs(np.max(eig)/np.min(eig))) model.obj_cons.egb_fim_block.outputs["ME-opt"].set_value(cond_number) # If the model has L, initialize it with the solved FIM From 15abd5d887fbd8a30f8cf445d96a98486969f3dc Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Wed, 1 Oct 2025 10:43:22 -0500 Subject: [PATCH 09/29] Fix jacobian formula for log condition number --- pyomo/contrib/doe/grey_box_utilities.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/pyomo/contrib/doe/grey_box_utilities.py b/pyomo/contrib/doe/grey_box_utilities.py index e10eb39f449..904e2e3527b 100644 --- a/pyomo/contrib/doe/grey_box_utilities.py +++ b/pyomo/contrib/doe/grey_box_utilities.py @@ -303,9 +303,7 @@ def evaluate_jacobian_outputs(self): # Combining the expression jac_M = ( - 1 - / (min_eig + np.sign(min_eig) * min_eig_epsilon) - * (max_eig_term - safe_cond_number * min_eig_term) + 1 / max_eig * max_eig_term - 1 / min_eig * min_eig_term ) else: ObjectiveLib(self.objective_option) From 459a5d9aaff7c30d15612efce8bc1ab96e27be67 Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Wed, 1 Oct 2025 10:43:47 -0500 Subject: [PATCH 10/29] Fix finite difference calculation to be log, add hessian test for log cond --- pyomo/contrib/doe/tests/test_greybox.py | 40 +++++++++++++++++++++---- 1 file changed, 34 insertions(+), 6 deletions(-) diff --git a/pyomo/contrib/doe/tests/test_greybox.py b/pyomo/contrib/doe/tests/test_greybox.py index d91217d7190..0b49daebf6b 100644 --- a/pyomo/contrib/doe/tests/test_greybox.py +++ b/pyomo/contrib/doe/tests/test_greybox.py @@ -88,7 +88,8 @@ def get_numerical_derivative(grey_box_object=None): vals_init, vecs_init = np.linalg.eig(current_FIM) unperturbed_value = np.min(vals_init) elif grey_box_object.objective_option == ObjectiveLib.condition_number: - unperturbed_value = np.linalg.cond(current_FIM) + vals_init, vecs_init = np.linalg.eig(current_FIM) + unperturbed_value = np.log(np.abs(np.max(vals_init) / np.min(vals_init))) # Calculate the numerical derivative, using forward difference numerical_derivative = np.zeros_like(current_FIM) @@ -109,7 +110,8 @@ def get_numerical_derivative(grey_box_object=None): vals, vecs = np.linalg.eig(FIM_perturbed) new_value_ij = np.min(vals) elif grey_box_object.objective_option == ObjectiveLib.condition_number: - new_value_ij = np.linalg.cond(FIM_perturbed) + vals, vecs = np.linalg.eig(FIM_perturbed) + new_value_ij = np.log(np.abs(np.max(vals) / np.min(vals))) # Calculate the derivative value from forward difference diff = (new_value_ij - unperturbed_value) / _FD_EPSILON_FIRST @@ -186,10 +188,14 @@ def get_numerical_second_derivative(grey_box_object=None, return_reduced=True): grey_box_object.objective_option == ObjectiveLib.condition_number ): - new_values[0] = np.linalg.cond(FIM_perturbed_1) - new_values[1] = np.linalg.cond(FIM_perturbed_2) - new_values[2] = np.linalg.cond(FIM_perturbed_3) - new_values[3] = np.linalg.cond(FIM_perturbed_4) + vals, vecs = np.linalg.eig(FIM_perturbed_1) + new_values[0] = np.log(np.abs(np.max(vals) / np.min(vals))) + vals, vecs = np.linalg.eig(FIM_perturbed_2) + new_values[1] = np.log(np.abs(np.max(vals) / np.min(vals))) + vals, vecs = np.linalg.eig(FIM_perturbed_3) + new_values[2] = np.log(np.abs(np.max(vals) / np.min(vals))) + vals, vecs = np.linalg.eig(FIM_perturbed_4) + new_values[3] = np.log(np.abs(np.max(vals) / np.min(vals))) # Calculate the derivative value from second order difference formula diff = ( @@ -675,6 +681,28 @@ def test_hessian_E_opt(self): # assert that each component is close self.assertTrue(np.all(np.isclose(hess_gb, hess_FD, rtol=1e-4, atol=1e-4))) + def test_hessian_ME_opt(self): + objective_option = "condition_number" + doe_obj, grey_box_object = make_greybox_and_doe_objects( + objective_option=objective_option + ) + + # Set input values to the random testing matrix + grey_box_object.set_input_values(testing_matrix[masking_matrix > 0]) + + # Grab the Jacobian values + hess_vals_from_gb = grey_box_object.evaluate_hessian_outputs().toarray() + + # Recover the Jacobian in Matrix Form + hess_gb = hess_vals_from_gb + hess_gb += hess_gb.transpose() - np.diag(np.diag(hess_gb)) + + # Get numerical derivative matrix + hess_FD = get_numerical_second_derivative(grey_box_object) + + # assert that each component is close + self.assertTrue(np.all(np.isclose(hess_gb, hess_FD, rtol=1e-4, atol=1e-4))) + def test_equality_constraint_names(self): objective_option = "condition_number" doe_obj, grey_box_object = make_greybox_and_doe_objects( From 941e49a73378f9507b4a8d960915e2f68235279a Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Wed, 1 Oct 2025 10:44:18 -0500 Subject: [PATCH 11/29] ran black --- pyomo/contrib/doe/doe.py | 2 +- pyomo/contrib/doe/grey_box_utilities.py | 4 +--- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py index c060afbc582..ece387dc033 100644 --- a/pyomo/contrib/doe/doe.py +++ b/pyomo/contrib/doe/doe.py @@ -371,7 +371,7 @@ def run_doe(self, model=None, results_file=None): model.obj_cons.egb_fim_block.outputs["E-opt"].set_value(np.min(eig)) elif self.objective_option == ObjectiveLib.condition_number: eig, _ = np.linalg.eig(np.array(self.get_FIM())) - cond_number = np.log(np.abs(np.max(eig)/np.min(eig))) + cond_number = np.log(np.abs(np.max(eig) / np.min(eig))) model.obj_cons.egb_fim_block.outputs["ME-opt"].set_value(cond_number) # If the model has L, initialize it with the solved FIM diff --git a/pyomo/contrib/doe/grey_box_utilities.py b/pyomo/contrib/doe/grey_box_utilities.py index 904e2e3527b..291ae6a3972 100644 --- a/pyomo/contrib/doe/grey_box_utilities.py +++ b/pyomo/contrib/doe/grey_box_utilities.py @@ -302,9 +302,7 @@ def evaluate_jacobian_outputs(self): safe_cond_number = max_eig / (min_eig + np.sign(min_eig) * min_eig_epsilon) # Combining the expression - jac_M = ( - 1 / max_eig * max_eig_term - 1 / min_eig * min_eig_term - ) + jac_M = 1 / max_eig * max_eig_term - 1 / min_eig * min_eig_term else: ObjectiveLib(self.objective_option) From 524a42e098a02ac3e12253ca9647a269038b92f1 Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Wed, 1 Oct 2025 15:55:58 -0500 Subject: [PATCH 12/29] Update pyomo/contrib/doe/grey_box_utilities.py Fixed typo in comment Co-authored-by: Bethany Nicholson --- pyomo/contrib/doe/grey_box_utilities.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyomo/contrib/doe/grey_box_utilities.py b/pyomo/contrib/doe/grey_box_utilities.py index e10eb39f449..d8708f488e8 100644 --- a/pyomo/contrib/doe/grey_box_utilities.py +++ b/pyomo/contrib/doe/grey_box_utilities.py @@ -562,7 +562,7 @@ def evaluate_hessian_outputs(self, FIM=None): # maximum eigenvalue) log_cond_term_1 = 0 for curr_eig in range(len(all_eig_vals)): - # Skip if we are at the minimum + # Skip if we are at the maximum # eigenvalue. Denominator is # zero. if curr_eig == max_eig_loc: From 4429592306f5c76c8a9ffa27b352cc1e72f45984 Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Thu, 2 Oct 2025 12:11:38 -0500 Subject: [PATCH 13/29] Changing A-optimality to test new hessian --- pyomo/contrib/doe/grey_box_utilities.py | 96 ++++++++++++++++++------- 1 file changed, 70 insertions(+), 26 deletions(-) diff --git a/pyomo/contrib/doe/grey_box_utilities.py b/pyomo/contrib/doe/grey_box_utilities.py index 9fe388c618c..dee865144f7 100644 --- a/pyomo/contrib/doe/grey_box_utilities.py +++ b/pyomo/contrib/doe/grey_box_utilities.py @@ -123,6 +123,25 @@ def _get_FIM(self): return current_FIM + def _reorder_pairs(self, i, j, k, l): + # Reorders the pairs (i, j) and + # (k, l) for considering only + # the symmetric portion of the FIM + # while calculating the Hessian + + # If the pairs ((i, j), (k, l)) are not + # in increasing order, we reorder + # the pairs. + if i > j: + if k > l: + return [j, i, l, k] + else: + return [j, i, k, l] + else: + if k > l: + return [i, j, l, k] + return [i, j, k, l] + def input_names(self): # Cartesian product gives us matrix indices flattened in row-first format # Can use itertools.combinations(self._param_names, 2) with added @@ -350,6 +369,12 @@ def evaluate_hessian_outputs(self, FIM=None): self._n_params, self._n_params ) + # Length of Hessian lists for the sparse + # matrix ar a function of number of parameters + hess_array_length = round( + (((self._n_params + 1) * self._n_params / 2) + 1) * (((self._n_params + 1) * self._n_params / 2)) / 2 + ) + # Hessian with correct size for using only the # lower (upper) triangle of the FIM hess_vals = [] @@ -365,6 +390,12 @@ def evaluate_hessian_outputs(self, FIM=None): from pyomo.contrib.doe import ObjectiveLib if self.objective_option == ObjectiveLib.trace: + input_differentials_2D = itertools.product( + self.input_names(), repeat=2 + ) + hess_vals = [0, ] * hess_array_length + hess_rows = [0, ] * hess_array_length + hess_cols = [0, ] * hess_array_length # Grab Inverse Minv = np.linalg.pinv(M) @@ -372,14 +403,7 @@ def evaluate_hessian_outputs(self, FIM=None): Minv_sq = Minv @ Minv for current_differential in input_differentials_2D: - # Row will be the location of the - # first ordered pair (d1) in input names - # - # Col will be the location of the - # second ordered pair (d2) in input names d1, d2 = current_differential - row = self.input_names().index(d2) - col = self.input_names().index(d1) # Grabbing the ordered quadruple (i, j, k, l) # `location` here refers to the index in the @@ -396,29 +420,46 @@ def evaluate_hessian_outputs(self, FIM=None): # New Formula (tested with finite differencing) # Will be cited from the Pyomo.DoE 2.0 paper - hess_vals.append( - (Minv[i, l] * Minv_sq[k, j]) + (Minv_sq[i, l] * Minv[k, j]) - ) - hess_rows.append(row) - hess_cols.append(col) + hess_contribution = (Minv[i, l] * Minv_sq[k, j]) + (Minv_sq[i, l] * Minv[k, j]) + + # Since we are considering the full matrix in + # this loop, we need to point the contribution + # to the correct index for the symmetric FIM + # Hessian. + reordered_ijkl = self._reorder_pairs(i, j, k, l) + d1_symmetric = (self._param_names[reordered_ijkl[0]], self._param_names[reordered_ijkl[1]]) + d2_symmetric = (self._param_names[reordered_ijkl[2]], self._param_names[reordered_ijkl[3]]) + + # Identify what index of the symmetric FIM + # Hessian arrays need to be updated + row = self.input_names().index(d1_symmetric) + col = self.input_names().index(d2_symmetric) + flattened_row_col_index = (row + 1) * row // 2 + col # Hessian needs to be handled carefully because of # the ``missing`` components when only passing - # a triangular version of the FIM. + # a symmetric version of the FIM. # - # The general format is that fully diagonal elements - # (i == j AND k == l) require no additional counting - # of the computed term. Off-diagonal elements - # (i != j OR k != l) add one instance. Having two - # off-diagonal elements (i != j AND k != l) will add - # two instances. Finally, if the off-diagonal elements - # are not the same ((i != j AND k != l) AND (i != k)) - # there is an additional element. A more detailed - # explanation is included in the documentation. - multiplier = ((i != j) + (k != l)) + ((i != j) and (k != l)) * (i != k) - hess_vals[-1] += multiplier * ( - (Minv[i, l] * Minv_sq[k, j]) + (Minv_sq[i, l] * Minv[k, j]) - ) + # When we reordered (i, j, k, l), we are correctly + # pointing to which index needs to be contributed to. + # However, when an element that is not included + # is being put into a diagonal element of the + # symmetric FIM hessian from the full FIM hessian, + # it needs to be counted twice. This only occurs + # when (i != j) and (k != l) and (i, j) and (k, l) + # are the conjugate of one another: + # (i == l) and (j == k). + # + # Otherwise, we only add the element once. + # Standard addition + hess_vals[flattened_row_col_index] += hess_contribution + + # Duplicate check and addition + if ((i != j) and (k != l)) and ((i == l) and ( j == k)): + hess_vals[flattened_row_col_index] += hess_contribution + + hess_rows[flattened_row_col_index] = row + hess_cols[flattened_row_col_index] = col elif self.objective_option == ObjectiveLib.determinant: # Grab inverse @@ -664,6 +705,9 @@ def evaluate_hessian_outputs(self, FIM=None): else: ObjectiveLib(self.objective_option) + print(hess_vals) + print(hess_rows) + print(hess_cols) # Returns coo_matrix of the correct shape return scipy.sparse.coo_matrix( (np.asarray(hess_vals), (hess_rows, hess_cols)), From 8b7287e0d5e9b279643078b45973afeed2a124e2 Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Thu, 2 Oct 2025 12:28:41 -0500 Subject: [PATCH 14/29] Fixing the other optimality criteria hessians --- pyomo/contrib/doe/grey_box_utilities.py | 152 +++++++++++++++--------- 1 file changed, 99 insertions(+), 53 deletions(-) diff --git a/pyomo/contrib/doe/grey_box_utilities.py b/pyomo/contrib/doe/grey_box_utilities.py index dee865144f7..a495eebdf77 100644 --- a/pyomo/contrib/doe/grey_box_utilities.py +++ b/pyomo/contrib/doe/grey_box_utilities.py @@ -377,25 +377,19 @@ def evaluate_hessian_outputs(self, FIM=None): # Hessian with correct size for using only the # lower (upper) triangle of the FIM - hess_vals = [] - hess_rows = [] - hess_cols = [] + hess_vals = [0, ] * hess_array_length + hess_rows = [0, ] * hess_array_length + hess_cols = [0, ] * hess_array_length # Need to iterate over the unique # differentials - input_differentials_2D = itertools.combinations_with_replacement( - self.input_names(), 2 + input_differentials_2D = itertools.product( + self.input_names(), repeat=2 ) from pyomo.contrib.doe import ObjectiveLib if self.objective_option == ObjectiveLib.trace: - input_differentials_2D = itertools.product( - self.input_names(), repeat=2 - ) - hess_vals = [0, ] * hess_array_length - hess_rows = [0, ] * hess_array_length - hess_cols = [0, ] * hess_array_length # Grab Inverse Minv = np.linalg.pinv(M) @@ -470,8 +464,6 @@ def evaluate_hessian_outputs(self, FIM=None): # obtained identically as in the trace # for loop above. d1, d2 = current_differential - row = self.input_names().index(d2) - col = self.input_names().index(d1) i = self._param_names.index(d1[0]) j = self._param_names.index(d1[1]) @@ -480,16 +472,36 @@ def evaluate_hessian_outputs(self, FIM=None): # New Formula (tested with finite differencing) # Will be cited from the Pyomo.DoE 2.0 paper - hess_vals.append(-(Minv[i, l] * Minv[k, j])) - hess_rows.append(row) - hess_cols.append(col) - - # See explanation above in trace hessian code - # for more brief information. Also, you may - # consult the documentation for a detailed - # description. - multiplier = ((i != j) + (k != l)) + ((i != j) and (k != l)) * (i != k) - hess_vals[-1] += -multiplier * (Minv[i, l] * Minv[k, j]) + hess_contribution = -(Minv[i, l] * Minv[k, j]) + + # Since we are considering the full matrix in + # this loop, we need to point the contribution + # to the correct index for the symmetric FIM + # Hessian. + reordered_ijkl = self._reorder_pairs(i, j, k, l) + d1_symmetric = (self._param_names[reordered_ijkl[0]], self._param_names[reordered_ijkl[1]]) + d2_symmetric = (self._param_names[reordered_ijkl[2]], self._param_names[reordered_ijkl[3]]) + + # Identify what index of the symmetric FIM + # Hessian arrays need to be updated + row = self.input_names().index(d1_symmetric) + col = self.input_names().index(d2_symmetric) + flattened_row_col_index = (row + 1) * row // 2 + col + + # Hessian needs to be handled carefully because of + # the ``missing`` components when only passing + # a symmetric version of the FIM. For a more + # detailed explanation, please see the trace + # for loop above + hess_vals[flattened_row_col_index] += hess_contribution + + # Duplicate check and addition + if ((i != j) and (k != l)) and ((i == l) and (j == k)): + hess_vals[flattened_row_col_index] += hess_contribution + + hess_rows[flattened_row_col_index] = row + hess_cols[flattened_row_col_index] = col + elif self.objective_option == ObjectiveLib.minimum_eigenvalue: # Grab eigenvalues and eigenvectors # Also need the min location @@ -506,8 +518,6 @@ def evaluate_hessian_outputs(self, FIM=None): # obtained identically as in the trace # for loop above. d1, d2 = current_differential - row = self.input_names().index(d2) - col = self.input_names().index(d1) i = self._param_names.index(d1[0]) j = self._param_names.index(d1[1]) @@ -516,7 +526,7 @@ def evaluate_hessian_outputs(self, FIM=None): # For loop to iterate over all # eigenvalues/vectors - curr_hess_val = 0 + hess_contribution = 0 for curr_eig in range(len(all_eig_vals)): # Skip if we are at the minimum # eigenvalue. Denominator is @@ -525,7 +535,7 @@ def evaluate_hessian_outputs(self, FIM=None): continue # Formula derived in Pyomo.DoE Paper - curr_hess_val += ( + hess_contribution += ( 1 * ( min_eig_vec[0, i] @@ -535,7 +545,7 @@ def evaluate_hessian_outputs(self, FIM=None): ) / (min_eig - all_eig_vals[curr_eig]) ) - curr_hess_val += ( + hess_contribution += ( 1 * ( min_eig_vec[0, k] @@ -545,16 +555,34 @@ def evaluate_hessian_outputs(self, FIM=None): ) / (min_eig - all_eig_vals[curr_eig]) ) - hess_vals.append(curr_hess_val) - hess_rows.append(row) - hess_cols.append(col) - - # See explanation above in trace hessian code - # for more brief information. Also, you may - # consult the documentation for a detailed - # description. - multiplier = ((i != j) + (k != l)) + ((i != j) and (k != l)) * (i != k) - hess_vals[-1] += multiplier * curr_hess_val + + # Since we are considering the full matrix in + # this loop, we need to point the contribution + # to the correct index for the symmetric FIM + # Hessian. + reordered_ijkl = self._reorder_pairs(i, j, k, l) + d1_symmetric = (self._param_names[reordered_ijkl[0]], self._param_names[reordered_ijkl[1]]) + d2_symmetric = (self._param_names[reordered_ijkl[2]], self._param_names[reordered_ijkl[3]]) + + # Identify what index of the symmetric FIM + # Hessian arrays need to be updated + row = self.input_names().index(d1_symmetric) + col = self.input_names().index(d2_symmetric) + flattened_row_col_index = (row + 1) * row // 2 + col + + # Hessian needs to be handled carefully because of + # the ``missing`` components when only passing + # a symmetric version of the FIM. See trace for loop + # for more detailed explanation + hess_vals[flattened_row_col_index] += hess_contribution + + # Duplicate check and addition + if ((i != j) and (k != l)) and ((i == l) and (j == k)): + hess_vals[flattened_row_col_index] += hess_contribution + + hess_rows[flattened_row_col_index] = row + hess_cols[flattened_row_col_index] = col + elif self.objective_option == ObjectiveLib.condition_number: # Hessian for log condition number has 4 # terms. Two are multiples of the second @@ -566,6 +594,13 @@ def evaluate_hessian_outputs(self, FIM=None): # # Grab eigenvalues and eigenvectors # Also need the max and min locations + input_differentials_2D = itertools.product( + self.input_names(), repeat=2 + ) + hess_vals = [0, ] * hess_array_length + hess_rows = [0, ] * hess_array_length + hess_cols = [0, ] * hess_array_length + all_eig_vals, all_eig_vecs = np.linalg.eig(M) min_eig_loc = np.argmin(all_eig_vals) max_eig_loc = np.argmax(all_eig_vals) @@ -585,8 +620,6 @@ def evaluate_hessian_outputs(self, FIM=None): # obtained identically as in the trace # for loop above. d1, d2 = current_differential - row = self.input_names().index(d2) - col = self.input_names().index(d1) i = self._param_names.index(d1[0]) j = self._param_names.index(d1[1]) @@ -685,29 +718,42 @@ def evaluate_hessian_outputs(self, FIM=None): ) # Combining all the components - curr_hess_val = ( + hess_contribution = ( log_cond_term_1 - log_cond_term_2 - log_cond_term_3 + log_cond_term_4 ) - hess_vals.append(curr_hess_val) - hess_rows.append(row) - hess_cols.append(col) + # Since we are considering the full matrix in + # this loop, we need to point the contribution + # to the correct index for the symmetric FIM + # Hessian. + reordered_ijkl = self._reorder_pairs(i, j, k, l) + d1_symmetric = (self._param_names[reordered_ijkl[0]], self._param_names[reordered_ijkl[1]]) + d2_symmetric = (self._param_names[reordered_ijkl[2]], self._param_names[reordered_ijkl[3]]) + + # Identify what index of the symmetric FIM + # Hessian arrays need to be updated + row = self.input_names().index(d1_symmetric) + col = self.input_names().index(d2_symmetric) + flattened_row_col_index = (row + 1) * row // 2 + col + + # Hessian needs to be handled carefully because of + # the ``missing`` components when only passing + # a symmetric version of the FIM. See trace for loop + # for more detailed explanation + hess_vals[flattened_row_col_index] += hess_contribution - # See explanation above in trace hessian code - # for more brief information. Also, you may - # consult the documentation for a detailed - # description. - multiplier = ((i != j) + (k != l)) + ((i != j) and (k != l)) * (i != k) - hess_vals[-1] += multiplier * curr_hess_val + # Duplicate check and addition + if ((i != j) and (k != l)) and ((i == l) and (j == k)): + hess_vals[flattened_row_col_index] += hess_contribution + + hess_rows[flattened_row_col_index] = row + hess_cols[flattened_row_col_index] = col else: ObjectiveLib(self.objective_option) - print(hess_vals) - print(hess_rows) - print(hess_cols) # Returns coo_matrix of the correct shape return scipy.sparse.coo_matrix( (np.asarray(hess_vals), (hess_rows, hess_cols)), From a6aff8ea900a15f1494d7b73aa2df68bc604ba97 Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Thu, 2 Oct 2025 12:29:01 -0500 Subject: [PATCH 15/29] ran black --- pyomo/contrib/doe/grey_box_utilities.py | 70 +++++++++++++++++-------- 1 file changed, 47 insertions(+), 23 deletions(-) diff --git a/pyomo/contrib/doe/grey_box_utilities.py b/pyomo/contrib/doe/grey_box_utilities.py index a495eebdf77..9fa8ac9aeb0 100644 --- a/pyomo/contrib/doe/grey_box_utilities.py +++ b/pyomo/contrib/doe/grey_box_utilities.py @@ -372,20 +372,20 @@ def evaluate_hessian_outputs(self, FIM=None): # Length of Hessian lists for the sparse # matrix ar a function of number of parameters hess_array_length = round( - (((self._n_params + 1) * self._n_params / 2) + 1) * (((self._n_params + 1) * self._n_params / 2)) / 2 + (((self._n_params + 1) * self._n_params / 2) + 1) + * (((self._n_params + 1) * self._n_params / 2)) + / 2 ) # Hessian with correct size for using only the # lower (upper) triangle of the FIM - hess_vals = [0, ] * hess_array_length - hess_rows = [0, ] * hess_array_length - hess_cols = [0, ] * hess_array_length + hess_vals = [0] * hess_array_length + hess_rows = [0] * hess_array_length + hess_cols = [0] * hess_array_length # Need to iterate over the unique # differentials - input_differentials_2D = itertools.product( - self.input_names(), repeat=2 - ) + input_differentials_2D = itertools.product(self.input_names(), repeat=2) from pyomo.contrib.doe import ObjectiveLib @@ -414,15 +414,23 @@ def evaluate_hessian_outputs(self, FIM=None): # New Formula (tested with finite differencing) # Will be cited from the Pyomo.DoE 2.0 paper - hess_contribution = (Minv[i, l] * Minv_sq[k, j]) + (Minv_sq[i, l] * Minv[k, j]) + hess_contribution = (Minv[i, l] * Minv_sq[k, j]) + ( + Minv_sq[i, l] * Minv[k, j] + ) # Since we are considering the full matrix in # this loop, we need to point the contribution # to the correct index for the symmetric FIM # Hessian. reordered_ijkl = self._reorder_pairs(i, j, k, l) - d1_symmetric = (self._param_names[reordered_ijkl[0]], self._param_names[reordered_ijkl[1]]) - d2_symmetric = (self._param_names[reordered_ijkl[2]], self._param_names[reordered_ijkl[3]]) + d1_symmetric = ( + self._param_names[reordered_ijkl[0]], + self._param_names[reordered_ijkl[1]], + ) + d2_symmetric = ( + self._param_names[reordered_ijkl[2]], + self._param_names[reordered_ijkl[3]], + ) # Identify what index of the symmetric FIM # Hessian arrays need to be updated @@ -449,7 +457,7 @@ def evaluate_hessian_outputs(self, FIM=None): hess_vals[flattened_row_col_index] += hess_contribution # Duplicate check and addition - if ((i != j) and (k != l)) and ((i == l) and ( j == k)): + if ((i != j) and (k != l)) and ((i == l) and (j == k)): hess_vals[flattened_row_col_index] += hess_contribution hess_rows[flattened_row_col_index] = row @@ -479,8 +487,14 @@ def evaluate_hessian_outputs(self, FIM=None): # to the correct index for the symmetric FIM # Hessian. reordered_ijkl = self._reorder_pairs(i, j, k, l) - d1_symmetric = (self._param_names[reordered_ijkl[0]], self._param_names[reordered_ijkl[1]]) - d2_symmetric = (self._param_names[reordered_ijkl[2]], self._param_names[reordered_ijkl[3]]) + d1_symmetric = ( + self._param_names[reordered_ijkl[0]], + self._param_names[reordered_ijkl[1]], + ) + d2_symmetric = ( + self._param_names[reordered_ijkl[2]], + self._param_names[reordered_ijkl[3]], + ) # Identify what index of the symmetric FIM # Hessian arrays need to be updated @@ -561,8 +575,14 @@ def evaluate_hessian_outputs(self, FIM=None): # to the correct index for the symmetric FIM # Hessian. reordered_ijkl = self._reorder_pairs(i, j, k, l) - d1_symmetric = (self._param_names[reordered_ijkl[0]], self._param_names[reordered_ijkl[1]]) - d2_symmetric = (self._param_names[reordered_ijkl[2]], self._param_names[reordered_ijkl[3]]) + d1_symmetric = ( + self._param_names[reordered_ijkl[0]], + self._param_names[reordered_ijkl[1]], + ) + d2_symmetric = ( + self._param_names[reordered_ijkl[2]], + self._param_names[reordered_ijkl[3]], + ) # Identify what index of the symmetric FIM # Hessian arrays need to be updated @@ -594,12 +614,10 @@ def evaluate_hessian_outputs(self, FIM=None): # # Grab eigenvalues and eigenvectors # Also need the max and min locations - input_differentials_2D = itertools.product( - self.input_names(), repeat=2 - ) - hess_vals = [0, ] * hess_array_length - hess_rows = [0, ] * hess_array_length - hess_cols = [0, ] * hess_array_length + input_differentials_2D = itertools.product(self.input_names(), repeat=2) + hess_vals = [0] * hess_array_length + hess_rows = [0] * hess_array_length + hess_cols = [0] * hess_array_length all_eig_vals, all_eig_vecs = np.linalg.eig(M) min_eig_loc = np.argmin(all_eig_vals) @@ -730,8 +748,14 @@ def evaluate_hessian_outputs(self, FIM=None): # to the correct index for the symmetric FIM # Hessian. reordered_ijkl = self._reorder_pairs(i, j, k, l) - d1_symmetric = (self._param_names[reordered_ijkl[0]], self._param_names[reordered_ijkl[1]]) - d2_symmetric = (self._param_names[reordered_ijkl[2]], self._param_names[reordered_ijkl[3]]) + d1_symmetric = ( + self._param_names[reordered_ijkl[0]], + self._param_names[reordered_ijkl[1]], + ) + d2_symmetric = ( + self._param_names[reordered_ijkl[2]], + self._param_names[reordered_ijkl[3]], + ) # Identify what index of the symmetric FIM # Hessian arrays need to be updated From b3f3aa9f9bb5826746005baf71b7923ce73063b8 Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Thu, 2 Oct 2025 12:31:16 -0500 Subject: [PATCH 16/29] Cleaned comments up to be accurate in tests --- pyomo/contrib/doe/tests/test_greybox.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/pyomo/contrib/doe/tests/test_greybox.py b/pyomo/contrib/doe/tests/test_greybox.py index 0b49daebf6b..5bebabebc00 100644 --- a/pyomo/contrib/doe/tests/test_greybox.py +++ b/pyomo/contrib/doe/tests/test_greybox.py @@ -624,10 +624,10 @@ def test_hessian_A_opt(self): # Set input values to the random testing matrix grey_box_object.set_input_values(testing_matrix[masking_matrix > 0]) - # Grab the Jacobian values + # Grab the Hessian values hess_vals_from_gb = grey_box_object.evaluate_hessian_outputs().toarray() - # Recover the Jacobian in Matrix Form + # Recover the Hessian in Matrix Form hess_gb = hess_vals_from_gb hess_gb += hess_gb.transpose() - np.diag(np.diag(hess_gb)) @@ -646,10 +646,10 @@ def test_hessian_D_opt(self): # Set input values to the random testing matrix grey_box_object.set_input_values(testing_matrix[masking_matrix > 0]) - # Grab the Jacobian values + # Grab the Hessian values hess_vals_from_gb = grey_box_object.evaluate_hessian_outputs().toarray() - # Recover the Jacobian in Matrix Form + # Recover the Hessian in Matrix Form hess_gb = hess_vals_from_gb hess_gb += hess_gb.transpose() - np.diag(np.diag(hess_gb)) @@ -668,10 +668,10 @@ def test_hessian_E_opt(self): # Set input values to the random testing matrix grey_box_object.set_input_values(testing_matrix[masking_matrix > 0]) - # Grab the Jacobian values + # Grab the Hessian values hess_vals_from_gb = grey_box_object.evaluate_hessian_outputs().toarray() - # Recover the Jacobian in Matrix Form + # Recover the Hessian in Matrix Form hess_gb = hess_vals_from_gb hess_gb += hess_gb.transpose() - np.diag(np.diag(hess_gb)) @@ -690,10 +690,10 @@ def test_hessian_ME_opt(self): # Set input values to the random testing matrix grey_box_object.set_input_values(testing_matrix[masking_matrix > 0]) - # Grab the Jacobian values + # Grab the Hessian values hess_vals_from_gb = grey_box_object.evaluate_hessian_outputs().toarray() - # Recover the Jacobian in Matrix Form + # Recover the Hessian in Matrix Form hess_gb = hess_vals_from_gb hess_gb += hess_gb.transpose() - np.diag(np.diag(hess_gb)) From 074fb613c6b949cd1b2d8192e5cd46f0e84c91e9 Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Thu, 2 Oct 2025 16:46:12 -0500 Subject: [PATCH 17/29] Fix indexing issue --- pyomo/contrib/doe/grey_box_utilities.py | 28 ++++++++++++------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/pyomo/contrib/doe/grey_box_utilities.py b/pyomo/contrib/doe/grey_box_utilities.py index 9fa8ac9aeb0..b2f58130ad3 100644 --- a/pyomo/contrib/doe/grey_box_utilities.py +++ b/pyomo/contrib/doe/grey_box_utilities.py @@ -385,7 +385,8 @@ def evaluate_hessian_outputs(self, FIM=None): # Need to iterate over the unique # differentials - input_differentials_2D = itertools.product(self.input_names(), repeat=2) + full_input_names = itertools.product(self._param_names, repeat=2) + input_differentials_2D = itertools.combinations_with_replacement(full_input_names, 2) from pyomo.contrib.doe import ObjectiveLib @@ -397,6 +398,7 @@ def evaluate_hessian_outputs(self, FIM=None): Minv_sq = Minv @ Minv for current_differential in input_differentials_2D: + # print(current_differential) d1, d2 = current_differential # Grabbing the ordered quadruple (i, j, k, l) @@ -434,8 +436,11 @@ def evaluate_hessian_outputs(self, FIM=None): # Identify what index of the symmetric FIM # Hessian arrays need to be updated - row = self.input_names().index(d1_symmetric) - col = self.input_names().index(d2_symmetric) + # Note, we have to account for lower triangular + # and only add to the element which is + # in that section. + row = max(self.input_names().index(d1_symmetric), self.input_names().index(d2_symmetric)) + col = min(self.input_names().index(d1_symmetric), self.input_names().index(d2_symmetric)) flattened_row_col_index = (row + 1) * row // 2 + col # Hessian needs to be handled carefully because of @@ -498,8 +503,8 @@ def evaluate_hessian_outputs(self, FIM=None): # Identify what index of the symmetric FIM # Hessian arrays need to be updated - row = self.input_names().index(d1_symmetric) - col = self.input_names().index(d2_symmetric) + row = max(self.input_names().index(d1_symmetric), self.input_names().index(d2_symmetric)) + col = min(self.input_names().index(d1_symmetric), self.input_names().index(d2_symmetric)) flattened_row_col_index = (row + 1) * row // 2 + col # Hessian needs to be handled carefully because of @@ -586,8 +591,8 @@ def evaluate_hessian_outputs(self, FIM=None): # Identify what index of the symmetric FIM # Hessian arrays need to be updated - row = self.input_names().index(d1_symmetric) - col = self.input_names().index(d2_symmetric) + row = max(self.input_names().index(d1_symmetric), self.input_names().index(d2_symmetric)) + col = min(self.input_names().index(d1_symmetric), self.input_names().index(d2_symmetric)) flattened_row_col_index = (row + 1) * row // 2 + col # Hessian needs to be handled carefully because of @@ -614,11 +619,6 @@ def evaluate_hessian_outputs(self, FIM=None): # # Grab eigenvalues and eigenvectors # Also need the max and min locations - input_differentials_2D = itertools.product(self.input_names(), repeat=2) - hess_vals = [0] * hess_array_length - hess_rows = [0] * hess_array_length - hess_cols = [0] * hess_array_length - all_eig_vals, all_eig_vecs = np.linalg.eig(M) min_eig_loc = np.argmin(all_eig_vals) max_eig_loc = np.argmax(all_eig_vals) @@ -759,8 +759,8 @@ def evaluate_hessian_outputs(self, FIM=None): # Identify what index of the symmetric FIM # Hessian arrays need to be updated - row = self.input_names().index(d1_symmetric) - col = self.input_names().index(d2_symmetric) + row = max(self.input_names().index(d1_symmetric), self.input_names().index(d2_symmetric)) + col = min(self.input_names().index(d1_symmetric), self.input_names().index(d2_symmetric)) flattened_row_col_index = (row + 1) * row // 2 + col # Hessian needs to be handled carefully because of From ab79164f731f94036cdf1eb71c9a018e886b98d3 Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Thu, 2 Oct 2025 16:46:27 -0500 Subject: [PATCH 18/29] Fix hessian computation for FD --- pyomo/contrib/doe/tests/test_greybox.py | 31 ++++++++++++++++++------- 1 file changed, 23 insertions(+), 8 deletions(-) diff --git a/pyomo/contrib/doe/tests/test_greybox.py b/pyomo/contrib/doe/tests/test_greybox.py index 5bebabebc00..1f2fd33ff84 100644 --- a/pyomo/contrib/doe/tests/test_greybox.py +++ b/pyomo/contrib/doe/tests/test_greybox.py @@ -212,21 +212,30 @@ def get_numerical_second_derivative(grey_box_object=None, return_reduced=True): # # Make ordered quads with no repeats # of the ordered pairs - ordered_pairs = itertools.combinations_with_replacement(range(4), 2) + ordered_pairs = itertools.product(range(4), repeat=2) ordered_pairs_list = list(itertools.combinations_with_replacement(range(4), 2)) ordered_quads = itertools.combinations_with_replacement(ordered_pairs, 2) numerical_derivative_reduced = np.zeros((10, 10)) - for i in ordered_quads: - row = ordered_pairs_list.index(i[0]) - col = ordered_pairs_list.index(i[1]) - i, j, k, l = i[0][0], i[0][1], i[1][0], i[1][1] - multiplier = 1 + ((i != j) + (k != l)) + ((i != j) and (k != l)) * (i != k) - numerical_derivative_reduced[row, col] = ( - multiplier * numerical_derivative[i, j, k, l] + for curr_quad in ordered_quads: + d1, d2 = curr_quad + i, j, k, l = d1[0], d1[1], d2[0], d2[1] + + reordered_ijkl = grey_box_object._reorder_pairs(i, j, k, l) + row = ordered_pairs_list.index((reordered_ijkl[2], reordered_ijkl[3])) + col = ordered_pairs_list.index((reordered_ijkl[0], reordered_ijkl[1])) + + numerical_derivative_reduced[row, col] += ( + numerical_derivative[i, j, k, l] ) + # Duplicate check and addition + if ((i != j) and (k != l)) and ((i == l) and (j == k)): + numerical_derivative_reduced[row, col] += ( + numerical_derivative[i, j, k, l] + ) + numerical_derivative_reduced += ( numerical_derivative_reduced.transpose() - np.diag(np.diag(numerical_derivative_reduced)) @@ -634,6 +643,9 @@ def test_hessian_A_opt(self): # Get numerical derivative matrix hess_FD = get_numerical_second_derivative(grey_box_object) + print(hess_gb) + print(hess_FD) + # assert that each component is close self.assertTrue(np.all(np.isclose(hess_gb, hess_FD, rtol=1e-4, atol=1e-4))) @@ -697,8 +709,11 @@ def test_hessian_ME_opt(self): hess_gb = hess_vals_from_gb hess_gb += hess_gb.transpose() - np.diag(np.diag(hess_gb)) + print(hess_gb) + # Get numerical derivative matrix hess_FD = get_numerical_second_derivative(grey_box_object) + print(hess_FD) # assert that each component is close self.assertTrue(np.all(np.isclose(hess_gb, hess_FD, rtol=1e-4, atol=1e-4))) From d709de83f2cfb26fa6bb2b8f8aab03b465292de2 Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Thu, 2 Oct 2025 16:47:02 -0500 Subject: [PATCH 19/29] Remove print statements --- pyomo/contrib/doe/tests/test_greybox.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/pyomo/contrib/doe/tests/test_greybox.py b/pyomo/contrib/doe/tests/test_greybox.py index 1f2fd33ff84..51f41cbe754 100644 --- a/pyomo/contrib/doe/tests/test_greybox.py +++ b/pyomo/contrib/doe/tests/test_greybox.py @@ -643,9 +643,6 @@ def test_hessian_A_opt(self): # Get numerical derivative matrix hess_FD = get_numerical_second_derivative(grey_box_object) - print(hess_gb) - print(hess_FD) - # assert that each component is close self.assertTrue(np.all(np.isclose(hess_gb, hess_FD, rtol=1e-4, atol=1e-4))) @@ -709,11 +706,8 @@ def test_hessian_ME_opt(self): hess_gb = hess_vals_from_gb hess_gb += hess_gb.transpose() - np.diag(np.diag(hess_gb)) - print(hess_gb) - # Get numerical derivative matrix hess_FD = get_numerical_second_derivative(grey_box_object) - print(hess_FD) # assert that each component is close self.assertTrue(np.all(np.isclose(hess_gb, hess_FD, rtol=1e-4, atol=1e-4))) From 87491a78b1693ee069faa04f6e0ef1ec17e8a1ef Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Thu, 2 Oct 2025 16:47:55 -0500 Subject: [PATCH 20/29] ran black --- pyomo/contrib/doe/grey_box_utilities.py | 44 ++++++++++++++++++++----- pyomo/contrib/doe/tests/test_greybox.py | 10 +++--- 2 files changed, 39 insertions(+), 15 deletions(-) diff --git a/pyomo/contrib/doe/grey_box_utilities.py b/pyomo/contrib/doe/grey_box_utilities.py index b2f58130ad3..93b52e967b2 100644 --- a/pyomo/contrib/doe/grey_box_utilities.py +++ b/pyomo/contrib/doe/grey_box_utilities.py @@ -386,7 +386,9 @@ def evaluate_hessian_outputs(self, FIM=None): # Need to iterate over the unique # differentials full_input_names = itertools.product(self._param_names, repeat=2) - input_differentials_2D = itertools.combinations_with_replacement(full_input_names, 2) + input_differentials_2D = itertools.combinations_with_replacement( + full_input_names, 2 + ) from pyomo.contrib.doe import ObjectiveLib @@ -439,8 +441,14 @@ def evaluate_hessian_outputs(self, FIM=None): # Note, we have to account for lower triangular # and only add to the element which is # in that section. - row = max(self.input_names().index(d1_symmetric), self.input_names().index(d2_symmetric)) - col = min(self.input_names().index(d1_symmetric), self.input_names().index(d2_symmetric)) + row = max( + self.input_names().index(d1_symmetric), + self.input_names().index(d2_symmetric), + ) + col = min( + self.input_names().index(d1_symmetric), + self.input_names().index(d2_symmetric), + ) flattened_row_col_index = (row + 1) * row // 2 + col # Hessian needs to be handled carefully because of @@ -503,8 +511,14 @@ def evaluate_hessian_outputs(self, FIM=None): # Identify what index of the symmetric FIM # Hessian arrays need to be updated - row = max(self.input_names().index(d1_symmetric), self.input_names().index(d2_symmetric)) - col = min(self.input_names().index(d1_symmetric), self.input_names().index(d2_symmetric)) + row = max( + self.input_names().index(d1_symmetric), + self.input_names().index(d2_symmetric), + ) + col = min( + self.input_names().index(d1_symmetric), + self.input_names().index(d2_symmetric), + ) flattened_row_col_index = (row + 1) * row // 2 + col # Hessian needs to be handled carefully because of @@ -591,8 +605,14 @@ def evaluate_hessian_outputs(self, FIM=None): # Identify what index of the symmetric FIM # Hessian arrays need to be updated - row = max(self.input_names().index(d1_symmetric), self.input_names().index(d2_symmetric)) - col = min(self.input_names().index(d1_symmetric), self.input_names().index(d2_symmetric)) + row = max( + self.input_names().index(d1_symmetric), + self.input_names().index(d2_symmetric), + ) + col = min( + self.input_names().index(d1_symmetric), + self.input_names().index(d2_symmetric), + ) flattened_row_col_index = (row + 1) * row // 2 + col # Hessian needs to be handled carefully because of @@ -759,8 +779,14 @@ def evaluate_hessian_outputs(self, FIM=None): # Identify what index of the symmetric FIM # Hessian arrays need to be updated - row = max(self.input_names().index(d1_symmetric), self.input_names().index(d2_symmetric)) - col = min(self.input_names().index(d1_symmetric), self.input_names().index(d2_symmetric)) + row = max( + self.input_names().index(d1_symmetric), + self.input_names().index(d2_symmetric), + ) + col = min( + self.input_names().index(d1_symmetric), + self.input_names().index(d2_symmetric), + ) flattened_row_col_index = (row + 1) * row // 2 + col # Hessian needs to be handled carefully because of diff --git a/pyomo/contrib/doe/tests/test_greybox.py b/pyomo/contrib/doe/tests/test_greybox.py index 51f41cbe754..1336b266685 100644 --- a/pyomo/contrib/doe/tests/test_greybox.py +++ b/pyomo/contrib/doe/tests/test_greybox.py @@ -226,15 +226,13 @@ def get_numerical_second_derivative(grey_box_object=None, return_reduced=True): row = ordered_pairs_list.index((reordered_ijkl[2], reordered_ijkl[3])) col = ordered_pairs_list.index((reordered_ijkl[0], reordered_ijkl[1])) - numerical_derivative_reduced[row, col] += ( - numerical_derivative[i, j, k, l] - ) + numerical_derivative_reduced[row, col] += numerical_derivative[i, j, k, l] # Duplicate check and addition if ((i != j) and (k != l)) and ((i == l) and (j == k)): - numerical_derivative_reduced[row, col] += ( - numerical_derivative[i, j, k, l] - ) + numerical_derivative_reduced[row, col] += numerical_derivative[ + i, j, k, l + ] numerical_derivative_reduced += ( numerical_derivative_reduced.transpose() From 6d526327e6d729da857253f02b1c08dce353f184 Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Wed, 8 Oct 2025 12:14:59 -0400 Subject: [PATCH 21/29] Housekeeping, remove old stuff --- pyomo/contrib/doe/grey_box_utilities.py | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/pyomo/contrib/doe/grey_box_utilities.py b/pyomo/contrib/doe/grey_box_utilities.py index 93b52e967b2..571681d8455 100644 --- a/pyomo/contrib/doe/grey_box_utilities.py +++ b/pyomo/contrib/doe/grey_box_utilities.py @@ -246,7 +246,7 @@ def evaluate_jacobian_equality_constraints(self): def evaluate_jacobian_outputs(self): # Compute the jacobian of the objective function with - # respect to the fisher information matrix. Then return + # respect to the fisher information matrix. Then, return # a coo_matrix that aligns with what IPOPT will expect. current_FIM = self._get_FIM() @@ -313,13 +313,6 @@ def evaluate_jacobian_outputs(self): min_eig_term = min_eig_vec * np.transpose(min_eig_vec) max_eig_term = max_eig_vec * np.transpose(max_eig_vec) - min_eig_epsilon = 2e-16 - - # Computing a (hopefully) nonsingular - # condition number for the jacobian - # expression. - safe_cond_number = max_eig / (min_eig + np.sign(min_eig) * min_eig_epsilon) - # Combining the expression jac_M = 1 / max_eig * max_eig_term - 1 / min_eig * min_eig_term else: From d08066004b540fe251f77ce17154f8658b1bc015 Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Wed, 8 Oct 2025 12:16:23 -0400 Subject: [PATCH 22/29] more descriptive explanations --- pyomo/contrib/doe/grey_box_utilities.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/pyomo/contrib/doe/grey_box_utilities.py b/pyomo/contrib/doe/grey_box_utilities.py index 571681d8455..92e4d67f25d 100644 --- a/pyomo/contrib/doe/grey_box_utilities.py +++ b/pyomo/contrib/doe/grey_box_utilities.py @@ -318,12 +318,13 @@ def evaluate_jacobian_outputs(self): else: ObjectiveLib(self.objective_option) - # We are only using the triangle, - # so we need to add the off-diagonals - # twice. + # We are only using a symmetric, triangular + # representation of the FIM, so we need + # to add the off-diagonal elements twice. jac_M = 2 * jac_M - np.diag(np.diag(jac_M)) - # Filter jac_M using the - # masking matrix + # Filter the Jacobian, jac_M, using the + # masking matrix to only select the + # symmetric, triangular components jac_M = jac_M[self._masking_matrix > 0] M_rows = np.zeros((len(jac_M.flatten()), 1)).flatten() M_cols = np.arange(len(jac_M.flatten())) From c8e8c0ad0318b9b3cd466b30d69b27a6876c2214 Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Wed, 8 Oct 2025 12:20:56 -0400 Subject: [PATCH 23/29] Update hessian commenting --- pyomo/contrib/doe/grey_box_utilities.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pyomo/contrib/doe/grey_box_utilities.py b/pyomo/contrib/doe/grey_box_utilities.py index 92e4d67f25d..daaac7474f6 100644 --- a/pyomo/contrib/doe/grey_box_utilities.py +++ b/pyomo/contrib/doe/grey_box_utilities.py @@ -355,8 +355,9 @@ def evaluate_hessian_equality_constraints(self): return None def evaluate_hessian_outputs(self, FIM=None): - # TODO: significant bookkeeping if the hessian's require vectorized - # operations. Just need mapping that works well and we are good. + # Compute the hessian of the objective function with + # respect to the fisher information matrix. Then, return + # a coo_matrix that aligns with what IPOPT will expect. current_FIM = self._get_FIM() M = np.asarray(current_FIM, dtype=np.float64).reshape( From b6a0c2e669694fdd6031ffb0474f219c04199c9d Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Wed, 8 Oct 2025 13:44:15 -0400 Subject: [PATCH 24/29] Remove input from old testing functions --- pyomo/contrib/doe/grey_box_utilities.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyomo/contrib/doe/grey_box_utilities.py b/pyomo/contrib/doe/grey_box_utilities.py index daaac7474f6..44037ba75a4 100644 --- a/pyomo/contrib/doe/grey_box_utilities.py +++ b/pyomo/contrib/doe/grey_box_utilities.py @@ -354,7 +354,7 @@ def evaluate_hessian_equality_constraints(self): # No constraints so this returns `None` return None - def evaluate_hessian_outputs(self, FIM=None): + def evaluate_hessian_outputs(self): # Compute the hessian of the objective function with # respect to the fisher information matrix. Then, return # a coo_matrix that aligns with what IPOPT will expect. From 769a1ba7b4a052e248b1d65cc84d1692c56f63f9 Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Wed, 8 Oct 2025 13:49:52 -0400 Subject: [PATCH 25/29] Increasing Hessian description verbosity in comments --- pyomo/contrib/doe/grey_box_utilities.py | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/pyomo/contrib/doe/grey_box_utilities.py b/pyomo/contrib/doe/grey_box_utilities.py index 44037ba75a4..d51132e0a76 100644 --- a/pyomo/contrib/doe/grey_box_utilities.py +++ b/pyomo/contrib/doe/grey_box_utilities.py @@ -364,23 +364,35 @@ def evaluate_hessian_outputs(self): self._n_params, self._n_params ) - # Length of Hessian lists for the sparse - # matrix ar a function of number of parameters + # We will store the Hessian values in + # vectorized (flattened) format. The length + # of the vectorized Hessian for the symmetric + # FIM representation scales by the number of + # unknown parameters. hess_array_length = round( (((self._n_params + 1) * self._n_params / 2) + 1) * (((self._n_params + 1) * self._n_params / 2)) / 2 ) - # Hessian with correct size for using only the - # lower (upper) triangle of the FIM + # Initializing lists of the correct length + # for the hessian values and the row and column + # of these data in the coo matrix to be returned hess_vals = [0] * hess_array_length hess_rows = [0] * hess_array_length hess_cols = [0] * hess_array_length - # Need to iterate over the unique - # differentials + # We are utilizing the symmetric Hessian, but we + # must consider the contribution from all elements. + # Therefore, we are required to use the full product + # space of the parameter names (full FIM) to compute + # to Hessian of the symmetric FIM. full_input_names = itertools.product(self._param_names, repeat=2) + + # Here, we use combination with replacement to only + # consider the upper triangle of the Hessian for the + # full FIM. We will map these second derivative values + # back onto the symmetric FIM Hessian. input_differentials_2D = itertools.combinations_with_replacement( full_input_names, 2 ) From 7ff2861a98e18095e6759c8de268514e04e4f0e2 Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Wed, 8 Oct 2025 13:50:03 -0400 Subject: [PATCH 26/29] remove old print statement --- pyomo/contrib/doe/grey_box_utilities.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pyomo/contrib/doe/grey_box_utilities.py b/pyomo/contrib/doe/grey_box_utilities.py index d51132e0a76..4ae0900b91b 100644 --- a/pyomo/contrib/doe/grey_box_utilities.py +++ b/pyomo/contrib/doe/grey_box_utilities.py @@ -407,7 +407,6 @@ def evaluate_hessian_outputs(self): Minv_sq = Minv @ Minv for current_differential in input_differentials_2D: - # print(current_differential) d1, d2 = current_differential # Grabbing the ordered quadruple (i, j, k, l) From c74f011fb443fa9b2ecbe7419bcca5208aeceee7 Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Wed, 8 Oct 2025 13:53:51 -0400 Subject: [PATCH 27/29] Improving verbosity and descriptions in comments --- pyomo/contrib/doe/grey_box_utilities.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/pyomo/contrib/doe/grey_box_utilities.py b/pyomo/contrib/doe/grey_box_utilities.py index 4ae0900b91b..c7c6501f96c 100644 --- a/pyomo/contrib/doe/grey_box_utilities.py +++ b/pyomo/contrib/doe/grey_box_utilities.py @@ -443,10 +443,9 @@ def evaluate_hessian_outputs(self): ) # Identify what index of the symmetric FIM - # Hessian arrays need to be updated - # Note, we have to account for lower triangular - # and only add to the element which is - # in that section. + # Hessian arrays need to be updated. + # Note: we are only interested in building + # the lower triangular portion of the Hessian. row = max( self.input_names().index(d1_symmetric), self.input_names().index(d2_symmetric), @@ -458,13 +457,13 @@ def evaluate_hessian_outputs(self): flattened_row_col_index = (row + 1) * row // 2 + col # Hessian needs to be handled carefully because of - # the ``missing`` components when only passing - # a symmetric version of the FIM. + # the ``missing`` components from the full FIM + # when only passing a symmetric version of the FIM. # # When we reordered (i, j, k, l), we are correctly # pointing to which index needs to be contributed to. # However, when an element that is not included - # is being put into a diagonal element of the + # is being mapped to a diagonal element of the # symmetric FIM hessian from the full FIM hessian, # it needs to be counted twice. This only occurs # when (i != j) and (k != l) and (i, j) and (k, l) @@ -472,10 +471,12 @@ def evaluate_hessian_outputs(self): # (i == l) and (j == k). # # Otherwise, we only add the element once. + # Standard addition hess_vals[flattened_row_col_index] += hess_contribution - # Duplicate check and addition + # Duplicate check and addition if + # criteria is satisfied. if ((i != j) and (k != l)) and ((i == l) and (j == k)): hess_vals[flattened_row_col_index] += hess_contribution From ac6ec630eb332ddd9b015dd2d8ea2786fccff142 Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Wed, 8 Oct 2025 14:22:35 -0400 Subject: [PATCH 28/29] More grammar updates --- pyomo/contrib/doe/grey_box_utilities.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/pyomo/contrib/doe/grey_box_utilities.py b/pyomo/contrib/doe/grey_box_utilities.py index c7c6501f96c..75f4090bfd5 100644 --- a/pyomo/contrib/doe/grey_box_utilities.py +++ b/pyomo/contrib/doe/grey_box_utilities.py @@ -637,11 +637,12 @@ def evaluate_hessian_outputs(self): elif self.objective_option == ObjectiveLib.condition_number: # Hessian for log condition number has 4 - # terms. Two are multiples of the second - # derivative of the maximum and minimum - # eigenvalues. The other two are tensor products + # terms. The first and third terms are + # multiples of the second derivative of the + # maximum and minimum eigenvalues, respectively + # The other two are tensor products # of the first derivative of the maximum - # eigenvalue with itself, and minimum + # eigenvalue with itself, and the minimum # eigenvalue with itself. # # Grab eigenvalues and eigenvectors From c4d9c9ee9b05c325b96413bbcdb8feb67937a156 Mon Sep 17 00:00:00 2001 From: Daniel Laky <29078718+djlaky@users.noreply.github.com> Date: Wed, 15 Oct 2025 17:13:52 -0400 Subject: [PATCH 29/29] Update pyomo/contrib/doe/grey_box_utilities.py Co-authored-by: Bethany Nicholson --- pyomo/contrib/doe/grey_box_utilities.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyomo/contrib/doe/grey_box_utilities.py b/pyomo/contrib/doe/grey_box_utilities.py index 75f4090bfd5..ca5a20250cb 100644 --- a/pyomo/contrib/doe/grey_box_utilities.py +++ b/pyomo/contrib/doe/grey_box_utilities.py @@ -386,7 +386,7 @@ def evaluate_hessian_outputs(self): # must consider the contribution from all elements. # Therefore, we are required to use the full product # space of the parameter names (full FIM) to compute - # to Hessian of the symmetric FIM. + # the Hessian of the symmetric FIM. full_input_names = itertools.product(self._param_names, repeat=2) # Here, we use combination with replacement to only