Skip to content

Commit 0a5e74b

Browse files
slilonfe5slilonfe5
committed
Updated covariance.rst, driver.rst, and parmest.py files
1 parent aa5a92a commit 0a5e74b

3 files changed

Lines changed: 152 additions & 19 deletions

File tree

‎doc/OnlineDocs/explanation/analysis/parmest/covariance.rst‎

Lines changed: 68 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -63,11 +63,77 @@ of data points, e.g.,
6363
)
6464
num_data = len(data)
6565

66+
# Create the Rooney-Biegler model
67+
def rooney_biegler_model():
68+
"""
69+
Formulates the Pyomo model of the Rooney-Biegler example
70+
71+
Returns:
72+
m: Pyomo model
73+
"""
74+
m = pyo.ConcreteModel()
75+
76+
m.asymptote = pyo.Var(within=pyo.NonNegativeReals, initialize=10)
77+
m.rate_constant = pyo.Var(within=pyo.NonNegativeReals, initialize=0.2)
78+
79+
m.hour = pyo.Var(within=pyo.PositiveReals, initialize=0.1)
80+
m.y = pyo.Var(within=pyo.NonNegativeReals)
81+
82+
@m.Constraint()
83+
def response_rule(m):
84+
return m.y == m.asymptote * (1 - pyo.exp(-m.rate_constant * m.hour))
85+
86+
return m
87+
88+
# Create the Experiment class
89+
from pyomo.contrib.parmest.experiment import Experiment
90+
class RooneyBieglerExperiment(Experiment):
91+
def __init__(self, hour, y):
92+
self.y = y
93+
self.hour = hour
94+
self.model = None
95+
96+
def get_labeled_model(self):
97+
self.create_model()
98+
self.finalize_model()
99+
self.label_model()
100+
101+
return self.model
102+
103+
def create_model(self):
104+
m = self.model = rooney_biegler_model()
105+
106+
return m
107+
108+
def finalize_model(self):
109+
m = self.model
110+
111+
# fix the input variable
112+
m.hour.fix(self.hour)
113+
114+
return m
115+
116+
def label_model(self):
117+
m = self.model
118+
119+
# add experiment outputs
120+
m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
121+
m.experiment_outputs.update([(m.y, self.y)])
122+
123+
# add unknown parameters
124+
m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL)
125+
m.unknown_parameters.update(
126+
(k, pyo.value(k)) for k in [m.asymptote, m.rate_constant]
127+
)
128+
129+
# create the measurement error
130+
m.measurement_error = pyo.Suffix(direction = pyo.Suffix.LOCAL)
131+
m.measurement_error.update([(m.y, None)])
132+
66133
# Create an experiment list
67-
from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import RooneyBieglerExperiment
68134
exp_list = []
69135
for i in range(data.shape[0]):
70-
exp_list.append(RooneyBieglerExperiment(data.loc[i, :]))
136+
exp_list.append(RooneyBieglerExperiment(data["hour"][i], data["y"][i]))
71137

72138
.. doctest::
73139
:skipif: not __import__('pyomo.contrib.parmest.parmest').contrib.parmest.parmest.parmest_available

‎doc/OnlineDocs/explanation/analysis/parmest/driver.rst‎

Lines changed: 68 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -67,11 +67,77 @@ Section.
6767
columns=['hour', 'y'],
6868
)
6969

