Merge pull request #3727 from djlaky/add-greybox

Pyomo.DoE: Bugfix for GreyBox

Author: mrmundt

pyomo/contrib/doe/doe.py

@@ -359,8 +359,7 @@ def run_doe(self, model=None, results_file=None):
                         )
             # Set objective value
             if self.objective_option == ObjectiveLib.trace:
-                # Do safe inverse here?
-                trace_val = 1 / np.trace(np.array(self.get_FIM()))
+                trace_val = np.trace(np.linalg.pinv(self.get_FIM()))
                 model.obj_cons.egb_fim_block.outputs["A-opt"].set_value(trace_val)
             elif self.objective_option == ObjectiveLib.determinant:
                 det_val = np.linalg.det(np.array(self.get_FIM()))

pyomo/contrib/doe/grey_box_utilities.py

@@ -310,6 +310,10 @@ 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.
+        jac_M = 2 * jac_M - np.diag(np.diag(jac_M))
         # Filter jac_M using the
         # masking matrix
         jac_M = jac_M[self._masking_matrix > 0]
@@ -402,6 +406,24 @@ def evaluate_hessian_outputs(self, FIM=None):
                 hess_rows.append(row)
                 hess_cols.append(col)
 
+                # Hessian needs to be handled carefully because of
+                # the ``missing`` components when only passing
+                # a triangular 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])
+                )
+
         elif self.objective_option == ObjectiveLib.determinant:
             # Grab inverse
             Minv = np.linalg.pinv(M)
@@ -424,6 +446,13 @@ def evaluate_hessian_outputs(self, FIM=None):
                 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])
         elif self.objective_option == ObjectiveLib.minimum_eigenvalue:
             # Grab eigenvalues and eigenvectors
             # Also need the min location
@@ -482,6 +511,13 @@ def evaluate_hessian_outputs(self, FIM=None):
                 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
         elif self.objective_option == ObjectiveLib.condition_number:
             pass
         else:

pyomo/contrib/doe/tests/test_greybox.py

@@ -215,9 +215,11 @@ def get_numerical_second_derivative(grey_box_object=None, return_reduced=True):
         for i in ordered_quads:
             row = ordered_pairs_list.index(i[0])
             col = ordered_pairs_list.index(i[1])
-            numerical_derivative_reduced[row, col] = numerical_derivative[
-                i[0][0], i[0][1], i[1][0], i[1][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]
+            )
 
         numerical_derivative_reduced += (
             numerical_derivative_reduced.transpose()
@@ -524,7 +526,8 @@ def test_jacobian_A_opt(self):
         # Recover the Jacobian in Matrix Form
         jac = np.zeros_like(grey_box_object._get_FIM())
         jac[np.triu_indices_from(jac)] = utri_vals_jac
-        jac += jac.transpose() - np.diag(np.diag(jac))
+        jac += jac.transpose()
+        jac = jac / 2
 
         # Get numerical derivative matrix
         jac_FD = get_numerical_derivative(grey_box_object)
@@ -547,7 +550,8 @@ def test_jacobian_D_opt(self):
         # Recover the Jacobian in Matrix Form
         jac = np.zeros_like(grey_box_object._get_FIM())
         jac[np.triu_indices_from(jac)] = utri_vals_jac
-        jac += jac.transpose() - np.diag(np.diag(jac))
+        jac += jac.transpose()
+        jac = jac / 2
 
         # Get numerical derivative matrix
         jac_FD = get_numerical_derivative(grey_box_object)
@@ -570,7 +574,8 @@ def test_jacobian_E_opt(self):
         # Recover the Jacobian in Matrix Form
         jac = np.zeros_like(grey_box_object._get_FIM())
         jac[np.triu_indices_from(jac)] = utri_vals_jac
-        jac += jac.transpose() - np.diag(np.diag(jac))
+        jac += jac.transpose()
+        jac = jac / 2
 
         # Get numerical derivative matrix
         jac_FD = get_numerical_derivative(grey_box_object)
@@ -593,7 +598,8 @@ def test_jacobian_ME_opt(self):
         # Recover the Jacobian in Matrix Form
         jac = np.zeros_like(grey_box_object._get_FIM())
         jac[np.triu_indices_from(jac)] = utri_vals_jac
-        jac += jac.transpose() - np.diag(np.diag(jac))
+        jac += jac.transpose()
+        jac = jac / 2
 
         # Get numerical derivative matrix
         jac_FD = get_numerical_derivative(grey_box_object)
@@ -990,7 +996,10 @@ def test_solve_D_optimality_log_determinant(self):
         # (time, optimal objective value)
         # Here, the objective value is
         # log-10(determinant) of the FIM
-        optimal_experimental_designs = [np.array([2.24, 4.33]), np.array([10.00, 4.35])]
+        optimal_experimental_designs = [
+            np.array([1.782, 4.34]),
+            np.array([10.00, 4.35]),
+        ]
         objective_option = "determinant"
         doe_object, grey_box_object = make_greybox_and_doe_objects_rooney_biegler(
             objective_option=objective_option
@@ -1002,12 +1011,6 @@ def test_solve_D_optimality_log_determinant(self):
         # Solve the model
         doe_object.run_doe()
 
-        print("Termination Message")
-        print(doe_object.results["Termination Message"])
-        print(cyipopt_call_working)
-        print(bad_message in doe_object.results["Termination Message"])
-        print("End Message")
-
         optimal_time_val = doe_object.results["Experiment Design"][0]
         optimal_obj_val = np.log10(np.exp(pyo.value(doe_object.model.objective)))
 
@@ -1035,7 +1038,7 @@ def test_solve_A_optimality_trace_of_inverse(self):
         # Here, the objective value is
         # trace(inverse(FIM))
         optimal_experimental_designs = [
-            np.array([1.94, 0.0295]),
+            np.array([1.33, 0.0277]),
             np.array([9.9, 0.0366]),
         ]
         objective_option = "trace"
@@ -1076,7 +1079,7 @@ def test_solve_E_optimality_minimum_eigenvalue(self):
         # Here, the objective value is
         # minimum eigenvalue of the FIM
         optimal_experimental_designs = [
-            np.array([1.92, 36.018]),
+            np.array([1.30, 38.70]),
             np.array([10.00, 28.349]),
         ]
         objective_option = "minimum_eigenvalue"
@@ -1117,7 +1120,7 @@ def test_solve_ME_optimality_condition_number(self):
         # Here, the objective value is
         # condition number of the FIM
         optimal_experimental_designs = [
-            np.array([1.59, 15.22]),
+            np.array([0.943, 13.524]),
             np.array([10.00, 27.675]),
         ]
         objective_option = "condition_number"