From cc2803483b05f6c9348b303aac4105d181346e49 Mon Sep 17 00:00:00 2001
From: Alex Dowling <alex@dowlinglab.org>
Date: Tue, 23 Apr 2024 21:14:32 -0400
Subject: [PATCH 01/11] Fixed error about if values being ambiguous

---
 pyomo/contrib/doe/measurements.py | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/pyomo/contrib/doe/measurements.py b/pyomo/contrib/doe/measurements.py
index 5a3c44a76e4..11b84b4231b 100644
--- a/pyomo/contrib/doe/measurements.py
+++ b/pyomo/contrib/doe/measurements.py
@@ -93,17 +93,17 @@ def add_variables(
             upper_bounds,
         )
 
-        if values:
+        if values is not None:
             # this dictionary keys are special set, values are its value
             self.variable_names_value.update(zip(added_names, values))
 
         # if a scalar (int or float) is given, set it as the lower bound for all variables
-        if lower_bounds:
+        if lower_bounds is not None:
             if type(lower_bounds) in [int, float]:
                 lower_bounds = [lower_bounds] * len(added_names)
             self.lower_bounds.update(zip(added_names, lower_bounds))
 
-        if upper_bounds:
+        if upper_bounds is not None:
             if type(upper_bounds) in [int, float]:
                 upper_bounds = [upper_bounds] * len(added_names)
             self.upper_bounds.update(zip(added_names, upper_bounds))