70+
# Create the Rooney-Biegler model
71+
def rooney_biegler_model():
72+
"""
73+
Formulates the Pyomo model of the Rooney-Biegler example
74+
75+
Returns:
76+
m: Pyomo model
77+
"""
78+
m = pyo.ConcreteModel()
79+
80+
m.asymptote = pyo.Var(within=pyo.NonNegativeReals, initialize=10)
81+
m.rate_constant = pyo.Var(within=pyo.NonNegativeReals, initialize=0.2)
82+
83+
m.hour = pyo.Var(within=pyo.PositiveReals, initialize=0.1)
84+
m.y = pyo.Var(within=pyo.NonNegativeReals)
85+
86+
@m.Constraint()
87+
def response_rule(m):
88+
return m.y == m.asymptote * (1 - pyo.exp(-m.rate_constant * m.hour))
89+
90+
return m
91+
92+
# Create the Experiment class
93+
from pyomo.contrib.parmest.experiment import Experiment
94+
class RooneyBieglerExperiment(Experiment):
95+
def __init__(self, hour, y):
96+
self.y = y
97+
self.hour = hour
98+
self.model = None
99+
100+
def get_labeled_model(self):
101+
self.create_model()
102+
self.finalize_model()
103+
self.label_model()
104+
105+
return self.model
106+
107+
def create_model(self):
108+
m = self.model = rooney_biegler_model()
109+
110+
return m
111+
112+
def finalize_model(self):
113+
m = self.model
114+
115+
# fix the input variable
116+
m.hour.fix(self.hour)
117+
118+
return m
119+
120+
def label_model(self):
121+
m = self.model
122+
123+
# add experiment outputs
124+
m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
125+
m.experiment_outputs.update([(m.y, self.y)])
126+
127+
# add unknown parameters
128+
m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL)
129+
m.unknown_parameters.update(
130+
(k, pyo.value(k)) for k in [m.asymptote, m.rate_constant]
131+
)
132+
133+
# create the measurement error
134+
m.measurement_error = pyo.Suffix(direction = pyo.Suffix.LOCAL)
135+
m.measurement_error.update([(m.y, None)])
136+
70137
# Create an experiment list
71-
from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import RooneyBieglerExperiment
72138
exp_list = []
73139
for i in range(data.shape[0]):
74-
exp_list.append(RooneyBieglerExperiment(data.loc[i, :]))
140+
exp_list.append(RooneyBieglerExperiment(data["hour"][i], data["y"][i]))
75141

76142
.. doctest::
77143
:skipif: not __import__('pyomo.contrib.parmest.parmest').contrib.parmest.parmest.parmest_available

‎pyomo/contrib/parmest/parmest.py‎

Lines changed: 16 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -1234,23 +1234,24 @@ def _cov_at_theta(self, method, solver, cov_n, step):
12341234
Returns:
12351235
cov: pd.DataFrame, covariance matrix of the estimated parameters
12361236
"""
1237-
# compute the inverse reduced hessian to be used
1238-
# in the "reduced_hessian" method
1239-
# parmest makes the fitted parameters stage 1 variables
1240-
ind_vars = []
1241-
for nd_name, Var, sol_val in ef_nonants(self.ef_instance):
1242-
ind_vars.append(Var)
1243-
# calculate the reduced hessian
1244-
(solve_result, inv_red_hes) = (
1245-
inverse_reduced_hessian.inv_reduced_hessian_barrier(
1246-
self.ef_instance,
1247-
independent_variables=ind_vars,
1248-
solver_options=self.solver_options,
1249-
tee=self.tee,
1237+
if method == CovarianceMethodLib.reduced_hessian.value:
1238+
# compute the inverse reduced hessian to be used
1239+
# in the "reduced_hessian" method
1240+
# parmest makes the fitted parameters stage 1 variables
1241+
ind_vars = []
1242+
for nd_name, Var, sol_val in ef_nonants(self.ef_instance):
1243+
ind_vars.append(Var)
1244+
# calculate the reduced hessian
1245+
(solve_result, inv_red_hes) = (
1246+
inverse_reduced_hessian.inv_reduced_hessian_barrier(
1247+
self.ef_instance,
1248+
independent_variables=ind_vars,
1249+
solver_options=self.solver_options,
1250+
tee=self.tee,
1251+
)
12501252
)
1251-
)
12521253

1253-
self.inv_red_hes = inv_red_hes
1254+
self.inv_red_hes = inv_red_hes
12541255

12551256
# Number of data points considered
12561257
n = cov_n

0 commit comments

Comments
 (0)