Pyomo.DoE: Correct A optimality #3803

smondal13 · 2025-12-16T14:08:17Z

Fixes

Summary/Motivation:

In the old definition of A-optimality (we will call it pseudo A-optimality from now on), DoE computed the trace of the Fisher Information Matrix (FIM), as few studies used that definition. However, that is not correct, and new literature mentions the trace of the inverse of the FIM as A-optimality. This PR has adopted this new definition of A-optimality in DoE.

Changes proposed in this PR:

Change the old A-optimality to pseudo A-optimality
Implement a new definition of A-optimality (=trace(inverse of FIM))

Legal Acknowledgement

By contributing to this software project, I have read the contribution guide and agree to the following terms and conditions for my contribution:

I agree my contributions are submitted under the BSD license.
I represent I am authorized to make the contributions and grant the license. If my employer has rights to intellectual property that includes these contributions, I represent that I have received permission to make contributions and grant the required license on behalf of that employer.

…t constraints

…ample did not solve to the desired point using A-optimality. However, when initialized to small value `1e-5`, it solved to the optimal point.

…scheme, the MM model solved without initialization from outside.

…`doe/tests`` files to reflect the changes.

…n `doe/test_utils.py`

…erstanding in ``doe/test_utils.py``

…and `_Cholesky_option = False`. Add flag for functions are the only used in `trace` calculation.

smondal13 · 2025-12-16T17:10:35Z

@adowling2 @blnicho This PR is ready for initial review.

pyomo/contrib/doe/doe.py

pyomo/contrib/doe/utils.py

pyomo/contrib/doe/tests/test_doe_build.py

pyomo/contrib/doe/tests/test_doe_errors.py

Replace `model` with `m` in `cholesky_LLinv_imp()` Co-authored-by: Bethany Nicholson <[email protected]>

blnicho

@smondal13 thanks for addressing my initial review comments. I think these changes look reasonable so far but would like to see tests for the corrected A-optimality before I'll approve the PR.

…ality

… and add `doe_example..py` based on the rooney_biegler.py

…g function.

…gler/doe_example.py`

…f glob

…oe/examples/rooney_biegler_doe_example.py`

adowling2

Started to leave review comments. I will add more comments soon.

adowling2 · 2026-01-20T18:49:52Z

pyomo/contrib/doe/doe.py

@@ -64,6 +64,7 @@
 class ObjectiveLib(Enum):
    determinant = "determinant"


Let's add the following comments after each line

# det(FIM), D-optimality # trace(inv(FIM)), A-optimalty # trace(FIM), "psuedo"-A-optimality # min(eig(FIM)), E-optimality # cond(FIM), ME-optimality # constant zero, useful for initialization and debugging

Added comments

adowling2

From a computational perspective, pseudo-A-optimality, i.e., trace(FIM), is the easiest objective to solve. I suggest we do a majority of the testing of permutations of DoE options with pseudo-A-optimality. I suggest that we then do a fewer number of tests with D-optimality (both Cholseky and GreyBox) and A-optimality (both Cholseky and GreyBox).

pyomo/contrib/doe/examples/rooney_biegler_doe_example.py

blnicho · 2026-01-21T06:22:18Z

pyomo/contrib/doe/examples/rooney_biegler_doe_example.py

+    rooney_biegler_doe = DesignOfExperiments(
+        experiment=rooney_biegler_experiment,
+        objective_option="determinant",
+        tee=tee,
+        prior_FIM=FIM,
+        improve_cholesky_roundoff_error=True,
+    )
+
+    if optimize_experiment_D:
+        # D-Optimality
+        rooney_biegler_doe_D = DesignOfExperiments(
+            experiment=rooney_biegler_experiment,
+            objective_option="determinant",
+            tee=tee,
+            prior_FIM=FIM,
+            improve_cholesky_roundoff_error=True,
+        )
+        rooney_biegler_doe_D.run_doe()


Why do you need both rooney_biegler_doe and rooney_biegler_doe_D? Seems like they are identical instances of the DesignOfExperiments class. Why not just call run_doe directly on rooney_biegler_doe when optimize_experiment_D is True?

That was repeated. It's corrected now.

blnicho · 2026-01-21T06:22:51Z

pyomo/contrib/doe/examples/rooney_biegler_doe_example.py

+        opt_D = results_container["optimization"]["D"]
+        opt_A = results_container["optimization"]["A"]


This will fail if either of optimize_experiment_A or optimize_experiment_D are set to False. You should either guard this code against different options sent to the run_rooney_biegler_doe function or remove the options and have this example always run DoE D/A optimization.

I have removed the plotting function. It is to visually confirm the result. It is not required.

blnicho · 2026-01-21T06:27:29Z

pyomo/contrib/doe/tests/test_doe_solve.py

+        A_opt_value = run_rooney_biegler_doe(optimize_experiment_A=True)[
+            "optimization"
+        ]["A"]["value"]
+        A_opt_design_value = run_rooney_biegler_doe(optimize_experiment_A=True)[
+            "optimization"
+        ]["A"]["design"][0]


Calling run_rooney_biegler_doe twice here seems inefficient. Why not save the return value in a local variable? Then you can extract any of the values in the saved results

Yeah, it was inefficient. I have saved the return as a variable and reused it.

blnicho · 2026-01-21T06:28:27Z

pyomo/contrib/doe/tests/test_doe_solve.py

+    @unittest.skipUnless(pandas_available, "test requires pandas")
+    @unittest.skipUnless(matplotlib_available, "test requires matplotlib")


The run_rooney_biegler_doe function has an option to plot or not. So why not have this test change the value of plot_optimal_design depending on the value of matplotlib_available? Then you can still test most of the example even when matplotlib isn't available.

Yes, I have done that.

blnicho · 2026-01-21T06:30:17Z

pyomo/contrib/doe/tests/test_doe_solve.py

+        # Switch backend to 'Agg' to prevent GUI windows from popping up
+        plt.switch_backend('Agg')


Is this still needed? Didn't you set the backend to 'Agg` up on L28?

