Skip to content

Commit e9e83b6

Browse files
committed
to test doe , added the "doe_test_example.py" script, which returns a doe object that uses the "reactor_experiment" example in doe directory. Also added the png files to gitignore that are stored when "reactor_example.py" is run.
1 parent 39f26e6 commit e9e83b6

File tree

3 files changed

+167
-2
lines changed

3 files changed

+167
-2
lines changed

.gitignore

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -34,3 +34,9 @@ cplex.log
3434

3535
# Mac tracking files
3636
*.DS_Store*
37+
38+
# DOE example .png files
39+
example_reactor_compute_FIM_E_opt.png
40+
example_reactor_compute_FIM_A_opt.png
41+
example_reactor_compute_FIM_D_opt.png
42+
example_reactor_compute_FIM_ME_opt.png
Lines changed: 150 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,150 @@
1+
# ___________________________________________________________________________
2+
#
3+
# Pyomo: Python Optimization Modeling Objects
4+
# Copyright (c) 2008-2025
5+
# National Technology and Engineering Solutions of Sandia, LLC
6+
# Under the terms of Contract DE-NA0003525 with National Technology and
7+
# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
8+
# rights in this software.
9+
# This software is distributed under the 3-clause BSD License.
10+
# ___________________________________________________________________________
11+
from pyomo.common.dependencies import numpy as np
12+
13+
from pyomo.contrib.doe.examples.reactor_experiment import ReactorExperiment
14+
from pyomo.contrib.doe import DesignOfExperiments
15+
16+
import pyomo.environ as pyo
17+
18+
import json
19+
from pathlib import Path
20+
import idaes.core.solvers.get_solver # import the idaes solvers
21+
22+
23+
# Example for sensitivity analysis on the reactor experiment
24+
# After sensitivity analysis is done, we perform optimal DoE
25+
def run_reactor_doe(
26+
run_doe=False,
27+
compute_FIM_full_factorial=False,
28+
draw_factorial_figure=False,
29+
save_draw_factorial_figure=False,
30+
):
31+
# Read in file
32+
DATA_DIR = Path(__file__).parent.parent / "examples"
33+
file_path = DATA_DIR / "result.json"
34+
35+
with open(file_path) as f:
36+
data_ex = json.load(f)
37+
38+
# Put temperature control time points into correct format for reactor experiment
39+
data_ex["control_points"] = {
40+
float(k): v for k, v in data_ex["control_points"].items()
41+
}
42+
43+
# Create a ReactorExperiment object; data and discretization information are part
44+
# of the constructor of this object
45+
experiment = ReactorExperiment(data=data_ex, nfe=10, ncp=3)
46+
47+
# Use a central difference, with step size 1e-3
48+
fd_formula = "central"
49+
step_size = 1e-3
50+
51+
# Use the determinant objective with scaled sensitivity matrix
52+
objective_option = "determinant"
53+
scale_nominal_param_value = True
54+
55+
# Create the DesignOfExperiments object
56+
# We will not be passing any prior information in this example
57+
# and allow the experiment object and the DesignOfExperiments
58+
# call of ``run_doe`` perform model initialization.
59+
doe_obj = DesignOfExperiments(
60+
experiment,
61+
fd_formula=fd_formula,
62+
step=step_size,
63+
objective_option=objective_option,
64+
scale_constant_value=1,
65+
scale_nominal_param_value=scale_nominal_param_value,
66+
prior_FIM=None,
67+
jac_initial=None,
68+
fim_initial=None,
69+
L_diagonal_lower_bound=1e-7,
70+
solver=None,
71+
tee=False,
72+
get_labeled_model_args=None,
73+
_Cholesky_option=True,
74+
_only_compute_fim_lower=True,
75+
)
76+
77+
# Make design ranges to compute the full factorial design
78+
design_ranges = {"CA[0]": [1, 5, 9], "T[0]": [300, 700, 9]}
79+
80+
# Compute the full factorial design with the sequential FIM calculation
81+
if compute_FIM_full_factorial:
82+
doe_obj.compute_FIM_full_factorial(
83+
design_ranges=design_ranges, method="sequential"
84+
)
85+
86+
# Plot the results
87+
if draw_factorial_figure:
88+
if save_draw_factorial_figure:
89+
figure_file_name = "example_reactor_compute_FIM"
90+
else:
91+
figure_file_name = None
92+
93+
doe_obj.draw_factorial_figure(
94+
sensitivity_design_variables=["CA[0]", "T[0]"],
95+
fixed_design_variables={
96+
"T[0.125]": 300,
97+
"T[0.25]": 300,
98+
"T[0.375]": 300,
99+
"T[0.5]": 300,
100+
"T[0.625]": 300,
101+
"T[0.75]": 300,
102+
"T[0.875]": 300,
103+
"T[1]": 300,
104+
},
105+
title_text="Reactor Example",
106+
xlabel_text="Concentration of A (M)",
107+
ylabel_text="Initial Temperature (K)",
108+
figure_file_name=figure_file_name,
109+
log_scale=False,
110+
)
111+
112+
###########################
113+
# End sensitivity analysis
114+
115+
# Begin optimal DoE
116+
####################
117+
if run_doe:
118+
doe_obj.run_doe()
119+
120+
# Print out a results summary
121+
print("Optimal experiment values: ")
122+
print(
123+
"\tInitial concentration: {:.2f}".format(
124+
doe_obj.results["Experiment Design"][0]
125+
)
126+
)
127+
print(
128+
("\tTemperature values: [" + "{:.2f}, " * 8 + "{:.2f}]").format(
129+
*doe_obj.results["Experiment Design"][1:]
130+
)
131+
)
132+
print("FIM at optimal design:\n {}".format(np.array(doe_obj.results["FIM"])))
133+
print(
134+
"Objective value at optimal design: {:.2f}".format(
135+
pyo.value(doe_obj.model.objective)
136+
)
137+
)
138+
139+
print(doe_obj.results["Experiment Design Names"])
140+
141+
###################
142+
# End optimal DoE
143+
144+
# retrun the doe object for further use
145+
return doe_obj
146+
147+
148+
# uncomment the following line to run the example
149+
# if __name__ == "__main__":
150+
# run_reactor_doe(run_doe=True)

pyomo/contrib/doe/tests/test_doe_FIM_metrics.py

Lines changed: 11 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -7,12 +7,21 @@
77

88
import pyomo.common.unittest as unittest
99
from pyomo.contrib.doe import DesignOfExperiments
10-
from pyomo.contrib.doe.doe import compute_FIM_full_factorial, _SMALL_TOLERANCE_IMG
10+
from pyomo.contrib.doe.tests import doe_test_example
11+
from pyomo.contrib.doe.doe import _SMALL_TOLERANCE_IMG
1112
from pyomo.opt import SolverFactory
13+
import pyomo.environ as pyo
14+
15+
import json
16+
from pathlib import Path
17+
import idaes
18+
19+
ipopt_available = SolverFactory("ipopt").available()
20+
k_aug_available = SolverFactory("k_aug", solver_io="nl", validate=False)
1221

1322

1423
@unittest.skipIf(not numpy_available, "Numpy is not available")
15-
class TestFullFactrialMetrics(unittest.TestCase):
24+
class TestFullFactorialMetrics(unittest.TestCase):
1625
def test_compute_FIM_full_factorial_metrics(self):
1726
# Create a sample Fisher Information Matrix (FIM)
1827
FIM = np.array([[4, 2], [2, 3]])

0 commit comments

Comments
 (0)