Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Changes to address Pyomo.DoE edits #16

Merged
merged 5 commits into from
May 23, 2023
Merged

Changes to address Pyomo.DoE edits #16

Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
0986950
Removing unused comment
adowling2 May 23, 2023
e4f46f0
Address model Pyomo.DoE review comments
adowling2 May 23, 2023
0262868
Added missing m: Pyomo model from doc string
adowling2 May 23, 2023
1f61bee
Fix black formatting error?
adowling2 May 23, 2023
5c0405a
Fixed variable name typo
adowling2 May 23, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
79 changes: 42 additions & 37 deletions pyomo/contrib/doe/doe.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -560,6 +560,10 @@ def _direct_kaug(self):
def _create_block(self):
"""
Create a pyomo Concrete model and add blocks with different parameter perturbation scenarios.

Returns
-------
mod: Concrete Pyomo model
"""

# create scenario information for block scenarios
Expand Down Expand Up @@ -770,7 +774,6 @@ def run_grid_search(
design_iter = self.design_vars.variable_names_value.copy()
# update the controlled value of certain time points for certain design variables
for i, names in enumerate(design_dimension_names):
# names = design_dimension_names[i]
# 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:
for n in names:
Expand Down Expand Up @@ -864,51 +867,51 @@ def _create_doe_model(self, no_obj=True):

Return
-------
m: the DOE model
model: the DOE model
"""
mod = self._create_block()
model = self._create_block()

# variables for jacobian and FIM
mod.regression_parameters = pyo.Set(initialize=list(self.param.keys()))
mod.measured_variables = pyo.Set(initialize=self.measure_name)
model.regression_parameters = pyo.Set(initialize=list(self.param.keys()))
model.measured_variables = pyo.Set(initialize=self.measure_name)

def identity_matrix(m, i, j):
if i == j:
return 1
else:
return 0

mod.sensitivity_jacobian = pyo.Var(
mod.regression_parameters, mod.measured_variables, initialize=0.1
model.sensitivity_jacobian = pyo.Var(
model.regression_parameters, model.measured_variables, initialize=0.1
)

if self.fim_initial:
dict_fim_initialize = {}
for i, bu in enumerate(mod.regression_parameters):
for j, un in enumerate(mod.regression_parameters):
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]

def initialize_fim(m, j, d):
return dict_fim_initialize[(j, d)]

if self.fim_initial:
mod.fim = pyo.Var(
mod.regression_parameters,
mod.regression_parameters,
model.fim = pyo.Var(
model.regression_parameters,
model.regression_parameters,
initialize=initialize_fim,
)
else:
mod.fim = pyo.Var(
mod.regression_parameters,
mod.regression_parameters,
model.fim = pyo.Var(
model.regression_parameters,
model.regression_parameters,
initialize=identity_matrix,
)

# move the L matrix initial point to a dictionary
if type(self.L_initial) != type(None):
dict_cho = {}
for i, bu in enumerate(mod.regression_parameters):
for j, un in enumerate(mod.regression_parameters):
for i, bu in enumerate(model.regression_parameters):
for j, un in enumerate(model.regression_parameters):
dict_cho[(bu, un)] = self.L_initial[i][j]

# use the L dictionary to initialize L matrix
Expand All @@ -920,33 +923,34 @@ def init_cho(m, i, j):
# Define elements of Cholesky decomposition matrix as Pyomo variables and either
# Initialize with L in L_initial
if type(self.L_initial) != type(None):
mod.L_ele = pyo.Var(
mod.regression_parameters,
mod.regression_parameters,
model.L_ele = pyo.Var(
model.regression_parameters,
model.regression_parameters,
initialize=init_cho,
)
# or initialize with the identity matrix
else:
mod.L_ele = pyo.Var(
mod.regression_parameters,
mod.regression_parameters,
model.L_ele = pyo.Var(
model.regression_parameters,
model.regression_parameters,
initialize=identity_matrix,
)

# loop over parameter name
for i, c in enumerate(mod.regression_parameters):
for j, d in enumerate(mod.regression_parameters):
for i, c in enumerate(model.regression_parameters):
for j, d in enumerate(model.regression_parameters):
# fix the 0 half of L matrix to be 0.0
if i < j:
mod.L_ele[c, d].fix(0.0)
model.L_ele[c, d].fix(0.0)
# Give LB to the diagonal entries
if self.L_LB:
if c == d:
mod.L_ele[c, d].setlb(self.L_LB)
model.L_ele[c, d].setlb(self.L_LB)

# jacobian rule
def jacobian_rule(m, p, n):
"""
m: Pyomo model
p: parameter
n: response
"""
Expand All @@ -970,19 +974,20 @@ def jacobian_rule(m, p, n):
# A constraint to calculate elements in Hessian matrix
# transfer prior FIM to be Expressions
fim_initial_dict = {}
for i, bu in enumerate(mod.regression_parameters):
for j, un in enumerate(mod.regression_parameters):
for i, bu in enumerate(model.regression_parameters):
for j, un in enumerate(model.regression_parameters):
fim_initial_dict[(bu, un)] = self.prior_FIM[i][j]

def read_prior(m, i, j):
return fim_initial_dict[(i, j)]

mod.priorFIM = pyo.Expression(
mod.regression_parameters, mod.regression_parameters, rule=read_prior
model.priorFIM = pyo.Expression(
model.regression_parameters, model.regression_parameters, rule=read_prior
)

def fim_rule(m, p, q):
"""
m: Pyomo model
p: parameter
q: parameter
"""
Expand All @@ -993,19 +998,19 @@ def fim_rule(m, p, q):
/ self.measurement_vars.variance[n]
* m.sensitivity_jacobian[p, n]
* m.sensitivity_jacobian[q, n]
for n in mod.measured_variables
for n in model.measured_variables
)
+ m.priorFIM[p, q] * self.fim_scale_constant_value
)

mod.jacobian_constraint = pyo.Constraint(
mod.regression_parameters, mod.measured_variables, rule=jacobian_rule
model.jacobian_constraint = pyo.Constraint(
model.regression_parameters, model.measured_variables, rule=jacobian_rule
)
mod.fim_constraint = pyo.Constraint(
mod.regression_parameters, mod.regression_parameters, rule=fim_rule
model.fim_constraint = pyo.Constraint(
model.regression_parameters, model.regression_parameters, rule=fim_rule
)

return mod
return model

def _add_objective(self, m):
def cholesky_imp(m, c, d):
Expand Down