Not needed, I have removed that

pyomo/contrib/doe/tests/test_doe_solve.py

blnicho · 2026-01-21T06:43:37Z

pyomo/contrib/doe/tests/test_doe_errors.py

+        fd_method = "central"
+        obj_used = "trace"
+
+        experiment = run_rooney_biegler_doe()["experiment"]


Could you rework this to use FullReactorExperiment like the other tests in this file?

Since we talked about using inexpensive tests during the IDAES summit in Pittsburgh, I have added simpler rooney_biegler example here for the test instead of FullReactorExperiment, which is more complex and sometimes fails to converge. Eventually, we want to make all tests inexpensive (or comparatively inexpensive) for Pyomo.DoE. However, it does not mean all tests that use FullReactorExperiment are expensive. I can add a test using FullReactorExperiment. Let me know whether I should add FullReactorExperiment.

Can we do a lightweight version of FullReactorExperiment, perhaps with pseudo A opt

@adowling2 @blnicho I have replaced all the tests in test_doe_errors.py with rooney_biegler experiment in #3828 . So, we are good here, right?

Note: This branch is merged with #3828

blnicho · 2026-01-21T07:03:45Z

pyomo/contrib/doe/tests/test_doe_build.py

+        params = list(model.parameter_names)
+
+        # Check cov_trace initialization
+        cov_trace_from_fim_inv = sum(pyo.value(model.fim_inv[j, j]) for j in params)


The expected value in a test should be a hard-coded number not a repeat of the implementation.

A hard-coded value is added.

blnicho · 2026-01-21T07:09:29Z

pyomo/contrib/doe/tests/test_doe_build.py

+                val = sum(
+                    pyo.value(model.L[c, params[k]])
+                    * pyo.value(model.L_inv[params[k], d])
+                    for k in range(len(params))
+                )


Suggested change

val = sum(

pyo.value(model.L[c, params[k]])

* pyo.value(model.L_inv[params[k], d])

for k in range(len(params))

)

val = pyo.value(sum(model.L[c, params[k]] * model.L_inv[params[k], d]

for k in range(len(params))

))

I'm pretty sure you can simplify this slightly by just calling pyo.value once. (Note: if you accept this change through the GitHub GUI the line indentation is off)

Added manually.

… and remove redundant/bad code.

mrmundt · 2026-01-21T18:59:18Z

@smondal13 - Gurobi just released a new version that supports Python 3.14, and we are seeing failures for our gurobi tests because of it. I cancelled your jobs because I have a fix already, just need to get it in.

TL;DR - if you see 3.14 failures in gurobi, not your fault, we're on it.

… change-a-optimality

adowling2

I have some concerns about how this testing framework can scale to all of the variants of Pyomo.DoE that we want to try, including:

