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

Pyomo / pyomo

Bug of running examples of Pyomo.DoE #3296

Summary

Steps to reproduce the issue from the website until step 3, and run main function

=== Required import ===

Error Message