pyomo
pyomo copied to clipboard
Pyomo
Bug of running examples of Pyomo.DoE
Summary
I am trying to run the examples on the website of Pyomo.DoE (https://pyomo.readthedocs.io/en/stable/contributed_packages/doe/doe.html). I copied the code until step 3, and I got errors saying that ModelOptionLib is not defined.
Steps to reproduce the issue from the website until step 3, and run main function
=== Required import ===
import pyomo.environ as pyo from pyomo.dae import ContinuousSet, DerivativeVar from pyomo.contrib.doe import DesignOfExperiments, MeasurementVariables, DesignVariables import numpy as np
def create_model( mod=None, model_option="stage2", control_time=[0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1], control_val=None, t_range=[0.0, 1], CA_init=1, C_init=0.1, ): """ This is an example user model provided to DoE library. It is a dynamic problem solved by Pyomo.DAE.
Arguments
---------
mod: Pyomo model. If None, a Pyomo concrete model is created
model_option: choose from the 3 options in model_option
if ModelOptionLib.parmest, create a process model.
if ModelOptionLib.stage1, create the global model.
if ModelOptionLib.stage2, add model variables and constraints for block.
control_time: a list of control timepoints
control_val: control design variable values T at corresponding timepoints
t_range: time range, h
CA_init: time-independent design (control) variable, an initial value for CA
C_init: An initial value for C
Return
------
m: a Pyomo.DAE model
"""
theta = {"A1": 84.79, "A2": 371.72, "E1": 7.78, "E2": 15.05}
model_option = ModelOptionLib(model_option)
if model_option == ModelOptionLib.parmest:
mod = pyo.ConcreteModel()
return_m = True
elif model_option == ModelOptionLib.stage1 or model_option == ModelOptionLib.stage2:
if not mod:
raise ValueError(
"If model option is stage1 or stage2, a created model needs to be provided."
)
return_m = False
else:
raise ValueError(
"model_option needs to be defined as parmest,stage1, or stage2."
)
if not control_val:
control_val = [300] * 9
controls = {}
for i, t in enumerate(control_time):
controls[t] = control_val[i]
mod.t0 = pyo.Set(initialize=[0])
mod.t_con = pyo.Set(initialize=control_time)
mod.CA0 = pyo.Var(
mod.t0, initialize=CA_init, bounds=(1.0, 5.0), within=pyo.NonNegativeReals
) # mol/L
# check if control_time is in time range
assert (
control_time[0] >= t_range[0] and control_time[-1] <= t_range[1]
), "control time is outside time range."
if model_option == ModelOptionLib.stage1:
mod.T = pyo.Var(
mod.t_con,
initialize=controls,
bounds=(300, 700),
within=pyo.NonNegativeReals,
)
return
else:
para_list = ["A1", "A2", "E1", "E2"]
### Add variables
mod.CA_init = CA_init
mod.para_list = para_list
# timepoints
mod.t = ContinuousSet(bounds=t_range, initialize=control_time)
# time-dependent design variable, initialized with the first control value
def T_initial(m, t):
if t in m.t_con:
return controls[t]
else:
# count how many control points are before the current t;
# locate the nearest neighbouring control point before this t
neighbour_t = max(tc for tc in control_time if tc < t)
return controls[neighbour_t]
mod.T = pyo.Var(
mod.t, initialize=T_initial, bounds=(300, 700), within=pyo.NonNegativeReals
)
mod.R = 8.31446261815324 # J / K / mole
# Define parameters as Param
mod.A1 = pyo.Var(initialize=theta["A1"])
mod.A2 = pyo.Var(initialize=theta["A2"])
mod.E1 = pyo.Var(initialize=theta["E1"])
mod.E2 = pyo.Var(initialize=theta["E2"])
# Concentration variables under perturbation
mod.C_set = pyo.Set(initialize=["CA", "CB", "CC"])
mod.C = pyo.Var(
mod.C_set, mod.t, initialize=C_init, within=pyo.NonNegativeReals
)
# time derivative of C
mod.dCdt = DerivativeVar(mod.C, wrt=mod.t)
# kinetic parameters
def kp1_init(m, t):
return m.A1 * pyo.exp(-m.E1 * 1000 / (m.R * m.T[t]))
def kp2_init(m, t):
return m.A2 * pyo.exp(-m.E2 * 1000 / (m.R * m.T[t]))
mod.kp1 = pyo.Var(mod.t, initialize=kp1_init)
mod.kp2 = pyo.Var(mod.t, initialize=kp2_init)
def T_control(m, t):
"""
T at interval timepoint equal to the T of the control time point at the beginning of this interval
Count how many control points are before the current t;
locate the nearest neighbouring control point before this t
"""
if t in m.t_con:
return pyo.Constraint.Skip
else:
neighbour_t = max(tc for tc in control_time if tc < t)
return m.T[t] == m.T[neighbour_t]
def cal_kp1(m, t):
"""
Create the perturbation parameter sets
m: model
t: time
"""
# LHS: 1/h
# RHS: 1/h*(kJ/mol *1000J/kJ / (J/mol/K) / K)
return m.kp1[t] == m.A1 * pyo.exp(-m.E1 * 1000 / (m.R * m.T[t]))
def cal_kp2(m, t):
"""
Create the perturbation parameter sets
m: model
t: time
"""
# LHS: 1/h
# RHS: 1/h*(kJ/mol *1000J/kJ / (J/mol/K) / K)
return m.kp2[t] == m.A2 * pyo.exp(-m.E2 * 1000 / (m.R * m.T[t]))
def dCdt_control(m, y, t):
"""
Calculate CA in Jacobian matrix analytically
y: CA, CB, CC
t: timepoints
"""
if y == "CA":
return m.dCdt[y, t] == -m.kp1[t] * m.C["CA", t]
elif y == "CB":
return m.dCdt[y, t] == m.kp1[t] * m.C["CA", t] - m.kp2[t] * m.C["CB", t]
elif y == "CC":
return pyo.Constraint.Skip
def alge(m, t):
"""
The algebraic equation for mole balance
z: m.pert
t: time
"""
return m.C["CA", t] + m.C["CB", t] + m.C["CC", t] == m.CA0[0]
# Control time
mod.T_rule = pyo.Constraint(mod.t, rule=T_control)
# calculating C, Jacobian, FIM
mod.k1_pert_rule = pyo.Constraint(mod.t, rule=cal_kp1)
mod.k2_pert_rule = pyo.Constraint(mod.t, rule=cal_kp2)
mod.dCdt_rule = pyo.Constraint(mod.C_set, mod.t, rule=dCdt_control)
mod.alge_rule = pyo.Constraint(mod.t, rule=alge)
# B.C.
mod.C["CB", 0.0].fix(0.0)
mod.C["CC", 0.0].fix(0.0)
if return_m:
return mod
def disc_for_measure(m, nfe=32, block=True): """Pyomo.DAE discretization
Arguments
---------
m: Pyomo model
nfe: number of finite elements b
block: if True, the input model has blocks
"""
discretizer = pyo.TransformationFactory("dae.collocation")
if block:
for s in range(len(m.block)):
discretizer.apply_to(m.block[s], nfe=nfe, ncp=3, wrt=m.block[s].t)
else:
discretizer.apply_to(m, nfe=nfe, ncp=3, wrt=m.t)
return m
# Control time set [h]
t_control = [0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1]
# Define parameter nominal value
parameter_dict = {"A1": 85, "A2": 370, "E1": 8, "E2": 15}
# Define measurement object
measurements = MeasurementVariables()
measurements.add_variables(
"C", # measurement variable name
indices={
0: ["CA", "CB", "CC"],
1: t_control,
}, # 0,1 are indices of the index sets
time_index_position=1,
)
# design object
exp_design = DesignVariables()
# add CAO as design variable
exp_design.add_variables(
"CA0", # design variable name
indices={0: [0]}, # index dictionary
time_index_position=0, # time index position
values=[5], # design variable values
lower_bounds=1, # design variable lower bounds
upper_bounds=5, # design variable upper bounds
)
# add T as design variable
exp_design.add_variables(
"T", # design variable name
indices={0: t_control}, # index dictionary
time_index_position=0, # time index position
values=[
570,
300,
300,
300,
300,
300,
300,
300,
300,
], # same length with t_control
lower_bounds=300, # design variable lower bounds
upper_bounds=700, # design variable upper bounds
)
main()
Error Message
INFO: =======Iteration Number: 1 ===== WARNING: :::::::::::Warning: Cannot converge this run.:::::::::::: INFO: =======Iteration Number: 2 ===== WARNING: :::::::::::Warning: Cannot converge this run.:::::::::::: INFO: =======Iteration Number: 3 ===== WARNING: :::::::::::Warning: Cannot converge this run.:::::::::::: INFO: =======Iteration Number: 4 ===== WARNING: :::::::::::Warning: Cannot converge this run.:::::::::::: INFO: =======Iteration Number: 5 ===== WARNING: :::::::::::Warning: Cannot converge this run.:::::::::::: INFO: =======Iteration Number: 6 ===== WARNING: :::::::::::Warning: Cannot converge this run.:::::::::::: INFO: =======Iteration Number: 7 ===== WARNING: :::::::::::Warning: Cannot converge this run.:::::::::::: INFO: =======Iteration Number: 8 ===== WARNING: :::::::::::Warning: Cannot converge this run.:::::::::::: INFO: =======Iteration Number: 9 ===== WARNING: :::::::::::Warning: Cannot converge this run.:::::::::::: INFO: Overall wall clock time [s]: 0 --- Logging error --- Traceback (most recent call last): File "c:\Users\LINL\miniconda3\envs\orca-env\lib\site-packages\pyomo\contrib\doe\doe.py", line 804, in run_grid_search result_iter = self.compute_FIM( File "c:\Users\LINL\miniconda3\envs\orca-env\lib\site-packages\pyomo\contrib\doe\doe.py", line 387, in compute_FIM FIM_analysis = self._direct_kaug() File "c:\Users\LINL\miniconda3\envs\orca-env\lib\site-packages\pyomo\contrib\doe\doe.py", line 475, in _direct_kaug mod = self.create_model(model_option=ModelOptionLib.parmest) File "C:\Users\LINL\AppData\Local\Temp\1\ipykernel_24524\3408820592.py", line 34, in create_model model_option = ModelOptionLib(model_option) NameError: name 'ModelOptionLib' is not defined