central finite difference versus symbolic differentiation
GreyBox: A, D, E, ME optimality
Cholesky: A, D

We do not need to add all of these variants in this PR. But I think we need to ensure this PR sets us up to easily add these extensions in other PRs. This is important because I could see multiple parallel PRs to add these features that will all be building on this testing foundation.

pyomo/contrib/doe/examples/reactor_example.py

adowling2 · 2026-01-30T15:31:13Z

pyomo/contrib/doe/examples/rooney_biegler_doe_example.py

+
+    Parameters
+    ----------
+    optimize_experiment_A : bool, optional


This does not seem very general... instead, do we want to support A, D, E, ME, pseudo A, ... objectives?

Do you mean that we pass an argument, e.g., objective="trace" / "pesudo_trace"... etc in run_rooney_biegler_doe()` ? I think that is a very good idea.

Yes, unless there are significant benefits to doing it as a single function (sequentially). If that is the case, you could make objectives the input and support either a single string or a list of strings.

I have added the objective and the gradient/sensitivity computation option. Maybe symbolic can be added there when symbolic is merged. And solving the example for only one type of objective is fine.

def run_rooney_biegler_doe( sensitivity_formula="central", optimization_objective="determinant", improve_cholesky_roundoff_error=False, compute_FIM_full_factorial=False, draw_factorial_figure=False, design_range={'hour': [0, 10, 40]}, tee=False, print_output=False, ):

adowling2 · 2026-01-30T15:33:35Z

pyomo/contrib/doe/tests/test_doe_errors.py

+        fd_method = "central"
+        obj_used = "trace"
+
+        experiment = run_rooney_biegler_doe()["experiment"]


Can we do a lightweight version of FullReactorExperiment, perhaps with pseudo A opt

…sitivity_formula parameters; update tests accordingly

smondal13 added 11 commits November 10, 2025 20:14

changed the definition of A-optimality in doe/utils.py

6a75654

used the correct definition of A-optimilty in both utils.py and doe.py

4ba1813

changed the label of the colorbar to trace(cov)

5c7462f

Merge branch 'Pyomo:main' into change-a-optimality

e3d4833

Adds A-optimlaity objective and constraints with round off improvemen…

4a50d2a

…t constraints

Add roundoff correction constraint. Without initialization, the MM ex…

b5a690b

…ample did not solve to the desired point using A-optimality. However, when initialized to small value `1e-5`, it solved to the optimal point.

Add initialization scheme for L_inv. After adding the initialization …

106781b

…scheme, the MM model solved without initialization from outside.

Add trace_FIM in `compute_FIM_full_factorial() and change related `…

91b0daf

…`doe/tests`` files to reflect the changes.

Update failing tests in and refactor the test with helper functions i…

682bf94

…n `doe/test_utils.py`

Reorder the helper functions' sequence and add comment for better und…

efe4505

…erstanding in ``doe/test_utils.py``

Raise ValueError in the edgecase when objective_option = 'trace' …

31347a9

…and `_Cholesky_option = False`. Add flag for functions are the only used in `trace` calculation.

smondal13 added this to ParmEst & Pyomo.DoE Development Dec 16, 2025

smondal13 moved this to Ready for design review in ParmEst & Pyomo.DoE Development Dec 16, 2025

blnicho requested review from adowling2 and blnicho December 16, 2025 19:40

Add pseudo A-optimality heatmap in the doe sensitivity plotting function

78f0a78

blnicho requested changes Dec 18, 2025

View reviewed changes

smondal13 and others added 2 commits December 18, 2025 15:12

Apply suggestions from code review

241d52c

Replace `model` with `m` in `cholesky_LLinv_imp()` Co-authored-by: Bethany Nicholson <[email protected]>

Add both trace_FIM (pseudo_A_optimality) and trace_cov(A-optimality)

2eba303

blnicho reviewed Dec 18, 2025

View reviewed changes

blnicho self-requested a review December 18, 2025 21:24

blnicho added the pyomo.DoE label Dec 18, 2025

smondal13 added 7 commits December 18, 2025 19:34

Merge branch 'parmest-consolidate-rooney-biegler' into change-a-optim…

51bd8c8

…ality

Add self.model is None in rooney_biegler.py get_labeled_model()…

76baff0

… and add `doe_example..py` based on the rooney_biegler.py

Add A-optimality test

cf2a7d2

Add tests for A-optimality in test_doe_build and test_doe_errors.py