@@ -177,20 +177,20 @@ def _check_valid_input(
             raise ValueError("time index cannot be found in indices.")
 
         # if given a list, check if bounds have the same length with flattened variable
-        if values and len(values) != len_indices:
+        if values is not None and len(values) != len_indices:
             raise ValueError("Values is of different length with indices.")
 
         if (
-            lower_bounds
-            and type(lower_bounds) == list
-            and len(lower_bounds) != len_indices
+            lower_bounds is not None # ensure not None
+            and type(lower_bounds) == list # ensure list
+            and len(lower_bounds) != len_indices # ensure same length
         ):
             raise ValueError("Lowerbounds is of different length with indices.")
 
         if (
-            upper_bounds
-            and type(upper_bounds) == list
-            and len(upper_bounds) != len_indices
+            upper_bounds is not None # ensure None
+            and type(upper_bounds) == list # ensure list
+            and len(upper_bounds) != len_indices # ensure same length
         ):
             raise ValueError("Upperbounds is of different length with indices.")
 

From 6d421a428c46a7621fb1f356406d392b3bc0e9e3 Mon Sep 17 00:00:00 2001
From: Alex Dowling <alex@dowlinglab.org>
Date: Tue, 23 Apr 2024 21:31:39 -0400
Subject: [PATCH 02/11] Added exception to prevent cryptic error messages
 later.

---
 pyomo/contrib/doe/doe.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py
index d2ba2f277d6..eba22b954cf 100644
--- a/pyomo/contrib/doe/doe.py
+++ b/pyomo/contrib/doe/doe.py
@@ -101,6 +101,8 @@ def __init__(
         """
 
         # parameters
+        if type(param_init) != dict:
+            raise ValueError("param_init should be a dictionary.")
         self.param = param_init
         # design variable name
         self.design_name = design_vars.variable_names

From 910cf2d1247f0d75e80b6bd776236e8e05ae6beb Mon Sep 17 00:00:00 2001
From: Alex Dowling <alex@dowlinglab.org>
Date: Wed, 24 Apr 2024 21:04:11 -0400
Subject: [PATCH 03/11] Added units.

---
 pyomo/contrib/doe/doe.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py
index eba22b954cf..8c83e3145ce 100644
--- a/pyomo/contrib/doe/doe.py
+++ b/pyomo/contrib/doe/doe.py
@@ -240,12 +240,12 @@ def stochastic_program(
         if self.optimize:
             analysis_optimize = self._optimize_stochastic_program(m)
             dT = sp_timer.toc(msg=None)
-            self.logger.info("elapsed time: %0.1f" % dT)
+            self.logger.info("elapsed time: %0.1f seconds" % dT)
             return analysis_square, analysis_optimize
 
         else:
             dT = sp_timer.toc(msg=None)
-            self.logger.info("elapsed time: %0.1f" % dT)
+            self.logger.info("elapsed time: %0.1f seconds" % dT)
             return analysis_square
 
     def _compute_stochastic_program(self, m, optimize_option):
@@ -389,7 +389,7 @@ def compute_FIM(
             FIM_analysis = self._direct_kaug()
 
         dT = square_timer.toc(msg=None)
-        self.logger.info("elapsed time: %0.1f" % dT)
+        self.logger.info("elapsed time: %0.1f seconds" % dT)
 
         return FIM_analysis
 

From 5df54523e46be0e66aac827e21286a9977592011 Mon Sep 17 00:00:00 2001
From: Alex Dowling <alex@dowlinglab.org>
Date: Wed, 24 Apr 2024 21:52:51 -0400
Subject: [PATCH 04/11] Switched to isinstance

---
 pyomo/contrib/doe/doe.py          | 6 +++---
 pyomo/contrib/doe/measurements.py | 7 ++++---
 2 files changed, 7 insertions(+), 6 deletions(-)

diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py
index 8c83e3145ce..b55f9aac0ac 100644
--- a/pyomo/contrib/doe/doe.py
+++ b/pyomo/contrib/doe/doe.py
@@ -38,7 +38,7 @@
 from pyomo.contrib.sensitivity_toolbox.sens import get_dsdp
 from pyomo.contrib.doe.scenario import ScenarioGenerator, FiniteDifferenceStep
 from pyomo.contrib.doe.result import FisherResults, GridSearchResult
-
+import collections
 
 class CalculationMode(Enum):
     sequential_finite = "sequential_finite"
@@ -101,7 +101,7 @@ def __init__(
         """
 
         # parameters
-        if type(param_init) != dict:
+        if not isinstance(param_init, collections.Mapping):
             raise ValueError("param_init should be a dictionary.")
         self.param = param_init
         # design variable name
@@ -777,7 +777,7 @@ def run_grid_search(
             # update the controlled value of certain time points for certain design variables
             for i, names in enumerate(design_dimension_names):
                 # if the element is a list, all design variables in this list share the same values
-                if type(names) is list or type(names) is tuple:
+                if isinstance(names, collections.Sequence):
                     for n in names:
                         design_iter[n] = list(design_set_iter)[i]
                 else:
diff --git a/pyomo/contrib/doe/measurements.py b/pyomo/contrib/doe/measurements.py
index 11b84b4231b..aa196ec9a49 100644
--- a/pyomo/contrib/doe/measurements.py
+++ b/pyomo/contrib/doe/measurements.py
@@ -26,6 +26,7 @@
 #  ___________________________________________________________________________
 
 import itertools
+import collections
 
 
 class VariablesWithIndices:
@@ -171,7 +172,7 @@ def _check_valid_input(
         """
         Check if the measurement information provided are valid to use.
         """
-        assert type(var_name) is str, "var_name should be a string."
+        assert isinstance(var_name, str), "var_name should be a string."
 
         if time_index_position not in indices:
             raise ValueError("time index cannot be found in indices.")
@@ -182,14 +183,14 @@ def _check_valid_input(
 
         if (
             lower_bounds is not None # ensure not None
-            and type(lower_bounds) == list # ensure list
+            and isinstance(lower_bounds, collections.Sequence) # ensure list-like
             and len(lower_bounds) != len_indices # ensure same length
         ):
             raise ValueError("Lowerbounds is of different length with indices.")
 
         if (
             upper_bounds is not None # ensure None
-            and type(upper_bounds) == list # ensure list
+            and isinstance(upper_bounds, collections.Sequence) # ensure list-like
             and len(upper_bounds) != len_indices # ensure same length
         ):
             raise ValueError("Upperbounds is of different length with indices.")

From 312f5722b1a703413133c887b415a0d249d1e032 Mon Sep 17 00:00:00 2001
From: Alex Dowling <alex@dowlinglab.org>
Date: Thu, 25 Apr 2024 21:46:59 -0400
Subject: [PATCH 05/11] Initialization improvements and degree of freedom fixes
 are important!

---
 pyomo/contrib/doe/doe.py          | 116 +++++++++++++++++++++++++++---
 pyomo/contrib/doe/measurements.py |   6 +-
 2 files changed, 108 insertions(+), 14 deletions(-)

diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py
index b55f9aac0ac..cac9bfd9271 100644
--- a/pyomo/contrib/doe/doe.py
+++ b/pyomo/contrib/doe/doe.py
@@ -38,7 +38,9 @@
 from pyomo.contrib.sensitivity_toolbox.sens import get_dsdp
 from pyomo.contrib.doe.scenario import ScenarioGenerator, FiniteDifferenceStep
 from pyomo.contrib.doe.result import FisherResults, GridSearchResult
-import collections
+import collections.abc
+
+import inspect
 
 class CalculationMode(Enum):
     sequential_finite = "sequential_finite"
@@ -101,7 +103,7 @@ def __init__(
         """
 
         # parameters
-        if not isinstance(param_init, collections.Mapping):
+        if not isinstance(param_init, collections.abc.Mapping):
             raise ValueError("param_init should be a dictionary.")
         self.param = param_init
         # design variable name
@@ -238,6 +240,10 @@ def stochastic_program(
         m, analysis_square = self._compute_stochastic_program(m, optimize_opt)
 
         if self.optimize:
+            # set max_iter to 0 to debug the initialization
+            # self.solver.options["max_iter"] = 0
+            # self.solver.options["bound_push"] = 1e-10
+            
             analysis_optimize = self._optimize_stochastic_program(m)
             dT = sp_timer.toc(msg=None)
             self.logger.info("elapsed time: %0.1f seconds" % dT)
@@ -588,13 +594,34 @@ def _create_block(self):
 
         # Set for block/scenarios
         mod.scenario = pyo.Set(initialize=self.scenario_data.scenario_indices)
+        
+        # Determine if create_model takes theta as an optional input
+        # print(inspect.getfullargspec(self.create_model))
+        pass_theta_to_initialize= ('theta' in inspect.getfullargspec(self.create_model).args)
+        #print("pass_theta_to_initialize =", pass_theta_to_initialize)
 
         # Allow user to self-define complex design variables
         self.create_model(mod=mod, model_option=ModelOptionLib.stage1)
 
+        # Fix parameter values in the copy of the stage1 model (if they exist)
+        for par in self.param:
+            cuid = pyo.ComponentUID(par)
+            var = cuid.find_component_on(mod)
+            if var is not None:
+                # Fix the parameter value
+                # Otherwise, the parameter does not exist on the stage 1 model
+                var.fix(self.param[par])
+
         def block_build(b, s):
             # create block scenarios
-            self.create_model(mod=b, model_option=ModelOptionLib.stage2)
+            # idea: check if create_model takes theta as an optional input, if so, pass parameter values to create_model
+
+            if pass_theta_to_initialize:
+                theta_initialize = self.scenario_data.scenario[s]
+                #print("Initializing with theta=", theta_initialize)
+                self.create_model(mod=b, model_option=ModelOptionLib.stage2, theta=theta_initialize)
+            else:
+                self.create_model(mod=b, model_option=ModelOptionLib.stage2)
 
             # fix parameter values to perturbed values
             for par in self.param:
@@ -605,7 +632,7 @@ def block_build(b, s):
         mod.block = pyo.Block(mod.scenario, rule=block_build)
 
         # discretize the model
-        if self.discretize_model:
+        if self.discretize_model is not None:
             mod = self.discretize_model(mod)
 
         # force design variables in blocks to be equal to global design values
@@ -777,7 +804,7 @@ def run_grid_search(
             # update the controlled value of certain time points for certain design variables
             for i, names in enumerate(design_dimension_names):
                 # if the element is a list, all design variables in this list share the same values
-                if isinstance(names, collections.Sequence):
+                if isinstance(names, collections.abc.Sequence):
                     for n in names:
                         design_iter[n] = list(design_set_iter)[i]
                 else:
@@ -881,20 +908,39 @@ def identity_matrix(m, i, j):
             else:
                 return 0
 
+        if self.jac_initial is not None:
+            dict_jac_initialize = {}
+            for i, bu in enumerate(model.regression_parameters):
+                for j, un in enumerate(model.measured_variables):
+                    if isinstance(self.jac_initial, dict):
+                        # Jacobian is a dictionary of arrays or lists where the key is the regression parameter name
+                        dict_jac_initialize[(bu, un)] = self.jac_initial[bu][j]
+                    elif isinstance(self.jac_initial, np.ndarray):
+                        # Jacobian is a numpy array, rows are regression parameters, columns are measured variables
+                        dict_jac_initialize[(bu, un)] = self.jac_initial[i][j]
+
+        def initialize_jac(m, i, j):
+            if self.jac_initial is not None:
+                return dict_jac_initialize[(i, j)]
+            else:
+                return 0.1
+
         model.sensitivity_jacobian = pyo.Var(
-            model.regression_parameters, model.measured_variables, initialize=0.1
+            model.regression_parameters, model.measured_variables, initialize=initialize_jac
         )
 
-        if self.fim_initial:
+        if self.fim_initial is not None:
             dict_fim_initialize = {}
             for i, bu in enumerate(model.regression_parameters):
                 for j, un in enumerate(model.regression_parameters):
                     dict_fim_initialize[(bu, un)] = self.fim_initial[i][j]
+            
+            #print(dict_fim_initialize)
 
         def initialize_fim(m, j, d):
             return dict_fim_initialize[(j, d)]
 
-        if self.fim_initial:
+        if self.fim_initial is not None:
             model.fim = pyo.Var(
                 model.regression_parameters,
                 model.regression_parameters,
@@ -1013,6 +1059,32 @@ def fim_rule(m, p, q):
         return model
 
     def _add_objective(self, m):
+
+        ### Initialize the Cholesky decomposition matrix
+        if self.Cholesky_option:
+
+            # Assemble the FIM matrix
+            fim = np.zeros((len(self.param), len(self.param)))
+            for i, bu in enumerate(m.regression_parameters):
+                for j, un in enumerate(m.regression_parameters):
+                    fim[i][j] = m.fim[bu, un].value
+
+            # Calculate the eigenvalues of the FIM matrix
+            eig = np.linalg.eigvals(fim)
+
+            # If the smallest eigenvalue is (pratcially) negative, add a diagonal matrix to make it positive definite
+            small_number = 1E-10
+            if min(eig) < small_number:
+                fim = fim + np.eye(len(self.param)) * (small_number - min(eig))
+
+            # Compute the Cholesky decomposition of the FIM matrix
+            L = np.linalg.cholesky(fim)
+
+        # Initialize the Cholesky matrix
+        for i, c in enumerate(m.regression_parameters):
+            for j, d in enumerate(m.regression_parameters):
+                m.L_ele[c, d].value = L[i, j]
+
         def cholesky_imp(m, c, d):
             """
             Calculate Cholesky L matrix using algebraic constraints
@@ -1103,14 +1175,20 @@ def _fix_design(self, m, design_val, fix_opt=True, optimize_option=None):
         m: model
         """
         for name in self.design_name:
+            # Loop over design variables
+            # Get Pyomo variable object
             cuid = pyo.ComponentUID(name)
             var = cuid.find_component_on(m)
             if fix_opt:
+                # If fix_opt is True, fix the design variable
                 var.fix(design_val[name])
             else:
+                # Otherwise check optimize_option
                 if optimize_option is None:
+                    # If optimize_option is None, unfix all design variables
                     var.unfix()
                 else:
+                    # Otherwise, unfix only the design variables listed in optimize_option with value True
                     if optimize_option[name]:
                         var.unfix()
         return m
@@ -1126,7 +1204,7 @@ def _get_default_ipopt_solver(self):
     def _solve_doe(self, m, fix=False, opt_option=None):
         """Solve DOE model.
         If it's a square problem, fix design variable and solve.
-        Else, fix design variable and solve square problem firstly, then unfix them and solve the optimization problem
+        Else, fix design variable and solve square problem first, then unfix them and solve the optimization problem
 
         Parameters
         ----------
@@ -1140,14 +1218,30 @@ def _solve_doe(self, m, fix=False, opt_option=None):
         -------
         solver_results: solver results
         """
-        ### Solve square problem
+        # if fix = False, solve the optimization problem
+        # if fix = True, solve the square problem
+
+        # either fix or unfix the design variables
         mod = self._fix_design(
             m, self.design_values, fix_opt=fix, optimize_option=opt_option
         )
 
+        '''
+        # This is for initialization diagnostics
+        # Remove before merging the PR
+        if not fix:
+            # halt at initial point
+            self.solver.options['max_iter'] = 0
+            self.solver.options['bound_push'] = 1E-10
+        else:
+            # resort to defaults
+            self.solver.options['max_iter'] = 3000
+            self.solver.options['bound_push'] = 0.01
+        '''
+
         # if user gives solver, use this solver. if not, use default IPOPT solver
         solver_result = self.solver.solve(mod, tee=self.tee_opt)
-
+            
         return solver_result
 
     def _sgn(self, p):
diff --git a/pyomo/contrib/doe/measurements.py b/pyomo/contrib/doe/measurements.py
index aa196ec9a49..ae5b3519498 100644
--- a/pyomo/contrib/doe/measurements.py
+++ b/pyomo/contrib/doe/measurements.py
@@ -26,7 +26,7 @@
 #  ___________________________________________________________________________
 
 import itertools
-import collections
+import collections.abc
 
 
 class VariablesWithIndices:
@@ -183,14 +183,14 @@ def _check_valid_input(
 
         if (
             lower_bounds is not None # ensure not None
-            and isinstance(lower_bounds, collections.Sequence) # ensure list-like
+            and isinstance(lower_bounds, collections.abc.Sequence) # ensure list-like
             and len(lower_bounds) != len_indices # ensure same length
         ):
             raise ValueError("Lowerbounds is of different length with indices.")
 
         if (
             upper_bounds is not None # ensure None
-            and isinstance(upper_bounds, collections.Sequence) # ensure list-like
+            and isinstance(upper_bounds, collections.abc.Sequence) # ensure list-like
             and len(upper_bounds) != len_indices # ensure same length
         ):
             raise ValueError("Upperbounds is of different length with indices.")

From 78c8558d3c0bc0beca6683ffe6b124223e501f6d Mon Sep 17 00:00:00 2001
From: Alex Dowling <alex@dowlinglab.org>
Date: Thu, 25 Apr 2024 21:53:56 -0400
Subject: [PATCH 06/11] Updated one type check and removed extra code from
 debugging.

---
 pyomo/contrib/doe/doe.py          | 29 +++++++----------------------
 pyomo/contrib/doe/measurements.py |  6 +++---
 2 files changed, 10 insertions(+), 25 deletions(-)

diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py
index cac9bfd9271..d10fc1b7fcc 100644
--- a/pyomo/contrib/doe/doe.py
+++ b/pyomo/contrib/doe/doe.py
@@ -230,6 +230,7 @@ def stochastic_program(
         # FIM = Jacobian.T@Jacobian, the FIM is scaled by squared value the Jacobian is scaled
         self.fim_scale_constant_value = self.scale_constant_value**2
 
+        # Start timer
         sp_timer = TicTocTimer()
         sp_timer.tic(msg=None)
 
@@ -240,18 +241,17 @@ def stochastic_program(
         m, analysis_square = self._compute_stochastic_program(m, optimize_opt)
 
         if self.optimize:
-            # set max_iter to 0 to debug the initialization
-            # self.solver.options["max_iter"] = 0
-            # self.solver.options["bound_push"] = 1e-10
-            
+            # If set to optimize, solve the optimization problem (with degrees of freedom)            
             analysis_optimize = self._optimize_stochastic_program(m)
             dT = sp_timer.toc(msg=None)
             self.logger.info("elapsed time: %0.1f seconds" % dT)
+            # Return both square problem and optimization problem results
             return analysis_square, analysis_optimize
 
         else:
             dT = sp_timer.toc(msg=None)
             self.logger.info("elapsed time: %0.1f seconds" % dT)
+            # Return only square problem results
             return analysis_square
 
     def _compute_stochastic_program(self, m, optimize_option):
@@ -596,9 +596,7 @@ def _create_block(self):
         mod.scenario = pyo.Set(initialize=self.scenario_data.scenario_indices)
         
         # Determine if create_model takes theta as an optional input
-        # print(inspect.getfullargspec(self.create_model))
         pass_theta_to_initialize= ('theta' in inspect.getfullargspec(self.create_model).args)
-        #print("pass_theta_to_initialize =", pass_theta_to_initialize)
 
         # Allow user to self-define complex design variables
         self.create_model(mod=mod, model_option=ModelOptionLib.stage1)
@@ -617,10 +615,12 @@ def block_build(b, s):
             # idea: check if create_model takes theta as an optional input, if so, pass parameter values to create_model
 
             if pass_theta_to_initialize:
+                # Grab the values of theta for this scenario/block
                 theta_initialize = self.scenario_data.scenario[s]
-                #print("Initializing with theta=", theta_initialize)
+                # Add model on block with theta values
                 self.create_model(mod=b, model_option=ModelOptionLib.stage2, theta=theta_initialize)
             else:
+                # Otherwise add model on block without theta values
                 self.create_model(mod=b, model_option=ModelOptionLib.stage2)
 
             # fix parameter values to perturbed values
@@ -934,8 +934,6 @@ def initialize_jac(m, i, j):
             for i, bu in enumerate(model.regression_parameters):
                 for j, un in enumerate(model.regression_parameters):
                     dict_fim_initialize[(bu, un)] = self.fim_initial[i][j]
-            
-            #print(dict_fim_initialize)
 
         def initialize_fim(m, j, d):
             return dict_fim_initialize[(j, d)]
@@ -1226,19 +1224,6 @@ def _solve_doe(self, m, fix=False, opt_option=None):
             m, self.design_values, fix_opt=fix, optimize_option=opt_option
         )
 
-        '''
-        # This is for initialization diagnostics
-        # Remove before merging the PR
-        if not fix:
-            # halt at initial point
-            self.solver.options['max_iter'] = 0
-            self.solver.options['bound_push'] = 1E-10
-        else:
-            # resort to defaults
-            self.solver.options['max_iter'] = 3000
-            self.solver.options['bound_push'] = 0.01
-        '''
-
         # if user gives solver, use this solver. if not, use default IPOPT solver
         solver_result = self.solver.solve(mod, tee=self.tee_opt)
             
diff --git a/pyomo/contrib/doe/measurements.py b/pyomo/contrib/doe/measurements.py
index ae5b3519498..dcaac1f14fd 100644
--- a/pyomo/contrib/doe/measurements.py
+++ b/pyomo/contrib/doe/measurements.py
@@ -27,7 +27,7 @@
 
 import itertools
 import collections.abc
-
+from pyomo.common.numeric_types import native_numeric_types
 
 class VariablesWithIndices:
     def __init__(self):
@@ -100,12 +100,12 @@ def add_variables(
 
         # if a scalar (int or float) is given, set it as the lower bound for all variables
         if lower_bounds is not None:
-            if type(lower_bounds) in [int, float]:
+            if type(lower_bounds) in native_numeric_types:
                 lower_bounds = [lower_bounds] * len(added_names)
             self.lower_bounds.update(zip(added_names, lower_bounds))
 
         if upper_bounds is not None:
-            if type(upper_bounds) in [int, float]:
+            if type(upper_bounds) in native_numeric_types:
                 upper_bounds = [upper_bounds] * len(added_names)
             self.upper_bounds.update(zip(added_names, upper_bounds))
 

From 845fc3fdcae84180cd7f777887a771f55074e2bf Mon Sep 17 00:00:00 2001
From: Alex Dowling <alex@dowlinglab.org>
Date: Thu, 25 Apr 2024 21:58:33 -0400
Subject: [PATCH 07/11] Ran black

---
 pyomo/contrib/doe/doe.py          | 27 +++++++++++++++++----------
 pyomo/contrib/doe/measurements.py | 13 +++++++------
 pyomo/contrib/doe/result.py       |  2 +-
 3 files changed, 25 insertions(+), 17 deletions(-)

diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py
index d10fc1b7fcc..ed5e81027dd 100644
--- a/pyomo/contrib/doe/doe.py
+++ b/pyomo/contrib/doe/doe.py
@@ -42,6 +42,7 @@
 
 import inspect
 
+
 class CalculationMode(Enum):
     sequential_finite = "sequential_finite"
     direct_kaug = "direct_kaug"
@@ -228,7 +229,7 @@ def stochastic_program(
 
         # calculate how much the FIM element is scaled by a constant number
         # FIM = Jacobian.T@Jacobian, the FIM is scaled by squared value the Jacobian is scaled
-        self.fim_scale_constant_value = self.scale_constant_value**2
+        self.fim_scale_constant_value = self.scale_constant_value ** 2
 
         # Start timer
         sp_timer = TicTocTimer()
@@ -241,7 +242,7 @@ def stochastic_program(
         m, analysis_square = self._compute_stochastic_program(m, optimize_opt)
 
         if self.optimize:
-            # If set to optimize, solve the optimization problem (with degrees of freedom)            
+            # If set to optimize, solve the optimization problem (with degrees of freedom)
             analysis_optimize = self._optimize_stochastic_program(m)
             dT = sp_timer.toc(msg=None)
             self.logger.info("elapsed time: %0.1f seconds" % dT)
@@ -382,7 +383,7 @@ def compute_FIM(
 
         # calculate how much the FIM element is scaled by a constant number
         # As FIM~Jacobian.T@Jacobian, FIM is scaled twice the number the Q is scaled
-        self.fim_scale_constant_value = self.scale_constant_value**2
+        self.fim_scale_constant_value = self.scale_constant_value ** 2
 
         square_timer = TicTocTimer()
         square_timer.tic(msg=None)
@@ -594,9 +595,11 @@ def _create_block(self):
 
         # Set for block/scenarios
         mod.scenario = pyo.Set(initialize=self.scenario_data.scenario_indices)
-        
+
         # Determine if create_model takes theta as an optional input
-        pass_theta_to_initialize= ('theta' in inspect.getfullargspec(self.create_model).args)
+        pass_theta_to_initialize = (
+            'theta' in inspect.getfullargspec(self.create_model).args
+        )
 
         # Allow user to self-define complex design variables
         self.create_model(mod=mod, model_option=ModelOptionLib.stage1)
@@ -618,7 +621,9 @@ def block_build(b, s):
                 # Grab the values of theta for this scenario/block
                 theta_initialize = self.scenario_data.scenario[s]
                 # Add model on block with theta values
-                self.create_model(mod=b, model_option=ModelOptionLib.stage2, theta=theta_initialize)
+                self.create_model(
+                    mod=b, model_option=ModelOptionLib.stage2, theta=theta_initialize
+                )
             else:
                 # Otherwise add model on block without theta values
                 self.create_model(mod=b, model_option=ModelOptionLib.stage2)
@@ -773,7 +778,7 @@ def run_grid_search(
         self.store_optimality_as_csv = store_optimality_as_csv
 
         # calculate how much the FIM element is scaled
-        self.fim_scale_constant_value = scale_constant_value**2
+        self.fim_scale_constant_value = scale_constant_value ** 2
 
         # to store all FIM results
         result_combine = {}
@@ -926,7 +931,9 @@ def initialize_jac(m, i, j):
                 return 0.1
 
         model.sensitivity_jacobian = pyo.Var(
-            model.regression_parameters, model.measured_variables, initialize=initialize_jac
+            model.regression_parameters,
+            model.measured_variables,
+            initialize=initialize_jac,
         )
 
         if self.fim_initial is not None:
@@ -1071,7 +1078,7 @@ def _add_objective(self, m):
             eig = np.linalg.eigvals(fim)
 
             # If the smallest eigenvalue is (pratcially) negative, add a diagonal matrix to make it positive definite
-            small_number = 1E-10
+            small_number = 1e-10
             if min(eig) < small_number:
                 fim = fim + np.eye(len(self.param)) * (small_number - min(eig))
 
@@ -1226,7 +1233,7 @@ def _solve_doe(self, m, fix=False, opt_option=None):
 
         # if user gives solver, use this solver. if not, use default IPOPT solver
         solver_result = self.solver.solve(mod, tee=self.tee_opt)
-            
+
         return solver_result
 
     def _sgn(self, p):
diff --git a/pyomo/contrib/doe/measurements.py b/pyomo/contrib/doe/measurements.py
index dcaac1f14fd..fd3962f7888 100644
--- a/pyomo/contrib/doe/measurements.py
+++ b/pyomo/contrib/doe/measurements.py
@@ -29,6 +29,7 @@
 import collections.abc
 from pyomo.common.numeric_types import native_numeric_types
 
+
 class VariablesWithIndices:
     def __init__(self):
         """This class provides utility methods for DesignVariables and MeasurementVariables to create
@@ -182,16 +183,16 @@ def _check_valid_input(
             raise ValueError("Values is of different length with indices.")
 
         if (
-            lower_bounds is not None # ensure not None
-            and isinstance(lower_bounds, collections.abc.Sequence) # ensure list-like
-            and len(lower_bounds) != len_indices # ensure same length
+            lower_bounds is not None  # ensure not None
+            and isinstance(lower_bounds, collections.abc.Sequence)  # ensure list-like
+            and len(lower_bounds) != len_indices  # ensure same length
         ):
             raise ValueError("Lowerbounds is of different length with indices.")
 
         if (
-            upper_bounds is not None # ensure None
-            and isinstance(upper_bounds, collections.abc.Sequence) # ensure list-like
-            and len(upper_bounds) != len_indices # ensure same length
+            upper_bounds is not None  # ensure None
+            and isinstance(upper_bounds, collections.abc.Sequence)  # ensure list-like
+            and len(upper_bounds) != len_indices  # ensure same length
         ):
             raise ValueError("Upperbounds is of different length with indices.")
 
diff --git a/pyomo/contrib/doe/result.py b/pyomo/contrib/doe/result.py
index 1593214c30a..8f98e74f159 100644
--- a/pyomo/contrib/doe/result.py
+++ b/pyomo/contrib/doe/result.py
@@ -81,7 +81,7 @@ def __init__(
         self.prior_FIM = prior_FIM
         self.store_FIM = store_FIM
         self.scale_constant_value = scale_constant_value
-        self.fim_scale_constant_value = scale_constant_value**2
+        self.fim_scale_constant_value = scale_constant_value ** 2
         self.max_condition_number = max_condition_number
         self.logger = logging.getLogger(__name__)
         self.logger.setLevel(level=logging.WARN)

From 5c7cab5bb4629f2dc5ed6dde08d006c78b3735bf Mon Sep 17 00:00:00 2001
From: Alex Dowling <alex@dowlinglab.org>
Date: Thu, 25 Apr 2024 21:59:45 -0400
Subject: [PATCH 08/11] Reran black.

---
 pyomo/contrib/doe/doe.py    | 6 +++---
 pyomo/contrib/doe/result.py | 2 +-
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py
index ed5e81027dd..0d2c7c2f982 100644
--- a/pyomo/contrib/doe/doe.py
+++ b/pyomo/contrib/doe/doe.py
@@ -229,7 +229,7 @@ def stochastic_program(
 
         # calculate how much the FIM element is scaled by a constant number
         # FIM = Jacobian.T@Jacobian, the FIM is scaled by squared value the Jacobian is scaled
-        self.fim_scale_constant_value = self.scale_constant_value ** 2
+        self.fim_scale_constant_value = self.scale_constant_value**2
 
         # Start timer
         sp_timer = TicTocTimer()
@@ -383,7 +383,7 @@ def compute_FIM(
 
         # calculate how much the FIM element is scaled by a constant number
         # As FIM~Jacobian.T@Jacobian, FIM is scaled twice the number the Q is scaled
-        self.fim_scale_constant_value = self.scale_constant_value ** 2
+        self.fim_scale_constant_value = self.scale_constant_value**2
 
         square_timer = TicTocTimer()
         square_timer.tic(msg=None)
@@ -778,7 +778,7 @@ def run_grid_search(
         self.store_optimality_as_csv = store_optimality_as_csv
 
         # calculate how much the FIM element is scaled
-        self.fim_scale_constant_value = scale_constant_value ** 2
+        self.fim_scale_constant_value = scale_constant_value**2
 
         # to store all FIM results
         result_combine = {}
diff --git a/pyomo/contrib/doe/result.py b/pyomo/contrib/doe/result.py
index 8f98e74f159..1593214c30a 100644
--- a/pyomo/contrib/doe/result.py
+++ b/pyomo/contrib/doe/result.py
@@ -81,7 +81,7 @@ def __init__(
         self.prior_FIM = prior_FIM
         self.store_FIM = store_FIM
         self.scale_constant_value = scale_constant_value
-        self.fim_scale_constant_value = scale_constant_value ** 2
+        self.fim_scale_constant_value = scale_constant_value**2
         self.max_condition_number = max_condition_number
         self.logger = logging.getLogger(__name__)
         self.logger.setLevel(level=logging.WARN)

From a8a5450fc5021b1c30daafcdb282f0e936a8d9ad Mon Sep 17 00:00:00 2001
From: Alex Dowling <alex@dowlinglab.org>
Date: Thu, 25 Apr 2024 22:06:38 -0400
Subject: [PATCH 09/11] Added more comments.

---
 pyomo/contrib/doe/doe.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py
index 0d2c7c2f982..28dad1f20c7 100644
--- a/pyomo/contrib/doe/doe.py
+++ b/pyomo/contrib/doe/doe.py
@@ -913,6 +913,9 @@ def identity_matrix(m, i, j):
             else:
                 return 0
 
+        ### Initialize the Jacobian if provided by the user
+        
+        # If the user provides an initial Jacobian, convert it to a dictionary
         if self.jac_initial is not None:
             dict_jac_initialize = {}
             for i, bu in enumerate(model.regression_parameters):
@@ -924,9 +927,12 @@ def identity_matrix(m, i, j):
                         # Jacobian is a numpy array, rows are regression parameters, columns are measured variables
                         dict_jac_initialize[(bu, un)] = self.jac_initial[i][j]
 
+        # Initialize the Jacobian matrix
         def initialize_jac(m, i, j):
+            # If provided by the user, use the values now stored in the dictionary
             if self.jac_initial is not None:
                 return dict_jac_initialize[(i, j)]
+            # Otherwise initialize to 0.1 (which is an arbitrary non-zero value)
             else:
                 return 0.1
 

From 1bac09e842eaaeeeaef8191e30e500f2e311fefb Mon Sep 17 00:00:00 2001
From: Alex Dowling <alex@dowlinglab.org>
Date: Thu, 25 Apr 2024 22:21:24 -0400
Subject: [PATCH 10/11] Reran black

---
 pyomo/contrib/doe/doe.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py
index 28dad1f20c7..ab9a5ad9f85 100644
--- a/pyomo/contrib/doe/doe.py
+++ b/pyomo/contrib/doe/doe.py
@@ -914,7 +914,7 @@ def identity_matrix(m, i, j):
                 return 0
 
         ### Initialize the Jacobian if provided by the user
-        
+
         # If the user provides an initial Jacobian, convert it to a dictionary
         if self.jac_initial is not None:
             dict_jac_initialize = {}

From 5439e9560b1d4b8602350d8de3c9006e1c5fd615 Mon Sep 17 00:00:00 2001
From: Bethany Nicholson <blnicho@users.noreply.github.com>
Date: Mon, 6 May 2024 15:34:15 -0600
Subject: [PATCH 11/11] Fix typo in doe.py

---
 pyomo/contrib/doe/doe.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py
index ab9a5ad9f85..0fc3e8770fe 100644
--- a/pyomo/contrib/doe/doe.py
+++ b/pyomo/contrib/doe/doe.py
@@ -1083,7 +1083,7 @@ def _add_objective(self, m):
             # Calculate the eigenvalues of the FIM matrix
             eig = np.linalg.eigvals(fim)
 
-            # If the smallest eigenvalue is (pratcially) negative, add a diagonal matrix to make it positive definite
+            # If the smallest eigenvalue is (practically) negative, add a diagonal matrix to make it positive definite
             small_number = 1e-10
             if min(eig) < small_number:
                 fim = fim + np.eye(len(self.param)) * (small_number - min(eig))