|
34 | 34 | testdir = this_file_dir() |
35 | 35 |
|
36 | 36 |
|
37 | | -@unittest.skipIf( |
38 | | - not parmest.parmest_available, |
39 | | - "Cannot test parmest: required dependencies are missing", |
40 | | -) |
41 | | -@unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available") |
42 | | -class TestExceptionRooneyBiegler(unittest.TestCase): |
43 | | - |
44 | | - def setUp(self): |
45 | | - self.data = pd.DataFrame( |
46 | | - data=[[1, 8.3], [2, 10.3], [3, 19.0], [4, 16.0], [5, 15.6], [7, 19.8]], |
47 | | - columns=["hour", "y"], |
48 | | - ) |
49 | | - |
50 | | - # create the Rooney-Biegler model |
51 | | - def rooney_biegler_model(): |
52 | | - """ |
53 | | - Formulates the Pyomo model of the Rooney-Biegler example |
54 | | -
|
55 | | - Returns: |
56 | | - m: Pyomo model |
57 | | - """ |
58 | | - m = pyo.ConcreteModel() |
59 | | - |
60 | | - m.asymptote = pyo.Var(within=pyo.NonNegativeReals, initialize=10) |
61 | | - m.rate_constant = pyo.Var(within=pyo.NonNegativeReals, initialize=0.2) |
62 | | - |
63 | | - m.hour = pyo.Var(within=pyo.PositiveReals, initialize=0.1) |
64 | | - m.y = pyo.Var(within=pyo.NonNegativeReals) |
65 | | - |
66 | | - @m.Constraint() |
67 | | - def response_rule(m): |
68 | | - return m.y == m.asymptote * (1 - pyo.exp(-m.rate_constant * m.hour)) |
69 | | - |
70 | | - return m |
71 | | - |
72 | | - # create the Experiment class |
73 | | - class RooneyBieglerExperiment(Experiment): |
74 | | - def __init__(self, experiment_number, hour, y): |
75 | | - self.y = y |
76 | | - self.hour = hour |
77 | | - self.experiment_number = experiment_number |
78 | | - self.model = None |
79 | | - |
80 | | - def get_labeled_model(self): |
81 | | - self.create_model() |
82 | | - self.finalize_model() |
83 | | - self.label_model() |
84 | | - |
85 | | - return self.model |
86 | | - |
87 | | - def create_model(self): |
88 | | - m = self.model = rooney_biegler_model() |
89 | | - |
90 | | - return m |
91 | | - |
92 | | - def finalize_model(self): |
93 | | - m = self.model |
94 | | - |
95 | | - # fix the input variable |
96 | | - m.hour.fix(self.hour) |
97 | | - |
98 | | - return m |
99 | | - |
100 | | - def label_model(self): |
101 | | - m = self.model |
102 | | - |
103 | | - # add unknown parameters |
104 | | - m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL) |
105 | | - m.unknown_parameters.update( |
106 | | - (k, pyo.value(k)) for k in [m.asymptote, m.rate_constant] |
107 | | - ) |
108 | | - |
109 | | - return m |
110 | | - |
111 | | - # extract the input and output variables |
112 | | - hour_data = self.data["hour"] |
113 | | - y_data = self.data["y"] |
114 | | - |
115 | | - # create the experiments list |
116 | | - rooney_biegler_exp_list = [] |
117 | | - for i in range(self.data.shape[0]): |
118 | | - rooney_biegler_exp_list.append( |
119 | | - RooneyBieglerExperiment(i, hour_data[i], y_data[i]) |
120 | | - ) |
121 | | - |
122 | | - self.exp_list = rooney_biegler_exp_list |
123 | | - |
124 | | - def test_parmest_exception(self): |
125 | | - """ |
126 | | - Test the exception raised by parmest as a result of not defining |
127 | | - the "experiment_outputs" attribute |
128 | | - """ |
129 | | - with self.assertRaises(RuntimeError) as context: |
130 | | - parmest.Estimator(self.exp_list, obj_function="SSE", tee=True) |
131 | | - |
132 | | - self.assertIn("experiment_outputs", str(context.exception)) |
133 | | - |
134 | | - |
135 | | -@unittest.skipIf( |
136 | | - not parmest.parmest_available, |
137 | | - "Cannot test parmest: required dependencies are missing", |
138 | | -) |
139 | | -@unittest.skipIf(not ipopt_available, "The 'ipopt' solver is not available") |
140 | | -class TestExceptionReactorDesign_DAE(unittest.TestCase): |
141 | | - # Based on a reactor example in `Chemical Reactor Analysis and Design Fundamentals`, |
142 | | - # https://sites.engineering.ucsb.edu/~jbraw/chemreacfun/ |
143 | | - # https://sites.engineering.ucsb.edu/~jbraw/chemreacfun/fig-html/appendix/fig-A-10.html |
144 | | - |
145 | | - def setUp(self): |
146 | | - def ABC_model(data): |
147 | | - |
148 | | - if isinstance(data, pd.DataFrame): |
149 | | - meas_t = data.index # time index |
150 | | - else: # dictionary |
151 | | - meas_t = list(data["ca"].keys()) # nested dictionary |
152 | | - |
153 | | - ca0 = 1.0 |
154 | | - cb0 = 0.0 |
155 | | - cc0 = 0.0 |
156 | | - |
157 | | - m = pyo.ConcreteModel() |
158 | | - |
159 | | - m.k1 = pyo.Var(initialize=0.5, bounds=(1e-4, 10)) |
160 | | - m.k2 = pyo.Var(initialize=3.0, bounds=(1e-4, 10)) |
161 | | - |
162 | | - m.time = dae.ContinuousSet(bounds=(0.0, 5.0), initialize=meas_t) |
163 | | - |
164 | | - # initialization and bounds |
165 | | - m.ca = pyo.Var(m.time, initialize=ca0, bounds=(-1e-3, ca0 + 1e-3)) |
166 | | - m.cb = pyo.Var(m.time, initialize=cb0, bounds=(-1e-3, ca0 + 1e-3)) |
167 | | - m.cc = pyo.Var(m.time, initialize=cc0, bounds=(-1e-3, ca0 + 1e-3)) |
168 | | - |
169 | | - m.dca = dae.DerivativeVar(m.ca, wrt=m.time) |
170 | | - m.dcb = dae.DerivativeVar(m.cb, wrt=m.time) |
171 | | - m.dcc = dae.DerivativeVar(m.cc, wrt=m.time) |
172 | | - |
173 | | - def _dcarate(m, t): |
174 | | - if t == 0: |
175 | | - return pyo.Constraint.Skip |
176 | | - else: |
177 | | - return m.dca[t] == -m.k1 * m.ca[t] |
178 | | - |
179 | | - m.dcarate = pyo.Constraint(m.time, rule=_dcarate) |
180 | | - |
181 | | - def _dcbrate(m, t): |
182 | | - if t == 0: |
183 | | - return pyo.Constraint.Skip |
184 | | - else: |
185 | | - return m.dcb[t] == m.k1 * m.ca[t] - m.k2 * m.cb[t] |
186 | | - |
187 | | - m.dcbrate = pyo.Constraint(m.time, rule=_dcbrate) |
188 | | - |
189 | | - def _dccrate(m, t): |
190 | | - if t == 0: |
191 | | - return pyo.Constraint.Skip |
192 | | - else: |
193 | | - return m.dcc[t] == m.k2 * m.cb[t] |
194 | | - |
195 | | - m.dccrate = pyo.Constraint(m.time, rule=_dccrate) |
196 | | - |
197 | | - return m |
198 | | - |
199 | | - class ReactorDesignExperimentDAE(Experiment): |
200 | | - |
201 | | - def __init__(self, data): |
202 | | - |
203 | | - self.data = data |
204 | | - self.model = None |
205 | | - |
206 | | - def create_model(self): |
207 | | - self.model = ABC_model(self.data) |
208 | | - |
209 | | - def label_model(self): |
210 | | - |
211 | | - m = self.model |
212 | | - |
213 | | - if isinstance(self.data, pd.DataFrame): |
214 | | - meas_time_points = self.data.index |
215 | | - else: |
216 | | - meas_time_points = list(self.data["ca"].keys()) |
217 | | - |
218 | | - m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL) |
219 | | - m.experiment_outputs.update( |
220 | | - (m.ca[t], self.data["ca"][t]) for t in meas_time_points |
221 | | - ) |
222 | | - m.experiment_outputs.update( |
223 | | - (m.cb[t], self.data["cb"][t]) for t in meas_time_points |
224 | | - ) |
225 | | - m.experiment_outputs.update( |
226 | | - (m.cc[t], self.data["cc"][t]) for t in meas_time_points |
227 | | - ) |
228 | | - |
229 | | - def get_labeled_model(self): |
230 | | - self.create_model() |
231 | | - self.label_model() |
232 | | - |
233 | | - return self.model |
234 | | - |
235 | | - # This example tests data formatted in 3 ways |
236 | | - # Each format holds 1 scenario |
237 | | - # 1. dataframe with time index |
238 | | - # 2. nested dictionary {ca: {t, val pairs}, ... } |
239 | | - data = [ |
240 | | - [0.000, 0.957, -0.031, -0.015], |
241 | | - [0.263, 0.557, 0.330, 0.044], |
242 | | - [0.526, 0.342, 0.512, 0.156], |
243 | | - [0.789, 0.224, 0.499, 0.310], |
244 | | - [1.053, 0.123, 0.428, 0.454], |
245 | | - [1.316, 0.079, 0.396, 0.556], |
246 | | - [1.579, 0.035, 0.303, 0.651], |
247 | | - [1.842, 0.029, 0.287, 0.658], |
248 | | - [2.105, 0.025, 0.221, 0.750], |
249 | | - [2.368, 0.017, 0.148, 0.854], |
250 | | - [2.632, -0.002, 0.182, 0.845], |
251 | | - [2.895, 0.009, 0.116, 0.893], |
252 | | - [3.158, -0.023, 0.079, 0.942], |
253 | | - [3.421, 0.006, 0.078, 0.899], |
254 | | - [3.684, 0.016, 0.059, 0.942], |
255 | | - [3.947, 0.014, 0.036, 0.991], |
256 | | - [4.211, -0.009, 0.014, 0.988], |
257 | | - [4.474, -0.030, 0.036, 0.941], |
258 | | - [4.737, 0.004, 0.036, 0.971], |
259 | | - [5.000, -0.024, 0.028, 0.985], |
260 | | - ] |
261 | | - data = pd.DataFrame(data, columns=["t", "ca", "cb", "cc"]) |
262 | | - data_df = data.set_index("t") |
263 | | - data_dict = { |
264 | | - "ca": {k: v for (k, v) in zip(data.t, data.ca)}, |
265 | | - "cb": {k: v for (k, v) in zip(data.t, data.cb)}, |
266 | | - "cc": {k: v for (k, v) in zip(data.t, data.cc)}, |
267 | | - } |
268 | | - |
269 | | - # Create an experiment list |
270 | | - self.exp_list_df = [ReactorDesignExperimentDAE(data_df)] |
271 | | - self.exp_list_dict = [ReactorDesignExperimentDAE(data_dict)] |
272 | | - |
273 | | - def test_parmest_exception(self): |
274 | | - """ |
275 | | - Test the exception raised by parmest as a result of not defining |
276 | | - the "unknown_parameters" attribute |
277 | | - """ |
278 | | - with self.assertRaises(RuntimeError) as context: |
279 | | - parmest.Estimator(self.exp_list_df, obj_function="SSE") |
280 | | - |
281 | | - self.assertIn("unknown_parameters", str(context.exception)) |
282 | | - |
283 | | - with self.assertRaises(RuntimeError) as context: |
284 | | - parmest.Estimator(self.exp_list_dict, obj_function="SSE") |
285 | | - |
286 | | - self.assertIn("unknown_parameters", str(context.exception)) |
287 | | - |
288 | | - |
289 | 37 | @unittest.skipIf( |
290 | 38 | not parmest.parmest_available, |
291 | 39 | "Cannot test parmest: required dependencies are missing", |
@@ -327,6 +75,39 @@ def SSE(model): |
327 | 75 | exp_list, obj_function=SSE, solver_options=solver_options, tee=True |
328 | 76 | ) |
329 | 77 |
|
| 78 | + def test_parmest_exception(self): |
| 79 | + """ |
| 80 | + Test the exception raised by parmest when the "experiment_outputs" |
| 81 | + attribute is not defined in the model |
| 82 | + """ |
| 83 | + from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import ( |
| 84 | + RooneyBieglerExperiment, |
| 85 | + ) |
| 86 | + |
| 87 | + # create an instance of the RooneyBieglerExperiment class |
| 88 | + # without the "experiment_outputs" attribute |
| 89 | + class RooneyBieglerExperimentException(RooneyBieglerExperiment): |
| 90 | + def label_model(self): |
| 91 | + m = self.model |
| 92 | + |
| 93 | + # add the unknown parameters |
| 94 | + m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL) |
| 95 | + m.unknown_parameters.update( |
| 96 | + (k, pyo.ComponentUID(k)) for k in [m.asymptote, m.rate_constant] |
| 97 | + ) |
| 98 | + |
| 99 | + # create an experiment list |
| 100 | + exp_list = [] |
| 101 | + for i in range(self.data.shape[0]): |
| 102 | + exp_list.append(RooneyBieglerExperimentException(self.data.loc[i, :])) |
| 103 | + |
| 104 | + # check the exception raised by parmest due to not defining |
| 105 | + # the "experiment_outputs" |
| 106 | + with self.assertRaises(RuntimeError) as context: |
| 107 | + parmest.Estimator(exp_list, obj_function="SSE", tee=True) |
| 108 | + |
| 109 | + self.assertIn("experiment_outputs", str(context.exception)) |
| 110 | + |
330 | 111 | def test_theta_est(self): |
331 | 112 | objval, thetavals = self.pest.theta_est() |
332 | 113 |
|
@@ -1137,6 +918,51 @@ def get_labeled_model(self): |
1137 | 918 | self.m_df = ABC_model(data_df) |
1138 | 919 | self.m_dict = ABC_model(data_dict) |
1139 | 920 |
|
| 921 | + # create an instance of the ReactorDesignExperimentDAE class |
| 922 | + # without the "unknown_parameters" attribute |
| 923 | + class ReactorDesignExperimentException(ReactorDesignExperimentDAE): |
| 924 | + def label_model(self): |
| 925 | + |
| 926 | + m = self.model |
| 927 | + |
| 928 | + if isinstance(self.data, pd.DataFrame): |
| 929 | + meas_time_points = self.data.index |
| 930 | + else: |
| 931 | + meas_time_points = list(self.data["ca"].keys()) |
| 932 | + |
| 933 | + m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL) |
| 934 | + m.experiment_outputs.update( |
| 935 | + (m.ca[t], self.data["ca"][t]) for t in meas_time_points |
| 936 | + ) |
| 937 | + m.experiment_outputs.update( |
| 938 | + (m.cb[t], self.data["cb"][t]) for t in meas_time_points |
| 939 | + ) |
| 940 | + m.experiment_outputs.update( |
| 941 | + (m.cc[t], self.data["cc"][t]) for t in meas_time_points |
| 942 | + ) |
| 943 | + |
| 944 | + # create an experiment list without the "unknown_parameters" attribute |
| 945 | + exp_list_df_no_params = [ReactorDesignExperimentException(data_df)] |
| 946 | + exp_list_dict_no_params = [ReactorDesignExperimentException(data_dict)] |
| 947 | + |
| 948 | + self.exp_list_df_no_params = exp_list_df_no_params |
| 949 | + self.exp_list_dict_no_params = exp_list_dict_no_params |
| 950 | + |
| 951 | + def test_parmest_exception(self): |
| 952 | + """ |
| 953 | + Test the exception raised by parmest when the "unknown_parameters" |
| 954 | + attribute is not defined in the model |
| 955 | + """ |
| 956 | + with self.assertRaises(RuntimeError) as context: |
| 957 | + parmest.Estimator(self.exp_list_df_no_params, obj_function="SSE") |
| 958 | + |
| 959 | + self.assertIn("unknown_parameters", str(context.exception)) |
| 960 | + |
| 961 | + with self.assertRaises(RuntimeError) as context: |
| 962 | + parmest.Estimator(self.exp_list_dict_no_params, obj_function="SSE") |
| 963 | + |
| 964 | + self.assertIn("unknown_parameters", str(context.exception)) |
| 965 | + |
1140 | 966 | def test_dataformats(self): |
1141 | 967 | obj1, theta1 = self.pest_df.theta_est() |
1142 | 968 | obj2, theta2 = self.pest_dict.theta_est() |
|
0 commit comments