50a81a6

Add test files for the objective trace and the draw_factorial plottin…

069aa50

…g function.

Activate timprove_cholesky_roundoff_error=True in the `rooney_bie…

f671418

…gler/doe_example.py`

Add attempt_import for matplotlib in test_doe_solve.py

785a449

smondal13 and others added 2 commits January 13, 2026 11:40

Add glob import at the beginning and removes inside function import o…

705d83a

…f glob

Merge branch 'main' into change-a-optimality

4e11069

smondal13 moved this from Ready for design review to Ready for final review in ParmEst & Pyomo.DoE Development Jan 15, 2026

smondal13 added 5 commits January 19, 2026 20:12

Move and renameparmest/examples/rooney_biegler/doe_example.py to `d…

352bb71

…oe/examples/rooney_biegler_doe_example.py`

Merge remote-tracking branch 'upstream/main' into change-a-optimality

9d1fece

Add correct import statements for rooney_biegler_doe_example

b9400a4

Move test forrooney_biegler_doe_example.py to test_doe_solve.py

240aba6

Correct name of a test and remove localmatplotlib.use("Agg") .

f445503

adowling2 reviewed Jan 20, 2026

View reviewed changes

adowling2 suggested changes Jan 20, 2026

View reviewed changes

pyomo/contrib/doe/examples/rooney_biegler_doe_example.py Show resolved Hide resolved

pyomo/contrib/doe/examples/rooney_biegler_doe_example.py Show resolved Hide resolved

blnicho added this to Pyomo 6.10 Jan 21, 2026

github-project-automation bot moved this to Todo in Pyomo 6.10 Jan 21, 2026

blnicho moved this from Todo to Review In Progress in Pyomo 6.10 Jan 21, 2026

blnicho requested changes Jan 21, 2026

View reviewed changes

smondal13 added 2 commits January 21, 2026 10:47

Address comments from Alex and Bethany. Add more tests, docstring…

0b9ce22

… and remove redundant/bad code.

Merge branch 'main' into change-a-optimality

b086222

smondal13 added 6 commits January 21, 2026 17:05

Merge branch 'main' into change-a-optimality

d376b72

Add minor docstring formatting change

bec7228

Merge branch 'change-a-optimality' of github.com:smondal13/pyomo into…

16392d6

… change-a-optimality

Merge branch 'main' into change-a-optimality

91da06f

Correct 1D plotting function y_label text

b036d81

Correct y_label in 1D plotting function

9e036da

mrmundt requested review from adowling2 and blnicho January 27, 2026 19:42

smondal13 mentioned this pull request Jan 28, 2026

[Pyomo.DoE] Replace computationally expensive test example with Rooney Biegler #3828

Open

adowling2 suggested changes Jan 30, 2026

View reviewed changes

smondal13 added 3 commits January 30, 2026 11:58

Remove save_plots argument from doe/reactor_exmple

89b9da5

Refactor run_rooney_biegler_doe to use optimization_objective and sen…

a7653b6

…sitivity_formula parameters; update tests accordingly

Merge branch 'main' into change-a-optimality

571780e

		@@ -64,6 +64,7 @@
		class ObjectiveLib(Enum):
		determinant = "determinant"

		opt_D = results_container["optimization"]["D"]
		opt_A = results_container["optimization"]["A"]

		@unittest.skipUnless(pandas_available, "test requires pandas")
		@unittest.skipUnless(matplotlib_available, "test requires matplotlib")

		# Switch backend to 'Agg' to prevent GUI windows from popping up
		plt.switch_backend('Agg')

Pyomo.DoE: Correct A optimality #3803

Are you sure you want to change the base?

Pyomo.DoE: Correct A optimality #3803

Conversation

smondal13 commented Dec 16, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Fixes

Summary/Motivation:

Changes proposed in this PR:

Legal Acknowledgement

Uh oh!

smondal13 commented Dec 16, 2025

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

blnicho left a comment

Choose a reason for hiding this comment

Uh oh!

adowling2 left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

adowling2 left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

smondal13 Jan 21, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

smondal13 Jan 30, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mrmundt commented Jan 21, 2026

Uh oh!

smondal13 commented Dec 16, 2025 •

edited

Loading

smondal13 Jan 21, 2026 •

edited

Loading

smondal13 Jan 30, 2026 •

edited

Loading

smondal13 Jan 30, 2026 •

edited

Loading